]> git.uio.no Git - u/mrichter/AliRoot.git/blob - EMCAL/AliEMCALPID.cxx
Update tracking class and macros, and pid class to read correctly the ESDs
[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 /* History of cvs commits:
18  *
19  * $Log$
20  * Revision 1.16  2007/11/23 13:39:05  gustavo
21  * Track matching and PID parameters added to AliEMCALRecParam
22  *
23  * Revision 1.15  2007/10/09 08:46:10  hristov
24  * The data members fEMCALClusterCluster and fPHOSCluster are removed from AliESDCaloCluster, the fClusterType is used to select PHOS or EMCAL clusters. Changes, needed to use correctly the new AliESDCaloCluster. (Christian)
25  *
26  * Revision 1.14  2007/07/26 16:54:53  morsch
27  * Changes in AliESDEvent fwd declarartions.
28  *
29  * Revision 1.13  2007/07/11 13:43:29  hristov
30  * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
31  *
32  * Revision 1.12  2007/06/11 20:43:06  hristov
33  * Changes required by the updated AliESDCaloCluster (Gustavo)
34  *
35  * Revision 1.11  2007/03/30 13:50:34  gustavo
36  * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
37  *
38  * Revision 1.10  2007/03/09 14:34:11  gustavo
39  * Correct probability calculation, added missing initialization of data members
40  *
41  * Revision 1.9  2007/02/20 20:17:43  hristov
42  * Corrected array size, removed warnings (icc)
43  *
44  * Revision 1.8  2006/12/19 08:49:35  gustavo
45  * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
46  *
47  *
48  */
49 //    to compute PID for all the clusters in ESDs.root file
50 //     the ESDs.root have to be in the same directory as the class
51 //
52 //
53 //
54 //
55 //
56 //      AliEMCALPID::CalculPID(Energy,Lambda0)
57 //        Calcul PID for all clusters in AliESDs.root file
58 //        keep this function for the moment for a simple verification, could be removed
59 //
60 //
61 //
62 //   AliEMCALPID::CalculPID(Energy,Lambda0)
63 //    calcul PID Weght for a cluster with Energy, Lambda0 .
64 //    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
65 //
66 //
67 //
68 //      kElectron :  fPIDFinal[0]
69 //      kMuon     :  fPIDFinal[1]
70 //      kPion     :  fPIDFinal[2]
71 //      kKaon     :  fPIDFinal[3]
72 //      kProton   :  fPIDFinal[4]
73 //      kPhoton   :  fPIDFinal[5]
74 //      kPi0      :  fPIDFinal[6]
75 //      kNeutron  :  fPIDFinal[7]
76 //      kKaon0    :  fPIDFinal[8]
77 //      kEleCon   :  fPIDFinal[9]
78 //      kUnknown  :  fPIDFinal[10]
79 //
80 //
81 //    PID[3] is a simple PID for
82 //      Electron & Photon  PID[0]
83 //                    Pi0  PID[1]
84 //                 Hadron  PID[2]
85 //
86 //  
87 //
88 //
89 //
90 // --- ROOT system ---
91
92 // standard C++ includes
93 #include <Riostream.h>
94
95 // ROOT includes
96 #include "TTree.h"
97 #include "TStyle.h"
98 #include "TVector3.h"
99 #include "TBranch.h"
100 #include "TClonesArray.h"
101 #include "TCanvas.h"
102 #include "TLorentzVector.h"
103 #include "TMath.h"
104 #include "TFile.h"
105 #include "TH1.h"
106 #include "TH2.h"
107 #include "TParticle.h"
108
109 // STEER includes
110 #include "AliLog.h"
111 #include "AliEMCALPID.h"
112 #include "AliESDCaloCluster.h"
113 #include "AliEMCALRecParam.h"
114 #include "AliEMCALReconstructor.h"
115
116   
117 ClassImp(AliEMCALPID)
118
119 //______________________________________________
120   AliEMCALPID::AliEMCALPID():
121     fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
122 {
123   //
124   // Constructor.
125   // Initialize all constant values which have to be used
126   // during PID algorithm execution
127   //
128  
129   fPIDWeight[0] = -1;
130   fPIDWeight[1] = -1;
131   fPIDWeight[2] = -1;
132
133   for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
134     fPIDFinal[i]= 0;
135
136   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
137   if(!recParam) {
138     AliFatal("Reconstruction parameters for EMCAL not set!");
139   }
140   else {
141     for(Int_t i=0; i<6; i++){
142       for(Int_t j=0; j<6; j++){
143         fGamma[i][j] = recParam->GetGamma(i,j);
144         fHadron[i][j] = recParam->GetHadron(i,j);
145         fPiZero5to10[i][j] = recParam->GetPiZero5to10(i,j);
146         fPiZero10to60[i][j] = recParam->GetPiZero10to60(i,j);
147         AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
148                         i,j, fGamma[i][j],fPiZero5to10[i][j],fHadron[i][j] ));
149       }
150     }
151     
152   }
153   
154 }
155
156 //______________________________________________
157 void AliEMCALPID::RunPID(AliESDEvent *esd)
158 {
159 //
160 // Make the PID for all the EMCAL clusters containedin the ESDs File
161 // but just gamma/PiO/Hadron
162 //
163         // trivial check against NULL object passed
164
165   if (esd == 0x0) {
166     AliInfo("NULL ESD object passed !!" );
167     return ;
168   }
169
170   Int_t nClusters = esd->GetNumberOfCaloClusters();
171   Int_t firstCluster = 0;
172   Double_t energy, lambda0;
173   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
174
175     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
176         if (!clust->IsEMCAL()) continue ; 
177     energy = clust->E();
178     lambda0 = clust->GetM02();
179     // verify cluster type
180     Int_t clusterType= clust->GetClusterType();
181     if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0  && energy < 1000) {
182
183
184       // reject clusters with lambda0 = 0
185
186
187       ComputePID(energy, lambda0);
188
189
190       if (fPrintInfo) {
191         AliInfo("___________________________________________________");
192         AliInfo(Form( "Particle Energy = %f",energy));
193         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
194         AliInfo("PIDWeight of the particle :" );
195         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
196         AliInfo(Form( " PiZero : %f",fPID[1] ));
197         AliInfo(Form( " HADRON : %f", fPID[2] ));
198         AliInfo("_________________________________________");
199         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
200         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
201         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
202         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
203         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
204         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
205         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
206         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
207         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
208         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
209         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
210         AliInfo("___________________________________________________");
211       }
212
213       if(fReconstructor) // In case it is called during reconstruction.
214         clust->SetPid(fPIDFinal);
215     } // end if (clusterType...)
216   } // end for (iCluster...)
217 }
218
219 //__________________________________________________________
220 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
221 {
222 //
223 // This is the main command, which uses the distributions computed and parametrised, 
224 // and gives the PID by the bayesian method.
225 //
226
227 if (energy<5){energy =6;}
228
229
230   TArrayD paramDistribGamma  = DistLambda0(energy, 1);
231   TArrayD paramDistribPiZero = DistLambda0(energy, 2);
232   TArrayD paramDistribHadron = DistLambda0(energy, 3);
233   
234   Bool_t norm = kFALSE;
235   
236   fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
237   fProbGamma  += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
238   fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
239   fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
240   fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
241   fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
242   
243   // compute PID Weight
244   fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
245   fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
246   fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
247   
248   SetPID(fPIDWeight[0], 0);
249   SetPID(fPIDWeight[1], 1);
250   SetPID(fPIDWeight[2], 2);
251   
252   // sortie ecran pid Weight only for control (= in english ???)
253   if (fPrintInfo) {
254     AliInfo(Form( "Energy in loop = %f", energy) );
255     AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
256     AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
257     //  AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
258     AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
259     AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
260     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
261     AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
262     AliInfo("********************************************************" );
263   }
264   
265   fPIDFinal[0]  = fPIDWeight[0]/2;
266   fPIDFinal[1]  = fPIDWeight[2]/8;
267   fPIDFinal[2]  = fPIDWeight[2]/8;
268   fPIDFinal[3]  = fPIDWeight[2]/8;
269   fPIDFinal[4]  = fPIDWeight[2]/8;
270   fPIDFinal[5]  = fPIDWeight[0]/2;
271   fPIDFinal[6]  = fPIDWeight[1]  ;
272   fPIDFinal[7]  = fPIDWeight[2]/8;
273   fPIDFinal[8]  = fPIDWeight[2]/8;
274   fPIDFinal[9]  = fPIDWeight[2]/8;
275   fPIDFinal[10] = fPIDWeight[2]/8;
276 }
277
278 //________________________________________________________
279 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
280 {
281   //
282   // Compute the values of the parametrised distributions using the data initialised before.
283   //
284   Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
285   Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
286   TArrayD  distributionParam(6);
287   
288   switch (type) {
289   case 1:
290     constGauss  = Polynomial(energy, fGamma[0]);
291     meanGauss   = Polynomial(energy, fGamma[1]);
292     sigmaGauss  = Polynomial(energy, fGamma[2]);
293     constLandau = Polynomial(energy, fGamma[3]);
294     mpvLandau   = Polynomial(energy, fGamma[4]);
295     sigmaLandau = Polynomial(energy, fGamma[5]);
296     
297
298     break;
299   case 2:
300     if (energy < 10) {
301       constGauss  = Polynomial(energy, fPiZero5to10[0]);
302       meanGauss   = Polynomial(energy, fPiZero5to10[1]);
303       sigmaGauss  = Polynomial(energy, fPiZero5to10[2]);
304       constLandau = Polynomial(energy, fPiZero5to10[3]);
305       mpvLandau   = Polynomial(energy, fPiZero5to10[4]);
306       sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
307
308     }
309     else {
310       constGauss  = Polynomial(energy, fPiZero10to60[0]);
311       meanGauss   = Polynomial(energy, fPiZero10to60[1]);
312       sigmaGauss  = Polynomial(energy, fPiZero10to60[2]);
313       constLandau = Polynomial(energy, fPiZero10to60[3]);
314       mpvLandau   = Polynomial(energy, fPiZero10to60[4]);
315       sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
316  
317     }
318     break;
319   case 3:
320     constGauss  = Polynomial(energy, fHadron[0]);
321     meanGauss   = Polynomial(energy, fHadron[1]);
322     sigmaGauss  = Polynomial(energy, fHadron[2]);
323     constLandau = Polynomial(energy, fHadron[3]);
324     mpvLandau   = Polynomial(energy, fHadron[4]);
325     sigmaLandau = Polynomial(energy, fHadron[5]);
326
327     break;
328   }
329   
330   distributionParam[0] = constGauss;
331   distributionParam[1] = meanGauss;
332   distributionParam[2] = sigmaGauss;
333   distributionParam[3] = constLandau;
334   distributionParam[4] = mpvLandau;
335   distributionParam[5] = sigmaLandau;
336   
337   return distributionParam;
338 }
339
340 //_______________________________________________________
341 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
342 {
343   //
344   // Compute a polynomial for a given value of 'x'
345   // with the array of parameters passed as the second arg
346   //
347   
348   Double_t y;
349   y  = params[0];
350   y += params[1] * x;
351   y += params[2] * x * x;
352   y += params[3] * x * x * x;
353   y += params[4] * x * x * x * x;
354   y += params[5] * x * x * x * x * x;
355   
356   return y;
357 }