IceRecoBase.cxx

Go to the documentation of this file.

 
#include "IceRecoBase.h"
#include "Riostream.h"
 
ClassImp(IceRecoBase); // Class implementation to enable ROOT I/O

IceRecoBase::IceRecoBase(const char* name,const char* title) : TTask(name,title)
{
 
 fFirst=1;
 fEvt=0;
 fCleanA=0;
 fMaxmodA=999999;
 fMinmodA=0;
 fMaxhitsA=0;
 fSingleA=0;
 fMinahitsA=0;
 fMinamodsA=0;
 fCleanI=0;
 fMaxmodI=999999;
 fMinmodI=0;
 fMaxhitsI=0;
 fSingleI=0;
 fMinahitsI=0;
 fMinamodsI=0;
 fCleanIC=0;
 fMaxmodIC=999999;
 fMinmodIC=0;
 fMaxhitsIC=0;
 fSingleIC=0;
 fMinahitsIC=0;
 fMinamodsIC=0;
 fCleanDC=0;
 fMaxmodDC=999999;
 fMinmodDC=0;
 fMaxhitsDC=0;
 fSingleDC=0;
 fMinahitsDC=0;
 fMinamodsDC=0;
 fSlcI=1;
 fSlcIC=1;
 fSlcDC=1;
 fThetatrk=-999;
 fThetahits=999;
 fLambdaA=33.3;
 fLambdaUD=30;
 fLambdaDL=5;
 fLambdaLD=40;
 fLabsA=50;
 fLabsUD=100;
 fLabsDL=10;
 fLabsLD=150;
 fTsigmaA=10;
 fTsigmaIC=5;
 fTsigmaDC=5;
 fTrackname="";
 fCharge=0;
 fUseNames=0;
 fUseNtk=0;
}

IceRecoBase::~IceRecoBase()
{
 
 if (fUseNames)
 {
  delete fUseNames;
  fUseNames=0;
 }
 if (fUseNtk)
 {
  delete fUseNtk;
  fUseNtk=0;
 }
}

void IceRecoBase::SetCleaned(Int_t flag,TString s)
{
 
 if (s=="A")
 {
  fCleanA=flag;
  fParams.AddNamedSlot("CleanA");
  fParams.SetSignal(fCleanA,"CleanA");
 }
 if (s=="I")
 {
  fCleanI=flag;
  fParams.AddNamedSlot("CleanI");
  fParams.SetSignal(fCleanI,"CleanI");
 }
 if (s=="IC")
 {
  fCleanIC=flag;
  fParams.AddNamedSlot("CleanIC");
  fParams.SetSignal(fCleanIC,"CleanIC");
 }
 if (s=="DC")
 {
  fCleanDC=flag;
  fParams.AddNamedSlot("CleanDC");
  fParams.SetSignal(fCleanDC,"CleanDC");
 }
}

void IceRecoBase::SetMaxMod(Int_t nmax,TString s)
{
 
 if (s=="A")
 {
  fMaxmodA=nmax;
  fParams.AddNamedSlot("MaxmodA");
  fParams.SetSignal(fMaxmodA,"MaxmodA");
 }
 if (s=="I")
 {
  fMaxmodI=nmax;
  fParams.AddNamedSlot("MaxmodI");
  fParams.SetSignal(fMaxmodI,"MaxmodI");
 }
 if (s=="IC")
 {
  fMaxmodIC=nmax;
  fParams.AddNamedSlot("MaxmodIC");
  fParams.SetSignal(fMaxmodIC,"MaxmodIC");
 }
 if (s=="DC")
 {
  fMaxmodDC=nmax;
  fParams.AddNamedSlot("MaxmodDC");
  fParams.SetSignal(fMaxmodDC,"MaxmodDC");
 }
}

void IceRecoBase::SetMinMod(Int_t nmin,TString s)
{
 
 if (s=="A")
 {
  fMinmodA=nmin;
  fParams.AddNamedSlot("MinmodA");
  fParams.SetSignal(fMinmodA,"MinmodA");
 }
 if (s=="I")
 {
  fMinmodI=nmin;
  fParams.AddNamedSlot("MinmodI");
  fParams.SetSignal(fMinmodI,"MinmodI");
 }
 if (s=="IC")
 {
  fMinmodIC=nmin;
  fParams.AddNamedSlot("MinmodIC");
  fParams.SetSignal(fMinmodIC,"MinmodIC");
 }
 if (s=="DC")
 {
  fMinmodDC=nmin;
  fParams.AddNamedSlot("MinmodDC");
  fParams.SetSignal(fMinmodDC,"MinmodDC");
 }
}

