Full PID delegated to AliPHOSPID
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.cxx
index 94abec4971ab7e8c8d9da0b9fc5c731fd84016e9..3f2d8ffad1d7519d2722535ad520f1fc29e22d02 100644 (file)
@@ -41,7 +41,7 @@
 #include "AliPHOSAnalyze.h"
 #include "AliPHOSClusterizerv1.h"
 #include "AliPHOSTrackSegmentMakerv1.h"
-#include "AliPHOSParticleGuesserv1.h"
+#include "AliPHOSPIDv1.h"
 #include "AliPHOSReconstructioner.h"
 #include "AliPHOSDigit.h"
 #include "AliPHOSTrackSegment.h"
@@ -90,8 +90,8 @@ AliPHOSAnalyze::~AliPHOSAnalyze()
   delete fClu ; 
   fClu = 0 ;
 
-  delete fPag ; 
-  fPag = 0 ;
+  delete fPID ; 
+  fPID = 0 ;
 
   delete fRec ; 
   fRec = 0 ;
@@ -148,7 +148,8 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
       AliPHOSDigit * digit ;
       AliPHOSEmcRecPoint * emc ;
       AliPHOSPpsdRecPoint * ppsd ;
-      AliPHOSTrackSegment * tracksegment ;
+      //      AliPHOSTrackSegment * tracksegment ;
+      AliPHOSRecParticle * recparticle;
       for ( ievent=0; ievent<Nevents; ievent++)
        {  
           if (ievent==0)  cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ; 
@@ -162,10 +163,12 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
          fClu->SetCalibrationParameters(0., 0.00000001) ;  
          //========== Creates the track segment maker
          fTrs = new AliPHOSTrackSegmentMakerv1()  ; 
-         //========== Creates the particle guesser
-         fPag = new AliPHOSParticleGuesserv1() ;
+         //========== Creates the particle identifier
+         fPID = new AliPHOSPIDv1() ;
+         fPID->SetShowerProfileCuts(0.5, 1.5, 0.5, 1.5 ) ; 
+         fPID->Print() ;           
          //========== Creates the Reconstructioner  
-         fRec = new AliPHOSReconstructioner(fClu, fTrs, fPag) ; 
+         fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
          //========== Event Number
          if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;  
          //=========== Connects the various Tree's for evt
@@ -214,45 +217,52 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
                }
            }
          //========== TRackSegments in the event
-         TIter nextTrackSegment(fPHOS->TrackSegments() ) ; 
-         while((tracksegment = (AliPHOSTrackSegment *)nextTrackSegment())) 
+         TIter nextRecParticle(fPHOS->RecParticles() ) ; 
+         while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
            {
-             if ( tracksegment->GetPHOSMod() == module )
+             if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
                { 
-                 switch(tracksegment->GetPartType())
+                 cout << "Particle type is " << recparticle->GetType() << endl ;  
+                 switch(recparticle->GetType())
                    {
-                   case 0:
-                     fhPhotonEnergy->Fill(tracksegment->GetEnergy() ) ; 
-                     //                      fhPhotonPositionX->Fill(tracksegment-> ) ;
-                     //fhPhotonPositionY->Fill(tracksegment-> ) ;                 
-                     //cout << "PHOTON" << endl;
+                   case kGAMMA:
+                     fhPhotonEnergy->Fill(recparticle->Energy() ) ; 
+                     //fhPhotonPositionX->Fill(recpart. ) ;
+                     //fhPhotonPositionY->Fill(recpart. ) ;                 
+                     cout << "PHOTON" << endl;
                      break;
-                   case :
-                     fhElectronEnergy->Fill(tracksegment->GetEnergy() ) ; 
-                     //fhElectronPositionX->Fill(tracksegment-> ) ;
-                     //fhElectronPositionY->Fill(tracksegment-> ) ; 
-                     //cout << "ELECTRON" << endl;
+                   case kELECTRON:
+                     fhElectronEnergy->Fill(recparticle->Energy() ) ; 
+                     //fhElectronPositionX->Fill(recpart. ) ;
+                     //fhElectronPositionY->Fill(recpart. ) ; 
+                     cout << "ELECTRON" << endl;
                      break;
-                   case :
-                     fhNeutralEnergy->Fill(tracksegment->GetEnergy() ) ; 
-                     //fhNeutralPositionX->Fill(tracksegment-> ) ;
-                     //fhNeutralPositionY->Fill(tracksegment-> ) ; 
-                     //cout << "NEUTRAL" << endl;
+                   case kNEUTRALHADRON:
+                     fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; 
+                     //fhNeutralHadronPositionX->Fill(recpart. ) ;
+                     //fhNeutralHadronPositionY->Fill(recpart. ) ; 
+                     cout << "NEUTRAl HADRON" << endl;
                      break ;
-                   case 3 :
-                     fhChargedEnergy->Fill(tracksegment->GetEnergy() ) ; 
-                     //fhChargedPositionX->Fill(tracksegment-> ) ;
-                     //fhChargedPositionY->Fill(tracksegment-> ) ; 
-                     //cout << "CHARGED" << endl;
+                   case kNEUTRALEM:
+                     fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; 
+                     //fhNeutralEMPositionX->Fill(recpart. ) ;
+                     //fhNeutralEMPositionY->Fill(recpart. ) ; 
+                     //cout << "NEUTRAL EM" << endl;
+                     break ;
+                   case kCHARGEDHADRON:
+                     fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; 
+                     //fhChargedHadronPositionX->Fill(recpart. ) ;
+                     //fhChargedHadronPositionY->Fill(recpart. ) ; 
+                     cout << "CHARGED HADRON" << endl;
                      break ;
                      
                    }
                }
            }
