]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSReconstructioner.cxx
Transition to NewIO
[u/mrichter/AliRoot.git] / PHOS / AliPHOSReconstructioner.cxx
index acd3ae5db2c77fb9d6de2d99cf91536255bfa1b8..86855cd6084a2f7cb43144dcfbcb40b94095d68a 100644 (file)
 /* $Id$ */
 
 //_________________________________________________________________________
-//  Algorithm class for the reconstruction: clusterizer
-//                                          track segment maker
-//                                          particle identifier   
 //*--
-//*-- Author: Gines Martinez & Yves Schutz (SUBATECH)
-
+//*-- Author: Gines Martinez & Yves Schutz (SUBATECH) 
+//*-- Compleetely redesigned by Dmitri Peressounko (SUBATECH & RRC KI) March 2001
+/////////////////////////////////////////////////////////////////////////////////////
+//  Wrapping class for reconstruction. Allows to produce reconstruction from 
+//  different steps: from previously produced hits,sdigits, etc. Each new reconstruction
+//  flow (e.g. digits, made from them RecPoints, subsequently made TrackSegments, 
+//  subsequently made RecParticles) are distinguished by the title of created branches. One can 
+//  use this title as a comment, see use case below. 
+//  Thanks to getters, one can set 
+//  parameters to reconstruction briks. The full set of parameters is saved in the 
+//  corresponding branch: e.g. parameters of clusterizer are stored in branch 
+//  TreeR::AliPHOSClusterizer with the same title as the branch containing the RecPoints. //  TTree does not support overwriting, therefore one can not produce several 
+//  branches with the same names and titles - use different titles.
+//
+//  Use case: 
+//
+//  root [0] AliPHOSReconstructioner * r = new AliPHOSReconstructioner("galice.root")
+//              //  Set the header file
+//  root [1] r->ExecuteTask() 
+//              //  Make full chain of reconstruction
+//
+//              // One can specify the title for each branch 
+//  root [2] r->SetBranchFileName("RecPoints","RecPoints1") ;
+//      
+//             // One can change parameters of reconstruction algorithms
+//  root [3] r->GetClusterizer()->SetEmcLocalMaxCut(0.02)
+//
+//             // One can specify the starting point of the reconstruction and title of all 
+//             // branches produced in this pass
+//  root [4] r->StartFrom("AliPHOSClusterizer","Local max cut 0.02") 
+//             // means that will use already generated Digits and produce only RecPoints, 
+//             // TS and RecParticles 
+//
+//             // And finally one can call ExecuteTask() with the following options
+//  root [5] r->ExecuteTask("debug all timing")
+//             // deb     - prints the numbers of produced SDigits, Digits etc.
+//             // deb all - prints in addition list of made SDigits, digits etc.
+//             // timing  - prints benchmarking results
+///////////////////////////////////////////////////////////////////////////////////////////////////
 
 // --- ROOT system ---
 
 #include "TClonesArray.h"
+#include "TROOT.h"
+#include "TTree.h"
+#include "TFile.h"
 
 // --- Standard library ---
 
-#include <iomanip.h>
-
 // --- AliRoot header files ---
-
+#include "AliRun.h"
+#include "AliRunLoader.h"
 #include "AliPHOSReconstructioner.h"
-#include "AliPHOSClusterizer.h"
+#include "AliPHOSClusterizerv1.h"
+#include "AliPHOSDigitizer.h"
+#include "AliPHOSSDigitizer.h"
+#include "AliPHOSTrackSegmentMakerv1.h"
+#include "AliPHOSPIDv1.h"
 #include "AliPHOSFastRecParticle.h"
 #include "AliPHOSCpvRecPoint.h"
+#include "AliPHOSLoader.h"
 
 ClassImp(AliPHOSReconstructioner)
 
 //____________________________________________________________________________
-AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer * Clusterizer, 
-                                                AliPHOSTrackSegmentMaker * Tracker,
-                                                AliPHOSPID * Pid)
+  AliPHOSReconstructioner::AliPHOSReconstructioner():TTask("AliPHOSReconstructioner","")
 {
   // ctor
-  
-  fClusterizer        = Clusterizer ;
-  fTrackSegmentMaker  = Tracker ;
-  fPID                = Pid ; 
-  fDebugReconstruction = kFALSE ;
-} 
+  fDigitizer   = 0 ;
+  fClusterizer = 0 ;
+  fTSMaker     = 0 ;
+  fPID         = 0 ; 
+  fSDigitizer  = 0 ;
 
+  fIsInitialized = kFALSE ;
+
+} 
 
 //____________________________________________________________________________
- void AliPHOSReconstructioner::Init(AliPHOSClusterizer * Clusterizer, 
-                                                AliPHOSTrackSegmentMaker * Tracker,
-                                                AliPHOSPID * Pid)
+AliPHOSReconstructioner::AliPHOSReconstructioner(const char* evFoldName,const char * branchName):
+TTask("AliPHOSReconstructioner",evFoldName)
 {
-  // Initialisation
+  // ctor
+  AliRunLoader* rl = AliRunLoader::GetRunLoader(evFoldName);
+  if (rl == 0x0)
+   {
+     Fatal("AliPHOSReconstructioner","Can not get Run Loader from folder %s.",evFoldName);
+   } 
+  if (rl->GetAliRun() == 0x0)
+   {
+     delete gAlice;
+     gAlice = 0x0;
+     rl->LoadgAlice();
+     gAlice = rl->GetAliRun();
+   }
+
+  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
+  if (gime == 0x0)
+   {
+     Error("AliPHOSReconstructioner","Can not get PHOS Loader");
+     return;  
+   }
+  
+  TString galicefn = rl->GetFileName();
+  TString method("AliPHOSReconstructioner::AliPHOSReconstructioner(");
+  method = (((method + evFoldName)+",")+branchName)+"): ";
+  
+  fSDigitsBranch= branchName; 
+  
+  //P.Skowronski remark
+  // Tasks has default fixed names
+  // other tasks can be added, even runtime
+  // with arbitrary name. See AliDataLoader::
+  cout<<"\n\n\n";
+  cout<<method<<"\n\nCreating SDigitizer\n";
+  fSDigitizer  = new AliPHOSSDigitizer(galicefn,GetTitle());
+  Add(fSDigitizer);
+  gime->PostSDigitizer(fSDigitizer);
+
+  fDigitsBranch=branchName ;
+  cout<<"\n\n\n";
+  cout<<method<<"\n\nCreating Digitizer\n";
+  fDigitizer   = new AliPHOSDigitizer(galicefn,GetTitle()) ;
+  Add(fDigitizer) ;
+  gime->PostDigitizer(fDigitizer);
+
+  fRecPointBranch=branchName ; 
+  cout<<"\n\n\n";
+  cout<<method<<"Creating Clusterizer\n";
+  fClusterizer = new AliPHOSClusterizerv1(galicefn,GetTitle());
+  Add(fClusterizer);
+  gime->PostReconstructioner(fClusterizer);
+  
+  fTSBranch=branchName ; 
+  fTSMaker     = new AliPHOSTrackSegmentMakerv1(galicefn,GetTitle());
+  Add(fTSMaker) ;
+  gime->PostTracker(fTSMaker);
 
-  fClusterizer        = Clusterizer ;
-  fTrackSegmentMaker  = Tracker ;
-  fPID                = Pid ; 
-  fDebugReconstruction = kFALSE ;
+  
+  fRecPartBranch=branchName ; 
+  cout<<"\n\n\n";
+  cout<<method<<"Creating PID\n";
+  fPID         = new AliPHOSPIDv1(galicefn,GetTitle());
+  Add(fPID);
+  cout<<"\nFINISHED \n\n"<<method;
+  
+  fIsInitialized = kTRUE ;
 } 
-
 //____________________________________________________________________________
