]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaElectron.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Major refactoring of the class, many parts of the code moved to different methods...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaElectron.cxx
CommitLineData
d9105d92 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 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/* $Id: AliAnaElectron.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18//
19// Class for the photon identification.
20// Clusters from calorimeters are identified as photons
21// and kept in the AOD. Few histograms produced.
22// Copy of AliAnaPhoton just add electron id.
23//
24// -- Author: Gustavo Conesa (LPSC-IN2P3-CRNS)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TH2F.h>
30#include <TH3D.h>
31#include <TClonesArray.h>
32#include <TObjString.h>
33//#include <Riostream.h>
34#include "TParticle.h"
35#include "TDatabasePDG.h"
36#include "AliVTrack.h"
37
38// --- Analysis system ---
39#include "AliAnaElectron.h"
40#include "AliCaloTrackReader.h"
41#include "AliStack.h"
42#include "AliCaloPID.h"
43#include "AliMCAnalysisUtils.h"
44#include "AliFiducialCut.h"
45#include "AliVCluster.h"
46#include "AliAODMCParticle.h"
47#include "AliMixedEvent.h"
48
49
50ClassImp(AliAnaElectron)
51
52//____________________________________________________________________________
53AliAnaElectron::AliAnaElectron() :
54 AliAnaPartCorrBaseClass(), fCalorimeter(""),
55 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
56 fTimeCutMin(-1), fTimeCutMax(999999),
57 fNCellsCut(0), fFillSSHistograms(kFALSE), fNOriginHistograms(8),
58 fdEdxMin(0.), fdEdxMax (200.),
59 fEOverPMin(0), fEOverPMax (2),
60 // Histograms
61 fhdEdxvsE(0), fhdEdxvsP(0),
62 fhEOverPvsE(0), fhEOverPvsP(0),
63 // MC histograms
64 // Electron SS MC histograms
65 fhMCElectronELambda0NoOverlap(0),
66 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
67
68 //Embedding
69 fhEmbeddedSignalFractionEnergy(0),
70 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
71 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
72{
73 //default ctor
74 for(Int_t index = 0; index < 2; index++){
75
76 fhNCellsE [index] = 0;
77 fhMaxCellDiffClusterE[index] = 0;
78 fhE [index] = 0;
79 fhPhi [index] = 0;
80 fhEta [index] = 0;
81 fhEtaPhi [index] = 0;
82 fhEtaPhi05[index] = 0;
83
84 // Shower shape histograms
85 fhDispE [index] = 0;
86 fhLam0E [index] = 0;
87 fhLam1E [index] = 0;
88 fhDispETRD[index] = 0;
89 fhLam0ETRD[index] = 0;
90 fhLam1ETRD[index] = 0;
91 fhNCellsLam0LowE [index] = 0;
92 fhNCellsLam0HighE[index] = 0;
93 fhEtaLam0LowE [index] = 0;
94 fhPhiLam0LowE [index] = 0;
95 fhEtaLam0HighE [index] = 0;
96 fhPhiLam0HighE [index] = 0;
97
98 for(Int_t i = 0; i < 10; i++){
99 fhMCPt [index][i] = 0;
100 fhMCE [index][i] = 0;
101 fhMCPhi [index][i] = 0;
102 fhMCEta [index][i] = 0;
103 fhMCDeltaE [index][i] = 0;
104 fhMC2E [index][i] = 0;
105 }
106
107 for(Int_t i = 0; i < 6; i++){
108 fhMCELambda0[index][i]= 0;
109 }
110
111 }
112 //Initialize parameters
113 InitParameters();
114
115}
116
117//__________________________________________________________________
118Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
119{
120 //Select clusters if they pass different cuts
121 if(GetDebug() > 2)
122 printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
123 GetReader()->GetEventNumber(),
124 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
125
126 //.......................................
127 //If too small or big energy, skip it
128 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
129 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
130
131 //.......................................
132 // TOF cut, BE CAREFUL WITH THIS CUT
133 Double_t tof = calo->GetTOF()*1e9;
134 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
135 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
136
137 //.......................................
138 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
139 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
140
141 //.......................................
142 //Check acceptance selection
143 if(IsFiducialCutOn()){
144 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
145 if(! in ) return kFALSE ;
146 }
147 if(GetDebug() > 2) printf("Fiducial cut passed \n");
148
149 //.......................................
150 //Skip not matched clusters with tracks
151 if(!IsTrackMatched(calo)) {
152 if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
153 return kFALSE ;
154 }
155 else if(GetDebug() > 2) printf(" Track-matching cut passed \n");
156
157 //.......................................
158 //Check Distance to Bad channel, set bit.
159 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
160 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
161 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
162 return kFALSE ;
163 }
164 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
165 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
166
167 if(GetDebug() > 0)
168 printf("AliAnaElectron::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
169 GetReader()->GetEventNumber(),
170 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
171
172 //All checks passed, cluster selected
173 return kTRUE;
174
175}
176
177//__________________________________________________________________
178void AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag, const Int_t pidTag){
179
180 //Fill cluster Shower Shape histograms
181
182 if(!fFillSSHistograms || GetMixedEvent()) return;
183
184 Int_t pidIndex = 0;// Electron
185 if (pidTag == AliCaloPID::kElectron) pidIndex = 0;
186 else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
187 else return;
188
189 Float_t energy = cluster->E();
190 Int_t ncells = cluster->GetNCells();
191 Float_t lambda0 = cluster->GetM02();
192 Float_t lambda1 = cluster->GetM20();
193 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
194
195 TLorentzVector mom;
196 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
197 cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
198 else{
199 Double_t vertex[]={0,0,0};
200 cluster->GetMomentum(mom,vertex) ;
201 }
202
203 Float_t eta = mom.Eta();
204 Float_t phi = mom.Phi();
205 if(phi < 0) phi+=TMath::TwoPi();
206
207 fhLam0E[pidIndex] ->Fill(energy,lambda0);
208 fhLam1E[pidIndex] ->Fill(energy,lambda1);
209 fhDispE[pidIndex] ->Fill(energy,disp);
210
211 if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
212 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
213 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
214 fhDispETRD[pidIndex]->Fill(energy,disp);
215 }
216
217 if(energy < 2){
218 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
219 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
220 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
221 }
222 else {
223 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
224 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
225 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
226 }
227
228 if(IsDataMC()){
229
230 AliVCaloCells* cells = 0;
231 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
232 else cells = GetPHOSCells();
233
234 //Fill histograms to check shape of embedded clusters
235 Float_t fraction = 0;
236 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
237
238 Float_t clusterE = 0; // recalculate in case corrections applied.
239 Float_t cellE = 0;
240 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
241 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
242 clusterE+=cellE;
243 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
244 }
245
246 //Fraction of total energy due to the embedded signal
247 fraction/=clusterE;
248
249 if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
250
251 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
252
253 } // embedded fraction
254
255 // Get the fraction of the cluster energy that carries the cell with highest energy
256 Int_t absID =-1 ;
257 Float_t maxCellFraction = 0.;
258
259 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
260
261 // Check the origin and fill histograms
262 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
263 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
264 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
265 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
266 fhMCELambda0[pidIndex][mcssPhoton] ->Fill(energy, lambda0);
267
268
269 }//photon no conversion
270 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
271 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))){
272 fhMCELambda0[pidIndex][mcssElectron] ->Fill(energy, lambda0);
273
274 if(!GetReader()->IsEmbeddedClusterSelectionOn()){
275 //Check particle overlaps in cluster
276
277 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
278 Int_t ancPDG = 0, ancStatus = -1;
279 TLorentzVector momentum; TVector3 prodVertex;
280 Int_t ancLabel = 0;
281 Int_t noverlaps = 1;
282 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
283 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
284 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
285 }
286
287 if(noverlaps == 1){
288 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
289 }
290 else if(noverlaps == 2){
291 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
292 }
293 else if(noverlaps > 2){
294 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
295 }
296 else {
297 printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
298 }
299 }//No embedding
300
301 //Fill histograms to check shape of embedded clusters
302 if(GetReader()->IsEmbeddedClusterSelectionOn()){
303
304 if (fraction > 0.9)
305 {
306 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
307 }
308 else if(fraction > 0.5)
309 {
310 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
311 }
312 else if(fraction > 0.1)
313 {
314 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
315 }
316 else
317 {
318 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
319 }
320 } // embedded
321 }//electron
322 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
323 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
324 fhMCELambda0[pidIndex][mcssConversion] ->Fill(energy, lambda0);
325 }//conversion photon
326 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ){
327 fhMCELambda0[pidIndex][mcssPi0] ->Fill(energy, lambda0);
328 }//pi0
329 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ){
330 fhMCELambda0[pidIndex][mcssEta] ->Fill(energy, lambda0);
331
332 }//eta
333 else {
334 fhMCELambda0[pidIndex][mcssOther] ->Fill(energy, lambda0);
335 }//other particles
336
337 }//MC data
338
339}
340
341//________________________________________________________________________
342TObjString * AliAnaElectron::GetAnalysisCuts()
343{
344 //Save parameters used for analysis
345 TString parList ; //this will be list of parameters used for this analysis.
346 const Int_t buffersize = 255;
347 char onePar[buffersize] ;
348
349 snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
350 parList+=onePar ;
351 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
352 parList+=onePar ;
353 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
354 parList+=onePar ;
355 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f \n",fEOverPMin, fEOverPMax) ;
356 parList+=onePar ;
357 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
358 parList+=onePar ;
359 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
360 parList+=onePar ;
361 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
362 parList+=onePar ;
363
364 //Get parameters set in base class.
365 parList += GetBaseParametersList() ;
366
367 //Get parameters set in PID class.
368 parList += GetCaloPID()->GetPIDParametersList() ;
369
370 //Get parameters set in FiducialCut class (not available yet)
371 //parlist += GetFidCut()->GetFidCutParametersList()
372
373 return new TObjString(parList) ;
374}
375
376//________________________________________________________________________
377TList * AliAnaElectron::GetCreateOutputObjects()
378{
379 // Create histograms to be saved in output file and
380 // store them in outputContainer
381 TList * outputContainer = new TList() ;
382 outputContainer->SetName("ElectronHistos") ;
383
384 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
385 Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
386 Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
387 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
388 Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
389 Int_t ndedxbins = GetHistodEdxBins(); Float_t dedxmax = GetHistodEdxMax(); Float_t dedxmin = GetHistodEdxMin();
390 Int_t nPoverEbins = GetHistoPOverEBins(); Float_t pOverEmax = GetHistoPOverEMax(); Float_t pOverEmin = GetHistoPOverEMin();
391
392 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
393 fhdEdxvsE->SetXTitle("E (GeV)");
394 fhdEdxvsE->SetYTitle("<dE/dx>");
395 outputContainer->Add(fhdEdxvsE);
396
397 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
398 fhdEdxvsP->SetXTitle("P (GeV/c)");
399 fhdEdxvsP->SetYTitle("<dE/dx>");
400 outputContainer->Add(fhdEdxvsP);
401
402 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
403 fhEOverPvsE->SetXTitle("E (GeV)");
404 fhEOverPvsE->SetYTitle("E/p");
405 outputContainer->Add(fhEOverPvsE);
406
407 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
408 fhEOverPvsP->SetXTitle("P (GeV/c)");
409 fhEOverPvsP->SetYTitle("E/p");
410 outputContainer->Add(fhEOverPvsP);
411
412
413 TString pidParticle[] = {"Electron","ChargedHadron"} ;
414
415 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++){
416
417 //Shower shape
418 if(fFillSSHistograms){
419 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
420 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
421 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
422 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
423 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
424 outputContainer->Add(fhLam0E[pidIndex]);
425
426 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
427 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
428 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
429 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
430 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
431 outputContainer->Add(fhLam1E[pidIndex]);
432
433 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
434 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
435 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
436 fhDispE[pidIndex]->SetYTitle("D^{2}");
437 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
438 outputContainer->Add(fhDispE[pidIndex]);
439
440 if(fCalorimeter == "EMCAL"){
441 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
442 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
443 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
444 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
445 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
446 outputContainer->Add(fhLam0ETRD[pidIndex]);
447
448 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
449 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
450 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
451 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
452 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
453 outputContainer->Add(fhLam1ETRD[pidIndex]);
454
455 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
456 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
457 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
458 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
459 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
460 outputContainer->Add(fhDispETRD[pidIndex]);
461 }
462
463 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
464 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
465 nbins,nmin, nmax, ssbins,ssmin,ssmax);
466 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
467 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
468 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
469
470 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
471 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
472 nbins,nmin, nmax, ssbins,ssmin,ssmax);
473 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
474 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
475 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
476
477
478 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
479 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
480 netabins,etamin,etamax, ssbins,ssmin,ssmax);
481 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
482 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
483 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
484
485 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
486 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
487 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
488 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
489 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
490 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
491
492 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
493 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
494 netabins,etamin,etamax, ssbins,ssmin,ssmax);
495 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
496 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
497 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
498
499 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
500 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
501 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
502 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
503 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
504 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
505
506 } // Shower shape
507
508
509 if(IsDataMC()){
510
511 if(fFillSSHistograms){
512
513 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
514
515 TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
516
517 for(Int_t i = 0; i < 6; i++){
518
519 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
520 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
521 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
522 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
523 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
524 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
525
526 }// loop
527 }
528 }
529
530 if(IsCaloPIDOn() && pidIndex > 0) continue;
531
532 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
533 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
534 nptbins,ptmin,ptmax, nbins,nmin,nmax);
535 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
536 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
537 outputContainer->Add(fhNCellsE[pidIndex]);
538
539 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
540 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
541 nptbins,ptmin,ptmax, 500,0,1.);
542 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
543 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
544 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
545
546 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
547 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
548 nptbins,ptmin,ptmax);
549 fhE[pidIndex]->SetYTitle("N");
550 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
551 outputContainer->Add(fhE[pidIndex]) ;
552
553 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
554 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
555 nptbins,ptmin,ptmax);
556 fhPt[pidIndex]->SetYTitle("N");
557 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
558 outputContainer->Add(fhPt[pidIndex]) ;
559
560 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
561 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
562 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
563 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
564 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
565 outputContainer->Add(fhPhi[pidIndex]) ;
566
567 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
568 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
569 nptbins,ptmin,ptmax,netabins,etamin,etamax);
570 fhEta[pidIndex]->SetYTitle("#eta");
571 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
572 outputContainer->Add(fhEta[pidIndex]) ;
573
574 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
575 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
576 netabins,etamin,etamax,nphibins,phimin,phimax);
577 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
578 fhEtaPhi[pidIndex]->SetXTitle("#eta");
579 outputContainer->Add(fhEtaPhi[pidIndex]) ;
580 if(GetMinPt() < 0.5){
581 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
582 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
583 netabins,etamin,etamax,nphibins,phimin,phimax);
584 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
585 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
586 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
587 }
588
589
590 if(IsDataMC()){
591
592
593 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
594 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
595
596 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
597 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
598
599 for(Int_t i = 0; i < fNOriginHistograms; i++){
600
601 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
602 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
603 nptbins,ptmin,ptmax);
604 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
605 outputContainer->Add(fhMCE[pidIndex][i]) ;
606
607 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
608 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
609 nptbins,ptmin,ptmax);
610 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
611 outputContainer->Add(fhMCPt[pidIndex][i]) ;
612
613 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
614 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
615 nptbins,ptmin,ptmax,netabins,etamin,etamax);
616 fhMCEta[pidIndex][i]->SetYTitle("#eta");
617 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
618 outputContainer->Add(fhMCEta[pidIndex][i]) ;
619
620 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
621 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
622 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
623 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
624 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
625 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
626
627
628 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
629 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
630 nptbins,ptmin,ptmax, 200,-50,50);
631 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
632 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
633
634 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
635 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
636 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
637 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
638 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
639 outputContainer->Add(fhMC2E[pidIndex][i]);
640
641 }
642 } // MC
643
644 }// pid Index
645
646
647 if(fFillSSHistograms){
648
649 if(IsDataMC()){
650
651 if(!GetReader()->IsEmbeddedClusterSelectionOn())
652 {
653 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
654 "cluster from Electron : E vs #lambda_{0}^{2}",
655 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
656 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
657 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
658 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
659
660 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
661 "cluster from Electron : E vs #lambda_{0}^{2}",
662 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
663 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
664 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
665 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
666
667 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
668 "cluster from Electron : E vs #lambda_{0}^{2}",
669 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
670 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
671 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
672 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
673
674 } //No embedding
675
676 //Fill histograms to check shape of embedded clusters
677 if(GetReader()->IsEmbeddedClusterSelectionOn())
678 {
679
680 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
681 "Energy Fraction of embedded signal versus cluster energy",
682 nptbins,ptmin,ptmax,100,0.,1.);
683 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
684 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
685 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
686
687 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
688 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
689 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
690 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
691 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
692 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
693
694 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
695 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
696 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
697 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
698 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
699 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
700
701 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
702 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
703 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
704 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
705 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
706 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
707
708 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
709 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
710 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
711 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
712 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
713 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
714
715
716 }// embedded histograms
717
718 }//Histos with MC
719
720 }// Fill SS MC histograms
721
722
723 //Store calo PID histograms
724 if(IsCaloPIDOn()){
725 TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
726 for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
727 outputContainer->Add(caloPIDHistos->At(i)) ;
728 }
729 delete caloPIDHistos;
730 }
731 return outputContainer ;
732
733}
734
735//____________________________________________________________________________
736void AliAnaElectron::Init()
737{
738
739 //Init
740 //Do some checks
741 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
742 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
743 abort();
744 }
745 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
746 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
747 abort();
748 }
749
750}
751
752//____________________________________________________________________________
753void AliAnaElectron::InitParameters()
754{
755
756 //Initialize the parameters of the analysis.
757 AddToHistogramsName("AnaElectron_");
758
759 fCalorimeter = "EMCAL" ;
760 fMinDist = 2.;
761 fMinDist2 = 4.;
762 fMinDist3 = 5.;
763
764 fTimeCutMin = -1;
765 fTimeCutMax = 9999999;
766 fNCellsCut = 0;
767
768 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
769 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
770
771 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
772 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
773
774}
775
776//__________________________________________________________________
777void AliAnaElectron::MakeAnalysisFillAOD()
778{
779 //Do photon analysis and fill aods
780
781 //Get the vertex
782 Double_t v[3] = {0,0,0}; //vertex ;
783 GetReader()->GetVertex(v);
784
785 //Select the Calorimeter of the photon
786 TObjArray * pl = 0x0;
787 if(fCalorimeter == "PHOS")
788 pl = GetPHOSClusters();
789 else if (fCalorimeter == "EMCAL")
790 pl = GetEMCALClusters();
791
792 if(!pl) {
793 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
794 return;
795 }
796
797 //Init arrays, variables, get number of clusters
798 TLorentzVector mom, mom2 ;
799 Int_t nCaloClusters = pl->GetEntriesFast();
800 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
801
802 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
803
804 //----------------------------------------------------
805 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
806 //----------------------------------------------------
807 // Loop on clusters
808 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
809
810 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
811 //printf("calo %d, %f\n",icalo,calo->E());
812
813 //Get the index where the cluster comes, to retrieve the corresponding vertex
814 Int_t evtIndex = 0 ;
815 if (GetMixedEvent()) {
816 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
817 //Get the vertex and check it is not too large in z
818 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
819 }
820
821 //Cluster selection, not charged, with photon id and in fiducial cut
822 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
823 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
824 else{
825 Double_t vertex[]={0,0,0};
826 calo->GetMomentum(mom,vertex) ;
827 }
828
829 //--------------------------------------
830 // Cluster selection
831 //--------------------------------------
832 if(!ClusterSelected(calo,mom)) continue;
833
834 //----------------------------
835 //Create AOD for analysis
836 //----------------------------
837 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
838
839 //...............................................
840 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
841 Int_t label = calo->GetLabel();
842 aodph.SetLabel(label);
843 aodph.SetCaloLabel(calo->GetID(),-1);
844 aodph.SetDetector(fCalorimeter);
845 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
846
847 //...............................................
848 //Set bad channel distance bit
849 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
850 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
851 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
852 else aodph.SetDistToBad(0) ;
853 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
854
855 //--------------------------------------------------------------------------------------
856 //Play with the MC stack if available
857 //--------------------------------------------------------------------------------------
858
859 //Check origin of the candidates
860 if(IsDataMC()){
861 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
862
863 if(GetDebug() > 0)
864 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
865 }//Work with stack also
866
867
868 //-------------------------------------
869 //PID selection via dedx
870 //-------------------------------------
871 Int_t pid = AliCaloPID::kChargedHadron;
872 AliVTrack *track = 0;
873 if(!strcmp("AliESDCaloCluster",Form("%s",calo->ClassName()))){
874 Int_t iESDtrack = calo->GetTrackMatchedIndex();
875 if(iESDtrack<0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Wrong track index\n");
876 AliVEvent * event = GetReader()->GetInputEvent();
877 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
878 }
879 else {
880 track = dynamic_cast<AliVTrack*>(calo->GetTrackMatched(0));
881 }
882
883 if(!track) {
884 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
885 continue;
886 }
887
e416be7d 888 Float_t dEdx = track->GetTPCsignal();
d9105d92 889 fhdEdxvsE->Fill(calo->E(), dEdx);
890 fhdEdxvsP->Fill(track->P(),dEdx);
891
d9105d92 892 if( dEdx < fdEdxMax && dEdx > fdEdxMin) {
893
e416be7d 894 Float_t eOverp = calo->E()/track->P();
895 fhEOverPvsE->Fill(calo->E(), eOverp);
896 fhEOverPvsP->Fill(track->P(), eOverp);
d9105d92 897
898 if( eOverp < fEOverPMax && eOverp > fEOverPMin) {
899
900 pid = AliCaloPID::kElectron;
901
902 } //E/p
903
904 }// dEdx
905
906 aodph.SetIdentifiedParticleType(pid);
907
908 Int_t pidIndex = 0;// Electron
909 if (pid == AliCaloPID::kElectron) pidIndex = 0;
910 else if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
911 else continue ;
912
913 //--------------------------------------------------------------------------------------
914 //Fill some shower shape histograms
915 //--------------------------------------------------------------------------------------
916
917 FillShowerShapeHistograms(calo,aodph.GetTag(),pid);
918
919 //-------------------------------------
920 //PID Shower Shape selection or bit setting
921 //-------------------------------------
922 // Data, PID check on
923 if(IsCaloPIDOn()){
924 //Get most probable PID, 2 options check PID weights
925 //or redo PID, recommended option for EMCal.
926 if(!IsCaloPIDRecalculationOn()){
927 if(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E())!=AliCaloPID::kPhoton) continue;
928 }
929 else{
930 if(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E())!=AliCaloPID::kPhoton) continue;
931 }
932
933 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
934
935 }
936
937 // //...............................................
938 // // Data, PID check off
939 // else{
940 // //Set PID bits for later selection (AliAnaPi0 for example)
941 // //GetPDG already called in SetPIDBits.
942 // GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
943 // if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - PID Bits set \n");
944 // }
945
946 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
947
948
949 //FIXME, this to MakeAnalysisFillHistograms ...
950 Int_t absID = 0;
951 Float_t maxCellFraction = 0;
952 AliVCaloCells* cells = 0;
953
954 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
955 else cells = GetPHOSCells();
956
957 absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
958
959 fhMaxCellDiffClusterE[pidIndex]->Fill(aodph.E(),maxCellFraction);
960 fhNCellsE[pidIndex] ->Fill(aodph.E(),calo->GetNCells());
961
962 //Add AOD with photon object to aod branch
963 AddAODParticle(aodph);
964
965 }//loop
966
967 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
968
969}
970
971//__________________________________________________________________
972void AliAnaElectron::MakeAnalysisFillHistograms()
973{
974 //Fill histograms
975
976 //-------------------------------------------------------------------
977 // Access MC information in stack if requested, check that it exists.
978 AliStack * stack = 0x0;
979 TParticle * primary = 0x0;
980 TClonesArray * mcparticles = 0x0;
981 AliAODMCParticle * aodprimary = 0x0;
982
983 if(IsDataMC()){
984
985 if(GetReader()->ReadStack()){
986 stack = GetMCStack() ;
987 if(!stack) {
988 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
989 abort();
990 }
991
992 }
993 else if(GetReader()->ReadAODMCParticles()){
994
995 //Get the list of MC particles
996 mcparticles = GetReader()->GetAODMCParticles(0);
997 if(!mcparticles && GetDebug() > 0) {
998 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
999 }
1000 }
1001 }// is data and MC
1002
1003
1004 // Get vertex
1005 Double_t v[3] = {0,0,0}; //vertex ;
1006 GetReader()->GetVertex(v);
1007 //fhVertex->Fill(v[0],v[1],v[2]);
1008 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1009
1010 //----------------------------------
1011 //Loop on stored AOD photons
1012 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1013 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1014
1015 for(Int_t iaod = 0; iaod < naod ; iaod++){
1016 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1017 Int_t pdg = ph->GetIdentifiedParticleType();
1018
1019 Int_t pidIndex = 0;// Electron
1020 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1021 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1022 else continue ;
1023
1024 if(ph->GetDetector() != fCalorimeter) continue;
1025
1026 if(GetDebug() > 2)
1027 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1028
1029 //................................
1030 //Fill photon histograms
1031 Float_t ptcluster = ph->Pt();
1032 Float_t phicluster = ph->Phi();
1033 Float_t etacluster = ph->Eta();
1034 Float_t ecluster = ph->E();
1035
1036 fhE[pidIndex] ->Fill(ecluster);
1037 fhPt[pidIndex] ->Fill(ptcluster);
1038 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1039 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1040 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1041 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1042
1043 //.......................................
1044 //Play with the MC data if available
1045 if(IsDataMC()){
1046
1047 //....................................................................
1048 // Access MC information in stack if requested, check that it exists.
1049 Int_t label =ph->GetLabel();
1050 if(label < 0) {
1051 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1052 continue;
1053 }
1054
1055 Float_t eprim = 0;
1056 Float_t ptprim = 0;
1057 if(GetReader()->ReadStack()){
1058
1059 if(label >= stack->GetNtrack()) {
1060 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1061 continue ;
1062 }
1063
1064 primary = stack->Particle(label);
1065 if(!primary){
1066 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1067 continue;
1068 }
1069
1070 eprim = primary->Energy();
1071 ptprim = primary->Pt();
1072
1073 }
1074 else if(GetReader()->ReadAODMCParticles()){
1075 //Check which is the input
1076 if(ph->GetInputFileIndex() == 0){
1077 if(!mcparticles) continue;
1078 if(label >= mcparticles->GetEntriesFast()) {
1079 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1080 label, mcparticles->GetEntriesFast());
1081 continue ;
1082 }
1083 //Get the particle
1084 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1085
1086 }
1087
1088 if(!aodprimary){
1089 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1090 continue;
1091 }
1092
1093 eprim = aodprimary->E();
1094 ptprim = aodprimary->Pt();
1095
1096 }
1097
1098 Int_t tag =ph->GetTag();
1099
1100 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][mcPhoton])
1101 {
1102 fhMCE [pidIndex][mcPhoton] ->Fill(ecluster);
1103 fhMCPt [pidIndex][mcPhoton] ->Fill(ptcluster);
1104 fhMCPhi[pidIndex][mcPhoton] ->Fill(ecluster,phicluster);
1105 fhMCEta[pidIndex][mcPhoton] ->Fill(ecluster,etacluster);
1106
1107 fhMC2E[pidIndex][mcPhoton] ->Fill(ecluster, eprim);
1108 fhMCDeltaE[pidIndex][mcPhoton] ->Fill(ecluster,eprim-ecluster);
1109
1110 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][mcConversion])
1111 {
1112 fhMCE [pidIndex][mcConversion] ->Fill(ecluster);
1113 fhMCPt [pidIndex][mcConversion] ->Fill(ptcluster);
1114 fhMCPhi[pidIndex][mcConversion] ->Fill(ecluster,phicluster);
1115 fhMCEta[pidIndex][mcConversion] ->Fill(ecluster,etacluster);
1116
1117 fhMC2E[pidIndex][mcConversion] ->Fill(ecluster, eprim);
1118 fhMCDeltaE[pidIndex][mcConversion] ->Fill(ecluster,eprim-ecluster);
1119
1120 }
1121 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1122 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][mcPi0Decay])
1123 {
1124 fhMCE [pidIndex][mcPi0Decay] ->Fill(ecluster);
1125 fhMCPt [pidIndex][mcPi0Decay] ->Fill(ptcluster);
1126 fhMCPhi[pidIndex][mcPi0Decay] ->Fill(ecluster,phicluster);
1127 fhMCEta[pidIndex][mcPi0Decay] ->Fill(ecluster,etacluster);
1128
1129 fhMC2E[pidIndex][mcPi0Decay] ->Fill(ecluster, eprim);
1130 fhMCDeltaE[pidIndex][mcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1131 }
1132 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1133 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][mcOtherDecay])
1134 {
1135 fhMCE [pidIndex][mcOtherDecay] ->Fill(ecluster);
1136 fhMCPt [pidIndex][mcOtherDecay] ->Fill(ptcluster);
1137 fhMCPhi[pidIndex][mcOtherDecay] ->Fill(ecluster,phicluster);
1138 fhMCEta[pidIndex][mcOtherDecay] ->Fill(ecluster,etacluster);
1139
1140 fhMC2E[pidIndex][mcOtherDecay] ->Fill(ecluster, eprim);
1141 fhMCDeltaE[pidIndex][mcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1142 }
1143 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][mcPi0])
1144 {
1145 fhMCE [pidIndex][mcPi0] ->Fill(ecluster);
1146 fhMCPt [pidIndex][mcPi0] ->Fill(ptcluster);
1147 fhMCPhi[pidIndex][mcPi0] ->Fill(ecluster,phicluster);
1148 fhMCEta[pidIndex][mcPi0] ->Fill(ecluster,etacluster);
1149
1150 fhMC2E[pidIndex][mcPi0] ->Fill(ecluster, eprim);
1151 fhMCDeltaE[pidIndex][mcPi0] ->Fill(ecluster,eprim-ecluster);
1152
1153 }
1154 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][mcEta])
1155 {
1156 fhMCE [pidIndex][mcEta] ->Fill(ecluster);
1157 fhMCPt [pidIndex][mcEta] ->Fill(ptcluster);
1158 fhMCPhi[pidIndex][mcEta] ->Fill(ecluster,phicluster);
1159 fhMCEta[pidIndex][mcEta] ->Fill(ecluster,etacluster);
1160
1161 fhMC2E[pidIndex][mcEta] ->Fill(ecluster, eprim);
1162 fhMCDeltaE[pidIndex][mcEta] ->Fill(ecluster,eprim-ecluster);
1163
1164 }
1165 }
1166 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][mcAntiNeutron])
1167 {
1168 fhMCE [pidIndex][mcAntiNeutron] ->Fill(ecluster);
1169 fhMCPt [pidIndex][mcAntiNeutron] ->Fill(ptcluster);
1170 fhMCPhi[pidIndex][mcAntiNeutron] ->Fill(ecluster,phicluster);
1171 fhMCEta[pidIndex][mcAntiNeutron] ->Fill(ecluster,etacluster);
1172
1173 fhMC2E[pidIndex][mcAntiNeutron] ->Fill(ecluster, eprim);
1174 fhMCDeltaE[pidIndex][mcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1175
1176 }
1177 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][mcAntiProton])
1178 {
1179 fhMCE [pidIndex][mcAntiProton] ->Fill(ecluster);
1180 fhMCPt [pidIndex][mcAntiProton] ->Fill(ptcluster);
1181 fhMCPhi[pidIndex][mcAntiProton] ->Fill(ecluster,phicluster);
1182 fhMCEta[pidIndex][mcAntiProton] ->Fill(ecluster,etacluster);
1183
1184 fhMC2E[pidIndex][mcAntiProton] ->Fill(ecluster, eprim);
1185 fhMCDeltaE[pidIndex][mcAntiProton] ->Fill(ecluster,eprim-ecluster);
1186
1187 }
1188 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][mcElectron])
1189 {
1190 fhMCE [pidIndex][mcElectron] ->Fill(ecluster);
1191 fhMCPt [pidIndex][mcElectron] ->Fill(ptcluster);
1192 fhMCPhi[pidIndex][mcElectron] ->Fill(ecluster,phicluster);
1193 fhMCEta[pidIndex][mcElectron] ->Fill(ecluster,etacluster);
1194
1195 fhMC2E[pidIndex][mcElectron] ->Fill(ecluster, eprim);
1196 fhMCDeltaE[pidIndex][mcElectron] ->Fill(ecluster,eprim-ecluster);
1197
1198 }
1199 else if( fhMCE[pidIndex][mcOther]){
1200 fhMCE [pidIndex][mcOther] ->Fill(ecluster);
1201 fhMCPt [pidIndex][mcOther] ->Fill(ptcluster);
1202 fhMCPhi[pidIndex][mcOther] ->Fill(ecluster,phicluster);
1203 fhMCEta[pidIndex][mcOther] ->Fill(ecluster,etacluster);
1204
1205 fhMC2E[pidIndex][mcOther] ->Fill(ecluster, eprim);
1206 fhMCDeltaE[pidIndex][mcOther] ->Fill(ecluster,eprim-ecluster);
1207
1208 }
1209
1210 }//Histograms with MC
1211
1212 }// aod loop
1213
1214}
1215
1216
1217//__________________________________________________________________
1218void AliAnaElectron::Print(const Option_t * opt) const
1219{
1220 //Print some relevant parameters set for the analysis
1221
1222 if(! opt)
1223 return;
1224
1225 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1226 AliAnaPartCorrBaseClass::Print(" ");
1227
1228 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1229 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1230 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1231 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1232 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1233 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1234 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1235 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1236 printf(" \n") ;
1237
1238}