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(), fCalorimeter(""),
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),
63 fhdEdxvsE(0), fhdEdxvsP(0),
64 fhEOverPvsE(0), fhEOverPvsP(0),
65 fhdEdxvsECutM02(0), fhdEdxvsPCutM02(0),
66 fhEOverPvsECutM02(0), fhEOverPvsPCutM02(0),
67 fhdEdxvsECutEOverP(0), fhdEdxvsPCutEOverP(0),
68 fhEOverPvsECutM02CutdEdx(0), fhEOverPvsPCutM02CutdEdx(0),
70 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
71 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
73 // Electron SS MC histograms
74 fhMCElectronELambda0NoOverlap(0),
75 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
78 fhEmbeddedSignalFractionEnergy(0),
79 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
80 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
83 for(Int_t index = 0; index < 2; index++)
85 fhNCellsE [index] = 0;
88 fhMaxCellDiffClusterE[index] = 0;
94 fhEtaPhi05[index] = 0;
96 // Shower shape histograms
100 fhDispETRD[index] = 0;
101 fhLam0ETRD[index] = 0;
102 fhLam1ETRD[index] = 0;
103 fhNCellsLam0LowE [index] = 0;
104 fhNCellsLam0HighE[index] = 0;
105 fhEtaLam0LowE [index] = 0;
106 fhPhiLam0LowE [index] = 0;
107 fhEtaLam0HighE [index] = 0;
108 fhPhiLam0HighE [index] = 0;
110 fhDispEtaE [index] = 0;
111 fhDispPhiE [index] = 0;
112 fhSumEtaE [index] = 0;
113 fhSumPhiE [index] = 0;
114 fhSumEtaPhiE [index] = 0;
115 fhDispEtaPhiDiffE[index] = 0;
116 fhSphericityE [index] = 0;
118 for(Int_t i = 0; i < 10; i++)
120 fhMCPt [index][i] = 0;
121 fhMCE [index][i] = 0;
122 fhMCPhi [index][i] = 0;
123 fhMCEta [index][i] = 0;
124 fhMCDeltaE [index][i] = 0;
125 fhMC2E [index][i] = 0;
128 fhMCEOverPvsE [i] = 0;
129 fhMCEOverPvsP [i] = 0;
132 for(Int_t i = 0; i < 6; i++)
134 fhMCELambda0 [index][i] = 0;
135 fhMCEDispEta [index][i] = 0;
136 fhMCEDispPhi [index][i] = 0;
137 fhMCESumEtaPhi [index][i] = 0;
138 fhMCEDispEtaPhiDiff[index][i] = 0;
139 fhMCESphericity [index][i] = 0;
142 for(Int_t i = 0; i < 5; i++)
144 fhDispEtaDispPhiEBin[index][i] = 0 ;
149 for(Int_t i =0; i < 14; i++)
151 fhLambda0ForW0[i] = 0;
152 //fhLambda1ForW0[i] = 0;
155 //Initialize parameters
160 //____________________________________________________________________________
161 Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
163 //Select clusters if they pass different cuts
165 printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
166 GetReader()->GetEventNumber(),
167 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
169 //.......................................
170 //If too small or big energy, skip it
171 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
172 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",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 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
180 //.......................................
181 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
182 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
184 //.......................................
185 //Check acceptance selection
186 if(IsFiducialCutOn()){
187 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
188 if(! in ) return kFALSE ;
190 if(GetDebug() > 2) printf("Fiducial cut passed \n");
192 //.......................................
193 //Skip not matched clusters with tracks
194 if(!IsTrackMatched(calo, GetReader()->GetInputEvent())) {
195 if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
198 else if(GetDebug() > 2) printf(" Track-matching cut passed \n");
200 //...........................................
201 // skip clusters with too many maxima
202 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
203 if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
205 //.......................................
206 //Check Distance to Bad channel, set bit.
207 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
208 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
209 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
212 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
213 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
216 printf("AliAnaElectron::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
217 GetReader()->GetEventNumber(),
218 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
220 //All checks passed, cluster selected
225 //______________________________________________________________________________________________
226 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.;
249 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
251 cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
253 Double_t vertex[]={0,0,0};
254 cluster->GetMomentum(mom,vertex) ;
257 Float_t eta = mom.Eta();
258 Float_t phi = mom.Phi();
259 if(phi < 0) phi+=TMath::TwoPi();
261 fhLam0E[pidIndex] ->Fill(energy,lambda0);
262 fhLam1E[pidIndex] ->Fill(energy,lambda1);
263 fhDispE[pidIndex] ->Fill(energy,disp);
265 if(fCalorimeter == "EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
266 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
268 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
269 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
270 fhDispETRD[pidIndex]->Fill(energy,disp);
273 if(!fFillOnlySimpleSSHisto)
277 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
278 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
279 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
283 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
284 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
285 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
288 if(fCalorimeter == "EMCAL")
290 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
291 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
292 fhDispEtaE [pidIndex]-> Fill(energy,dEta);
293 fhDispPhiE [pidIndex]-> Fill(energy,dPhi);
294 fhSumEtaE [pidIndex]-> Fill(energy,sEta);
295 fhSumPhiE [pidIndex]-> Fill(energy,sPhi);
296 fhSumEtaPhiE [pidIndex]-> Fill(energy,sEtaPhi);
297 fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
298 if(dEta+dPhi>0)fhSphericityE [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
300 if (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
301 else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
302 else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
303 else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
304 else fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
311 AliVCaloCells* cells = 0;
312 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
313 else cells = GetPHOSCells();
315 //Fill histograms to check shape of embedded clusters
316 Float_t fraction = 0;
317 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
319 Float_t clusterE = 0; // recalculate in case corrections applied.
321 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
322 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
324 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
327 //Fraction of total energy due to the embedded signal
330 if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
332 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
334 } // embedded fraction
336 // Check the origin and fill histograms
339 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
340 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
341 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
342 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
346 }//photon no conversion
347 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
348 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
350 index = kmcssElectron;
352 if(!GetReader()->IsEmbeddedClusterSelectionOn())
354 //Check particle overlaps in cluster
356 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
357 Int_t ancPDG = 0, ancStatus = -1;
358 TLorentzVector momentum; TVector3 prodVertex;
361 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
362 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
363 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
367 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
369 else if(noverlaps == 2){
370 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
372 else if(noverlaps > 2){
373 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
376 printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d for ancestor %d!!", noverlaps, ancLabel);
380 //Fill histograms to check shape of embedded clusters
381 if(GetReader()->IsEmbeddedClusterSelectionOn())
385 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
387 else if(fraction > 0.5)
389 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
391 else if(fraction > 0.1)
393 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
397 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
401 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
402 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
404 index = kmcssConversion;
406 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
410 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
419 fhMCELambda0[pidIndex][index] ->Fill(energy, lambda0);
421 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
423 fhMCEDispEta [pidIndex][index]-> Fill(energy,dEta);
424 fhMCEDispPhi [pidIndex][index]-> Fill(energy,dPhi);
425 fhMCESumEtaPhi [pidIndex][index]-> Fill(energy,sEtaPhi);
426 fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
427 if(dEta+dPhi>0)fhMCESphericity [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
434 //_____________________________________________
435 TObjString * AliAnaElectron::GetAnalysisCuts()
437 //Save parameters used for analysis
438 TString parList ; //this will be list of parameters used for this analysis.
439 const Int_t buffersize = 255;
440 char onePar[buffersize] ;
442 snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
444 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
446 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
448 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f \n",fEOverPMin, fEOverPMax) ;
450 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
452 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
454 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
457 //Get parameters set in base class.
458 parList += GetBaseParametersList() ;
460 //Get parameters set in PID class.
461 parList += GetCaloPID()->GetPIDParametersList() ;
463 //Get parameters set in FiducialCut class (not available yet)
464 //parlist += GetFidCut()->GetFidCutParametersList()
466 return new TObjString(parList) ;
469 //_______________________________________________
470 TList * AliAnaElectron::GetCreateOutputObjects()
472 // Create histograms to be saved in output file and
473 // store them in outputContainer
474 TList * outputContainer = new TList() ;
475 outputContainer->SetName("ElectronHistos") ;
477 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
478 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
479 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
480 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
481 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
482 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
483 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
484 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
487 // MC labels, titles, for originator particles
488 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
489 TString pnamess[] = { "Photon","Hadron" ,"Pi0" ,"Eta" ,"Conversion" ,"Electron"} ;
490 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
491 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
493 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
494 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
497 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
498 fhdEdxvsE->SetXTitle("E (GeV)");
499 fhdEdxvsE->SetYTitle("<dE/dx>");
500 outputContainer->Add(fhdEdxvsE);
502 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
503 fhdEdxvsP->SetXTitle("P (GeV/c)");
504 fhdEdxvsP->SetYTitle("<dE/dx>");
505 outputContainer->Add(fhdEdxvsP);
507 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
508 fhEOverPvsE->SetXTitle("E (GeV)");
509 fhEOverPvsE->SetYTitle("E/p");
510 outputContainer->Add(fhEOverPvsE);
512 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
513 fhEOverPvsP->SetXTitle("P (GeV/c)");
514 fhEOverPvsP->SetYTitle("E/p");
515 outputContainer->Add(fhEOverPvsP);
518 fhdEdxvsECutM02 = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
519 fhdEdxvsECutM02->SetXTitle("E (GeV)");
520 fhdEdxvsECutM02->SetYTitle("<dE/dx>");
521 outputContainer->Add(fhdEdxvsECutM02);
523 fhdEdxvsPCutM02 = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
524 fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
525 fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
526 outputContainer->Add(fhdEdxvsPCutM02);
528 fhEOverPvsECutM02 = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
529 fhEOverPvsECutM02->SetXTitle("E (GeV)");
530 fhEOverPvsECutM02->SetYTitle("E/p");
531 outputContainer->Add(fhEOverPvsECutM02);
533 fhEOverPvsPCutM02 = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
534 fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
535 fhEOverPvsPCutM02->SetYTitle("E/p");
536 outputContainer->Add(fhEOverPvsPCutM02);
539 fhdEdxvsECutEOverP = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
540 fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
541 fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
542 outputContainer->Add(fhdEdxvsECutEOverP);
544 fhdEdxvsPCutEOverP = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
545 fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
546 fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
547 outputContainer->Add(fhdEdxvsPCutEOverP);
549 fhEOverPvsECutM02CutdEdx = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
550 fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
551 fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
552 outputContainer->Add(fhEOverPvsECutM02CutdEdx);
554 fhEOverPvsPCutM02CutdEdx = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
555 fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
556 fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
557 outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
561 for(Int_t i = 0; i < fNOriginHistograms; i++)
563 fhMCdEdxvsE[i] = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
564 Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
565 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
566 fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
567 fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
568 outputContainer->Add(fhMCdEdxvsE[i]) ;
570 fhMCdEdxvsP[i] = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
571 Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
572 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
573 fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
574 fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
575 outputContainer->Add(fhMCdEdxvsP[i]) ;
578 fhMCEOverPvsE[i] = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
579 Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
580 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
581 fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
582 fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
583 outputContainer->Add(fhMCEOverPvsE[i]) ;
585 fhMCEOverPvsP[i] = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
586 Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
587 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
588 fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
589 fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
590 outputContainer->Add(fhMCEOverPvsP[i]) ;
595 TString pidParticle[] = {"Electron","ChargedHadron"} ;
597 if(fFillWeightHistograms)
600 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
601 nptbins,ptmin,ptmax, 100,0,1.);
602 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
603 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
604 outputContainer->Add(fhECellClusterRatio);
606 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
607 nptbins,ptmin,ptmax, 100,-10,0);
608 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
609 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
610 outputContainer->Add(fhECellClusterLogRatio);
612 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
613 nptbins,ptmin,ptmax, 100,0,1.);
614 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
615 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
616 outputContainer->Add(fhEMaxCellClusterRatio);
618 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
619 nptbins,ptmin,ptmax, 100,-10,0);
620 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
621 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
622 outputContainer->Add(fhEMaxCellClusterLogRatio);
624 for(Int_t iw = 0; iw < 14; iw++)
626 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),
627 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
628 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
629 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
630 outputContainer->Add(fhLambda0ForW0[iw]);
632 // 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),
633 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
634 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
635 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
636 // outputContainer->Add(fhLambda1ForW0[iw]);
641 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
644 if(fFillSSHistograms)
646 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
647 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
648 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
649 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
650 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
651 outputContainer->Add(fhLam0E[pidIndex]);
653 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
654 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
655 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
656 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
657 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
658 outputContainer->Add(fhLam1E[pidIndex]);
660 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
661 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
662 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
663 fhDispE[pidIndex]->SetYTitle("D^{2}");
664 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
665 outputContainer->Add(fhDispE[pidIndex]);
667 if(fCalorimeter == "EMCAL" && GetFirstSMCoveredByTRD() >=0 )
669 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
670 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
671 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
672 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
673 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
674 outputContainer->Add(fhLam0ETRD[pidIndex]);
676 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
677 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
678 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
679 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
680 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
681 outputContainer->Add(fhLam1ETRD[pidIndex]);
683 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
684 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
685 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
686 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
687 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
688 outputContainer->Add(fhDispETRD[pidIndex]);
691 if(!fFillOnlySimpleSSHisto)
694 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",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 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
698 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
699 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
701 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
702 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
703 nbins,nmin, nmax, ssbins,ssmin,ssmax);
704 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
705 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
706 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
709 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
710 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
711 netabins,etamin,etamax, ssbins,ssmin,ssmax);
712 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
713 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
714 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
716 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
717 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
718 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
719 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
720 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
721 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
723 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
724 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
725 netabins,etamin,etamax, ssbins,ssmin,ssmax);
726 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
727 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
728 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
730 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
731 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
732 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
733 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
734 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
735 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
737 if(fCalorimeter == "EMCAL")
739 fhDispEtaE[pidIndex] = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
740 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
741 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
742 fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
743 fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
744 outputContainer->Add(fhDispEtaE[pidIndex]);
746 fhDispPhiE[pidIndex] = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
747 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
748 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
749 fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
750 fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
751 outputContainer->Add(fhDispPhiE[pidIndex]);
753 fhSumEtaE[pidIndex] = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
754 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
755 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
756 fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
757 fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
758 outputContainer->Add(fhSumEtaE[pidIndex]);
760 fhSumPhiE[pidIndex] = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
761 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),
762 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
763 fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
764 fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
765 outputContainer->Add(fhSumPhiE[pidIndex]);
767 fhSumEtaPhiE[pidIndex] = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
768 Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),
769 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
770 fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
771 fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
772 outputContainer->Add(fhSumEtaPhiE[pidIndex]);
774 fhDispEtaPhiDiffE[pidIndex] = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
775 Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
776 nptbins,ptmin,ptmax,200, -10,10);
777 fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
778 fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
779 outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
781 fhSphericityE[pidIndex] = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
782 Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),
783 nptbins,ptmin,ptmax, 200, -1,1);
784 fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
785 fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
786 outputContainer->Add(fhSphericityE[pidIndex]);
788 Int_t bin[] = {0,2,4,6,10,1000};
789 for(Int_t i = 0; i < 5; i++)
791 fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
792 Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
793 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
794 fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
795 fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
796 outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
804 if(fFillSSHistograms)
806 for(Int_t i = 0; i < 6; i++)
808 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
809 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
810 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
811 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
812 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
813 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
815 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
817 fhMCEDispEta[pidIndex][i] = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
818 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()),
819 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
820 fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
821 fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
822 outputContainer->Add(fhMCEDispEta[pidIndex][i]);
824 fhMCEDispPhi[pidIndex][i] = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
825 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()),
826 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
827 fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
828 fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
829 outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
831 fhMCESumEtaPhi[pidIndex][i] = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
832 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()),
833 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
834 fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
835 fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
836 outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
838 fhMCEDispEtaPhiDiff[pidIndex][i] = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
839 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
840 nptbins,ptmin,ptmax,200,-10,10);
841 fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
842 fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
843 outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
845 fhMCESphericity[pidIndex][i] = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
846 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()),
847 nptbins,ptmin,ptmax, 200,-1,1);
848 fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
849 fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
850 outputContainer->Add(fhMCESphericity[pidIndex][i]);
857 //if(IsCaloPIDOn() && pidIndex > 0) continue;
859 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
860 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
861 nptbins,ptmin,ptmax, nbins,nmin,nmax);
862 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
863 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
864 outputContainer->Add(fhNCellsE[pidIndex]);
866 fhNLME[pidIndex] = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
867 Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
868 nptbins,ptmin,ptmax, 10,0,10);
869 fhNLME[pidIndex]->SetXTitle("E (GeV)");
870 fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
871 outputContainer->Add(fhNLME[pidIndex]);
873 fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
874 Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
875 ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
876 fhTimeE[pidIndex]->SetXTitle("E (GeV)");
877 fhTimeE[pidIndex]->SetYTitle(" t (ns)");
878 outputContainer->Add(fhTimeE[pidIndex]);
880 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
881 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
882 nptbins,ptmin,ptmax, 500,0,1.);
883 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
884 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
885 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
887 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
888 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
889 nptbins,ptmin,ptmax);
890 fhE[pidIndex]->SetYTitle("N");
891 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
892 outputContainer->Add(fhE[pidIndex]) ;
894 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
895 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
896 nptbins,ptmin,ptmax);
897 fhPt[pidIndex]->SetYTitle("N");
898 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
899 outputContainer->Add(fhPt[pidIndex]) ;
901 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
902 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
903 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
904 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
905 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
906 outputContainer->Add(fhPhi[pidIndex]) ;
908 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
909 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
910 nptbins,ptmin,ptmax,netabins,etamin,etamax);
911 fhEta[pidIndex]->SetYTitle("#eta");
912 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
913 outputContainer->Add(fhEta[pidIndex]) ;
915 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
916 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
917 netabins,etamin,etamax,nphibins,phimin,phimax);
918 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
919 fhEtaPhi[pidIndex]->SetXTitle("#eta");
920 outputContainer->Add(fhEtaPhi[pidIndex]) ;
923 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
924 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
925 netabins,etamin,etamax,nphibins,phimin,phimax);
926 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
927 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
928 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
934 for(Int_t i = 0; i < fNOriginHistograms; i++)
936 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
937 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
938 nptbins,ptmin,ptmax);
939 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
940 outputContainer->Add(fhMCE[pidIndex][i]) ;
942 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
943 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
944 nptbins,ptmin,ptmax);
945 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
946 outputContainer->Add(fhMCPt[pidIndex][i]) ;
948 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
949 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
950 nptbins,ptmin,ptmax,netabins,etamin,etamax);
951 fhMCEta[pidIndex][i]->SetYTitle("#eta");
952 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
953 outputContainer->Add(fhMCEta[pidIndex][i]) ;
955 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
956 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
957 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
958 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
959 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
960 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
963 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
964 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
965 nptbins,ptmin,ptmax, 200,-50,50);
966 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
967 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
969 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
970 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
971 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
972 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
973 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
974 outputContainer->Add(fhMC2E[pidIndex][i]);
982 if(fFillSSHistograms)
986 if(!GetReader()->IsEmbeddedClusterSelectionOn())
988 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
989 "cluster from Electron : E vs #lambda_{0}^{2}",
990 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
991 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
992 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
993 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
995 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
996 "cluster from Electron : E vs #lambda_{0}^{2}",
997 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
998 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
999 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
1000 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
1002 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
1003 "cluster from Electron : E vs #lambda_{0}^{2}",
1004 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1005 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1006 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
1007 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
1011 //Fill histograms to check shape of embedded clusters
1012 if(GetReader()->IsEmbeddedClusterSelectionOn())
1015 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1016 "Energy Fraction of embedded signal versus cluster energy",
1017 nptbins,ptmin,ptmax,100,0.,1.);
1018 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1019 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1020 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1022 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
1023 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1024 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1025 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1026 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1027 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1029 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1030 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1031 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1032 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1033 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1034 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1036 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1037 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1038 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1039 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1040 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1041 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1043 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
1044 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1045 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1046 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1047 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1048 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1051 }// embedded histograms
1055 }// Fill SS MC histograms
1057 return outputContainer ;
1061 //_________________________
1062 void AliAnaElectron::Init()
1067 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1068 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1071 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1072 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1078 //___________________________________
1079 void AliAnaElectron::InitParameters()
1082 //Initialize the parameters of the analysis.
1083 AddToHistogramsName("AnaElectron_");
1085 fCalorimeter = "EMCAL" ;
1091 fTimeCutMax = 9999999;
1094 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1095 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1097 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1098 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1102 //_________________________________________
1103 void AliAnaElectron::MakeAnalysisFillAOD()
1105 //Do photon analysis and fill aods
1108 Double_t v[3] = {0,0,0}; //vertex ;
1109 GetReader()->GetVertex(v);
1111 //Select the Calorimeter of the photon
1112 TObjArray * pl = 0x0;
1113 if(fCalorimeter == "PHOS")
1114 pl = GetPHOSClusters();
1115 else if (fCalorimeter == "EMCAL")
1116 pl = GetEMCALClusters();
1120 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1124 //Init arrays, variables, get number of clusters
1125 TLorentzVector mom, mom2 ;
1126 Int_t nCaloClusters = pl->GetEntriesFast();
1127 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1129 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
1131 //----------------------------------------------------
1132 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1133 //----------------------------------------------------
1135 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
1137 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1138 //printf("calo %d, %f\n",icalo,calo->E());
1140 //Get the index where the cluster comes, to retrieve the corresponding vertex
1141 Int_t evtIndex = 0 ;
1142 if (GetMixedEvent())
1144 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1145 //Get the vertex and check it is not too large in z
1146 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1149 //Cluster selection, not charged, with photon id and in fiducial cut
1150 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1152 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1154 Double_t vertex[]={0,0,0};
1155 calo->GetMomentum(mom,vertex) ;
1158 //--------------------------------------
1159 // Cluster selection
1160 //--------------------------------------
1161 AliVCaloCells* cells = 0;
1162 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1163 else cells = GetPHOSCells();
1165 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1166 if(!ClusterSelected(calo,mom,nMaxima)) continue;
1168 //-------------------------------------
1169 // PID selection via dE/dx
1170 //-------------------------------------
1172 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1176 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
1180 //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1181 Float_t dEdx = track->GetTPCsignal();
1182 Float_t eOverp = calo->E()/track->P();
1184 fhdEdxvsE->Fill(calo->E(), dEdx);
1185 fhdEdxvsP->Fill(track->P(),dEdx);
1187 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1189 fhdEdxvsECutEOverP ->Fill(calo->E(), dEdx);
1190 fhdEdxvsPCutEOverP ->Fill(track->P(),dEdx);
1193 //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1194 Float_t m02 = calo->GetM02();
1195 if(m02 > 0.1 && m02 < 0.4)
1197 fhdEdxvsECutM02 ->Fill(calo->E(), dEdx);
1198 fhdEdxvsPCutM02 ->Fill(track->P(),dEdx);
1199 fhEOverPvsECutM02->Fill(calo->E(), eOverp);
1200 fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1203 Int_t pid = AliCaloPID::kChargedHadron;
1205 if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1207 fhEOverPvsE->Fill(calo->E(), eOverp);
1208 fhEOverPvsP->Fill(track->P(), eOverp);
1210 if(m02 > 0.1 && m02 < 0.4)
1212 fhEOverPvsECutM02CutdEdx->Fill(calo->E(), eOverp);
1213 fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1216 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1218 pid = AliCaloPID::kElectron;
1223 Int_t pidIndex = 0;// Electron
1224 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1226 //--------------------------------------------------------------------------------------
1227 //Play with the MC stack if available
1228 //--------------------------------------------------------------------------------------
1230 //Check origin of the candidates
1234 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
1237 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",tag);
1239 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1241 fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx);
1242 fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx);
1243 fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1244 fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1246 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1248 fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx);
1249 fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx);
1250 fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1251 fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1253 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1254 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1256 fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1257 fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx);
1258 fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1259 fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1261 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1262 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && 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);
1269 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1271 fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx);
1272 fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx);
1273 fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1274 fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1276 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1278 fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx);
1279 fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx);
1280 fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1281 fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1284 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1286 fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1287 fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx);
1288 fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1289 fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1291 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1293 fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx);
1294 fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx);
1295 fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1296 fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1298 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1300 fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx);
1301 fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx);
1302 fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1303 fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1305 else if( fhMCE[pidIndex][kmcOther])
1307 fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx);
1308 fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx);
1309 fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1310 fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1312 }// set MC tag and fill Histograms with MC
1314 //---------------------------------
1315 //Fill some shower shape histograms
1316 //---------------------------------
1318 FillShowerShapeHistograms(calo,tag,pid);
1320 if(pid == AliCaloPID::kElectron)
1321 WeightHistograms(calo);
1323 //-----------------------------------------
1324 //PID Shower Shape selection or bit setting
1325 //-----------------------------------------
1327 // Data, PID check on
1330 // Get most probable PID, 2 options check bayesian PID weights or redo PID
1331 // By default, redo PID
1333 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1335 if(fAODParticle == AliCaloPID::kElectron)
1338 if(fAODParticle == 0 )
1339 pid = AliCaloPID::kChargedHadron ;
1341 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - PDG of identified particle %d\n",pid);
1344 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
1347 Float_t maxCellFraction = 0;
1348 Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1349 if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(mom.E(),maxCellFraction);
1351 fhNCellsE[pidIndex] ->Fill(mom.E(),calo->GetNCells());
1352 fhNLME [pidIndex] ->Fill(mom.E(),nMaxima);
1353 fhTimeE [pidIndex] ->Fill(mom.E(),calo->GetTOF()*1.e9);
1356 //----------------------------
1357 //Create AOD for analysis
1358 //----------------------------
1360 //Add AOD with electron/hadron object to aod branch
1361 if ( pid == fAODParticle || fAODParticle == 0 )
1363 AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
1365 //...............................................
1366 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1367 Int_t label = calo->GetLabel();
1368 aodpart.SetLabel(label);
1369 aodpart.SetCaloLabel (calo ->GetID(),-1);
1370 aodpart.SetTrackLabel(track->GetID(),-1);
1372 aodpart.SetDetector(fCalorimeter);
1373 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1375 //...............................................
1376 //Set bad channel distance bit
1377 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1378 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1379 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1380 else aodpart.SetDistToBad(0) ;
1381 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1384 aodpart.SetTag(tag);
1387 aodpart.SetIdentifiedParticleType(pid);
1389 AddAODParticle(aodpart);
1394 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1398 //________________________________________________
1399 void AliAnaElectron::MakeAnalysisFillHistograms()
1403 //-------------------------------------------------------------------
1404 // Access MC information in stack if requested, check that it exists.
1405 AliStack * stack = 0x0;
1406 TParticle * primary = 0x0;
1407 TClonesArray * mcparticles = 0x0;
1408 AliAODMCParticle * aodprimary = 0x0;
1412 if(GetReader()->ReadStack())
1414 stack = GetMCStack() ;
1417 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1421 else if(GetReader()->ReadAODMCParticles())
1423 //Get the list of MC particles
1424 mcparticles = GetReader()->GetAODMCParticles();
1427 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1435 Double_t v[3] = {0,0,0}; //vertex ;
1436 GetReader()->GetVertex(v);
1437 //fhVertex->Fill(v[0],v[1],v[2]);
1438 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1440 //----------------------------------
1441 //Loop on stored AOD photons
1442 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1443 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1445 for(Int_t iaod = 0; iaod < naod ; iaod++)
1447 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1448 Int_t pdg = ph->GetIdentifiedParticleType();
1450 Int_t pidIndex = 0;// Electron
1451 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1452 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1455 if(ph->GetDetector() != fCalorimeter) continue;
1458 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1460 //................................
1461 //Fill photon histograms
1462 Float_t ptcluster = ph->Pt();
1463 Float_t phicluster = ph->Phi();
1464 Float_t etacluster = ph->Eta();
1465 Float_t ecluster = ph->E();
1467 fhE[pidIndex] ->Fill(ecluster);
1468 fhPt[pidIndex] ->Fill(ptcluster);
1469 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1470 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1471 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1472 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1474 //.......................................
1475 //Play with the MC data if available
1478 //....................................................................
1479 // Access MC information in stack if requested, check that it exists.
1480 Int_t label =ph->GetLabel();
1482 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1487 //Float_t ptprim = 0;
1488 if(GetReader()->ReadStack())
1490 if(label >= stack->GetNtrack())
1492 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1496 primary = stack->Particle(label);
1499 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1503 eprim = primary->Energy();
1504 //ptprim = primary->Pt();
1507 else if(GetReader()->ReadAODMCParticles())
1509 //Check which is the input
1510 if(ph->GetInputFileIndex() == 0)
1512 if(!mcparticles) continue;
1514 if(label >= mcparticles->GetEntriesFast())
1516 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1517 label, mcparticles->GetEntriesFast());
1521 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1527 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1531 eprim = aodprimary->E();
1532 //ptprim = aodprimary->Pt();
1535 Int_t tag =ph->GetTag();
1537 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1539 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1540 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1541 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1542 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1544 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1545 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1547 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1549 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1550 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1551 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1552 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1554 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1555 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1558 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1559 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1561 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1562 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1563 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1564 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1566 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1567 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1569 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1570 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1572 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1573 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1574 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1575 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1577 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1578 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1580 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1582 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1583 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1584 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1585 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1587 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1588 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1591 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1593 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1594 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1595 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1596 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1598 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1599 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1603 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1605 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1606 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1607 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1608 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1610 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1611 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1614 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1616 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1617 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1618 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1619 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1621 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1622 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1625 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1627 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1628 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1629 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1630 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1632 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1633 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1636 else if( fhMCE[pidIndex][kmcOther])
1638 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1639 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1640 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1641 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1643 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1644 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1648 }//Histograms with MC
1654 //____________________________________________________
1655 void AliAnaElectron::Print(const Option_t * opt) const
1657 //Print some relevant parameters set for the analysis
1662 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1663 AliAnaCaloTrackCorrBaseClass::Print(" ");
1665 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1666 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1667 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1668 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1669 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1670 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1671 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1672 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1677 //______________________________________________________
1678 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1680 // Calculate weights and fill histograms
1682 if(!fFillWeightHistograms || GetMixedEvent()) return;
1684 AliVCaloCells* cells = 0;
1685 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1686 else cells = GetPHOSCells();
1688 // First recalculate energy in case non linearity was applied
1691 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1693 Int_t id = clus->GetCellsAbsId()[ipos];
1695 //Recalibrate cell energy if needed
1696 Float_t amp = cells->GetCellAmplitude(id);
1697 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
1707 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1711 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1712 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1713 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1715 //Get the ratio and log ratio to all cells in cluster
1716 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1717 Int_t id = clus->GetCellsAbsId()[ipos];
1719 //Recalibrate cell energy if needed
1720 Float_t amp = cells->GetCellAmplitude(id);
1721 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
1723 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1724 fhECellClusterRatio ->Fill(energy,amp/energy);
1725 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1728 //Recalculate shower shape for different W0
1729 if(fCalorimeter=="EMCAL"){
1731 Float_t l0org = clus->GetM02();
1732 Float_t l1org = clus->GetM20();
1733 Float_t dorg = clus->GetDispersion();
1735 for(Int_t iw = 0; iw < 14; iw++){
1737 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1738 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1740 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1741 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1743 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1747 // Set the original values back
1748 clus->SetM02(l0org);
1749 clus->SetM20(l1org);
1750 clus->SetDispersion(dorg);