-         // Deleting fClu, fTrs, fPag et fRec
+         // Deleting fClu, fTrs, fPID et fRec
          fClu->Delete();
          fTrs->Delete();
-         fPag->Delete();
+         fPID->Delete();
          fRec->Delete();
 
        }   // endfor
@@ -264,32 +274,43 @@ void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
 //____________________________________________________________________________
 void  AliPHOSAnalyze::BookingHistograms()
 {
-  if (fhEmcDigit )         delete fhEmcDigit  ;
-  if (fhVetoDigit )        delete fhVetoDigit  ;
-  if (fhConvertorDigit )   delete fhConvertorDigit   ;
-  if (fhEmcCluster   )     delete  fhEmcCluster   ;
-  if (fhVetoCluster )      delete fhVetoCluster   ;
-  if (fhConvertorCluster ) delete fhConvertorCluster  ;
-  if (fhConvertorEmc )     delete fhConvertorEmc  ;
-  fhEmcDigit         = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
-  fhVetoDigit        = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
-  fhConvertorDigit   = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
-  fhEmcCluster       = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
-  fhVetoCluster      = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
-  fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
-  fhConvertorEmc     = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
-  fhPhotonEnergy     = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
-  fhElectronEnergy   = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
-  fhNeutralEnergy    = new TH1F("hNeutralEnergy", "hNeutralEnergy",    1000,  0. ,  30.);
-  fhChargedEnergy    = new TH1F("hChargedEnergy", "hChargedEnergy",    1000,  0. ,  30.);
-  fhPhotonPositionX  = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
-  fhElectronPositionX= new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
-  fhNeutralPositionX  = new TH1F("hNeutralPositionX","hNeutralPositionX",500,-80. , 80.);
-  fhChargedPositionX  = new TH1F("hChargedPositionX","hChargedPositionX",500,-80. , 80.);
-  fhPhotonPositionY  = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
-  fhElectronPositionY= new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
-  fhNeutralPositionY  = new TH1F("hNeutralPositionY","hNeutralPositionY",500,-80. , 80.);
-  fhChargedPositionY  = new TH1F("hChargedPositionY","hChargedPositionY",500,-80. , 80.);
+  if (fhEmcDigit )         
+    delete fhEmcDigit  ;
+  if (fhVetoDigit )     
+    delete fhVetoDigit  ;
+  if (fhConvertorDigit ) 
+    delete fhConvertorDigit   ;
+  if (fhEmcCluster   )  
+    delete  fhEmcCluster   ;
+  if (fhVetoCluster )     
+    delete fhVetoCluster   ;
+  if (fhConvertorCluster )
+    delete fhConvertorCluster  ;
+  if (fhConvertorEmc )    
+    delete fhConvertorEmc  ;
+  fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
+  fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
+  fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
+  fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
+  fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
+  fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
+  fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
+  fhPhotonEnergy            = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
+  fhElectronEnergy          = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
+  fhNeutralHadronEnergy     = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy",    1000,  0. ,  30.);
+  fhNeutralEMEnergy         = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy",    1000,  0. ,  30.);
+  fhChargedHadronEnergy     = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy",    1000,  0. ,  30.);
+  fhPhotonPositionX         = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
+  fhElectronPositionX       = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
+  fhNeutralHadronPositionX  = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
+  fhNeutralEMPositionX      = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
+  fhChargedHadronPositionX  = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
+  fhPhotonPositionY         = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
+  fhElectronPositionY       = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
+  fhNeutralHadronPositionY  = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
+  fhNeutralEMPositionY      = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
+  fhChargedHadronPositionY  = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
 
 }
 //____________________________________________________________________________
@@ -338,14 +359,14 @@ Bool_t AliPHOSAnalyze::Init(Int_t evt)
     fTrs = new AliPHOSTrackSegmentMakerv1() ;
     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; 
     
-    //========== Creates the particle guesser
+    //========== Creates the particle identifier
     
