1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
18 // Class for the photon identification.
19 // Clusters from calorimeters are identified as photons
20 // and kept in the AOD. Few histograms produced.
21 // Copy of AliAnaPhoton just add electron id.
23 // -- Author: Gustavo Conesa (LPSC-IN2P3-CRNS)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
30 #include <TClonesArray.h>
31 #include <TObjString.h>
32 //#include <Riostream.h>
33 #include "TParticle.h"
34 #include "TDatabasePDG.h"
35 #include "AliVTrack.h"
37 // --- Analysis system ---
38 #include "AliAnaElectron.h"
39 #include "AliCaloTrackReader.h"
41 #include "AliCaloPID.h"
42 #include "AliMCAnalysisUtils.h"
43 #include "AliFiducialCut.h"
44 #include "AliVCluster.h"
45 #include "AliAODMCParticle.h"
46 #include "AliMixedEvent.h"
49 ClassImp(AliAnaElectron)
51 //________________________________
52 AliAnaElectron::AliAnaElectron() :
53 AliAnaCaloTrackCorrBaseClass(),
54 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
55 fTimeCutMin(-1), fTimeCutMax(999999),
56 fNCellsCut(0), fNLMCutMin(-1), fNLMCutMax(10),
57 fFillSSHistograms(kFALSE), fFillOnlySimpleSSHisto(1),
58 fFillWeightHistograms(kFALSE), fNOriginHistograms(8),
59 fdEdxMin(0.), fdEdxMax (200.),
60 fEOverPMin(0), fEOverPMax (2),
62 fMomentum(), fMomentumMC(), fProdVertex(),
64 fhdEdxvsE(0), fhdEdxvsP(0),
65 fhEOverPvsE(0), fhEOverPvsP(0),
66 fhdEdxvsECutM02(0), fhdEdxvsPCutM02(0),
67 fhEOverPvsECutM02(0), fhEOverPvsPCutM02(0),
68 fhdEdxvsECutEOverP(0), fhdEdxvsPCutEOverP(0),
69 fhEOverPvsECutM02CutdEdx(0), fhEOverPvsPCutM02CutdEdx(0),
71 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
72 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
74 // Electron SS MC histograms
75 fhMCElectronELambda0NoOverlap(0),
76 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
79 fhEmbeddedSignalFractionEnergy(0),
80 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
81 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
84 for(Int_t index = 0; index < 2; index++)
86 fhNCellsE [index] = 0;
89 fhMaxCellDiffClusterE[index] = 0;
95 fhEtaPhi05[index] = 0;
97 // Shower shape histograms
101 fhDispETRD[index] = 0;
102 fhLam0ETRD[index] = 0;
103 fhLam1ETRD[index] = 0;
104 fhNCellsLam0LowE [index] = 0;
105 fhNCellsLam0HighE[index] = 0;
106 fhEtaLam0LowE [index] = 0;
107 fhPhiLam0LowE [index] = 0;
108 fhEtaLam0HighE [index] = 0;
109 fhPhiLam0HighE [index] = 0;
111 fhDispEtaE [index] = 0;
112 fhDispPhiE [index] = 0;
113 fhSumEtaE [index] = 0;
114 fhSumPhiE [index] = 0;
115 fhSumEtaPhiE [index] = 0;
116 fhDispEtaPhiDiffE[index] = 0;
117 fhSphericityE [index] = 0;
119 for(Int_t i = 0; i < 10; i++)
121 fhMCPt [index][i] = 0;
122 fhMCE [index][i] = 0;
123 fhMCPhi [index][i] = 0;
124 fhMCEta [index][i] = 0;
125 fhMCDeltaE [index][i] = 0;
126 fhMC2E [index][i] = 0;
129 fhMCEOverPvsE [i] = 0;
130 fhMCEOverPvsP [i] = 0;
133 for(Int_t i = 0; i < 6; i++)
135 fhMCELambda0 [index][i] = 0;
136 fhMCEDispEta [index][i] = 0;
137 fhMCEDispPhi [index][i] = 0;
138 fhMCESumEtaPhi [index][i] = 0;
139 fhMCEDispEtaPhiDiff[index][i] = 0;
140 fhMCESphericity [index][i] = 0;
143 for(Int_t i = 0; i < 5; i++)
145 fhDispEtaDispPhiEBin[index][i] = 0 ;
150 for(Int_t i =0; i < 14; i++)
152 fhLambda0ForW0[i] = 0;
153 //fhLambda1ForW0[i] = 0;
156 //Initialize parameters
161 //____________________________________________________________________________
162 Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, Int_t nMaxima)
164 // Select clusters if they pass different cuts
166 AliDebug(1,Form("Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
167 GetReader()->GetEventNumber(),fMomentum.E(),fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
169 //.......................................
170 //If too small or big energy, skip it
171 if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) return kFALSE ;
172 AliDebug(2,Form("\t Cluster %d Pass E Cut",calo->GetID()));
174 //.......................................
175 // TOF cut, BE CAREFUL WITH THIS CUT
176 Double_t tof = calo->GetTOF()*1e9;
177 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
178 AliDebug(2,Form("\t Cluster %d Pass Time Cut",calo->GetID()));
180 //.......................................
181 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
182 AliDebug(2,Form("\t Cluster %d Pass NCell Cut",calo->GetID()));
184 //.......................................
185 //Check acceptance selection
186 if(IsFiducialCutOn())
188 Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
189 if(! in ) return kFALSE ;
191 AliDebug(2,"\t Fiducial cut passed");
193 //.......................................
194 //Skip not matched clusters with tracks
195 if(!IsTrackMatched(calo, GetReader()->GetInputEvent()))
197 AliDebug(1,"\t Reject non track-matched clusters");
200 else AliDebug(2,"\t Track-matching cut passed");
202 //...........................................
203 // skip clusters with too many maxima
204 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
205 AliDebug(2,Form("\t Cluster %d pass NLM %d of out of range",calo->GetID(), nMaxima));
207 //.......................................
208 //Check Distance to Bad channel, set bit.
209 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
210 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
211 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
214 else AliDebug(2,Form("\t Bad channel cut passed %4.2f > %2.2f",distBad, fMinDist));
215 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
217 AliDebug(1,Form("Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
218 GetReader()->GetEventNumber(),
219 fMomentum.E(), fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
221 //All checks passed, cluster selected
226 //______________________________________________________________________________________________
227 void AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
229 // Fill cluster Shower Shape histograms
231 if(!fFillSSHistograms || GetMixedEvent()) return;
233 Int_t pidIndex = 0;// Electron
234 if (pidTag == AliCaloPID::kElectron) pidIndex = 0;
235 else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
238 Float_t energy = cluster->E();
239 Int_t ncells = cluster->GetNCells();
240 Float_t lambda0 = cluster->GetM02();
241 Float_t lambda1 = cluster->GetM20();
242 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
244 Float_t l0 = 0., l1 = 0.;
245 Float_t dispp= 0., dEta = 0., dPhi = 0.;
246 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
248 Float_t eta = fMomentum.Eta();
249 Float_t phi = fMomentum.Phi();
250 if(phi < 0) phi+=TMath::TwoPi();
252 fhLam0E[pidIndex] ->Fill(energy,lambda0);
253 fhLam1E[pidIndex] ->Fill(energy,lambda1);
254 fhDispE[pidIndex] ->Fill(energy,disp);
256 if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
257 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
259 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
260 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
261 fhDispETRD[pidIndex]->Fill(energy,disp);
264 if(!fFillOnlySimpleSSHisto)
268 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
269 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
270 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
274 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
275 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
276 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
279 if(GetCalorimeter() == kEMCAL)
281 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
282 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
283 fhDispEtaE [pidIndex]-> Fill(energy,dEta);
284 fhDispPhiE [pidIndex]-> Fill(energy,dPhi);
285 fhSumEtaE [pidIndex]-> Fill(energy,sEta);
286 fhSumPhiE [pidIndex]-> Fill(energy,sPhi);
287 fhSumEtaPhiE [pidIndex]-> Fill(energy,sEtaPhi);
288 fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
289 if(dEta+dPhi>0)fhSphericityE [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
291 if (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
292 else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
293 else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
294 else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
295 else fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
302 AliVCaloCells* cells = 0;
303 if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
304 else cells = GetPHOSCells();
306 //Fill histograms to check shape of embedded clusters
307 Float_t fraction = 0;
308 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
310 Float_t clusterE = 0; // recalculate in case corrections applied.
312 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
313 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
315 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
318 //Fraction of total energy due to the embedded signal
321 AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
323 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
325 } // embedded fraction
327 // Check the origin and fill histograms
330 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
331 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
332 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
333 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
337 }//photon no conversion
338 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
339 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
341 index = kmcssElectron;
343 if(!GetReader()->IsEmbeddedClusterSelectionOn())
345 //Check particle overlaps in cluster
347 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
348 Int_t ancPDG = 0, ancStatus = -1;
351 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
353 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),
354 ancPDG,ancStatus,fMomentumMC,fProdVertex);
355 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
359 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
361 else if(noverlaps == 2){
362 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
364 else if(noverlaps > 2){
365 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
369 AliWarning(Form("N overlaps = %d for ancestor %d!!", noverlaps, ancLabel));
373 //Fill histograms to check shape of embedded clusters
374 if(GetReader()->IsEmbeddedClusterSelectionOn())
378 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
380 else if(fraction > 0.5)
382 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
384 else if(fraction > 0.1)
386 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
390 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
394 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
395 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
397 index = kmcssConversion;
399 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
403 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
412 fhMCELambda0[pidIndex][index] ->Fill(energy, lambda0);
414 if(GetCalorimeter() == kEMCAL && !fFillOnlySimpleSSHisto)
416 fhMCEDispEta [pidIndex][index]-> Fill(energy,dEta);
417 fhMCEDispPhi [pidIndex][index]-> Fill(energy,dPhi);
418 fhMCESumEtaPhi [pidIndex][index]-> Fill(energy,sEtaPhi);
419 fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
420 if(dEta+dPhi>0) fhMCESphericity[pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
427 //_____________________________________________
428 TObjString * AliAnaElectron::GetAnalysisCuts()
430 //Save parameters used for analysis
431 TString parList ; //this will be list of parameters used for this analysis.
432 const Int_t buffersize = 255;
433 char onePar[buffersize] ;
435 snprintf(onePar,buffersize,"--- AliAnaElectron ---: ") ;
437 snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
439 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f;",fdEdxMin,fdEdxMax) ;
441 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f;",fEOverPMin, fEOverPMax) ;
443 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
445 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
447 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
450 //Get parameters set in base class.
451 parList += GetBaseParametersList() ;
453 //Get parameters set in PID class.
454 parList += GetCaloPID()->GetPIDParametersList() ;
456 //Get parameters set in FiducialCut class (not available yet)
457 //parlist += GetFidCut()->GetFidCutParametersList()
459 return new TObjString(parList) ;
462 //_______________________________________________
463 TList * AliAnaElectron::GetCreateOutputObjects()
465 // Create histograms to be saved in output file and
466 // store them in outputContainer
467 TList * outputContainer = new TList() ;
468 outputContainer->SetName("ElectronHistos") ;
470 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
471 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
472 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
473 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
474 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
475 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
476 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
477 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
480 // MC labels, titles, for originator particles
481 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
482 TString pnamess[] = { "Photon","Hadron" ,"Pi0" ,"Eta" ,"Conversion" ,"Electron"} ;
483 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
484 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
486 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
487 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
490 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
491 fhdEdxvsE->SetXTitle("E (GeV)");
492 fhdEdxvsE->SetYTitle("<dE/dx>");
493 outputContainer->Add(fhdEdxvsE);
495 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
496 fhdEdxvsP->SetXTitle("P (GeV/c)");
497 fhdEdxvsP->SetYTitle("<dE/dx>");
498 outputContainer->Add(fhdEdxvsP);
500 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
501 fhEOverPvsE->SetXTitle("E (GeV)");
502 fhEOverPvsE->SetYTitle("E/p");
503 outputContainer->Add(fhEOverPvsE);
505 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
506 fhEOverPvsP->SetXTitle("P (GeV/c)");
507 fhEOverPvsP->SetYTitle("E/p");
508 outputContainer->Add(fhEOverPvsP);
511 fhdEdxvsECutM02 = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
512 fhdEdxvsECutM02->SetXTitle("E (GeV)");
513 fhdEdxvsECutM02->SetYTitle("<dE/dx>");
514 outputContainer->Add(fhdEdxvsECutM02);
516 fhdEdxvsPCutM02 = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
517 fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
518 fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
519 outputContainer->Add(fhdEdxvsPCutM02);
521 fhEOverPvsECutM02 = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
522 fhEOverPvsECutM02->SetXTitle("E (GeV)");
523 fhEOverPvsECutM02->SetYTitle("E/p");
524 outputContainer->Add(fhEOverPvsECutM02);
526 fhEOverPvsPCutM02 = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
527 fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
528 fhEOverPvsPCutM02->SetYTitle("E/p");
529 outputContainer->Add(fhEOverPvsPCutM02);
532 fhdEdxvsECutEOverP = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
533 fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
534 fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
535 outputContainer->Add(fhdEdxvsECutEOverP);
537 fhdEdxvsPCutEOverP = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
538 fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
539 fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
540 outputContainer->Add(fhdEdxvsPCutEOverP);
542 fhEOverPvsECutM02CutdEdx = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
543 fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
544 fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
545 outputContainer->Add(fhEOverPvsECutM02CutdEdx);
547 fhEOverPvsPCutM02CutdEdx = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
548 fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
549 fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
550 outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
554 for(Int_t i = 0; i < fNOriginHistograms; i++)
556 fhMCdEdxvsE[i] = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
557 Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
558 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
559 fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
560 fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
561 outputContainer->Add(fhMCdEdxvsE[i]) ;
563 fhMCdEdxvsP[i] = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
564 Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
565 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
566 fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
567 fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
568 outputContainer->Add(fhMCdEdxvsP[i]) ;
571 fhMCEOverPvsE[i] = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
572 Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
573 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
574 fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
575 fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
576 outputContainer->Add(fhMCEOverPvsE[i]) ;
578 fhMCEOverPvsP[i] = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
579 Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
580 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
581 fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
582 fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
583 outputContainer->Add(fhMCEOverPvsP[i]) ;
588 TString pidParticle[] = {"Electron","ChargedHadron"} ;
590 if(fFillWeightHistograms)
593 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
594 nptbins,ptmin,ptmax, 100,0,1.);
595 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
596 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
597 outputContainer->Add(fhECellClusterRatio);
599 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
600 nptbins,ptmin,ptmax, 100,-10,0);
601 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
602 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
603 outputContainer->Add(fhECellClusterLogRatio);
605 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
606 nptbins,ptmin,ptmax, 100,0,1.);
607 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
608 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
609 outputContainer->Add(fhEMaxCellClusterRatio);
611 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
612 nptbins,ptmin,ptmax, 100,-10,0);
613 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
614 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
615 outputContainer->Add(fhEMaxCellClusterLogRatio);
617 for(Int_t iw = 0; iw < 14; iw++)
619 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
620 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
621 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
622 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
623 outputContainer->Add(fhLambda0ForW0[iw]);
625 // fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
626 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
627 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
628 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
629 // outputContainer->Add(fhLambda1ForW0[iw]);
634 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
637 if(fFillSSHistograms)
639 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
640 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
641 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
642 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
643 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
644 outputContainer->Add(fhLam0E[pidIndex]);
646 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
647 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
648 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
649 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
650 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
651 outputContainer->Add(fhLam1E[pidIndex]);
653 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
654 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
655 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
656 fhDispE[pidIndex]->SetYTitle("D^{2}");
657 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
658 outputContainer->Add(fhDispE[pidIndex]);
660 if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >=0 )
662 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
663 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
664 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
665 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
666 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
667 outputContainer->Add(fhLam0ETRD[pidIndex]);
669 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
670 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
671 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
672 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
673 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
674 outputContainer->Add(fhLam1ETRD[pidIndex]);
676 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
677 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
678 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
679 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
680 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
681 outputContainer->Add(fhDispETRD[pidIndex]);
684 if(!fFillOnlySimpleSSHisto)
687 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
688 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
689 nbins,nmin, nmax, ssbins,ssmin,ssmax);
690 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
691 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
692 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
694 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
695 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
696 nbins,nmin, nmax, ssbins,ssmin,ssmax);
697 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
698 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
699 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
702 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
703 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
704 netabins,etamin,etamax, ssbins,ssmin,ssmax);
705 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
706 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
707 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
709 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
710 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
711 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
712 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
713 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
714 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
716 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
717 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
718 netabins,etamin,etamax, ssbins,ssmin,ssmax);
719 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
720 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
721 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
723 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
724 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
725 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
726 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
727 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
728 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
730 if(GetCalorimeter() == kEMCAL)
732 fhDispEtaE[pidIndex] = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
733 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
734 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
735 fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
736 fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
737 outputContainer->Add(fhDispEtaE[pidIndex]);
739 fhDispPhiE[pidIndex] = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
740 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
741 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
742 fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
743 fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
744 outputContainer->Add(fhDispPhiE[pidIndex]);
746 fhSumEtaE[pidIndex] = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
747 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
748 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
749 fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
750 fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
751 outputContainer->Add(fhSumEtaE[pidIndex]);
753 fhSumPhiE[pidIndex] = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
754 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),
755 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
756 fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
757 fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
758 outputContainer->Add(fhSumPhiE[pidIndex]);
760 fhSumEtaPhiE[pidIndex] = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
761 Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),
762 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
763 fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
764 fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
765 outputContainer->Add(fhSumEtaPhiE[pidIndex]);
767 fhDispEtaPhiDiffE[pidIndex] = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
768 Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
769 nptbins,ptmin,ptmax,200, -10,10);
770 fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
771 fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
772 outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
774 fhSphericityE[pidIndex] = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
775 Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),
776 nptbins,ptmin,ptmax, 200, -1,1);
777 fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
778 fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
779 outputContainer->Add(fhSphericityE[pidIndex]);
781 Int_t bin[] = {0,2,4,6,10,1000};
782 for(Int_t i = 0; i < 5; i++)
784 fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
785 Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
786 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
787 fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
788 fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
789 outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
797 if(fFillSSHistograms)
799 for(Int_t i = 0; i < 6; i++)
801 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
802 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
803 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
804 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
805 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
806 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
808 if(GetCalorimeter()==kEMCAL && !fFillOnlySimpleSSHisto)
810 fhMCEDispEta[pidIndex][i] = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
811 Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
812 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
813 fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
814 fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
815 outputContainer->Add(fhMCEDispEta[pidIndex][i]);
817 fhMCEDispPhi[pidIndex][i] = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
818 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
819 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
820 fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
821 fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
822 outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
824 fhMCESumEtaPhi[pidIndex][i] = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
825 Form("cluster from %s : %s like, #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
826 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
827 fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
828 fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
829 outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
831 fhMCEDispEtaPhiDiff[pidIndex][i] = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
832 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
833 nptbins,ptmin,ptmax,200,-10,10);
834 fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
835 fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
836 outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
838 fhMCESphericity[pidIndex][i] = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
839 Form("cluster from %s : %s like, (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
840 nptbins,ptmin,ptmax, 200,-1,1);
841 fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
842 fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
843 outputContainer->Add(fhMCESphericity[pidIndex][i]);
850 //if(IsCaloPIDOn() && pidIndex > 0) continue;
852 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
853 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
854 nptbins,ptmin,ptmax, nbins,nmin,nmax);
855 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
856 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
857 outputContainer->Add(fhNCellsE[pidIndex]);
859 fhNLME[pidIndex] = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
860 Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
861 nptbins,ptmin,ptmax, 10,0,10);
862 fhNLME[pidIndex]->SetXTitle("E (GeV)");
863 fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
864 outputContainer->Add(fhNLME[pidIndex]);
866 fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
867 Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
868 ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
869 fhTimeE[pidIndex]->SetXTitle("E (GeV)");
870 fhTimeE[pidIndex]->SetYTitle(" t (ns)");
871 outputContainer->Add(fhTimeE[pidIndex]);
873 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
874 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
875 nptbins,ptmin,ptmax, 500,0,1.);
876 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
877 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
878 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
880 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
881 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
882 nptbins,ptmin,ptmax);
883 fhE[pidIndex]->SetYTitle("N");
884 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
885 outputContainer->Add(fhE[pidIndex]) ;
887 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
888 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
889 nptbins,ptmin,ptmax);
890 fhPt[pidIndex]->SetYTitle("N");
891 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
892 outputContainer->Add(fhPt[pidIndex]) ;
894 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
895 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
896 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
897 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
898 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
899 outputContainer->Add(fhPhi[pidIndex]) ;
901 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
902 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
903 nptbins,ptmin,ptmax,netabins,etamin,etamax);
904 fhEta[pidIndex]->SetYTitle("#eta");
905 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
906 outputContainer->Add(fhEta[pidIndex]) ;
908 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
909 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
910 netabins,etamin,etamax,nphibins,phimin,phimax);
911 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
912 fhEtaPhi[pidIndex]->SetXTitle("#eta");
913 outputContainer->Add(fhEtaPhi[pidIndex]) ;
916 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
917 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
918 netabins,etamin,etamax,nphibins,phimin,phimax);
919 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
920 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
921 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
927 for(Int_t i = 0; i < fNOriginHistograms; i++)
929 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
930 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
931 nptbins,ptmin,ptmax);
932 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
933 outputContainer->Add(fhMCE[pidIndex][i]) ;
935 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
936 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
937 nptbins,ptmin,ptmax);
938 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
939 outputContainer->Add(fhMCPt[pidIndex][i]) ;
941 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
942 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
943 nptbins,ptmin,ptmax,netabins,etamin,etamax);
944 fhMCEta[pidIndex][i]->SetYTitle("#eta");
945 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
946 outputContainer->Add(fhMCEta[pidIndex][i]) ;
948 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
949 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
950 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
951 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
952 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
953 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
956 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
957 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
958 nptbins,ptmin,ptmax, 200,-50,50);
959 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
960 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
962 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
963 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
964 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
965 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
966 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
967 outputContainer->Add(fhMC2E[pidIndex][i]);
975 if(fFillSSHistograms)
979 if(!GetReader()->IsEmbeddedClusterSelectionOn())
981 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
982 "cluster from Electron : E vs #lambda_{0}^{2}",
983 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
984 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
985 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
986 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
988 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
989 "cluster from Electron : E vs #lambda_{0}^{2}",
990 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
991 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
992 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
993 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
995 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
996 "cluster from Electron : E vs #lambda_{0}^{2}",
997 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
998 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
999 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
1000 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
1004 //Fill histograms to check shape of embedded clusters
1005 if(GetReader()->IsEmbeddedClusterSelectionOn())
1008 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1009 "Energy Fraction of embedded signal versus cluster energy",
1010 nptbins,ptmin,ptmax,100,0.,1.);
1011 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1012 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1013 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1015 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
1016 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1017 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1018 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1019 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1020 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1022 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1023 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1024 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1025 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1026 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1027 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1029 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1030 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1031 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1032 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1033 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1034 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1036 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
1037 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1038 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1039 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1040 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1041 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1044 }// embedded histograms
1048 }// Fill SS MC histograms
1050 return outputContainer ;
1054 //_________________________
1055 void AliAnaElectron::Init()
1061 if ( GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD() )
1062 AliFatal("STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
1063 else if ( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
1064 AliFatal("STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
1068 //___________________________________
1069 void AliAnaElectron::InitParameters()
1072 //Initialize the parameters of the analysis.
1073 AddToHistogramsName("AnaElectron_");
1080 fTimeCutMax = 9999999;
1083 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1084 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1086 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1087 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1091 //_________________________________________
1092 void AliAnaElectron::MakeAnalysisFillAOD()
1094 //Do photon analysis and fill aods
1097 Double_t v[3] = {0,0,0}; //vertex ;
1098 GetReader()->GetVertex(v);
1100 //Select the Calorimeter of the photon
1101 TObjArray * pl = 0x0;
1102 if (GetCalorimeter() == kPHOS ) pl = GetPHOSClusters ();
1103 else if (GetCalorimeter() == kEMCAL) pl = GetEMCALClusters();
1107 AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
1111 //Init arrays, variables, get number of clusters
1112 Int_t nCaloClusters = pl->GetEntriesFast();
1113 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1115 AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
1117 //----------------------------------------------------
1118 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1119 //----------------------------------------------------
1121 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
1123 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1124 //printf("calo %d, %f\n",icalo,calo->E());
1126 //Get the index where the cluster comes, to retrieve the corresponding vertex
1127 Int_t evtIndex = 0 ;
1128 if (GetMixedEvent())
1130 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1131 //Get the vertex and check it is not too large in z
1132 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1135 //Cluster selection, not charged, with photon id and in fiducial cut
1136 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1138 calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
1139 }//Assume that come from vertex in straight line
1142 Double_t vertex[]={0,0,0};
1143 calo->GetMomentum(fMomentum,vertex) ;
1146 //--------------------------------------
1147 // Cluster selection
1148 //--------------------------------------
1149 AliVCaloCells* cells = 0;
1150 if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1151 else cells = GetPHOSCells();
1153 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1154 if(!ClusterSelected(calo,nMaxima)) continue;
1156 //-------------------------------------
1157 // PID selection via dE/dx
1158 //-------------------------------------
1160 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1164 AliWarning("Null track");
1168 //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1169 Float_t dEdx = track->GetTPCsignal();
1170 Float_t eOverp = calo->E()/track->P();
1172 fhdEdxvsE->Fill(calo->E(), dEdx);
1173 fhdEdxvsP->Fill(track->P(),dEdx);
1175 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1177 fhdEdxvsECutEOverP ->Fill(calo->E(), dEdx);
1178 fhdEdxvsPCutEOverP ->Fill(track->P(),dEdx);
1181 // Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1182 Float_t m02 = calo->GetM02();
1183 if(m02 > 0.1 && m02 < 0.4)
1185 fhdEdxvsECutM02 ->Fill(calo->E(), dEdx);
1186 fhdEdxvsPCutM02 ->Fill(track->P(),dEdx);
1187 fhEOverPvsECutM02->Fill(calo->E(), eOverp);
1188 fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1191 Int_t pid = AliCaloPID::kChargedHadron;
1193 if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1195 fhEOverPvsE->Fill(calo->E(), eOverp);
1196 fhEOverPvsP->Fill(track->P(), eOverp);
1198 if(m02 > 0.1 && m02 < 0.4)
1200 fhEOverPvsECutM02CutdEdx->Fill(calo->E(), eOverp);
1201 fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1204 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1206 pid = AliCaloPID::kElectron;
1211 Int_t pidIndex = 0;// Electron
1212 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1214 //--------------------------------------------------------------------------------------
1215 // Play with the MC stack if available
1216 //--------------------------------------------------------------------------------------
1218 //Check origin of the candidates
1222 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
1224 AliDebug(1,Form("Origin of candidate, bit map %d",tag));
1226 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1228 fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx);
1229 fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx);
1230 fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1231 fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1233 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1235 fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx);
1236 fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx);
1237 fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1238 fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1240 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1241 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1243 fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1244 fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx);
1245 fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1246 fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1248 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1250 fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx);
1251 fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx);
1252 fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1253 fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1255 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1257 fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx);
1258 fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx);
1259 fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1260 fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1262 else if( fhMCE[pidIndex][kmcOtherDecay] )
1264 fhMCdEdxvsE [kmcOtherDecay]->Fill(calo ->E(), dEdx);
1265 fhMCdEdxvsP [kmcOtherDecay]->Fill(track->P(), dEdx);
1266 fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
1267 fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
1270 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1272 fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1273 fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx);
1274 fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1275 fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1277 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1279 fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx);
1280 fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx);
1281 fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1282 fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1284 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1286 fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx);
1287 fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx);
1288 fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1289 fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1291 else if( fhMCE[pidIndex][kmcOther])
1293 fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx);
1294 fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx);
1295 fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1296 fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1298 }// set MC tag and fill Histograms with MC
1300 //---------------------------------
1301 //Fill some shower shape histograms
1302 //---------------------------------
1304 FillShowerShapeHistograms(calo,tag,pid);
1306 if(pid == AliCaloPID::kElectron)
1307 WeightHistograms(calo);
1309 //-----------------------------------------
1310 // PID Shower Shape selection or bit setting
1311 //-----------------------------------------
1313 // Data, PID check on
1316 // Get most probable PID, 2 options check bayesian PID weights or redo PID
1317 // By default, redo PID
1319 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1321 if(fAODParticle == AliCaloPID::kElectron)
1324 if(fAODParticle == 0 )
1325 pid = AliCaloPID::kChargedHadron ;
1328 AliDebug(1,Form("PDG of identified particle %d",pid));
1331 AliDebug(1,Form("Photon selection cuts passed: pT %3.2f, pdg %d",fMomentum.Pt(),pid));
1333 Float_t maxCellFraction = 0;
1334 Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1335 if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(fMomentum.E(),maxCellFraction);
1337 fhNCellsE[pidIndex] ->Fill(fMomentum.E(),calo->GetNCells());
1338 fhNLME [pidIndex] ->Fill(fMomentum.E(),nMaxima);
1339 fhTimeE [pidIndex] ->Fill(fMomentum.E(),calo->GetTOF()*1.e9);
1341 //----------------------------
1342 // Create AOD for analysis
1343 //----------------------------
1345 //Add AOD with electron/hadron object to aod branch
1346 if ( pid == fAODParticle || fAODParticle == 0 )
1348 AliAODPWG4Particle aodpart = AliAODPWG4Particle(fMomentum);
1350 //...............................................
1351 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1352 Int_t label = calo->GetLabel();
1353 aodpart.SetLabel(label);
1354 aodpart.SetCaloLabel (calo ->GetID(),-1);
1355 aodpart.SetTrackLabel(track->GetID(),-1);
1357 aodpart.SetDetectorTag(GetCalorimeter());
1358 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1360 aodpart.SetM02(calo->GetM02());
1361 aodpart.SetNLM(nMaxima);
1362 aodpart.SetTime(calo->GetTOF()*1e9);
1363 aodpart.SetNCells(calo->GetNCells());
1364 Int_t nSM = GetModuleNumber(calo);
1365 aodpart.SetSModNumber(nSM);
1367 //...............................................
1368 //Set bad channel distance bit
1369 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1370 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1371 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1372 else aodpart.SetDistToBad(0) ;
1373 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1376 aodpart.SetTag(tag);
1379 aodpart.SetIdentifiedParticleType(pid);
1381 AddAODParticle(aodpart);
1386 AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
1390 //________________________________________________
1391 void AliAnaElectron::MakeAnalysisFillHistograms()
1395 //-------------------------------------------------------------------
1396 // Access MC information in stack if requested, check that it exists.
1397 AliStack * stack = 0x0;
1398 TParticle * primary = 0x0;
1399 TClonesArray * mcparticles = 0x0;
1400 AliAODMCParticle * aodprimary = 0x0;
1404 if(GetReader()->ReadStack())
1406 stack = GetMCStack() ;
1409 AliFatal("Stack not available, is the MC handler called? STOP");
1413 else if(GetReader()->ReadAODMCParticles())
1415 //Get the list of MC particles
1416 mcparticles = GetReader()->GetAODMCParticles();
1419 AliFatal("Standard MCParticles not available! STOP");
1427 Double_t v[3] = {0,0,0}; //vertex ;
1428 GetReader()->GetVertex(v);
1429 //fhVertex->Fill(v[0],v[1],v[2]);
1430 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1432 //----------------------------------
1433 //Loop on stored AOD photons
1434 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1435 AliDebug(1,Form("AOD branch entries %d", naod));
1437 for(Int_t iaod = 0; iaod < naod ; iaod++)
1439 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1440 Int_t pdg = ph->GetIdentifiedParticleType();
1442 Int_t pidIndex = 0;// Electron
1443 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1444 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1447 if(((Int_t) ph->GetDetectorTag()) != GetCalorimeter()) continue;
1449 AliDebug(1,Form("ID Electron: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
1451 //................................
1452 //Fill photon histograms
1453 Float_t ptcluster = ph->Pt();
1454 Float_t phicluster = ph->Phi();
1455 Float_t etacluster = ph->Eta();
1456 Float_t ecluster = ph->E();
1458 fhE[pidIndex] ->Fill(ecluster);
1459 fhPt[pidIndex] ->Fill(ptcluster);
1460 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1461 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1462 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1463 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1465 //.......................................
1466 //Play with the MC data if available
1469 //....................................................................
1470 // Access MC information in stack if requested, check that it exists.
1471 Int_t label =ph->GetLabel();
1474 AliDebug(1,Form("*** bad label ***: label %d", label));
1479 //Float_t ptprim = 0;
1480 if( GetReader()->ReadStack() )
1482 if(label >= stack->GetNtrack())
1484 AliDebug(1,Form("*** large label ***: label %d, n tracks %d", label, stack->GetNtrack()));
1488 primary = stack->Particle(label);
1491 AliWarning(Form("*** no primary ***: label %d", label));
1495 eprim = primary->Energy();
1496 //ptprim = primary->Pt();
1499 else if( GetReader()->ReadAODMCParticles() )
1501 if(label >= mcparticles->GetEntriesFast())
1503 AliDebug(1,Form("*** large label ***: label %d, n tracks %d",label, mcparticles->GetEntriesFast()));
1507 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1511 AliWarning(Form("*** no primary ***: label %d", label));
1515 eprim = aodprimary->E();
1516 //ptprim = aodprimary->Pt();
1519 Int_t tag =ph->GetTag();
1521 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1523 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1524 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1525 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1526 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1528 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1529 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1531 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1533 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1534 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1535 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1536 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1538 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1539 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1542 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1543 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1545 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1546 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1547 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1548 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1550 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1551 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1553 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1554 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1556 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1557 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1558 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1559 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1561 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1562 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1564 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1566 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1567 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1568 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1569 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1571 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1572 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1575 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1577 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1578 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1579 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1580 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1582 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1583 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1587 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1589 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1590 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1591 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1592 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1594 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1595 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1598 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1600 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1601 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1602 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1603 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1605 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1606 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1609 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1611 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1612 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1613 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1614 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1616 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1617 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1620 else if( fhMCE[pidIndex][kmcOther])
1622 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1623 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1624 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1625 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1627 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1628 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1632 }//Histograms with MC
1638 //____________________________________________________
1639 void AliAnaElectron::Print(const Option_t * opt) const
1641 //Print some relevant parameters set for the analysis
1646 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1647 AliAnaCaloTrackCorrBaseClass::Print(" ");
1649 printf("Calorimeter = %s\n", GetCalorimeterString().Data()) ;
1650 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1651 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1652 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1653 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1654 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1655 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1656 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1661 //______________________________________________________
1662 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1664 // Calculate weights and fill histograms
1666 if(!fFillWeightHistograms || GetMixedEvent()) return;
1668 AliVCaloCells* cells = 0;
1669 if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1670 else cells = GetPHOSCells();
1672 // First recalculate energy in case non linearity was applied
1675 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1677 Int_t id = clus->GetCellsAbsId()[ipos];
1679 //Recalibrate cell energy if needed
1680 Float_t amp = cells->GetCellAmplitude(id);
1681 GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
1692 AliWarning(Form("Wrong calculated energy %f",energy));
1696 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1697 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1698 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1700 //Get the ratio and log ratio to all cells in cluster
1701 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1702 Int_t id = clus->GetCellsAbsId()[ipos];
1704 //Recalibrate cell energy if needed
1705 Float_t amp = cells->GetCellAmplitude(id);
1706 GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
1708 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1709 fhECellClusterRatio ->Fill(energy,amp/energy);
1710 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1713 //Recalculate shower shape for different W0
1714 if(GetCalorimeter()==kEMCAL)
1716 Float_t l0org = clus->GetM02();
1717 Float_t l1org = clus->GetM20();
1718 Float_t dorg = clus->GetDispersion();
1720 for(Int_t iw = 0; iw < 14; iw++){
1722 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1723 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1725 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1726 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1728 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1732 // Set the original values back
1733 clus->SetM02(l0org);
1734 clus->SetM20(l1org);
1735 clus->SetDispersion(dorg);