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),
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,GetCalorimeter()) ;
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(GetCalorimeter() == "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(GetCalorimeter() == "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(GetCalorimeter() == "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(GetCalorimeter() == "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",GetCalorimeter().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(GetCalorimeter() == "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(GetCalorimeter() == "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(GetCalorimeter()=="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(GetCalorimeter() == "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(GetCalorimeter() == "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_");
1090 fTimeCutMax = 9999999;
1093 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1094 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1096 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1097 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1101 //_________________________________________
1102 void AliAnaElectron::MakeAnalysisFillAOD()
1104 //Do photon analysis and fill aods
1107 Double_t v[3] = {0,0,0}; //vertex ;
1108 GetReader()->GetVertex(v);
1110 //Select the Calorimeter of the photon
1111 TObjArray * pl = 0x0;
1112 if(GetCalorimeter() == "PHOS")
1113 pl = GetPHOSClusters();
1114 else if (GetCalorimeter() == "EMCAL")
1115 pl = GetEMCALClusters();
1119 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",GetCalorimeter().Data());
1123 //Init arrays, variables, get number of clusters
1124 TLorentzVector mom, mom2 ;
1125 Int_t nCaloClusters = pl->GetEntriesFast();
1126 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1128 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", GetCalorimeter().Data(), nCaloClusters);
1130 //----------------------------------------------------
1131 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1132 //----------------------------------------------------
1134 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
1136 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1137 //printf("calo %d, %f\n",icalo,calo->E());
1139 //Get the index where the cluster comes, to retrieve the corresponding vertex
1140 Int_t evtIndex = 0 ;
1141 if (GetMixedEvent())
1143 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1144 //Get the vertex and check it is not too large in z
1145 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1148 //Cluster selection, not charged, with photon id and in fiducial cut
1149 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1151 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1153 Double_t vertex[]={0,0,0};
1154 calo->GetMomentum(mom,vertex) ;
1157 //--------------------------------------
1158 // Cluster selection
1159 //--------------------------------------
1160 AliVCaloCells* cells = 0;
1161 if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
1162 else cells = GetPHOSCells();
1164 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1165 if(!ClusterSelected(calo,mom,nMaxima)) continue;
1167 //-------------------------------------
1168 // PID selection via dE/dx
1169 //-------------------------------------
1171 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1175 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
1179 //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1180 Float_t dEdx = track->GetTPCsignal();
1181 Float_t eOverp = calo->E()/track->P();
1183 fhdEdxvsE->Fill(calo->E(), dEdx);
1184 fhdEdxvsP->Fill(track->P(),dEdx);
1186 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1188 fhdEdxvsECutEOverP ->Fill(calo->E(), dEdx);
1189 fhdEdxvsPCutEOverP ->Fill(track->P(),dEdx);
1192 // Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1193 Float_t m02 = calo->GetM02();
1194 if(m02 > 0.1 && m02 < 0.4)
1196 fhdEdxvsECutM02 ->Fill(calo->E(), dEdx);
1197 fhdEdxvsPCutM02 ->Fill(track->P(),dEdx);
1198 fhEOverPvsECutM02->Fill(calo->E(), eOverp);
1199 fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1202 Int_t pid = AliCaloPID::kChargedHadron;
1204 if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1206 fhEOverPvsE->Fill(calo->E(), eOverp);
1207 fhEOverPvsP->Fill(track->P(), eOverp);
1209 if(m02 > 0.1 && m02 < 0.4)
1211 fhEOverPvsECutM02CutdEdx->Fill(calo->E(), eOverp);
1212 fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1215 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1217 pid = AliCaloPID::kElectron;
1222 Int_t pidIndex = 0;// Electron
1223 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1225 //--------------------------------------------------------------------------------------
1226 // Play with the MC stack if available
1227 //--------------------------------------------------------------------------------------
1229 //Check origin of the candidates
1233 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
1236 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",tag);
1238 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1240 fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx);
1241 fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx);
1242 fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1243 fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1245 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1247 fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx);
1248 fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx);
1249 fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1250 fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1252 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1253 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1255 fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1256 fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx);
1257 fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1258 fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1260 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1261 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1263 fhMCdEdxvsE [kmcOtherDecay]->Fill(calo ->E(), dEdx);
1264 fhMCdEdxvsP [kmcOtherDecay]->Fill(track->P(), dEdx);
1265 fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
1266 fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
1268 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1270 fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx);
1271 fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx);
1272 fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1273 fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1275 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1277 fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx);
1278 fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx);
1279 fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1280 fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1283 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1285 fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1286 fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx);
1287 fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1288 fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1290 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1292 fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx);
1293 fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx);
1294 fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1295 fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1297 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1299 fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx);
1300 fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx);
1301 fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1302 fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1304 else if( fhMCE[pidIndex][kmcOther])
1306 fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx);
1307 fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx);
1308 fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1309 fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1311 }// set MC tag and fill Histograms with MC
1313 //---------------------------------
1314 //Fill some shower shape histograms
1315 //---------------------------------
1317 FillShowerShapeHistograms(calo,tag,pid);
1319 if(pid == AliCaloPID::kElectron)
1320 WeightHistograms(calo);
1322 //-----------------------------------------
1323 // PID Shower Shape selection or bit setting
1324 //-----------------------------------------
1326 // Data, PID check on
1329 // Get most probable PID, 2 options check bayesian PID weights or redo PID
1330 // By default, redo PID
1332 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1334 if(fAODParticle == AliCaloPID::kElectron)
1337 if(fAODParticle == 0 )
1338 pid = AliCaloPID::kChargedHadron ;
1340 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - PDG of identified particle %d\n",pid);
1343 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
1346 Float_t maxCellFraction = 0;
1347 Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1348 if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(mom.E(),maxCellFraction);
1350 fhNCellsE[pidIndex] ->Fill(mom.E(),calo->GetNCells());
1351 fhNLME [pidIndex] ->Fill(mom.E(),nMaxima);
1352 fhTimeE [pidIndex] ->Fill(mom.E(),calo->GetTOF()*1.e9);
1355 //----------------------------
1356 //Create AOD for analysis
1357 //----------------------------
1359 //Add AOD with electron/hadron object to aod branch
1360 if ( pid == fAODParticle || fAODParticle == 0 )
1362 AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
1364 //...............................................
1365 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1366 Int_t label = calo->GetLabel();
1367 aodpart.SetLabel(label);
1368 aodpart.SetCaloLabel (calo ->GetID(),-1);
1369 aodpart.SetTrackLabel(track->GetID(),-1);
1371 aodpart.SetDetector(GetCalorimeter());
1372 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1374 aodpart.SetM02(calo->GetM02());
1375 aodpart.SetNLM(nMaxima);
1377 //...............................................
1378 //Set bad channel distance bit
1379 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1380 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1381 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1382 else aodpart.SetDistToBad(0) ;
1383 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1386 aodpart.SetTag(tag);
1389 aodpart.SetIdentifiedParticleType(pid);
1391 AddAODParticle(aodpart);
1396 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1400 //________________________________________________
1401 void AliAnaElectron::MakeAnalysisFillHistograms()
1405 //-------------------------------------------------------------------
1406 // Access MC information in stack if requested, check that it exists.
1407 AliStack * stack = 0x0;
1408 TParticle * primary = 0x0;
1409 TClonesArray * mcparticles = 0x0;
1410 AliAODMCParticle * aodprimary = 0x0;
1414 if(GetReader()->ReadStack())
1416 stack = GetMCStack() ;
1419 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1423 else if(GetReader()->ReadAODMCParticles())
1425 //Get the list of MC particles
1426 mcparticles = GetReader()->GetAODMCParticles();
1429 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1437 Double_t v[3] = {0,0,0}; //vertex ;
1438 GetReader()->GetVertex(v);
1439 //fhVertex->Fill(v[0],v[1],v[2]);
1440 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1442 //----------------------------------
1443 //Loop on stored AOD photons
1444 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1445 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1447 for(Int_t iaod = 0; iaod < naod ; iaod++)
1449 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1450 Int_t pdg = ph->GetIdentifiedParticleType();
1452 Int_t pidIndex = 0;// Electron
1453 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1454 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1457 if(ph->GetDetector() != GetCalorimeter()) continue;
1460 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1462 //................................
1463 //Fill photon histograms
1464 Float_t ptcluster = ph->Pt();
1465 Float_t phicluster = ph->Phi();
1466 Float_t etacluster = ph->Eta();
1467 Float_t ecluster = ph->E();
1469 fhE[pidIndex] ->Fill(ecluster);
1470 fhPt[pidIndex] ->Fill(ptcluster);
1471 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1472 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1473 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1474 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1476 //.......................................
1477 //Play with the MC data if available
1480 //....................................................................
1481 // Access MC information in stack if requested, check that it exists.
1482 Int_t label =ph->GetLabel();
1484 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1489 //Float_t ptprim = 0;
1490 if(GetReader()->ReadStack())
1492 if(label >= stack->GetNtrack())
1494 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1498 primary = stack->Particle(label);
1501 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1505 eprim = primary->Energy();
1506 //ptprim = primary->Pt();
1509 else if(GetReader()->ReadAODMCParticles())
1511 //Check which is the input
1512 if(ph->GetInputFileIndex() == 0)
1514 if(!mcparticles) continue;
1516 if(label >= mcparticles->GetEntriesFast())
1518 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1519 label, mcparticles->GetEntriesFast());
1523 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1529 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1533 eprim = aodprimary->E();
1534 //ptprim = aodprimary->Pt();
1537 Int_t tag =ph->GetTag();
1539 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1541 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1542 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1543 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1544 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1546 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1547 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1549 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1551 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1552 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1553 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1554 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1556 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1557 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1560 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1561 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1563 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1564 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1565 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1566 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1568 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1569 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1571 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1572 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1574 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1575 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1576 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1577 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1579 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1580 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1582 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1584 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1585 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1586 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1587 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1589 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1590 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1593 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1595 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1596 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1597 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1598 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1600 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1601 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1605 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1607 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1608 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1609 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1610 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1612 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1613 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1616 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1618 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1619 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1620 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1621 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1623 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1624 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1627 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1629 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1630 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1631 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1632 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1634 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1635 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1638 else if( fhMCE[pidIndex][kmcOther])
1640 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1641 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1642 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1643 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1645 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1646 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1650 }//Histograms with MC
1656 //____________________________________________________
1657 void AliAnaElectron::Print(const Option_t * opt) const
1659 //Print some relevant parameters set for the analysis
1664 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1665 AliAnaCaloTrackCorrBaseClass::Print(" ");
1667 printf("Calorimeter = %s\n", GetCalorimeter().Data()) ;
1668 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1669 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1670 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1671 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1672 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1673 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1674 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1679 //______________________________________________________
1680 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1682 // Calculate weights and fill histograms
1684 if(!fFillWeightHistograms || GetMixedEvent()) return;
1686 AliVCaloCells* cells = 0;
1687 if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
1688 else cells = GetPHOSCells();
1690 // First recalculate energy in case non linearity was applied
1693 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1695 Int_t id = clus->GetCellsAbsId()[ipos];
1697 //Recalibrate cell energy if needed
1698 Float_t amp = cells->GetCellAmplitude(id);
1699 GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
1709 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1713 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1714 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1715 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1717 //Get the ratio and log ratio to all cells in cluster
1718 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1719 Int_t id = clus->GetCellsAbsId()[ipos];
1721 //Recalibrate cell energy if needed
1722 Float_t amp = cells->GetCellAmplitude(id);
1723 GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
1725 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1726 fhECellClusterRatio ->Fill(energy,amp/energy);
1727 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1730 //Recalculate shower shape for different W0
1731 if(GetCalorimeter()=="EMCAL"){
1733 Float_t l0org = clus->GetM02();
1734 Float_t l1org = clus->GetM20();
1735 Float_t dorg = clus->GetDispersion();
1737 for(Int_t iw = 0; iw < 14; iw++){
1739 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1740 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1742 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1743 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1745 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1749 // Set the original values back
1750 clus->SetM02(l0org);
1751 clus->SetM20(l1org);
1752 clus->SetDispersion(dorg);