Loaders removed from the reconstruction code (C.Cheshkov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPIDv0.cxx
index ff914db..e281da8 100644 (file)
 
 /* $Id$ */
 
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.15  2007/03/06 06:57:05  kharlov
+ * DP:calculation of distance to CPV done in TSM
+ *
+ * Revision 1.14  2006/09/07 18:31:08  kharlov
+ * Effective c++ corrections (T.Pocheptsov)
+ *
+ * Revision 1.13  2005/05/28 14:19:04  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
 // Implementation version v0 of the PHOS particle identifier 
 // Particle identification based on the 
 //            Completely redesined by Dmitri Peressounko, March 2001
 
 // --- ROOT system ---
-#include "TROOT.h"
 #include "TTree.h"
-#include "TFile.h"
 #include "TF2.h"
 #include "TFormula.h"
 #include "TCanvas.h"
-#include "TFolder.h"
-#include "TSystem.h"
+#include "TClonesArray.h"
+
 #include "TBenchmark.h"
 // --- Standard library ---
 
 // --- AliRoot header files ---
-
-#include "AliRun.h"
-#include "AliGenerator.h"
-#include "AliPHOS.h"
+#include "AliLog.h"
 #include "AliPHOSPIDv0.h"
 #include "AliPHOSEmcRecPoint.h"
-#include "AliPHOSClusterizerv1.h"
 #include "AliPHOSTrackSegment.h"
-#include "AliPHOSTrackSegmentMakerv1.h"
 #include "AliPHOSRecParticle.h"
 #include "AliPHOSGeometry.h"
-#include "AliPHOSLoader.h"
 
 ClassImp( AliPHOSPIDv0) 
 
 //____________________________________________________________________________
-AliPHOSPIDv0::AliPHOSPIDv0():AliPHOSPID()
+AliPHOSPIDv0::AliPHOSPIDv0():
+  fIDOptions("dis time"),
+  fClusterizer(0),
+  fTSMaker(0),
+  fFormula(0),
+  fDispersion(0.f),
+  fCpvEmcDistance(0.f),
+  fTimeGate(2.e-9f)
 { 
   // default ctor
-  fFormula           = 0 ;
-  fDispersion        = 0. ; 
-  fCpvEmcDistance    = 0 ; 
-  fTimeGate          = 2.e-9 ;
-  fEventFolderName    = "" ; 
-  fTrackSegmentsTitle= "" ; 
-  fRecPointsTitle    = "" ; 
-  fRecParticlesTitle = "" ; 
-  fIDOptions         = "dis time" ; 
-  fRecParticlesInRun = 0 ;
-  fClusterizer = 0;
-  fTSMaker = 0;
 }
 
 //____________________________________________________________________________
-AliPHOSPIDv0::AliPHOSPIDv0(const char * evFolderName,const char * name) : AliPHOSPID(evFolderName, name)
+AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) : 
+  AliPHOSPID(geom),
+  fIDOptions("dis time"),
+  fClusterizer(0),
+  fTSMaker(0),
+  fFormula(new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)")),
+  fDispersion(2.f),
+  fCpvEmcDistance(3.f),
+  fTimeGate(2.e-9f)
 { 
-  //ctor with the indication on where to look for the track segments
-
-  fFormula        = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)") ;   
-  fDispersion     = 2.0 ; 
-  fCpvEmcDistance = 3.0 ;
-  fTimeGate          = 2.e-9 ;
-  fEventFolderName     = GetTitle() ; 
-  fTrackSegmentsTitle = GetName();
-  fRecPointsTitle     = GetName();
-  fRecParticlesTitle  = GetName();
-  fIDOptions          = "dis time" ;
-    
-  fRecParticlesInRun = 0 ; 
-
-  Init() ;
-
+  //ctor
 }
 
 //____________________________________________________________________________
-AliPHOSPIDv0::~AliPHOSPIDv0()
-{ 
+AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & rhs) :
+  AliPHOSPID(rhs),
+  fIDOptions(rhs.fIDOptions),
+  fClusterizer(rhs.fClusterizer),
+  fTSMaker(rhs.fTSMaker),
+  fFormula(rhs.fFormula),
+  fDispersion(rhs.fDispersion),
+  fCpvEmcDistance(rhs.fCpvEmcDistance),
+  fTimeGate(rhs.fTimeGate)
+{
+  //Copy ctor, the same as compiler-generated, possibly wrong if
+  //someone implements dtor correctly.
+}
+  
+//____________________________________________________________________________
+AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs)
+{
+  //Copy-assignment, emulates compiler generated, possibly wrong.
+  AliPHOSPID::operator = (rhs);
+  fIDOptions = rhs.fIDOptions;
+  fClusterizer = rhs.fClusterizer;
+  fTSMaker = rhs.fTSMaker;
+  fFormula = rhs.fFormula;
+  fDispersion = rhs.fDispersion;
+  fCpvEmcDistance = rhs.fCpvEmcDistance;
+  fTimeGate = rhs.fTimeGate;
+
+  return *this;
 }
 
 //____________________________________________________________________________