void IceRecoBase::SetMaxHits(Int_t nmax,TString s)
{
 
 if (s=="A")
 {
  fMaxhitsA=nmax;
  fParams.AddNamedSlot("MaxhitsA");
  fParams.SetSignal(fMaxhitsA,"MaxhitsA");
 }
 if (s=="I")
 {
  fMaxhitsI=nmax;
  fParams.AddNamedSlot("MaxhitsI");
  fParams.SetSignal(fMaxhitsI,"MaxhitsI");
 }
 if (s=="IC")
 {
  fMaxhitsIC=nmax;
  fParams.AddNamedSlot("MaxhitsIC");
  fParams.SetSignal(fMaxhitsIC,"MaxhitsIC");
 }
 if (s=="DC")
 {
  fMaxhitsDC=nmax;
  fParams.AddNamedSlot("MaxhitsDC");
  fParams.SetSignal(fMaxhitsDC,"MaxhitsDC");
 }
}

void IceRecoBase::SetSingleHit(Int_t ndoms,TString s,Int_t ndoms1)
{
 
 if (s=="A")
 {
  if (ndoms>=0)
  {
   fSingleA=ndoms;
   fParams.AddNamedSlot("SingleA");
   fParams.SetSignal(fSingleA,"SingleA");
  }
  if (ndoms1>=0)
  {
   fSingle1A=ndoms1;
   fParams.AddNamedSlot("Single1A");
   fParams.SetSignal(fSingle1A,"Single1A");
  }
 }
 if (s=="IC")
 {
  if (ndoms>=0)
  {
   fSingleIC=ndoms;
   fParams.AddNamedSlot("SingleIC");
   fParams.SetSignal(fSingleIC,"SingleIC");
  }
  if (ndoms1>=0)
  {
   fSingle1IC=ndoms1;
   fParams.AddNamedSlot("Single1IC");
   fParams.SetSignal(fSingle1IC,"Single1IC");
  }
 }
 if (s=="DC")
 {
  if (ndoms>=0)
  {
   fSingleDC=ndoms;
   fParams.AddNamedSlot("SingleDC");
   fParams.SetSignal(fSingleDC,"SingleDC");
  }
  if (ndoms1>=0)
  {
   fSingle1DC=ndoms1;
   fParams.AddNamedSlot("Single1DC");
   fParams.SetSignal(fSingle1DC,"Single1DC");
  }
 }
 if (s=="I")
 {
  if (ndoms>=0)
  {
   fSingleI=ndoms;
   fParams.AddNamedSlot("SingleI");
   fParams.SetSignal(fSingleI,"SingleI");
  }
  if (ndoms1>=0)
  {
   fSingle1I=ndoms1;
   fParams.AddNamedSlot("Single1I");
   fParams.SetSignal(fSingle1I,"Single1I");
  }
 }
}

void IceRecoBase::SetMinAhits(Int_t nmin,TString s)
{
 
 if (s=="A")
 {
  fMinahitsA=nmin;
  fParams.AddNamedSlot("MinahitsA");
  fParams.SetSignal(fMinahitsA,"MinahitsA");
 }
 if (s=="I")
 {
  fMinahitsI=nmin;
  fParams.AddNamedSlot("MinahitsI");
  fParams.SetSignal(fMinahitsI,"MinahitsI");
 }
 if (s=="IC")
 {
  fMinahitsIC=nmin;
  fParams.AddNamedSlot("MinahitsIC");
  fParams.SetSignal(fMinahitsIC,"MinahitsIC");
 }
 if (s=="DC")
 {
  fMinahitsDC=nmin;
  fParams.AddNamedSlot("MinahitsDC");
  fParams.SetSignal(fMinahitsDC,"MinahitsDC");
 }
}

void IceRecoBase::SetMinAmods(Int_t nmin,TString s)
{
 
 if (s=="A")
 {
  fMinamodsA=nmin;
  fParams.AddNamedSlot("MinamodsA");
  fParams.SetSignal(fMinamodsA,"MinamodsA");
 }
 if (s=="I")
 {
  fMinamodsI=nmin;
  fParams.AddNamedSlot("MinamodsI");
  fParams.SetSignal(fMinamodsI,"MinamodsI");
 }
 if (s=="IC")
 {
  fMinamodsIC=nmin;
  fParams.AddNamedSlot("MinamodsIC");
  fParams.SetSignal(fMinamodsIC,"MinamodsIC");
 }
 if (s=="DC")
 {
  fMinamodsDC=nmin;
  fParams.AddNamedSlot("MinamodsDC");
  fParams.SetSignal(fMinamodsDC,"MinamodsDC");
 }
}