- void AliPHOSReconstructioner::Make(DigitsList * dl, 
-                                   AliPHOSRecPoint::RecPointsList * emccl, 
-                                   AliPHOSRecPoint::RecPointsList * ppsdl, 
-                                   AliPHOSTrackSegment::TrackSegmentsList * trsl, 
-                                   AliPHOSRecParticle::RecParticlesList * rpl)
+void AliPHOSReconstructioner::Exec(Option_t *option)
 {
-  // Launches the Reconstruction process in the sequence: Make the reconstructed poins (clusterize)
-  //                                                      Make the track segments 
-  //                                                      Make the reconstructed particles
-  Int_t index ;   
-  if  (fDebugReconstruction)
-    cout << "\n\nDebugReconstruction>>> " << "Start making reconstructed points (clusterizing!!)" << endl;
-  
-  fClusterizer->MakeClusters(dl, emccl, ppsdl);
+  //chesk, if the names of branches, which should be made conicide with already
+  //existing
+  if(!fIsInitialized)
+    Init() ;
+}
+//____________________________________________________________________________
+ void AliPHOSReconstructioner::Init()
+{
+  // initiliaze Reconstructioner if necessary: we can not do this in default constructor
 
-  if  (fDebugReconstruction){
-    cout << "DebugReconstruction>>> " << "AliPHOSReconstructioner: Digit list entries is " << dl->GetEntries() << endl ;
-    cout << "AliPHOSReconstructioner: Emc  list entries is " << emccl->GetEntries() << endl ;
-    cout << "AliPHOSReconstructioner: Ppsd list entries is " << ppsdl->GetEntries() << endl ;
-  }
+  if(!fIsInitialized){
+    // Initialisation
+
+    fSDigitsBranch="Default" ; 
+    fSDigitizer  = new AliPHOSSDigitizer(GetTitle(),fSDigitsBranch.Data()) ; 
+    Add(fSDigitizer) ;
 
-  // Digit Debuging
-  if  (fDebugReconstruction)     {
-    cout << ">>>>>>>>>>>>>>>>>>>>>> DebugReconstruction  <<<<<<<<<<<<<<<<<<<<<<<<<<"  << endl ;
-    cout << "DebugReconstruction>>> Digit list entries is " <<    dl->GetEntries() << endl ;
-    AliPHOSDigit * digit;
-    Bool_t calorimeter ;
-    Float_t factor;
-    cout << "DebugReconstruction>>>    Vol Id " << 
-      " Ene(MeV, KeV) "              <<                         
-      " Index "                      << 
-      " Nprim "                      << 
-      " Primaries list "             <<  endl;      
-    for (index = 0 ; index < dl->GetEntries() ; index++) {
-      digit = (AliPHOSDigit * )  dl->At(index) ;
-      calorimeter = fClusterizer->IsInEmc(digit);
-      if (calorimeter) factor =1000. ; else factor=1000000.;
-      cout << "DebugReconstruction>>>  " << 
-        setw(8)  <<  digit->GetId() << " "  <<
-       setw(3)  <<  (Int_t) calorimeter <<  
-       setw(10) <<  factor*fClusterizer->Calibrate(digit->GetAmp()) << "  "  <<                   
-       setw(6)  <<  digit->GetIndexInList() << "  "  << 
-       setw(5)  <<  digit->GetNprimary() <<"  ";
-      for (Int_t iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
-       cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << " ";
-      cout << endl;     
-    }
+    fDigitsBranch="Default" ; 
+    fDigitizer   = new AliPHOSDigitizer(GetTitle(),fDigitsBranch.Data());
+    Add(fDigitizer) ;
+
+    fRecPointBranch="Default" ; 
+    fClusterizer = new AliPHOSClusterizerv1(GetTitle(),fRecPointBranch.Data());
+    Add(fClusterizer) ;
+
+    fTSBranch="Default" ; 
+    fTSMaker     = new AliPHOSTrackSegmentMakerv1(GetTitle(),fTSBranch.Data());
+    Add(fTSMaker) ;
+
+
+    fRecPartBranch="Default"; 
+    fPID         = new AliPHOSPIDv1(GetTitle(),fRecPartBranch.Data()) ;
+    Add(fPID) ;
+    
+    fIsInitialized = kTRUE ;
     
   }
+} 
+//____________________________________________________________________________
+AliPHOSReconstructioner::~AliPHOSReconstructioner()
+{
+  // Delete data members if any
+} 
 
-  // Making Clusters
-  if  (fDebugReconstruction)  cout << "DebugReconstruction>>> Start making reconstructed points (clusterizing)" << endl;
+void AliPHOSReconstructioner::Print(Option_t * option)const {
+  // Print reconstructioner data  
 
-  // mark the position of the RecPoints in the array
-  AliPHOSEmcRecPoint * emcrp ; 
-  for (index = 0 ; index < emccl->GetEntries() ; index++) {
-    emcrp = (AliPHOSEmcRecPoint * )emccl->At(index) ; 
-    emcrp->SetIndexInList(index) ; 
+  TString message ; 
+  message  = "-----------------AliPHOSReconstructioner---------------\n" ;
+  message += " Reconstruction of the header file %s\n" ;
+  message += " with the following modules:\n" ;
+
+  if(fSDigitizer->IsActive()){
+    message += "   (+)   %s to branch %s\n" ; 
   }
-  AliPHOSPpsdRecPoint * ppsdrp ; 
-  for (index = 0 ; index < ppsdl->GetEntries() ; index++) {
-    ppsdrp = (AliPHOSPpsdRecPoint * )ppsdl->At(index) ; 
-    ppsdrp->SetIndexInList(index) ; 
+  if(fDigitizer->IsActive()){
+    message += "   (+)   %s to branch %s\n" ; 
   }
   
-  if  (fDebugReconstruction) {
-    cout << "DebugReconstruction>>> Cluster emc list entries is " <<    emccl->GetEntries() << endl ;
-    AliPHOSEmcRecPoint * recpoint;
-    cout << "DebugReconstruction>>> Module "  << 
-      "Ene(MeV) "             <<                         
-      "Index "                << 
-      "Multi "                << 
-      "   X     "             << 
-      "   Y     "             << 
-      "   Z    "              << 
-      " Lambda 1   "          <<  
-      " Lambda 2   "          <<
-      "MaxEnergy(MeV) "       <<
-      "Nprim "                <<
-      " Primaries list "      <<  endl;      
-    for (index = 0 ; index < emccl->GetEntries() ; index++) {
-      recpoint = (AliPHOSEmcRecPoint * )emccl->At(index) ; 
-      TVector3  locpos;  recpoint->GetLocalPosition(locpos);
-      Float_t lambda[2]; recpoint->GetElipsAxis(lambda);
-      Int_t * primaries; 
-      Int_t nprimaries;
-      primaries = recpoint->GetPrimaries(nprimaries);
-      cout << "DebugReconstruction>>>  " << 
-       setw(2) <<recpoint->GetPHOSMod() << " "  << 
-       setw(9) << 1000.*recpoint->GetTotalEnergy() <<       " "  <<                   
-       setw(6) <<  recpoint->GetIndexInList() << " "  << 
-       setw(5) <<  recpoint->GetMultiplicity() <<" "  << 
-       setw(8) <<  locpos.X() <<" "  << 
-       setw(8) <<  locpos.Y() <<" "  << 
-       setw(8) <<  locpos.Z() << " " <<
-       setw(10) << lambda[0] << "  " <<
-       setw(10) << lambda[1] << "  " <<
-       setw(9) << 1000*recpoint->GetMaximalEnergy() << "  " << 
-       setw(9) << nprimaries << "  ";
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-       cout << setw(4)  <<  primaries[iprimary] << " ";
-      cout << endl;     
-    }
-      
-    cout << "DebugReconstruction>>> Cluster ppsd list entries is " <<    ppsdl->GetEntries() << endl ;
-    AliPHOSPpsdRecPoint * ppsdrecpoint;
-    Text_t detector[4];
-    cout << "DebugReconstruction>>> Module "  << 
-      "Det     "             <<      
-      "Ene(KeV) "            <<                         
-      "Index "               << 
-      "Multi "               << 
-      "   X     "            << 
-      "   Y     "            << 
-      "   Z         "        << 
-      "Nprim "               <<
-      " Primaries list "     <<  endl;      
-    for (index = 0 ; index < ppsdl->GetEntries() ; index++) {
-      ppsdrecpoint = (AliPHOSPpsdRecPoint * ) ppsdl->At(index) ; 
-      TVector3  locpos; ppsdrecpoint->GetLocalPosition(locpos);
-      Int_t * primaries; 
-      Int_t nprimaries;
-      if (ppsdrecpoint->GetUp()) 
-       strcpy(detector, "CPV"); 
-      else 
-       strcpy(detector, "PC ");
-      primaries = ppsdrecpoint->GetPrimaries(nprimaries);
-      cout << "DebugReconstruction>>> " << 
-       setw(4) << ppsdrecpoint->GetPHOSMod() << "  "  << 
-       setw(4)  << detector << " "      <<
-       setw(9) << 1000000.*ppsdrecpoint->GetTotalEnergy() <<       " "  <<             
-       setw(6) <<  ppsdrecpoint->GetIndexInList() << " "  << 
-       setw(5) <<  ppsdrecpoint->GetMultiplicity() <<" "  << 
-       setw(8) <<  locpos.X() <<" "  << 
-       setw(8) <<  locpos.Y() <<" "  << 
-       setw(8) <<  locpos.Z() << " " <<
-       setw(9) <<  nprimaries << "  ";
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-       cout << setw(4)  <<  primaries[iprimary] << " ";
-      cout << endl;     
-    }
-  }  
-  
-  
-  if  (fDebugReconstruction)  cout << "DebugReconstruction>>>> Start making track segments(unfolding+tracksegments)" << endl;
-  fTrackSegmentMaker->MakeTrackSegments(dl, emccl, ppsdl, trsl) ;   
-  
-  // mark the position of the TrackSegments in the array
-  AliPHOSTrackSegment * trs ; 
-  for (index = 0 ; index < trsl->GetEntries() ; index++) {
-    trs = (AliPHOSTrackSegment * )trsl->At(index) ; 
-    trs->SetIndexInList(index) ; 
-  }
-  if  (fDebugReconstruction){
-    cout << "DebugReconstruction>>> Track segment list entries is " <<    trsl->GetEntries() << endl ;
-    cout << "DebugReconstruction>>> Module "  << 
-      "Ene(KeV) "             <<                         
-      "Index "                << 
-      "   X      "            << 
-      "   Y      "            << 
-      "   Z       "           <<
-      " rX        "           << 
-      " rY        "           << 
-      " rZ      "             << 
-      "Nprim "                <<
-      " Primaries list "      <<  endl;      
-    
-    for (index = 0 ; index < trsl->GetEntries() ; index++) {
-      trs = (AliPHOSTrackSegment * )trsl->At(index) ; 
-      TVector3 locpos; trs->GetPosition(locpos);
-      Int_t * primaries; 
-      Int_t nprimaries;
-      primaries = trs->GetPrimariesEmc(nprimaries);
-      cout << "DebugReconstruction>>> " << 
-       setw(4) << trs->GetPHOSMod() << "  "  << 
-       setw(9) << 1000.*trs->GetEnergy() <<       " "  <<             
-       setw(3) <<  trs->GetIndexInList() << " "  <<  
-       setw(9) <<  locpos.X() <<" "  << 
-       setw(9) <<  locpos.Y() <<" "  << 
-       setw(9) <<  locpos.Z() << " " <<
-       setw(10) <<  (trs->GetMomentumDirection()).X() << " " <<
-       setw(10) <<  (trs->GetMomentumDirection()).Y() << " " <<
-       setw(10) <<  (trs->GetMomentumDirection()).Z() << " " <<
-       setw(4) << nprimaries << "  ";
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-       cout << setw(4)  <<  primaries[iprimary] << " ";
-      cout << endl;     
-    }
-    
+  if(fClusterizer->IsActive()){
+    message += "   (+)   %s to branch %s\n" ;
   }
-  if  (fDebugReconstruction)  cout << "DebugReconstruction>>>> Start making reconstructed particles" << endl;
-  
-  if (fPID) {
-    fPID->MakeParticles(trsl, rpl) ; 
-  
-    // mark the position of the RecParticles in the array
-    AliPHOSRecParticle * rp ; 
-    for (index = 0 ; index < rpl->GetEntries() ; index++) {
-      rp = (AliPHOSRecParticle * )rpl->At(index) ; 
-      rp->SetIndexInList(index) ; 
-    }
-    //Debugger of RecParticles
-    if  (fDebugReconstruction){
-      cout << "DebugReconstruction>>>  Reconstructed particle list entries is " <<    rpl->GetEntries() << endl ;
-      cout << "DebugReconstruction>>> Module "  << 
-       "    PARTICLE     "   <<
-       "Ene(KeV) "           <<                         
-       "Index "              << 
-       "   X      "          << 
-       "   Y      "          << 
-       "   Z       "         <<
-       "Nprim "              <<
-       " Primaries list "    <<  endl;      
-      for (index = 0 ; index < rpl->GetEntries() ; index++) {
-       rp = (AliPHOSRecParticle * ) rpl->At(index) ;       
-       TVector3 locpos; (rp->GetPHOSTrackSegment())->GetPosition(locpos);
-       Int_t * primaries; 
-       Int_t nprimaries;
-       Text_t particle[11];
-       primaries = (rp->GetPHOSTrackSegment())->GetPrimariesEmc(nprimaries);
-       switch(rp->GetType()) {
-       case  AliPHOSFastRecParticle::kNEUTRALEM:
-         strcpy( particle, "NEUTRAL_EM");
-         break;
-       case  AliPHOSFastRecParticle::kNEUTRALHA:
-         strcpy(particle, "NEUTRAL_HA");
-         break;
-       case  AliPHOSFastRecParticle::kGAMMA:
-         strcpy(particle, "GAMMA");
-         break ;
-       case  AliPHOSFastRecParticle::kGAMMAHA: 
-         strcpy(particle, "GAMMA_H");
-         break ;
-       case  AliPHOSFastRecParticle::kABSURDEM:
-         strcpy(particle, "ABSURD_EM") ;
-         break ;
-       case  AliPHOSFastRecParticle::kABSURDHA:
-         strcpy(particle, "ABSURD_HA") ;
-         break ;       
-       case  AliPHOSFastRecParticle::kELECTRON:
-         strcpy(particle, "ELECTRON") ;
-         break ;
-       case  AliPHOSFastRecParticle::kCHARGEDHA:
-         strcpy(particle, "CHARGED_HA") ;
-         break ; 
-       }
-       
-       cout << "DebugReconstruction>>> " << 
-         setw(4) << (rp->GetPHOSTrackSegment())->GetPHOSMod() << "  "  <<
-         setw(15) << particle << "  " <<
-         setw(9) << 1000.*(rp->GetPHOSTrackSegment())->GetEnergy() <<       " "  <<             
-         setw(3) <<  rp->GetIndexInList() << " "  <<  
-         setw(9) <<  locpos.X() <<" "  << 
-         setw(9) <<  locpos.Y() <<" "  << 
-         setw(9) <<  locpos.Z() << " " <<
-         setw(4) << nprimaries << "  ";
-       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-         cout << setw(4)  <<  primaries[iprimary] << " ";
-       cout << endl;    
-      }
-    }
-    
+
+  if(fTSMaker->IsActive()){
+    message += "   (+)   %s to branch %s\n" ; 
   }
 
+  if(fPID->IsActive()){
+    message += "   (+)   %s to branch %s\n" ;  
+  }
+  Info("Print", message.Data(), 
+       GetTitle(), 
+       fSDigitizer->GetName(), fSDigitsBranch.Data(), 
+       fDigitizer->GetName(), fDigitsBranch.Data() , 
+       fClusterizer->GetName(), fRecPointBranch.Data(), 
+       fTSMaker->GetName(), fTSBranch.Data() , 
+       fPID->GetName(), fRecPartBranch.Data() ) ; 
 }