-Float_t  AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
+AliPHOSPIDv0::~AliPHOSPIDv0()
 {
-  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
-  const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ; 
-  TVector3 vecEmc ;
-  TVector3 vecCpv ;
-  
-  emc->GetLocalPosition(vecEmc) ;
-  cpv->GetLocalPosition(vecCpv) ; 
-  if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ 
-    
-    // Correct to difference in CPV and EMC position due to different distance to center.
-    // we assume, that particle moves from center
-    Float_t dCPV = geom->GetIPtoOuterCoverDistance();
-    Float_t dEMC = geom->GetIPtoCrystalSurface() ;
-    dEMC         = dEMC / dCPV ;
-    vecCpv = dEMC * vecCpv  - vecEmc ; 
-    if (Axis == "X") return vecCpv.X();
-    if (Axis == "Y") return vecCpv.Y();
-    if (Axis == "Z") return vecCpv.Z();
-    if (Axis == "R") return vecCpv.Mag();
-  } 
-  return 100000000 ;
+  //Empty dtor, fFormula leaks 
 }
 
+//DP
+////____________________________________________________________________________
+//Float_t  AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
+//{
+//  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
+// 
+//  const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ; 
+//  TVector3 vecEmc ;
+//  TVector3 vecCpv ;
+//  
+//  emc->GetLocalPosition(vecEmc) ;
+//  cpv->GetLocalPosition(vecCpv) ; 
+//  if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ 
+//    
+//    // Correct to difference in CPV and EMC position due to different distance to center.
+//    // we assume, that particle moves from center
+//    Float_t dCPV = geom->GetIPtoOuterCoverDistance();
+//    Float_t dEMC = geom->GetIPtoCrystalSurface() ;
+//    dEMC         = dEMC / dCPV ;
+//    vecCpv = dEMC * vecCpv  - vecEmc ; 
+//    if (Axis == "X") return vecCpv.X();
+//    if (Axis == "Y") return vecCpv.Y();
+//    if (Axis == "Z") return vecCpv.Z();
+//    if (Axis == "R") return vecCpv.Mag();
+//  } 
+// 
+//  return 100000000 ;
+//}
+
 //____________________________________________________________________________
-void  AliPHOSPIDv0::Exec(Option_t * option) 
+void  AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option) 
 {
   //Steering method
-  
-  if( strcmp(GetName(), "")== 0 ) 
-    Init() ;
-  
+
   if(strstr(option,"tim"))
     gBenchmark->Start("PHOSPID");
   
@@ -173,103 +196,31 @@ void  AliPHOSPIDv0::Exec(Option_t * option)
     return ; 
   }
 
-  AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
-  if(runget == 0x0) 
-   {
-     Error("Exec","Can not find run getter in event folder \"%s\"",GetTitle());
-     return;
-   }
-  
-  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
-  if ( gime == 0 ) 
-   {
-     Error("Exec","Could not obtain the Loader object !"); 
-     return ;
-   } 
-
-  if(gime->BranchExists("RecParticles") )
-    return ;
-
-  
-  Int_t nevents = runget->GetNumberOfEvents() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
+  AliInfo(Form("%d emc clusters, %d track segments", 
+              fEMCRecPoints->GetEntriesFast(), 
+              fTrackSegments->GetEntriesFast())) ;
 
-  Int_t ievent ;
-  
-  for(ievent = 0; ievent < nevents; ievent++){
-    runget->GetEvent(ievent);
-    Info("Exec", "event %d %d %d", ievent, gime->EmcRecPoints(), gime->TrackSegments()) ;
-    MakeRecParticles() ;
-    
-    WriteRecParticles();
+  MakeRecParticles() ;
     
-    if(strstr(option,"deb"))
-      PrintRecParticles(option) ;
-
-    //increment the total number of rec particles per run 
-    fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
+  if(strstr(option,"deb"))
+    PrintRecParticles(option) ;
 
-  }
-  
   if(strstr(option,"tim")){
     gBenchmark->Stop("PHOSPID");
-    Info("Exec", "took %f seconds for PID %f seconds per event", 
-        gBenchmark->GetCpuTime("PHOSPID"), gBenchmark->GetCpuTime("PHOSPID")/nevents) ; 
+    AliInfo(Form("took %f seconds for PID", 
+                gBenchmark->GetCpuTime("PHOSPID"))); 
   }
   
 }