void IceRecoBase::SetSLChitUsage(Int_t flag,TString s)
{
 
 if (s=="I")
 {
  fSlcI=flag;
  fParams.AddNamedSlot("SlcI");
  fParams.SetSignal(fSlcI,"SlcI");
 }
 if (s=="IC")
 {
  fSlcIC=flag;
  fParams.AddNamedSlot("SlcIC");
  fParams.SetSignal(fSlcIC,"SlcIC");
 }
 if (s=="DC")
 {
  fSlcDC=flag;
  fParams.AddNamedSlot("SlcDC");
  fParams.SetSignal(fSlcDC,"SlcDC");
 }
}

void IceRecoBase::SetFlipAngles(Float_t thetatrk,Float_t thetahits)
{
 
 fThetatrk=thetatrk;
 fThetahits=thetahits;
 
 fParams.AddNamedSlot("Thetatrk");
 fParams.AddNamedSlot("Thetahits");
 fParams.SetSignal(fThetatrk,"Thetatrk");
 fParams.SetSignal(fThetahits,"Thetahits");
}

void IceRecoBase::SetScatteringLength(Float_t lambda,TString s)
{
 
 if (s=="A")
 {
  fLambdaA=lambda;
  fParams.AddNamedSlot("LambdaA");
  fParams.SetSignal(fLambdaA,"LambdaA");
 }
 if (s=="UD")
 {
  fLambdaUD=lambda;
  fParams.AddNamedSlot("LambdaUD");
  fParams.SetSignal(fLambdaUD,"LambdaUD");
 }
 if (s=="DL")
 {
  fLambdaDL=lambda;
  fParams.AddNamedSlot("LambdaDL");
  fParams.SetSignal(fLambdaDL,"LambdaDL");
 }
 if (s=="LD")
 {
  fLambdaLD=lambda;
  fParams.AddNamedSlot("LambdaLD");
  fParams.SetSignal(fLambdaLD,"LambdaLD");
 }
}

void IceRecoBase::SetAbsorptionLength(Float_t lambda,TString s)
{
 
 if (s=="A")
 {
  fLabsA=lambda;
  fParams.AddNamedSlot("LabsA");
  fParams.SetSignal(fLabsA,"LabsA");
 }
 if (s=="UD")
 {
  fLabsUD=lambda;
  fParams.AddNamedSlot("LabsUD");
  fParams.SetSignal(fLabsUD,"LabsUD");
 }
 if (s=="DL")
 {
  fLabsDL=lambda;
  fParams.AddNamedSlot("LabsDL");
  fParams.SetSignal(fLabsDL,"LabsDL");
 }
 if (s=="LD")
 {
  fLabsLD=lambda;
  fParams.AddNamedSlot("LabsLD");
  fParams.SetSignal(fLabsLD,"LabsLD");
 }
}

void IceRecoBase::SetTimeJitter(Float_t sigma,TString s)
{
 
 if (s=="A")
 {
  fTsigmaA=sigma;
  fParams.AddNamedSlot("TsigmaA");
  fParams.SetSignal(fTsigmaA,"TsigmaA");
 }
 if (s=="IC")
 {
  fTsigmaIC=sigma;
  fParams.AddNamedSlot("TsigmaIC");
  fParams.SetSignal(fTsigmaIC,"TsigmaIC");
 }
 if (s=="DC")
 {
  fTsigmaDC=sigma;
  fParams.AddNamedSlot("TsigmaDC");
  fParams.SetSignal(fTsigmaDC,"TsigmaDC");
 }
}

void IceRecoBase::SetVgroupUsage(Int_t flag,TString s)
{
 
 if (s=="A")
 {
  fVgroupA=flag;
  fParams.AddNamedSlot("VgroupA");
  fParams.SetSignal(fVgroupA,"VgroupA");
 }
 if (s=="IC")
 {
  fVgroupIC=flag;
  fParams.AddNamedSlot("VgroupIC");
  fParams.SetSignal(fVgroupIC,"VgroupIC");
 }
 if (s=="DC")
 {
  fVgroupDC=flag;
  fParams.AddNamedSlot("VgroupDC");
  fParams.SetSignal(fVgroupDC,"VgroupDC");
 }
 if (s=="I")
 {
  fVgroupI=flag;
  fParams.AddNamedSlot("VgroupI");
  fParams.SetSignal(fVgroupI,"VgroupI");
 }
}

void IceRecoBase::SetTrackName(TString s)
{
 
 fTrackname=s;
 
 fParams.AddNamedSlot(fTrackname);
 fParams.SetSignal(1,fTrackname);
}

void IceRecoBase::SetCharge(Float_t charge)
{
 
 fCharge=charge;
}

