]>
Commit | Line | Data |
---|---|---|
dc293ae9 | 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$ | |
46363b15 | 20 | * Revision 1.16 2007/11/23 13:39:05 gustavo |
21 | * Track matching and PID parameters added to AliEMCALRecParam | |
22 | * | |
8ba062b1 | 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 | * | |
8ada0ffe | 26 | * Revision 1.14 2007/07/26 16:54:53 morsch |
27 | * Changes in AliESDEvent fwd declarartions. | |
28 | * | |
89ffc0b0 | 29 | * Revision 1.13 2007/07/11 13:43:29 hristov |
30 | * New class AliESDEvent, backward compatibility with the old AliESD (Christian) | |
31 | * | |
af885e0f | 32 | * Revision 1.12 2007/06/11 20:43:06 hristov |
33 | * Changes required by the updated AliESDCaloCluster (Gustavo) | |
34 | * | |
7592dfc4 | 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 | * | |
4ce73a00 | 38 | * Revision 1.10 2007/03/09 14:34:11 gustavo |
39 | * Correct probability calculation, added missing initialization of data members | |
40 | * | |
4cfd8ef9 | 41 | * Revision 1.9 2007/02/20 20:17:43 hristov |
42 | * Corrected array size, removed warnings (icc) | |
43 | * | |
d69ab345 | 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 | * | |
dc293ae9 | 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> | |
dc293ae9 | 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 | ||
46363b15 | 109 | // STEER includes |
dc293ae9 | 110 | #include "AliLog.h" |
111 | #include "AliEMCALPID.h" | |
89ffc0b0 | 112 | #include "AliESDCaloCluster.h" |
8ba062b1 | 113 | #include "AliEMCALRecParam.h" |
114 | #include "AliEMCALReconstructor.h" | |
115 | ||
dc293ae9 | 116 | |
117 | ClassImp(AliEMCALPID) | |
118 | ||
4cfd8ef9 | 119 | //______________________________________________ |
120 | AliEMCALPID::AliEMCALPID(): | |
121 | fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE) | |
dc293ae9 | 122 | { |
4cfd8ef9 | 123 | // |
124 | // Constructor. | |
125 | // Initialize all constant values which have to be used | |
126 | // during PID algorithm execution | |
127 | // | |
46363b15 | 128 | |
4cfd8ef9 | 129 | fPIDWeight[0] = -1; |
130 | fPIDWeight[1] = -1; | |
131 | fPIDWeight[2] = -1; | |
8ba062b1 | 132 | |
46363b15 | 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) { | |
8ba062b1 | 138 | AliFatal("Reconstruction parameters for EMCAL not set!"); |
139 | } | |
140 | else { | |
46363b15 | 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 | ||
8ba062b1 | 152 | } |
46363b15 | 153 | |
dc293ae9 | 154 | } |
4cfd8ef9 | 155 | |
156 | //______________________________________________ | |
af885e0f | 157 | void AliEMCALPID::RunPID(AliESDEvent *esd) |
dc293ae9 | 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 | |
8ba062b1 | 164 | |
dc293ae9 | 165 | if (esd == 0x0) { |
4ce73a00 | 166 | AliInfo("NULL ESD object passed !!" ); |
dc293ae9 | 167 | return ; |
168 | } | |
4cfd8ef9 | 169 | |
dc293ae9 | 170 | Int_t nClusters = esd->GetNumberOfEMCALClusters(); |
171 | Int_t firstCluster = esd->GetFirstEMCALCluster(); | |
172 | Double_t energy, lambda0; | |
173 | for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) { | |
8ba062b1 | 174 | |
dc293ae9 | 175 | AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster); |
7592dfc4 | 176 | energy = clust->E(); |
dc293ae9 | 177 | lambda0 = clust->GetM02(); |
178 | // verify cluster type | |
179 | Int_t clusterType= clust->GetClusterType(); | |
8ada0ffe | 180 | if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0 && energy < 1000) { |
4ce73a00 | 181 | |
182 | ||
dc293ae9 | 183 | // reject clusters with lambda0 = 0 |
4ce73a00 | 184 | |
185 | ||
dc293ae9 | 186 | ComputePID(energy, lambda0); |
4ce73a00 | 187 | |
188 | ||
dc293ae9 | 189 | if (fPrintInfo) { |
190 | AliInfo("___________________________________________________"); | |
191 | AliInfo(Form( "Particle Energy = %f",energy)); | |
192 | AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) ); | |
193 | AliInfo("PIDWeight of the particle :" ); | |
194 | AliInfo(Form( " GAMMA : %f",fPID[0] )); | |
195 | AliInfo(Form( " PiZero : %f",fPID[1] )); | |
196 | AliInfo(Form( " HADRON : %f", fPID[2] )); | |
197 | AliInfo("_________________________________________"); | |
198 | AliInfo(Form( " kElectron : %f", fPIDFinal[0]) ); | |
199 | AliInfo(Form( " kMuon : %f", fPIDFinal[1] )); | |
200 | AliInfo(Form( " kPion : %f", fPIDFinal[2] )); | |
201 | AliInfo(Form( " kKaon : %f", fPIDFinal[3] )); | |
202 | AliInfo(Form( " kProton : %f", fPIDFinal[4] )); | |
203 | AliInfo(Form( " kPhoton : %f", fPIDFinal[5] )); | |
204 | AliInfo(Form( " kPi0 : %f", fPIDFinal[6] )); | |
205 | AliInfo(Form( " kNeutron : %f", fPIDFinal[7] )); | |
206 | AliInfo(Form( " kKaon0 : %f", fPIDFinal[8] )); | |
207 | AliInfo(Form( " kEleCon : %f", fPIDFinal[9] )); | |
208 | AliInfo(Form( " kUnknown : %f", fPIDFinal[10] )); | |
209 | AliInfo("___________________________________________________"); | |
210 | } | |
4cfd8ef9 | 211 | |
dc293ae9 | 212 | if(fReconstructor) // In case it is called during reconstruction. |
213 | clust->SetPid(fPIDFinal); | |
214 | } // end if (clusterType...) | |
215 | } // end for (iCluster...) | |
216 | } | |
4cfd8ef9 | 217 | |
218 | //__________________________________________________________ | |
dc293ae9 | 219 | void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0) |
220 | { | |
221 | // | |
222 | // This is the main command, which uses the distributions computed and parametrised, | |
223 | // and gives the PID by the bayesian method. | |
224 | // | |
4ce73a00 | 225 | |
226 | if (energy<5){energy =6;} | |
227 | ||
228 | ||
4cfd8ef9 | 229 | TArrayD paramDistribGamma = DistLambda0(energy, 1); |
230 | TArrayD paramDistribPiZero = DistLambda0(energy, 2); | |
231 | TArrayD paramDistribHadron = DistLambda0(energy, 3); | |
232 | ||
233 | Bool_t norm = kFALSE; | |
234 | ||
235 | fProbGamma = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0]; | |
236 | fProbGamma += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3]; | |
237 | fProbPiZero = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0]; | |
238 | fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3]; | |
239 | fProbHadron = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0]; | |
240 | fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3]; | |
241 | ||
242 | // compute PID Weight | |
243 | fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron); | |
244 | fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron); | |
245 | fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron); | |
246 | ||
247 | SetPID(fPIDWeight[0], 0); | |
248 | SetPID(fPIDWeight[1], 1); | |
249 | SetPID(fPIDWeight[2], 2); | |
250 | ||
251 | // sortie ecran pid Weight only for control (= in english ???) | |
252 | if (fPrintInfo) { | |
253 | AliInfo(Form( "Energy in loop = %f", energy) ); | |
254 | AliInfo(Form( "Lambda0 in loop = %f", lambda0) ); | |
255 | AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) ); | |
256 | // AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) ); | |
257 | AliInfo(Form( "fProbaPiZero = %f", fProbPiZero )); | |
258 | AliInfo(Form( "fProbaHadron = %f", fProbHadron) ); | |
259 | AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f", fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) ); | |
260 | AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] )); | |
261 | AliInfo("********************************************************" ); | |
262 | } | |
263 | ||
264 | fPIDFinal[0] = fPIDWeight[0]/2; | |
265 | fPIDFinal[1] = fPIDWeight[2]/8; | |
266 | fPIDFinal[2] = fPIDWeight[2]/8; | |
267 | fPIDFinal[3] = fPIDWeight[2]/8; | |
268 | fPIDFinal[4] = fPIDWeight[2]/8; | |
269 | fPIDFinal[5] = fPIDWeight[0]/2; | |
270 | fPIDFinal[6] = fPIDWeight[1] ; | |
271 | fPIDFinal[7] = fPIDWeight[2]/8; | |
272 | fPIDFinal[8] = fPIDWeight[2]/8; | |
273 | fPIDFinal[9] = fPIDWeight[2]/8; | |
274 | fPIDFinal[10] = fPIDWeight[2]/8; | |
dc293ae9 | 275 | } |
4cfd8ef9 | 276 | |
277 | //________________________________________________________ | |
dc293ae9 | 278 | TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type) |
279 | { | |
4cfd8ef9 | 280 | // |
281 | // Compute the values of the parametrised distributions using the data initialised before. | |
282 | // | |
283 | Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.; | |
284 | Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.; | |
285 | TArrayD distributionParam(6); | |
286 | ||
287 | switch (type) { | |
288 | case 1: | |
289 | constGauss = Polynomial(energy, fGamma[0]); | |
290 | meanGauss = Polynomial(energy, fGamma[1]); | |
291 | sigmaGauss = Polynomial(energy, fGamma[2]); | |
292 | constLandau = Polynomial(energy, fGamma[3]); | |
293 | mpvLandau = Polynomial(energy, fGamma[4]); | |
294 | sigmaLandau = Polynomial(energy, fGamma[5]); | |
46363b15 | 295 | |
296 | ||
4cfd8ef9 | 297 | break; |
298 | case 2: | |
299 | if (energy < 10) { | |
300 | constGauss = Polynomial(energy, fPiZero5to10[0]); | |
301 | meanGauss = Polynomial(energy, fPiZero5to10[1]); | |
302 | sigmaGauss = Polynomial(energy, fPiZero5to10[2]); | |
303 | constLandau = Polynomial(energy, fPiZero5to10[3]); | |
304 | mpvLandau = Polynomial(energy, fPiZero5to10[4]); | |
305 | sigmaLandau = Polynomial(energy, fPiZero5to10[5]); | |
46363b15 | 306 | |
4cfd8ef9 | 307 | } |
308 | else { | |
309 | constGauss = Polynomial(energy, fPiZero10to60[0]); | |
310 | meanGauss = Polynomial(energy, fPiZero10to60[1]); | |
311 | sigmaGauss = Polynomial(energy, fPiZero10to60[2]); | |
312 | constLandau = Polynomial(energy, fPiZero10to60[3]); | |
313 | mpvLandau = Polynomial(energy, fPiZero10to60[4]); | |
314 | sigmaLandau = Polynomial(energy, fPiZero10to60[5]); | |
46363b15 | 315 | |
4cfd8ef9 | 316 | } |
317 | break; | |
318 | case 3: | |
319 | constGauss = Polynomial(energy, fHadron[0]); | |
320 | meanGauss = Polynomial(energy, fHadron[1]); | |
321 | sigmaGauss = Polynomial(energy, fHadron[2]); | |
322 | constLandau = Polynomial(energy, fHadron[3]); | |
323 | mpvLandau = Polynomial(energy, fHadron[4]); | |
324 | sigmaLandau = Polynomial(energy, fHadron[5]); | |
46363b15 | 325 | |
4cfd8ef9 | 326 | break; |
327 | } | |
328 | ||
329 | distributionParam[0] = constGauss; | |
330 | distributionParam[1] = meanGauss; | |
331 | distributionParam[2] = sigmaGauss; | |
332 | distributionParam[3] = constLandau; | |
333 | distributionParam[4] = mpvLandau; | |
334 | distributionParam[5] = sigmaLandau; | |
335 | ||
336 | return distributionParam; | |
dc293ae9 | 337 | } |
4cfd8ef9 | 338 | |
339 | //_______________________________________________________ | |
dc293ae9 | 340 | Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params) |
341 | { | |
4cfd8ef9 | 342 | // |
343 | // Compute a polynomial for a given value of 'x' | |
344 | // with the array of parameters passed as the second arg | |
345 | // | |
46363b15 | 346 | |
4cfd8ef9 | 347 | Double_t y; |
348 | y = params[0]; | |
349 | y += params[1] * x; | |
350 | y += params[2] * x * x; | |
351 | y += params[3] * x * x * x; | |
352 | y += params[4] * x * x * x * x; | |
353 | y += params[5] * x * x * x * x * x; | |
354 | ||
355 | return y; | |
dc293ae9 | 356 | } |