MCEvent and MCEventHandler can be used as containers. This allows us to compose MC...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPID.cxx
CommitLineData
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
117ClassImp(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 157void 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
3a2a23e1 170 Int_t nClusters = esd->GetNumberOfCaloClusters();
171 Int_t firstCluster = 0;
dc293ae9 172 Double_t energy, lambda0;
173 for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
8ba062b1 174
dc293ae9 175 AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
3a2a23e1 176 if (!clust->IsEMCAL()) continue ;
7592dfc4 177 energy = clust->E();
dc293ae9 178 lambda0 = clust->GetM02();
179 // verify cluster type
180 Int_t clusterType= clust->GetClusterType();
8ada0ffe 181 if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0 && energy < 1000) {
4ce73a00 182
183
dc293ae9 184 // reject clusters with lambda0 = 0
4ce73a00 185
186
dc293ae9 187 ComputePID(energy, lambda0);
4ce73a00 188
189
dc293ae9 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 }
4cfd8ef9 212
dc293ae9 213 if(fReconstructor) // In case it is called during reconstruction.
214 clust->SetPid(fPIDFinal);
215 } // end if (clusterType...)
216 } // end for (iCluster...)
217}
4cfd8ef9 218
219//__________________________________________________________
dc293ae9 220void 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//
4ce73a00 226
227if (energy<5){energy =6;}
228
229
4cfd8ef9 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;
dc293ae9 276}
4cfd8ef9 277
278//________________________________________________________
dc293ae9 279TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
280{
4cfd8ef9 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]);
46363b15 296
297
4cfd8ef9 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]);
46363b15 307
4cfd8ef9 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]);
46363b15 316
4cfd8ef9 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]);
46363b15 326
4cfd8ef9 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;
dc293ae9 338}
4cfd8ef9 339
340//_______________________________________________________
dc293ae9 341Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
342{
4cfd8ef9 343 //
344 // Compute a polynomial for a given value of 'x'
345 // with the array of parameters passed as the second arg
346 //
46363b15 347
4cfd8ef9 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;
dc293ae9 357}