void IceRecoBase::UseTracks(TString classname,Int_t n)
{
 
 if (!fUseNames)
 {
  fUseNames=new TObjArray();
  fUseNames->SetOwner();
 }
 
 if (!fUseNtk) fUseNtk=new TArrayI();
 
 // Check if this classname has already been specified before 
 TString s;
 Int_t nen=fUseNames->GetEntries();
 for (Int_t i=0; i<nen; i++)
 {
  TObjString* sx=(TObjString*)fUseNames->At(i);
  if (!sx) continue;
  s=sx->GetString();
  if (s==classname) return;
 }
 
 // New classname to be added into the storage
 if (nen >= fUseNames->GetSize()) fUseNames->Expand(nen+1);
 if (nen >= fUseNtk->GetSize()) fUseNtk->Set(nen+1);
 
 TObjString* name=new TObjString();
 name->SetString(classname);
 fUseNames->Add(name);
 fUseNtk->AddAt(n,nen);
 
 fParams.AddNamedSlot("UseTracks");
 fParams.SetSignal(1,"UseTracks");
}

void IceRecoBase::Exec(Option_t* opt)
{
// Template for the implementation of a reconstruction processor. //
~~~
**/
 
 // Obtain a pointer to the parent NcJob of this reconstruction task
 TString name=opt;
 NcJob* parent=(NcJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 
 if (!parent) return;
 
 // Obtain a pointer to the IceCube event data structure
 fEvt=(IceEvent*)parent->GetObject("IceEvent");
 if (!fEvt) return;
 
 // Only process accepted events
 NcDevice* seldev=(NcDevice*)fEvt->GetDevice("NcEventSelector");
 if (seldev)
 {
  if (seldev->GetSignal("Select") < 0.1) return;
 }
 
 // Provide a name for the fParams device in the event
 if (!fUseNames) // Reconstruction procedure on complete event
 {
  fParams.SetNameTitle("IceRecoBase","IceRecoBase complete event reco parameters");
 }
 else // Reconstruction procedure on track associated data
 {
  fParams.SetNameTitle("IceRecoBase4Track","IceRecoBase track based reco parameters");
 }
 
 // Add the fParams device to the IceEvent structure
 fEvt->AddDevice(fParams);
 
 // Printout information on used tracks (if any) at first startup of the processor task
 if (fUseNames && fFirst)
 {
  Int_t nclasses=0;
  if (fUseNames) nclasses=fUseNames->GetEntries(); // Number of first guess classes to be processed
  Int_t ntkmax=0; // Max. number of tracks for a certain class
  TObjString* strx=0;
  TString str;
  cout << " *IceRecoBase* First guess selections to be processed (-1=all)." << endl;
  for (Int_t i=0; i<nclasses; i++)
  {
   strx=(TObjString*)fUseNames->At(i);
   if (!strx) continue;
   str=strx->GetString();
   ntkmax=fUseNtk->At(i);
   cout << " Maximally " << ntkmax << " track(s) per event for procedure : " << str.Data() << endl;
  }
  cout << endl;
 
  fFirst=0;
 }

 // Enter here your code for the actual reconstruction process. //
 
}

void IceRecoBase::FlipTrack(NcTrack* t) const
{
 
 if (!t || fThetatrk<0 || fThetahits>180) return;
 
 Nc3Vector p;
 Float_t theta=0; 
 p=t->Get3Momentum();
 theta=p.GetX(2,"sph","deg");
 
 if (theta>=fThetatrk) return;
 
 // Check the HitPath of the time ordered hits associated to this track
 TObjArray* hits=t->GetSignals("IceGOM",1);
 if (!hits) return;
 NcDevice sorter;
 TObjArray* sorted=sorter.SortHits("LE",1,hits,7);
 if (!sorted) return;
 
 // Make sure that the starting hit is an HLC hit
 for (Int_t i=0; i<sorted->GetEntries(); i++)
 {
  NcSignal* sx=(NcSignal*)sorted->At(i);
  if (!sx) continue;
 
  // For Amanda the concept of HLC and SLC hits didn't exist
  NcDevice* omx=sx->GetDevice();
  if (!omx) continue;
  if (omx->InheritsFrom("IceAOM")) continue;
 
  if (sx->GetSignal("SLC")<0.5) break;
  sorted->Remove(sx);
 }
 sorted->Compress();
 
 Nc3Vector hitpath;
 hitpath=fEvt->GetHitPath(sorted,1);
 theta=hitpath.GetX(2,"sph","deg");
 
 if (theta<=fThetahits) return;
 
 // Flip the track direction
 p*=-1.;
 t->Set3Momentum(p);
 NcSignal* fitstats=(NcSignal*)t->GetFitDetails();
 if (fitstats)
 {
  fitstats->AddNamedSlot("TrackFlip");
  fitstats->SetSignal(1,"TrackFlip");
 }
}