-    fPag = new AliPHOSParticleGuesserv1() ;
-    cout <<  "AnalyzeOneEvent > using particle guess " << fPag->GetName() << endl ; 
+    fPID = new AliPHOSPIDv1() ;
+    cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; 
     
     //========== Creates the Reconstructioner  
     
-    fRec = new AliPHOSReconstructioner(fClu, fTrs, fPag) ;     
+    fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;     
     
     //=========== Connect the various Tree's for evt
     
@@ -594,13 +615,16 @@ void AliPHOSAnalyze::DisplayRecPoints()
       AliPHOSEmcRecPoint * emc ;
       while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
        {
+         Int_t numberofprimaries ;
+         Int_t * primariesarray = new Int_t[10] ;
+         emc->GetPrimaries(numberofprimaries, primariesarray) ;
          totalnClusters++ ;
          if ( emc->GetPHOSMod() == module )
            { 
              nClusters++ ; 
              energy = emc->GetTotalEnergy() ;   
              etot+= energy ;  
-             emc->Draw("P") ;
+             emc->Draw("M") ;
            }
        }
       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; 
@@ -740,25 +764,50 @@ void AliPHOSAnalyze::SavingHistograms()
   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
   TFile output(outputname,"RECREATE");
   output.cd();
-  if (fhEmcDigit )         fhEmcDigit->Write()  ;
-  if (fhVetoDigit )        fhVetoDigit->Write()  ;
-  if (fhConvertorDigit )   fhConvertorDigit->Write()   ;
-  if (fhEmcCluster   )     fhEmcCluster->Write()   ;
-  if (fhVetoCluster )      fhVetoCluster->Write()   ;
-  if (fhConvertorCluster ) fhConvertorCluster->Write()  ;
-  if (fhConvertorEmc )     fhConvertorEmc->Write()  ;
-  if (fhPhotonEnergy)      fhPhotonEnergy->Write() ;
-  if (fhPhotonPositionX)   fhPhotonPositionX->Write() ;
-  if (fhPhotonPositionY)   fhPhotonPositionX->Write() ;
-  if (fhElectronEnergy)      fhElectronEnergy->Write() ;
-  if (fhElectronPositionX) fhElectronPositionX->Write() ;
-  if (fhElectronPositionY) fhElectronPositionX->Write() ;
-  if (fhNeutralEnergy)     fhNeutralEnergy->Write() ;
-  if (fhNeutralPositionX)  fhNeutralPositionX->Write() ;
-  if (fhNeutralPositionY)  fhNeutralPositionX->Write() ;
-  if (fhChargedEnergy)     fhChargedEnergy->Write() ;
-  if (fhChargedPositionX)  fhChargedPositionX->Write() ;
-  if (fhChargedPositionY)  fhChargedPositionX->Write() ;
+  if (fhEmcDigit )         
+    fhEmcDigit->Write()  ;
+  if (fhVetoDigit )  
+    fhVetoDigit->Write()  ;
+  if (fhConvertorDigit ) 
+    fhConvertorDigit->Write()   ;
+  if (fhEmcCluster   )
+    fhEmcCluster->Write()   ;
+  if (fhVetoCluster ) 
+    fhVetoCluster->Write()   ;
+  if (fhConvertorCluster )
+    fhConvertorCluster->Write()  ;
+  if (fhConvertorEmc ) 
+    fhConvertorEmc->Write()  ;
+  if (fhPhotonEnergy)    
+    fhPhotonEnergy->Write() ;
+  if (fhPhotonPositionX)  
+    fhPhotonPositionX->Write() ;
+  if (fhPhotonPositionY)  
+    fhPhotonPositionX->Write() ;
+  if (fhElectronEnergy)  
+    fhElectronEnergy->Write() ;
+  if (fhElectronPositionX)
+    fhElectronPositionX->Write() ;
+  if (fhElectronPositionY) 
+    fhElectronPositionX->Write() ;
+  if (fhNeutralHadronEnergy) 
+    fhNeutralHadronEnergy->Write() ;
+  if (fhNeutralHadronPositionX)
+    fhNeutralHadronPositionX->Write() ;
+  if (fhNeutralHadronPositionY) 
+    fhNeutralHadronPositionX->Write() ;
+  if (fhNeutralEMEnergy)   
+    fhNeutralEMEnergy->Write() ;
+  if (fhNeutralEMPositionX)
+    fhNeutralEMPositionX->Write() ;
+  if (fhNeutralEMPositionY) 
+    fhNeutralEMPositionX->Write() ;
+  if (fhChargedHadronEnergy) 
+    fhChargedHadronEnergy->Write() ;
+  if (fhChargedHadronPositionX) 
+    fhChargedHadronPositionX->Write() ;
+  if (fhChargedHadronPositionY)
+    fhChargedHadronPositionX->Write() ;
 
   output.Write();
   output.Close();