IceLinefit.cxx

Go to the documentation of this file.

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

IceLinefit::IceLinefit(const char* name,const char* title) : IceRecoBase(name,title)
{
 
 SetCleaned(0,"A");
 SetCleaned(0,"I");
 SetCleaned(0,"IC");
 SetCleaned(0,"DC");
 
 SetMaxMod(999999,"A");
 SetMaxMod(999999,"I");
 SetMaxMod(999999,"IC");
 SetMaxMod(999999,"DC");
 
 SetMinMod(0,"A");
 SetMinMod(0,"I");
 SetMinMod(0,"IC");
 SetMinMod(0,"DC");
 
 SetMaxHits(0,"A");
 SetMaxHits(0,"I");
 SetMaxHits(-1,"IC");
 SetMaxHits(-1,"DC");
 
 SetMinAhits(0,"A");
 SetMinAhits(0,"I");
 SetMinAhits(0,"IC");
 SetMinAhits(0,"DC");
 
 SetMinAmods(0,"A");
 SetMinAmods(0,"I");
 SetMinAmods(0,"IC");
 SetMinAmods(0,"DC");
 
 SetSLChitUsage(0,"I");
 SetSLChitUsage(0,"IC");
 SetSLChitUsage(0,"DC");
 
 SetFlipAngles(-999,999);
}

IceLinefit::~IceLinefit()
{
}