-//____________________________________________________________________________
-void AliPHOSPIDv0::Init()
-{
-  // Make all memory allocations that are not possible in default constructor
-  // Add the PID task to the list of PHOS tasks
-  
-  AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
-  if(runget == 0x0) 
-   {
-     Error("Exec","Can not find run getter in event folder \"%s\"",GetTitle());
-     return;
-   }
-  
-  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
-  if ( gime == 0 ) 
-   {
-     Error("Exec","Could not obtain the Loader object !"); 
-     return ;
-   } 
-   
-  gime->PostPID(this);
-  gime->LoadRecParticles("UPDATE");
-  
-}
 
 //____________________________________________________________________________
-void  AliPHOSPIDv0::MakeRecParticles(){
-
-  TString taskName(GetName()) ; 
-  taskName.Remove(taskName.Index(Version())-1) ;
+void  AliPHOSPIDv0::MakeRecParticles()
+{
+  // Reconstructs the particles from the tracksegments
 
-  AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
-  if(runget == 0x0) 
-   {
-     Error("Exec","Can not find run getter in event folder \"%s\"",GetTitle());
-     return;
-   }
-  
-  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
-  if ( gime == 0 ) 
-   {
-     Error("Exec","Could not obtain the Loader object !"); 
-     return ;
-   } 
-
-  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
-  TClonesArray * trackSegments = gime->TrackSegments() ; 
-  TClonesArray * recParticles  = gime->RecParticles() ; 
-
-  recParticles->Clear();
+  fRecParticles->Clear();
   
-  TIter next(trackSegments) ; 
+  TIter next(fTrackSegments) ; 
   AliPHOSTrackSegment * ts ; 
   Int_t index = 0 ; 
   AliPHOSRecParticle * rp ; 
@@ -280,18 +231,18 @@ void  AliPHOSPIDv0::MakeRecParticles(){
   
   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
     
-    new( (*recParticles)[index] ) AliPHOSRecParticle() ;
-    rp = (AliPHOSRecParticle *)recParticles->At(index) ; 
+    new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
+    rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; 
     rp->SetTrackSegment(index) ;
     rp->SetIndexInList(index) ;
     
     AliPHOSEmcRecPoint * emc = 0 ;
     if(ts->GetEmcIndex()>=0)
-      emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
+      emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
     
     AliPHOSRecPoint    * cpv = 0 ;
     if(ts->GetCpvIndex()>=0)
-      cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
+      cpv = (AliPHOSRecPoint *)   fCPVRecPoints->At(ts->GetCpvIndex()) ;
     
     //set momentum and energy first
     Float_t    e = emc->GetEnergy() ;
@@ -323,7 +274,7 @@ void  AliPHOSPIDv0::MakeRecParticles(){
     // Looking at the CPV detector
     Int_t cpvdetector= 0 ;  //1 hit and 0 no hit     
     if(cpv)
-      if(GetDistance(emc, cpv,  "R") < fCpvEmcDistance) 
+      if(ts->GetCpvDistance("R") < fCpvEmcDistance) 
        cpvdetector = 1 ;  
     
     Int_t type = showerprofile + 2 * slow  + 4 * cpvdetector ;
@@ -340,31 +291,24 @@ void  AliPHOSPIDv0::MakeRecParticles(){
 }
 
 //____________________________________________________________________________
-void  AliPHOSPIDv0:: Print() const
+void  AliPHOSPIDv0:: Print(const Option_t *) const
 {
   // Print the parameters used for the particle type identification
   TString message ; 
   message  = "=============== AliPHOSPIDv0 ================\n" ;
   message += "Making PID\n" ;
-  message += "    Headers file:               %s\n" ; 
-  message += "    RecPoints branch title:     %s\n" ;
-  message += "    TrackSegments Branch title: %s\n" ; 
-  message += "    RecParticles Branch title   %s\n" ;  
   message += "with parameters:\n"  ;
   message += "    Maximal EMC - CPV  distance (cm) %f\n" ;
-  Info("Print", message.Data(),  
-       GetTitle(), 
-       fRecPointsTitle.Data(), 
-       fTrackSegmentsTitle.Data(), 
-       fRecParticlesTitle.Data(), 
-       fCpvEmcDistance );
+  AliInfo(Form( message.Data(),  
+       fCpvEmcDistance ));
 
   if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
-    Info("Print", "                    dispersion cut %f",  fDispersion ) ;
+    AliInfo(Form("                    dispersion cut %f",  fDispersion )) ;
   if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
-    Info("Print", "             Eliptic cuts function: %s",  fFormula->GetTitle() ) ;
+    AliInfo(Form("             Eliptic cuts function: %s",  
+                fFormula->GetTitle() )) ;
   if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
-    Info("Print",  "             Time Gate used: %f",  fTimeGate) ;
+    AliInfo(Form("             Time Gate used: %f",  fTimeGate)) ;
 }
 
 //____________________________________________________________________________
@@ -376,52 +320,6 @@ void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
     delete fFormula; 
   fFormula = new TFormula("Lambda Cut",formula) ;
 }
-//____________________________________________________________________________
-void  AliPHOSPIDv0::WriteRecParticles()
-{
-  AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
-  if(runget == 0x0) 
-   {
-     Error("Exec","Can not find run getter in event folder \"%s\"",GetTitle());
-     return;
-   }
-  
-  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
-  if ( gime == 0 ) 
-   {
-     Error("Exec","Could not obtain the Loader object !"); 
-     return ;
-   } 
-
-  TClonesArray * recParticles = gime->RecParticles() ; 
-  recParticles->Expand(recParticles->GetEntriesFast() ) ;
-
-  TTree * treeR = gime->TreeR();
-  
-  if(!treeR){
-    gime->MakeTree("R");
-    treeR = gime->TreeR() ;
-  }
-  
-  //First rp
-  Int_t bufferSize = 32000 ;    
-  TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
-  rpBranch->SetTitle(fRecParticlesTitle);
-  
-  //second, pid
-  Int_t splitlevel = 0 ; 
-  AliPHOSPIDv0 * pid = this ;
-  TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv0",&pid,bufferSize,splitlevel);
-  pidBranch->SetTitle(fRecParticlesTitle.Data());
-  
-  rpBranch->Fill() ;
-  pidBranch->Fill() ;
-
-  gime->WriteRecParticles("OVERWRITE");
-  gime->WritePID("OVERWRITE");
-
-}
 
 //____________________________________________________________________________
 void  AliPHOSPIDv0::PlotDispersionCuts()const
@@ -493,11 +391,14 @@ TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRec
   dir.SetZ( -dir.Z() ) ;   // why ?  
   dir.SetMag(1.) ;
 
+  // One can not access MC information in the reconstruction!!
+  // PLEASE FIT IT, EITHER BY TAKING 0,0,0 OR ACCESSING THE
+  // VERTEX DIAMOND FROM CDB GRP FOLDER.
   //account correction to the position of IP
-  Float_t xo,yo,zo ; //Coordinates of the origin
-  gAlice->Generator()->GetOrigin(xo,yo,zo) ;
-  TVector3 origin(xo,yo,zo);
-  dir = dir - origin ;
+  //  Float_t xo,yo,zo ; //Coordinates of the origin
+  //  gAlice->Generator()->GetOrigin(xo,yo,zo) ;
+  //  TVector3 origin(xo,yo,zo);
+  //  dir = dir - origin ;
 
   return dir ;  
 }
@@ -506,35 +407,17 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
 {
   // Print table of reconstructed particles
 
-  AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
-  if(runget == 0x0) 
-   {
-     Error("WriteRecParticles","Can not find run getter in event folder \"%s\"",GetTitle());
-     return;
-   }
-  
-  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
-  if ( gime == 0 ) 
-   {
-     Error("WriteRecParticles","Could not obtain the Loader object !"); 
-     return ;
-   } 
-
-  TString taskName(GetName()) ; 
-  taskName.Remove(taskName.Index(Version())-1) ;
-  TClonesArray * recParticles = gime->RecParticles() ; 
-  
   TString message ; 
-  message  = "event %d\n" ; 
-  message += "       found %d RecParticles\n" ; 
-  Info("PrintRecParticles", message.Data(), gAlice->GetEvNumber(), recParticles->GetEntriesFast() ) ;   
+  message = "Found %d RecParticles\n" ; 
+  AliInfo(Form(message.Data(), 
+              fRecParticles->GetEntriesFast() )) ;   
 
   if(strstr(option,"all")) {  // printing found TS
-    Info("PrintRecParticles","  PARTICLE   Index    \n"  ) ; 
+    AliInfo("  PARTICLE   Index"  ) ; 
    
     Int_t index ;
-    for (index = 0 ; index < recParticles->GetEntries() ; index++) {
-      AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;       
+    for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
+      AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;       
       
       Text_t particle[11];
       switch(rp->GetType()) {
@@ -568,10 +451,8 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
       //    Int_t nprimaries;
       //    primaries = rp->GetPrimaries(nprimaries);
       
-      Info("PrintRecParticles", "          %s     %d\n",  particle, rp->GetIndexInList()) ;
+      AliInfo(Form("          %s     %d",  
+                  particle, rp->GetIndexInList())) ;
     }
   }  
 }
-
-
-