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 //...............................................
1375 //Set bad channel distance bit
1376 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1377 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1378 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1379 else aodpart.SetDistToBad(0) ;
1380 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1383 aodpart.SetTag(tag);
1386 aodpart.SetIdentifiedParticleType(pid);
1388 AddAODParticle(aodpart);
1393 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1397 //________________________________________________
1398 void AliAnaElectron::MakeAnalysisFillHistograms()
1402 //-------------------------------------------------------------------
1403 // Access MC information in stack if requested, check that it exists.
1404 AliStack * stack = 0x0;
1405 TParticle * primary = 0x0;
1406 TClonesArray * mcparticles = 0x0;
1407 AliAODMCParticle * aodprimary = 0x0;
1411 if(GetReader()->ReadStack())
1413 stack = GetMCStack() ;
1416 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1420 else if(GetReader()->ReadAODMCParticles())
1422 //Get the list of MC particles
1423 mcparticles = GetReader()->GetAODMCParticles();
1426 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1434 Double_t v[3] = {0,0,0}; //vertex ;
1435 GetReader()->GetVertex(v);
1436 //fhVertex->Fill(v[0],v[1],v[2]);
1437 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1439 //----------------------------------
1440 //Loop on stored AOD photons
1441 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1442 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1444 for(Int_t iaod = 0; iaod < naod ; iaod++)
1446 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1447 Int_t pdg = ph->GetIdentifiedParticleType();
1449 Int_t pidIndex = 0;// Electron
1450 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1451 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1454 if(ph->GetDetector() != GetCalorimeter()) continue;
1457 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1459 //................................
1460 //Fill photon histograms
1461 Float_t ptcluster = ph->Pt();
1462 Float_t phicluster = ph->Phi();
1463 Float_t etacluster = ph->Eta();
1464 Float_t ecluster = ph->E();
1466 fhE[pidIndex] ->Fill(ecluster);
1467 fhPt[pidIndex] ->Fill(ptcluster);
1468 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1469 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1470 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1471 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1473 //.......................................
1474 //Play with the MC data if available
1477 //....................................................................
1478 // Access MC information in stack if requested, check that it exists.
1479 Int_t label =ph->GetLabel();
1481 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1486 //Float_t ptprim = 0;
1487 if(GetReader()->ReadStack())
1489 if(label >= stack->GetNtrack())
1491 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1495 primary = stack->Particle(label);
1498 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1502 eprim = primary->Energy();
1503 //ptprim = primary->Pt();
1506 else if(GetReader()->ReadAODMCParticles())
1508 //Check which is the input
1509 if(ph->GetInputFileIndex() == 0)
1511 if(!mcparticles) continue;
1513 if(label >= mcparticles->GetEntriesFast())
1515 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1516 label, mcparticles->GetEntriesFast());
1520 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1526 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1530 eprim = aodprimary->E();
1531 //ptprim = aodprimary->Pt();
1534 Int_t tag =ph->GetTag();
1536 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1538 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1539 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1540 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1541 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1543 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1544 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1546 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1548 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1549 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1550 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1551 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1553 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1554 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1557 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1558 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1560 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1561 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1562 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1563 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1565 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1566 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1568 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1569 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1571 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1572 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1573 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1574 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1576 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1577 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1579 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1581 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1582 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1583 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1584 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1586 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1587 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1590 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1592 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1593 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1594 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1595 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1597 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1598 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1602 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1604 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1605 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1606 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1607 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1609 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1610 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1613 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1615 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1616 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1617 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1618 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1620 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1621 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1624 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1626 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1627 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1628 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1629 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1631 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1632 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1635 else if( fhMCE[pidIndex][kmcOther])
1637 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1638 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1639 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1640 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1642 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1643 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1647 }//Histograms with MC
1653 //____________________________________________________
1654 void AliAnaElectron::Print(const Option_t * opt) const
1656 //Print some relevant parameters set for the analysis
1661 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1662 AliAnaCaloTrackCorrBaseClass::Print(" ");
1664 printf("Calorimeter = %s\n", GetCalorimeter().Data()) ;
1665 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1666 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1667 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1668 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1669 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1670 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1671 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1676 //______________________________________________________
1677 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1679 // Calculate weights and fill histograms
1681 if(!fFillWeightHistograms || GetMixedEvent()) return;
1683 AliVCaloCells* cells = 0;
1684 if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
1685 else cells = GetPHOSCells();
1687 // First recalculate energy in case non linearity was applied
1690 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1692 Int_t id = clus->GetCellsAbsId()[ipos];
1694 //Recalibrate cell energy if needed
1695 Float_t amp = cells->GetCellAmplitude(id);
1696 GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
1706 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1710 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1711 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1712 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1714 //Get the ratio and log ratio to all cells in cluster
1715 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1716 Int_t id = clus->GetCellsAbsId()[ipos];
1718 //Recalibrate cell energy if needed
1719 Float_t amp = cells->GetCellAmplitude(id);
1720 GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
1722 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1723 fhECellClusterRatio ->Fill(energy,amp/energy);
1724 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1727 //Recalculate shower shape for different W0
1728 if(GetCalorimeter()=="EMCAL"){
1730 Float_t l0org = clus->GetM02();
1731 Float_t l1org = clus->GetM20();
1732 Float_t dorg = clus->GetDispersion();
1734 for(Int_t iw = 0; iw < 14; iw++){
1736 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1737 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1739 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1740 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1742 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1746 // Set the original values back
1747 clus->SetM02(l0org);
1748 clus->SetM20(l1org);
1749 clus->SetDispersion(dorg);