void IceLinefit::Exec(Option_t* opt)
{
 
 // 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) // Linefit procedure on complete event
 {
  fParams.SetNameTitle("IceLinefit","IceLinefit complete event reco parameters");
 }
 else // Linefit procedure on track associated hits
 {
  fParams.SetNameTitle("IceLinefit4Track","IceLinefit 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 << " *IceLinefit* 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;
 }
 
 // Perform linefit reconstruction for the various (associated) hits
 
 if (!fUseNames) // Linefit procedure on complete event
 {
  Amanda();
  InIce();
  IceCube();
  DeepCore();
 }
 else // Linefit procedure on track associated hits
 {
  Tracks();
 }
}

Int_t IceLinefit::Amanda()
{
 
 if (fMaxhitsA<0) return 0;
 
 // Fetch all fired Amanda OMs for this event
 TObjArray* aoms=fEvt->GetDevices("IceAOM");
 
 // Perform the reconstruction
 NcTrack* trk=Reconstruct(aoms,fCleanA,fMinmodA,fMaxmodA,fMaxhitsA,fMinahitsA,fMinamodsA,1);
 
 if (!trk) return 0;
 
 // Give the reconstructed track the proper name and title
 TString name=fTrackname;
 if (name=="") name=ClassName();
 name+="A";
 TString title=ClassName();
 title+=" Amanda track";
 trk->SetNameTitle(name.Data(),title.Data());
 
 return 1;
}

Int_t IceLinefit::InIce()
{
 
 if (fMaxhitsI<0) return 0;
 
 // Fetch all fired InIce DOMs for this event
 TObjArray* doms=fEvt->GetDevices("IceIDOM");
 
 // Perform the reconstruction
 NcTrack* trk=Reconstruct(doms,fCleanI,fMinmodI,fMaxmodI,fMaxhitsI,fMinahitsI,fMinamodsI,fSlcI);
 
 if (!trk) return 0;
 
 // Give the reconstructed track the proper name and title
 TString name=fTrackname;
 if (name=="") name=ClassName();
 name+="I";
 TString title=ClassName();
 title+=" InIce track";
 trk->SetNameTitle(name.Data(),title.Data());
 
 return 1;
}

Int_t IceLinefit::IceCube()
{
 
 if (fMaxhitsIC<0) return 0;
 
 // Fetch all fired standard IceCube InIce DOMs for this event
 TObjArray* doms=fEvt->GetDevices("IceICDOM");
 
 // Perform the reconstruction
 NcTrack* trk=Reconstruct(doms,fCleanIC,fMinmodIC,fMaxmodIC,fMaxhitsIC,fMinahitsIC,fMinamodsIC,fSlcIC);
 
 if (!trk) return 0;
 
 // Give the reconstructed track the proper name and title
 TString name=fTrackname;
 if (name=="") name=ClassName();
 name+="IC";
 TString title=ClassName();
 title+=" Standard IceCube InIce track";
 trk->SetNameTitle(name.Data(),title.Data());
 
 return 1;
}

Int_t IceLinefit::DeepCore()
{
 
 if (fMaxhitsDC<0) return 0;
 
 // Fetch all fired DeepCore DOMs for this event
 TObjArray* doms=fEvt->GetDevices("IceDCDOM");
 
 // Perform the reconstruction
 NcTrack* trk=Reconstruct(doms,fCleanDC,fMinmodDC,fMaxmodDC,fMaxhitsDC,fMinahitsDC,fMinamodsDC,fSlcDC);
 
 if (!trk) return 0;
 
 // Give the reconstructed track the proper name and title
 TString name=fTrackname;
 if (name=="") name=ClassName();
 name+="DC";
 TString title=ClassName();
 title+=" DeepCore track";
 trk->SetNameTitle(name.Data(),title.Data());
 
 return 1;
}

Int_t IceLinefit::Tracks()
{
 
 if (!fUseNames) return 0;
 
 Int_t 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;
 
 // Track by track processing of the selected first guess classes
 Int_t ntk=0;
 NcTrack* track=0;
 TObjArray* hits=0;
 Int_t trkflag=0;
 TString trackname;
 TString newname;
 TObjArray mytracks; // Temp. storage for the extracted tracks per class
 for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
 {
  strx=(TObjString*)fUseNames->At(iclass);
  if (!strx) continue;
  str=strx->GetString();
  ntkmax=fUseNtk->At(iclass);
  TObjArray* tracks=fEvt->GetTracks(str);
  ntk=0;
  if (tracks) ntk=tracks->GetEntries();
  if (ntkmax>0 && ntk>ntkmax) ntk=ntkmax;
 
  // Store track pointers in private array to prevent overwriting of TObjArray
  mytracks.Clear();
  for (Int_t i=0; i<ntk; i++)
  {
   mytracks.Add(tracks->At(i));
  }
 
  for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
  {
   track=(NcTrack*)mytracks.At(jtk);
   if (!track) continue;
 
   // The name of the first guess track without the intial "Ice"
   trackname=track->GetName();
   trackname.ReplaceAll("Ice","");
 
   // Retrieval of the associated hits to be used in the reconstruction procedure
   hits=track->GetSignals("IceGOM",2);
 
   // Perform the reconstruction for this track
   NcTrack* trk=Reconstruct(hits,0,0,0,0,fMinahitsI,fMinamodsI,1);
 
   if (!trk) continue;
 
   // Give the reconstructed track the proper name and title
   trkflag=1;
   newname=fTrackname;
   if (newname=="") newname=ClassName();
   newname+="4";
   newname+=trackname;
   trk->SetNameTitle(newname.Data(),"Linefit reco for all assoc. hits");
   // Link this newly created track as a hypothesis to the parent first guess track
   track->SetHypCopy(0);
   track->AddTrackHypothesis(*trk);
  } // End of loop over the tracks of the current first guess class
 } // End of loop over first guess classes
 
 return trkflag;
}

NcTrack* IceLinefit::Reconstruct(TObjArray* arr,Int_t cln,Int_t minmod,Int_t maxmod,Int_t maxhits,Int_t minahits,Int_t minamods,Int_t slc)
{
 
 if (maxhits<0) return 0;
 
 if (!arr) return 0;
 Int_t narr=arr->GetEntries();
 if (narr<=0) return 0;
 
 if (!fUseNames) // Complete event reco via provided array of (D)OMs
 {
  // Check for the minimum and/or maximum number of good fired DOMs
  Int_t ngood=0;
  for (Int_t iarr=0; iarr<narr; iarr++)
  {
   IceGOM* omx=(IceGOM*)arr->At(iarr);
   if (!omx) continue;
   if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
   ngood++;
  } 
  if (ngood<minmod || ngood>maxmod) return 0;
 }
 
 const Float_t c=0.299792; // Light speed in vacuum in meters per ns
 
 IceGOM* omx=0;
 NcSignal* sx=0;
 Nc3Vector rom,sumr;
 Nc3Vector rt,sumrt;
 Float_t thit;
 Float_t sumt=0,sumt2=0;
 TObjArray hits;
 TObjArray* ordered=0;
 Int_t nhitloop=0;
 Int_t nh;
 
 // Loop over all DOMs and/or hits to determine the linefit parameters.
 // Also all the used hits are recorded for association with newly created reco track.
 for (Int_t iarr=0; iarr<narr; iarr++)
 {
  if (!fUseNames) // Complete event reco via provided array of (D)OMs
  {
   omx=(IceGOM*)arr->At(iarr);
   if (!omx) continue;
   if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
   // Use the specified good hits of this (D)OM
   ordered=0;
   if (maxhits>0 && omx->GetNhits()>maxhits)
   {
    if (omx->InheritsFrom("IceAOM")) ordered=omx->SortHits("LE",1,0,7);
    if (omx->InheritsFrom("IceDOM")) ordered=omx->SortHits("ADC",-1,0,7);
   }
   nhitloop=omx->GetNhits();
  }
  else // Track based reco via provided array of associated hits
  {
   nhitloop=1;
  }
 
  // Loop over all the hits of the current (D)OM.
  // In case an associated hit array was provided, only a single pass is made.
  nh=0;
  for (Int_t ih=1; ih<=nhitloop; ih++)
  {
   if (ordered) // Only take the first "maxhits" ordered hits for this (D)OM 
   {
    if (nh>=maxhits) break;
    sx=(NcSignal*)ordered->At(ih-1);
   }
   else // Take all available hits
   {
    if (!fUseNames) // Obtain the hits from the provided (D)OMs
    {
     sx=omx->GetHit(ih);
    }
    else // Obtain the hits directly from the provided array
    {
     sx=(NcSignal*)arr->At(iarr);
     omx=(IceGOM*)sx->GetDevice();
    }
   }
 
   if (!omx || !sx) continue;
 
   if (cln && (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT"))) continue;
   if (!slc && sx->GetSignal("SLC")) continue; 
 
   thit=sx->GetSignal("LE",7);
   rom=(Nc3Vector)omx->GetPosition();
   rt=rom*thit;
   sumr+=rom;
   sumrt+=rt;
   sumt+=thit;
   sumt2+=thit*thit;
 
   // Record this hit for association with the track
   hits.Add(sx);
   nh++;
  }
 }
 
 Int_t nused=hits.GetEntries();
 if (!nused || nused<minahits) return 0;
 
 Int_t nasmods=fEvt->GetNdevices("IceGOM",&hits);
 if (nasmods<minamods) return 0;
 
 sumr/=float(nused);
 sumrt/=float(nused);
 sumt/=float(nused);
 sumt2/=float(nused);
 
 Nc3Vector v;
 Nc3Vector temp;
 temp=sumr*sumt;
 v=sumrt-temp;
 Float_t dum=sumt2-(sumt*sumt);
 if (dum) v/=dum;
 
 Float_t beta=v.GetNorm()/c;
 NcSignal fitstats;
 fitstats.SetNameTitle("Fitstats","Fit stats for IceLinefit");
 fitstats.AddNamedSlot("Beta");
 fitstats.SetSignal(beta,"Beta");
 
 Nc3Vector r;
 temp=v*sumt;
 r=sumr-temp;
 
 NcTrack t; 
 t.SetCharge(fCharge);
 fEvt->AddTrack(t);
 NcTrack* trk=fEvt->GetTrack(fEvt->GetNtracks());
 if (!trk) return 0;
 
 trk->SetId(fEvt->GetNtracks(1)+1);
 
 Nc3Vector p;
 Float_t vec[3];
 v.GetVector(vec,"sph");
 vec[0]=1;
 p.SetVector(vec,"sph");
 
 NcPosition r0;
 r0.SetPosition(r);
 r0.SetTimestamp((NcTimestamp&)*fEvt);
 NcTimestamp* t0=r0.GetTimestamp();
 t0->Add(0,0,(int)sumt);
 
 trk->Set3Momentum(p);
 trk->SetReferencePoint(r0);
 trk->SetTimestamp(*t0);
 trk->SetFitDetails(fitstats);
 
 // Link the used hits to the track (and vice versa)
 for (Int_t i=0; i<nused; i++)
 {
  sx=(NcSignal*)hits.At(i);
  if (sx) sx->AddTrack(*trk);
 }
 
 // Check whether the track direction should be flipped
 FlipTrack(trk);
 
 return trk; 
}