bug in loop, wrong varialble in increment
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPID.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //   Compute PID weights for all the clusters that are in AliESDs.root file
19 //   the AliESDs.root have to be in the same directory as the class
20 //
21 //   and do:    
22 //   AliEMCALPID *pid = new AliEMCALPID(kFALSE); // this calls the constructor which avoids the call to recparam 
23 //   pid->SetReconstructor(kFALSE);
24 //   pid->SetPrintInfo(kTRUE);
25 //   pid->SetHighFluxParam(); //   pid->SetLowFluxParam(); 
26 //   
27 //   then in cluster loop do
28 //   pid->ComputePID(energy, lambda0);
29 //        
30 //        Compute PID Weight for all clusters in AliESDs.root file
31 //        keep this function for the moment for a simple verification, could be removed
32 //
33 //   pid->GetPIDFinal(idx) gives the probabilities
34 //
35 //   Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
36 //
37 //      kElectron :  fPIDFinal[0]
38 //      kMuon     :  fPIDFinal[1]
39 //      kPion       :  fPIDFinal[2]
40 //      kKaon       :  fPIDFinal[3]
41 //      kProton   :  fPIDFinal[4]
42 //      kPhoton   :  fPIDFinal[5]
43 //      kPi0        :  fPIDFinal[6]
44 //      kNeutron  :  fPIDFinal[7]
45 //      kKaon0    :  fPIDFinal[8]
46 //      kEleCon   :  fPIDFinal[9]
47 //      kUnknown  :  fPIDFinal[10]
48 //
49 //
50 //    PID[3] is a simple PID for
51 //      Electron & Photon  PID[0]
52 //                    Pi0  PID[1]
53 //                 Hadron  PID[2]
54 //
55 // Author: Genole Bourdaud 2007 (SUBATECH)
56 //         Marie Germain 07/2009 (SUBATECH), new parametrization for low and high flux environment
57 //         Gustavo Conesa 08/2009 (LNF), divide class in AliEMCALPID and AliEMCALPIDUtils, PIDUtils belong to library EMCALUtils 
58 // --- standard c ---
59
60 // standard C++ includes
61 //#include <Riostream.h>
62
63 // ROOT includes
64
65 // STEER includes
66 #include "AliESDEvent.h"
67 #include "AliEMCALPID.h"
68 #include "AliESDCaloCluster.h"
69 #include "AliEMCALReconstructor.h"
70
71   
72 ClassImp(AliEMCALPID)
73   
74 //______________________________________________
75   AliEMCALPID::AliEMCALPID()
76         : AliEMCALPIDUtils(), fReconstructor(kTRUE)
77 {
78   //
79   // Constructor.
80   // Initialize all constant values which have to be used
81   // during PID algorithm execution
82   //
83   
84   InitParameters(); 
85   
86   
87 }
88
89 //______________________________________________
90 AliEMCALPID::AliEMCALPID(Bool_t reconstructor)
91 : AliEMCALPIDUtils(), fReconstructor(reconstructor)
92 {
93   //
94   // Constructor.
95   // Initialize all constant values which have to be used
96   // during PID algorithm execution called when used in standalone mode 
97   //
98   
99   InitParameters(); 
100   
101 }
102
103 //______________________________________________
104 void AliEMCALPID::RunPID(AliESDEvent *esd)
105 {
106   //
107   // Make the PID for all the EMCAL clusters containedin the ESDs File
108   // but just gamma/PiO/Hadron
109   //
110   // trivial check against NULL object passed
111   
112   if (esd == 0x0) {
113     AliInfo("NULL ESD object passed !!" );
114     return ;
115   }
116   
117   Int_t nClusters = esd->GetNumberOfCaloClusters();
118   Int_t firstCluster = 0;
119   Double_t energy=0., lambda0=0.;
120   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
121     
122     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
123     if (!clust->IsEMCAL()) continue ; 
124     
125     energy = clust->E();
126     lambda0 = clust->GetM02();
127    
128     if (lambda0 != 0  && energy < 1000) {
129       
130       // reject clusters with lambda0 = 0
131       
132       
133       ComputePID(energy, lambda0);
134       
135       
136       if (fPrintInfo) {
137         AliInfo("___________________________________________________");
138         AliInfo(Form( "Particle Energy = %f",energy));
139         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
140         AliInfo("PIDWeight of the particle :" );
141         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
142         AliInfo(Form( " PiZero : %f",fPID[1] ));
143         AliInfo(Form( " HADRON : %f", fPID[2] ));
144         AliInfo("_________________________________________");
145         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
146         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
147         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
148         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
149         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
150         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
151         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
152         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
153         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
154         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
155         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
156         AliInfo("___________________________________________________");
157       }
158       
159       if(fReconstructor){ // In case it is called during reconstruction.
160         //      cout << "############# Fill ESDs with PIDWeight ##########" << endl;
161         clust->SetPID(fPIDFinal);}
162     } // end if (lambda0 != 0  && energy < 1000)
163   } // end for (iCluster...)
164 }
165
166
167 //_______________________________________________________
168 void AliEMCALPID::InitParameters()
169 {
170   // Initialize PID parameters, depending on the use or not of the reconstructor
171   // and the kind of event type if the reconstructor is not used.
172   //  fWeightHadronEnergy=0.;
173   //  fWeightPiZeroEnergy=0.;
174   //  fWeightGammaEnergy=0.;
175   
176   fPIDWeight[0] = -1;
177   fPIDWeight[1] = -1;
178   fPIDWeight[2] = -1;
179   
180   for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
181     fPIDFinal[i]= 0;
182   
183   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
184   
185   if(fReconstructor){
186     
187     if(!recParam) {
188       AliFatal("Reconstruction parameters for EMCAL not set!");
189     }
190     else {
191       
192       for(Int_t i=0; i<6; i++){
193         for(Int_t j=0; j<6; j++){
194           fGamma[i][j]       = recParam->GetGamma(i,j);
195           fGamma1to10[i][j]  = recParam->GetGamma1to10(i,j);
196           fHadron[i][j]      = recParam->GetHadron(i,j);
197           fHadron1to10[i][j] = recParam->GetHadron1to10(i,j);
198           fPiZero[i][j]      = recParam->GetPiZero(i,j);
199           
200           
201           //    AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
202           //                    i,j, fGamma[i][j],fPiZero[i][j],fHadron[i][j] ));
203           //    cout << "PID parameters (" << i << " ,"<<j<<") fGamma= "<<  fGamma[i][j]<<" fPi0 ="<<  fPiZero[i][j]<< endl;
204           
205         } // end loop j
206         fHadronEnergyProb[i] = recParam->GetHadronEnergyProb(i);
207         fPiZeroEnergyProb[i] = recParam->GetPiZeroEnergyProb(i);
208         fGammaEnergyProb[i]  = recParam->GetGammaEnergyProb(i);
209       } //end loop i
210       
211       
212     } // end if !recparam 
213     
214   } 
215   
216   else{
217     //   init the parameters here instead of from loading from recparam
218     //   default parameters are PbPb parameters.
219     SetHighFluxParam();
220     
221   }
222   
223 }
224