1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
18 // Class for the photon identification.
19 // Clusters from calorimeters are identified as photons
20 // and kept in the AOD. Few histograms produced.
21 // Copy of AliAnaPhoton just add electron id.
23 // -- Author: Gustavo Conesa (LPSC-IN2P3-CRNS)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
30 #include <TClonesArray.h>
31 #include <TObjString.h>
32 //#include <Riostream.h>
33 #include "TParticle.h"
34 #include "TDatabasePDG.h"
35 #include "AliVTrack.h"
37 // --- Analysis system ---
38 #include "AliAnaElectron.h"
39 #include "AliCaloTrackReader.h"
41 #include "AliCaloPID.h"
42 #include "AliMCAnalysisUtils.h"
43 #include "AliFiducialCut.h"
44 #include "AliVCluster.h"
45 #include "AliAODMCParticle.h"
46 #include "AliMixedEvent.h"
49 ClassImp(AliAnaElectron)
51 //________________________________
52 AliAnaElectron::AliAnaElectron() :
53 AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
54 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
55 fTimeCutMin(-1), fTimeCutMax(999999),
56 fNCellsCut(0), fNLMCutMin(-1), fNLMCutMax(10),
57 fFillSSHistograms(kFALSE), fFillOnlySimpleSSHisto(1),
58 fFillWeightHistograms(kFALSE), fNOriginHistograms(8),
59 fdEdxMin(0.), fdEdxMax (200.),
60 fEOverPMin(0), fEOverPMax (2),
63 fhdEdxvsE(0), fhdEdxvsP(0),
64 fhEOverPvsE(0), fhEOverPvsP(0),
65 fhdEdxvsECutM02(0), fhdEdxvsPCutM02(0),
66 fhEOverPvsECutM02(0), fhEOverPvsPCutM02(0),
67 fhdEdxvsECutEOverP(0), fhdEdxvsPCutEOverP(0),
68 fhEOverPvsECutM02CutdEdx(0), fhEOverPvsPCutM02CutdEdx(0),
70 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
71 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
73 // Electron SS MC histograms
74 fhMCElectronELambda0NoOverlap(0),
75 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
78 fhEmbeddedSignalFractionEnergy(0),
79 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
80 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
83 for(Int_t index = 0; index < 2; index++)
85 fhNCellsE [index] = 0;
88 fhMaxCellDiffClusterE[index] = 0;
94 fhEtaPhi05[index] = 0;
96 // Shower shape histograms
100 fhDispETRD[index] = 0;
101 fhLam0ETRD[index] = 0;
102 fhLam1ETRD[index] = 0;
103 fhNCellsLam0LowE [index] = 0;
104 fhNCellsLam0HighE[index] = 0;
105 fhEtaLam0LowE [index] = 0;
106 fhPhiLam0LowE [index] = 0;
107 fhEtaLam0HighE [index] = 0;
108 fhPhiLam0HighE [index] = 0;
110 fhDispEtaE [index] = 0;
111 fhDispPhiE [index] = 0;
112 fhSumEtaE [index] = 0;
113 fhSumPhiE [index] = 0;
114 fhSumEtaPhiE [index] = 0;
115 fhDispEtaPhiDiffE[index] = 0;
116 fhSphericityE [index] = 0;
118 for(Int_t i = 0; i < 10; i++)
120 fhMCPt [index][i] = 0;
121 fhMCE [index][i] = 0;
122 fhMCPhi [index][i] = 0;
123 fhMCEta [index][i] = 0;
124 fhMCDeltaE [index][i] = 0;
125 fhMC2E [index][i] = 0;
128 fhMCEOverPvsE [i] = 0;
129 fhMCEOverPvsP [i] = 0;
132 for(Int_t i = 0; i < 6; i++)
134 fhMCELambda0 [index][i] = 0;
135 fhMCEDispEta [index][i] = 0;
136 fhMCEDispPhi [index][i] = 0;
137 fhMCESumEtaPhi [index][i] = 0;
138 fhMCEDispEtaPhiDiff[index][i] = 0;
139 fhMCESphericity [index][i] = 0;
142 for(Int_t i = 0; i < 5; i++)
144 fhDispEtaDispPhiEBin[index][i] = 0 ;
149 for(Int_t i =0; i < 14; i++)
151 fhLambda0ForW0[i] = 0;
152 //fhLambda1ForW0[i] = 0;
155 //Initialize parameters
160 //____________________________________________________________________________
161 Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
163 //Select clusters if they pass different cuts
165 printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
166 GetReader()->GetEventNumber(),
167 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
169 //.......................................
170 //If too small or big energy, skip it
171 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
172 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
174 //.......................................
175 // TOF cut, BE CAREFUL WITH THIS CUT
176 Double_t tof = calo->GetTOF()*1e9;
177 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
178 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
180 //.......................................
181 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
182 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
184 //.......................................
185 //Check acceptance selection
186 if(IsFiducialCutOn()){
187 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
188 if(! in ) return kFALSE ;
190 if(GetDebug() > 2) printf("Fiducial cut passed \n");
192 //.......................................
193 //Skip not matched clusters with tracks
194 if(!IsTrackMatched(calo, GetReader()->GetInputEvent())) {
195 if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
198 else if(GetDebug() > 2) printf(" Track-matching cut passed \n");
200 //...........................................
201 // skip clusters with too many maxima
202 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
203 if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
205 //.......................................
206 //Check Distance to Bad channel, set bit.
207 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
208 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
209 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
212 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
213 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
216 printf("AliAnaElectron::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
217 GetReader()->GetEventNumber(),
218 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
220 //All checks passed, cluster selected
225 //______________________________________________________________________________________________
226 void AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
229 //Fill cluster Shower Shape histograms
231 if(!fFillSSHistograms || GetMixedEvent()) return;
233 Int_t pidIndex = 0;// Electron
234 if (pidTag == AliCaloPID::kElectron) pidIndex = 0;
235 else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
238 Float_t energy = cluster->E();
239 Int_t ncells = cluster->GetNCells();
240 Float_t lambda0 = cluster->GetM02();
241 Float_t lambda1 = cluster->GetM20();
242 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
244 Float_t l0 = 0., l1 = 0.;
245 Float_t dispp= 0., dEta = 0., dPhi = 0.;
246 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
249 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
251 cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
253 Double_t vertex[]={0,0,0};
254 cluster->GetMomentum(mom,vertex) ;
257 Float_t eta = mom.Eta();
258 Float_t phi = mom.Phi();
259 if(phi < 0) phi+=TMath::TwoPi();
261 fhLam0E[pidIndex] ->Fill(energy,lambda0);
262 fhLam1E[pidIndex] ->Fill(energy,lambda1);
263 fhDispE[pidIndex] ->Fill(energy,disp);
265 if(fCalorimeter == "EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
266 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
268 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
269 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
270 fhDispETRD[pidIndex]->Fill(energy,disp);
273 if(!fFillOnlySimpleSSHisto)
277 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
278 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
279 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
283 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
284 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
285 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
288 if(fCalorimeter == "EMCAL")
290 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
291 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
292 fhDispEtaE [pidIndex]-> Fill(energy,dEta);
293 fhDispPhiE [pidIndex]-> Fill(energy,dPhi);
294 fhSumEtaE [pidIndex]-> Fill(energy,sEta);
295 fhSumPhiE [pidIndex]-> Fill(energy,sPhi);
296 fhSumEtaPhiE [pidIndex]-> Fill(energy,sEtaPhi);
297 fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
298 if(dEta+dPhi>0)fhSphericityE [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
300 if (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
301 else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
302 else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
303 else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
304 else fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
311 AliVCaloCells* cells = 0;
312 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
313 else cells = GetPHOSCells();
315 //Fill histograms to check shape of embedded clusters
316 Float_t fraction = 0;
317 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
319 Float_t clusterE = 0; // recalculate in case corrections applied.
321 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
322 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
324 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
327 //Fraction of total energy due to the embedded signal
330 if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
332 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
334 } // embedded fraction
336 // Check the origin and fill histograms
339 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
340 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
341 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
342 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
346 }//photon no conversion
347 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
348 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
350 index = kmcssElectron;
352 if(!GetReader()->IsEmbeddedClusterSelectionOn())
354 //Check particle overlaps in cluster
356 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
357 Int_t ancPDG = 0, ancStatus = -1;
358 TLorentzVector momentum; TVector3 prodVertex;
361 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
362 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
363 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
367 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
369 else if(noverlaps == 2){
370 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
372 else if(noverlaps > 2){
373 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
376 printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d for ancestor %d!!", noverlaps, ancLabel);
380 //Fill histograms to check shape of embedded clusters
381 if(GetReader()->IsEmbeddedClusterSelectionOn())
385 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
387 else if(fraction > 0.5)
389 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
391 else if(fraction > 0.1)
393 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
397 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
401 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
402 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
404 index = kmcssConversion;
406 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
410 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
419 fhMCELambda0[pidIndex][index] ->Fill(energy, lambda0);
421 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
423 fhMCEDispEta [pidIndex][index]-> Fill(energy,dEta);
424 fhMCEDispPhi [pidIndex][index]-> Fill(energy,dPhi);
425 fhMCESumEtaPhi [pidIndex][index]-> Fill(energy,sEtaPhi);
426 fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
427 if(dEta+dPhi>0)fhMCESphericity [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
434 //_____________________________________________
435 TObjString * AliAnaElectron::GetAnalysisCuts()
437 //Save parameters used for analysis
438 TString parList ; //this will be list of parameters used for this analysis.
439 const Int_t buffersize = 255;
440 char onePar[buffersize] ;
442 snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
444 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
446 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
448 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f \n",fEOverPMin, fEOverPMax) ;
450 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
452 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
454 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
457 //Get parameters set in base class.
458 parList += GetBaseParametersList() ;
460 //Get parameters set in PID class.
461 parList += GetCaloPID()->GetPIDParametersList() ;
463 //Get parameters set in FiducialCut class (not available yet)
464 //parlist += GetFidCut()->GetFidCutParametersList()
466 return new TObjString(parList) ;
469 //_______________________________________________
470 TList * AliAnaElectron::GetCreateOutputObjects()
472 // Create histograms to be saved in output file and
473 // store them in outputContainer
474 TList * outputContainer = new TList() ;
475 outputContainer->SetName("ElectronHistos") ;
477 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
478 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
479 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
480 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
481 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
482 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
483 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
484 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
487 // MC labels, titles, for originator particles
488 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
489 TString pnamess[] = { "Photon","Hadron" ,"Pi0" ,"Eta" ,"Conversion" ,"Electron"} ;
490 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
491 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
493 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
494 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
497 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
498 fhdEdxvsE->SetXTitle("E (GeV)");
499 fhdEdxvsE->SetYTitle("<dE/dx>");
500 outputContainer->Add(fhdEdxvsE);
502 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
503 fhdEdxvsP->SetXTitle("P (GeV/c)");
504 fhdEdxvsP->SetYTitle("<dE/dx>");
505 outputContainer->Add(fhdEdxvsP);
507 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
508 fhEOverPvsE->SetXTitle("E (GeV)");
509 fhEOverPvsE->SetYTitle("E/p");
510 outputContainer->Add(fhEOverPvsE);
512 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
513 fhEOverPvsP->SetXTitle("P (GeV/c)");
514 fhEOverPvsP->SetYTitle("E/p");
515 outputContainer->Add(fhEOverPvsP);
518 fhdEdxvsECutM02 = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
519 fhdEdxvsECutM02->SetXTitle("E (GeV)");
520 fhdEdxvsECutM02->SetYTitle("<dE/dx>");
521 outputContainer->Add(fhdEdxvsECutM02);
523 fhdEdxvsPCutM02 = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
524 fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
525 fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
526 outputContainer->Add(fhdEdxvsPCutM02);
528 fhEOverPvsECutM02 = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
529 fhEOverPvsECutM02->SetXTitle("E (GeV)");
530 fhEOverPvsECutM02->SetYTitle("E/p");
531 outputContainer->Add(fhEOverPvsECutM02);
533 fhEOverPvsPCutM02 = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
534 fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
535 fhEOverPvsPCutM02->SetYTitle("E/p");
536 outputContainer->Add(fhEOverPvsPCutM02);
539 fhdEdxvsECutEOverP = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
540 fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
541 fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
542 outputContainer->Add(fhdEdxvsECutEOverP);
544 fhdEdxvsPCutEOverP = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
545 fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
546 fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
547 outputContainer->Add(fhdEdxvsPCutEOverP);
549 fhEOverPvsECutM02CutdEdx = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
550 fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
551 fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
552 outputContainer->Add(fhEOverPvsECutM02CutdEdx);
554 fhEOverPvsPCutM02CutdEdx = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
555 fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
556 fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
557 outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
561 for(Int_t i = 0; i < fNOriginHistograms; i++)
563 fhMCdEdxvsE[i] = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
564 Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
565 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
566 fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
567 fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
568 outputContainer->Add(fhMCdEdxvsE[i]) ;
570 fhMCdEdxvsP[i] = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
571 Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
572 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
573 fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
574 fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
575 outputContainer->Add(fhMCdEdxvsP[i]) ;
578 fhMCEOverPvsE[i] = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
579 Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
580 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
581 fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
582 fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
583 outputContainer->Add(fhMCEOverPvsE[i]) ;
585 fhMCEOverPvsP[i] = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
586 Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
587 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
588 fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
589 fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
590 outputContainer->Add(fhMCEOverPvsP[i]) ;
595 TString pidParticle[] = {"Electron","ChargedHadron"} ;
597 if(fFillWeightHistograms)
600 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
601 nptbins,ptmin,ptmax, 100,0,1.);
602 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
603 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
604 outputContainer->Add(fhECellClusterRatio);
606 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
607 nptbins,ptmin,ptmax, 100,-10,0);
608 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
609 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
610 outputContainer->Add(fhECellClusterLogRatio);
612 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
613 nptbins,ptmin,ptmax, 100,0,1.);
614 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
615 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
616 outputContainer->Add(fhEMaxCellClusterRatio);
618 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
619 nptbins,ptmin,ptmax, 100,-10,0);
620 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
621 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
622 outputContainer->Add(fhEMaxCellClusterLogRatio);
624 for(Int_t iw = 0; iw < 14; iw++)
626 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
627 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
628 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
629 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
630 outputContainer->Add(fhLambda0ForW0[iw]);
632 // fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
633 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
634 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
635 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
636 // outputContainer->Add(fhLambda1ForW0[iw]);
641 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
644 if(fFillSSHistograms)
646 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
647 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
648 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
649 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
650 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
651 outputContainer->Add(fhLam0E[pidIndex]);
653 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
654 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
655 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
656 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
657 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
658 outputContainer->Add(fhLam1E[pidIndex]);
660 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
661 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
662 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
663 fhDispE[pidIndex]->SetYTitle("D^{2}");
664 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
665 outputContainer->Add(fhDispE[pidIndex]);
667 if(fCalorimeter == "EMCAL" && GetFirstSMCoveredByTRD() >=0 )
669 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
670 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
671 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
672 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
673 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
674 outputContainer->Add(fhLam0ETRD[pidIndex]);
676 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
677 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
678 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
679 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
680 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
681 outputContainer->Add(fhLam1ETRD[pidIndex]);
683 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
684 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
685 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
686 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
687 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
688 outputContainer->Add(fhDispETRD[pidIndex]);
691 if(!fFillOnlySimpleSSHisto)
694 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
695 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
696 nbins,nmin, nmax, ssbins,ssmin,ssmax);
697 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
698 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
699 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
701 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
702 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
703 nbins,nmin, nmax, ssbins,ssmin,ssmax);
704 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
705 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
706 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
709 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
710 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
711 netabins,etamin,etamax, ssbins,ssmin,ssmax);
712 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
713 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
714 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
716 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
717 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
718 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
719 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
720 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
721 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
723 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
724 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
725 netabins,etamin,etamax, ssbins,ssmin,ssmax);
726 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
727 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
728 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
730 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
731 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
732 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
733 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
734 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
735 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
737 if(fCalorimeter == "EMCAL")
739 fhDispEtaE[pidIndex] = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
740 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
741 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
742 fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
743 fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
744 outputContainer->Add(fhDispEtaE[pidIndex]);
746 fhDispPhiE[pidIndex] = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
747 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
748 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
749 fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
750 fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
751 outputContainer->Add(fhDispPhiE[pidIndex]);
753 fhSumEtaE[pidIndex] = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
754 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
755 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
756 fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
757 fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
758 outputContainer->Add(fhSumEtaE[pidIndex]);
760 fhSumPhiE[pidIndex] = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
761 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),
762 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
763 fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
764 fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
765 outputContainer->Add(fhSumPhiE[pidIndex]);
767 fhSumEtaPhiE[pidIndex] = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
768 Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),
769 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
770 fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
771 fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
772 outputContainer->Add(fhSumEtaPhiE[pidIndex]);
774 fhDispEtaPhiDiffE[pidIndex] = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
775 Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
776 nptbins,ptmin,ptmax,200, -10,10);
777 fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
778 fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
779 outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
781 fhSphericityE[pidIndex] = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
782 Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),
783 nptbins,ptmin,ptmax, 200, -1,1);
784 fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
785 fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
786 outputContainer->Add(fhSphericityE[pidIndex]);
788 Int_t bin[] = {0,2,4,6,10,1000};
789 for(Int_t i = 0; i < 5; i++)
791 fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
792 Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
793 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
794 fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
795 fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
796 outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
804 if(fFillSSHistograms)
806 for(Int_t i = 0; i < 6; i++)
808 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
809 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
810 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
811 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
812 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
813 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
815 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
817 fhMCEDispEta[pidIndex][i] = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
818 Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
819 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
820 fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
821 fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
822 outputContainer->Add(fhMCEDispEta[pidIndex][i]);
824 fhMCEDispPhi[pidIndex][i] = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
825 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
826 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
827 fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
828 fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
829 outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
831 fhMCESumEtaPhi[pidIndex][i] = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
832 Form("cluster from %s : %s like, #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
833 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
834 fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
835 fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
836 outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
838 fhMCEDispEtaPhiDiff[pidIndex][i] = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
839 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
840 nptbins,ptmin,ptmax,200,-10,10);
841 fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
842 fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
843 outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
845 fhMCESphericity[pidIndex][i] = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
846 Form("cluster from %s : %s like, (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
847 nptbins,ptmin,ptmax, 200,-1,1);
848 fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
849 fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
850 outputContainer->Add(fhMCESphericity[pidIndex][i]);
857 //if(IsCaloPIDOn() && pidIndex > 0) continue;
859 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
860 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
861 nptbins,ptmin,ptmax, nbins,nmin,nmax);
862 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
863 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
864 outputContainer->Add(fhNCellsE[pidIndex]);
866 fhNLME[pidIndex] = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
867 Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
868 nptbins,ptmin,ptmax, 10,0,10);
869 fhNLME[pidIndex]->SetXTitle("E (GeV)");
870 fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
871 outputContainer->Add(fhNLME[pidIndex]);
873 fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
874 Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
875 ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
876 fhTimeE[pidIndex]->SetXTitle("E (GeV)");
877 fhTimeE[pidIndex]->SetYTitle(" t (ns)");
878 outputContainer->Add(fhTimeE[pidIndex]);
880 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
881 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
882 nptbins,ptmin,ptmax, 500,0,1.);
883 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
884 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
885 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
887 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
888 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
889 nptbins,ptmin,ptmax);
890 fhE[pidIndex]->SetYTitle("N");
891 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
892 outputContainer->Add(fhE[pidIndex]) ;
894 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
895 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
896 nptbins,ptmin,ptmax);
897 fhPt[pidIndex]->SetYTitle("N");
898 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
899 outputContainer->Add(fhPt[pidIndex]) ;
901 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
902 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
903 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
904 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
905 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
906 outputContainer->Add(fhPhi[pidIndex]) ;
908 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
909 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
910 nptbins,ptmin,ptmax,netabins,etamin,etamax);
911 fhEta[pidIndex]->SetYTitle("#eta");
912 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
913 outputContainer->Add(fhEta[pidIndex]) ;
915 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
916 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
917 netabins,etamin,etamax,nphibins,phimin,phimax);
918 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
919 fhEtaPhi[pidIndex]->SetXTitle("#eta");
920 outputContainer->Add(fhEtaPhi[pidIndex]) ;
923 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
924 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
925 netabins,etamin,etamax,nphibins,phimin,phimax);
926 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
927 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
928 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
934 for(Int_t i = 0; i < fNOriginHistograms; i++)
936 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
937 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
938 nptbins,ptmin,ptmax);
939 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
940 outputContainer->Add(fhMCE[pidIndex][i]) ;
942 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
943 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
944 nptbins,ptmin,ptmax);
945 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
946 outputContainer->Add(fhMCPt[pidIndex][i]) ;
948 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
949 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
950 nptbins,ptmin,ptmax,netabins,etamin,etamax);
951 fhMCEta[pidIndex][i]->SetYTitle("#eta");
952 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
953 outputContainer->Add(fhMCEta[pidIndex][i]) ;
955 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
956 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
957 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
958 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
959 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
960 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
963 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
964 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
965 nptbins,ptmin,ptmax, 200,-50,50);
966 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
967 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
969 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
970 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
971 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
972 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
973 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
974 outputContainer->Add(fhMC2E[pidIndex][i]);
982 if(fFillSSHistograms)
986 if(!GetReader()->IsEmbeddedClusterSelectionOn())
988 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
989 "cluster from Electron : E vs #lambda_{0}^{2}",
990 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
991 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
992 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
993 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
995 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
996 "cluster from Electron : E vs #lambda_{0}^{2}",
997 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
998 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
999 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
1000 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
1002 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
1003 "cluster from Electron : E vs #lambda_{0}^{2}",
1004 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1005 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1006 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
1007 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
1011 //Fill histograms to check shape of embedded clusters
1012 if(GetReader()->IsEmbeddedClusterSelectionOn())
1015 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1016 "Energy Fraction of embedded signal versus cluster energy",
1017 nptbins,ptmin,ptmax,100,0.,1.);
1018 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1019 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1020 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1022 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
1023 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1024 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1025 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1026 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1027 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1029 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1030 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1031 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1032 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1033 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1034 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1036 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1037 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1038 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1039 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1040 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1041 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1043 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
1044 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1045 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1046 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1047 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1048 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1051 }// embedded histograms
1055 }// Fill SS MC histograms
1057 return outputContainer ;
1061 //_________________________
1062 void AliAnaElectron::Init()
1067 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1068 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1071 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1072 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1078 //___________________________________
1079 void AliAnaElectron::InitParameters()
1082 //Initialize the parameters of the analysis.
1083 AddToHistogramsName("AnaElectron_");
1085 fCalorimeter = "EMCAL" ;
1091 fTimeCutMax = 9999999;
1094 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1095 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1097 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1098 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1102 //_________________________________________
1103 void AliAnaElectron::MakeAnalysisFillAOD()
1105 //Do photon analysis and fill aods
1108 Double_t v[3] = {0,0,0}; //vertex ;
1109 GetReader()->GetVertex(v);
1111 //Select the Calorimeter of the photon
1112 TObjArray * pl = 0x0;
1113 if(fCalorimeter == "PHOS")
1114 pl = GetPHOSClusters();
1115 else if (fCalorimeter == "EMCAL")
1116 pl = GetEMCALClusters();
1119 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.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", fCalorimeter.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()) {
1142 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1143 //Get the vertex and check it is not too large in z
1144 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1147 //Cluster selection, not charged, with photon id and in fiducial cut
1148 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1149 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1151 Double_t vertex[]={0,0,0};
1152 calo->GetMomentum(mom,vertex) ;
1155 //--------------------------------------
1156 // Cluster selection
1157 //--------------------------------------
1158 AliVCaloCells* cells = 0;
1159 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1160 else cells = GetPHOSCells();
1162 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1163 if(!ClusterSelected(calo,mom,nMaxima)) continue;
1165 //----------------------------
1166 //Create AOD for analysis
1167 //----------------------------
1168 AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
1170 //...............................................
1171 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1172 Int_t label = calo->GetLabel();
1173 aodpart.SetLabel(label);
1174 aodpart.SetCaloLabel(calo->GetID(),-1);
1175 aodpart.SetDetector(fCalorimeter);
1176 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1178 //...............................................
1179 //Set bad channel distance bit
1180 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1181 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1182 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1183 else aodpart.SetDistToBad(0) ;
1184 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1186 //-------------------------------------
1187 //PID selection via dEdx
1188 //-------------------------------------
1190 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1194 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
1198 //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1199 Float_t dEdx = track->GetTPCsignal();
1200 Float_t eOverp = calo->E()/track->P();
1202 fhdEdxvsE->Fill(calo->E(), dEdx);
1203 fhdEdxvsP->Fill(track->P(),dEdx);
1205 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1207 fhdEdxvsECutEOverP ->Fill(calo->E(), dEdx);
1208 fhdEdxvsPCutEOverP ->Fill(track->P(),dEdx);
1211 //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1212 Float_t m02 = calo->GetM02();
1213 if(m02 > 0.1 && m02 < 0.4)
1215 fhdEdxvsECutM02 ->Fill(calo->E(), dEdx);
1216 fhdEdxvsPCutM02 ->Fill(track->P(),dEdx);
1217 fhEOverPvsECutM02->Fill(calo->E(), eOverp);
1218 fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1221 Int_t pid = AliCaloPID::kChargedHadron;
1223 if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1225 fhEOverPvsE->Fill(calo->E(), eOverp);
1226 fhEOverPvsP->Fill(track->P(), eOverp);
1228 if(m02 > 0.1 && m02 < 0.4)
1230 fhEOverPvsECutM02CutdEdx->Fill(calo->E(), eOverp);
1231 fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1234 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1236 pid = AliCaloPID::kElectron;
1241 aodpart.SetIdentifiedParticleType(pid);
1243 Int_t pidIndex = 0;// Electron
1244 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1246 //--------------------------------------------------------------------------------------
1247 //Play with the MC stack if available
1248 //--------------------------------------------------------------------------------------
1250 //Check origin of the candidates
1253 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
1254 aodpart.SetTag(tag);
1257 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodpart.GetTag());
1259 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1261 fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx);
1262 fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx);
1263 fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1264 fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1266 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1268 fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx);
1269 fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx);
1270 fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1271 fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1273 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1274 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1276 fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1277 fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx);
1278 fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1279 fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1281 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1282 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1284 fhMCdEdxvsE [kmcOtherDecay]->Fill(calo ->E(), dEdx);
1285 fhMCdEdxvsP [kmcOtherDecay]->Fill(track->P(), dEdx);
1286 fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
1287 fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
1289 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1291 fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx);
1292 fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx);
1293 fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1294 fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1296 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1298 fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx);
1299 fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx);
1300 fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1301 fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1304 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1306 fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1307 fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx);
1308 fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1309 fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1311 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1313 fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx);
1314 fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx);
1315 fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1316 fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1318 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1320 fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx);
1321 fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx);
1322 fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1323 fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1325 else if( fhMCE[pidIndex][kmcOther])
1327 fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx);
1328 fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx);
1329 fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1330 fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1332 }// set MC tag and fill Histograms with MC
1334 //---------------------------------
1335 //Fill some shower shape histograms
1336 //---------------------------------
1338 FillShowerShapeHistograms(calo,aodpart.GetTag(),pid);
1340 if(pid == AliCaloPID::kElectron)
1341 WeightHistograms(calo);
1343 //-----------------------------------------
1344 //PID Shower Shape selection or bit setting
1345 //-----------------------------------------
1347 // Data, PID check on
1350 // Get most probable PID, 2 options check bayesian PID weights or redo PID
1351 // By default, redo PID
1353 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1357 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodpart.GetIdentifiedParticleType());
1361 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
1362 aodpart.Pt(), aodpart.GetIdentifiedParticleType());
1364 //FIXME, this to MakeAnalysisFillHistograms ...
1365 Float_t maxCellFraction = 0;
1367 Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1368 if(absID>=0)fhMaxCellDiffClusterE[pidIndex]->Fill(aodpart.E(),maxCellFraction);
1369 fhNCellsE[pidIndex] ->Fill(aodpart.E(),calo->GetNCells());
1370 fhNLME[pidIndex] ->Fill(aodpart.E(),nMaxima);
1371 fhTimeE[pidIndex] ->Fill(aodpart.E(),calo->GetTOF()*1.e9);
1373 //Add AOD with electron/hadron object to aod branch
1374 if ( pid == fAODParticle || fAODParticle == 0 )
1376 AddAODParticle(aodpart);
1381 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1385 //________________________________________________
1386 void AliAnaElectron::MakeAnalysisFillHistograms()
1390 //-------------------------------------------------------------------
1391 // Access MC information in stack if requested, check that it exists.
1392 AliStack * stack = 0x0;
1393 TParticle * primary = 0x0;
1394 TClonesArray * mcparticles = 0x0;
1395 AliAODMCParticle * aodprimary = 0x0;
1399 if(GetReader()->ReadStack()){
1400 stack = GetMCStack() ;
1402 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1407 else if(GetReader()->ReadAODMCParticles()){
1409 //Get the list of MC particles
1410 mcparticles = GetReader()->GetAODMCParticles();
1411 if(!mcparticles && GetDebug() > 0) {
1412 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1419 Double_t v[3] = {0,0,0}; //vertex ;
1420 GetReader()->GetVertex(v);
1421 //fhVertex->Fill(v[0],v[1],v[2]);
1422 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1424 //----------------------------------
1425 //Loop on stored AOD photons
1426 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1427 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1429 for(Int_t iaod = 0; iaod < naod ; iaod++)
1431 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1432 Int_t pdg = ph->GetIdentifiedParticleType();
1434 Int_t pidIndex = 0;// Electron
1435 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1436 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1439 if(ph->GetDetector() != fCalorimeter) continue;
1442 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1444 //................................
1445 //Fill photon histograms
1446 Float_t ptcluster = ph->Pt();
1447 Float_t phicluster = ph->Phi();
1448 Float_t etacluster = ph->Eta();
1449 Float_t ecluster = ph->E();
1451 fhE[pidIndex] ->Fill(ecluster);
1452 fhPt[pidIndex] ->Fill(ptcluster);
1453 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1454 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1455 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1456 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1458 //.......................................
1459 //Play with the MC data if available
1462 //....................................................................
1463 // Access MC information in stack if requested, check that it exists.
1464 Int_t label =ph->GetLabel();
1466 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1471 //Float_t ptprim = 0;
1472 if(GetReader()->ReadStack()){
1474 if(label >= stack->GetNtrack()) {
1475 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1479 primary = stack->Particle(label);
1481 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1485 eprim = primary->Energy();
1486 //ptprim = primary->Pt();
1489 else if(GetReader()->ReadAODMCParticles()){
1490 //Check which is the input
1491 if(ph->GetInputFileIndex() == 0){
1492 if(!mcparticles) continue;
1493 if(label >= mcparticles->GetEntriesFast()) {
1494 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1495 label, mcparticles->GetEntriesFast());
1499 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1504 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1508 eprim = aodprimary->E();
1509 //ptprim = aodprimary->Pt();
1513 Int_t tag =ph->GetTag();
1515 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1517 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1518 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1519 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1520 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1522 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1523 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1525 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1527 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1528 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1529 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1530 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1532 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1533 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1536 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1537 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1539 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1540 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1541 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1542 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1544 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1545 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1547 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1548 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1550 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1551 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1552 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1553 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1555 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1556 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1558 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1560 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1561 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1562 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1563 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1565 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1566 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1569 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1571 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1572 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1573 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1574 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1576 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1577 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1581 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1583 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1584 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1585 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1586 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1588 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1589 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1592 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1594 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1595 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1596 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1597 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1599 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1600 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1603 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1605 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1606 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1607 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1608 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1610 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1611 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1614 else if( fhMCE[pidIndex][kmcOther]){
1615 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1616 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1617 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1618 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1620 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1621 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1625 }//Histograms with MC
1631 //____________________________________________________
1632 void AliAnaElectron::Print(const Option_t * opt) const
1634 //Print some relevant parameters set for the analysis
1639 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1640 AliAnaCaloTrackCorrBaseClass::Print(" ");
1642 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1643 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1644 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1645 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1646 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1647 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1648 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1649 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1654 //______________________________________________________
1655 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1657 // Calculate weights and fill histograms
1659 if(!fFillWeightHistograms || GetMixedEvent()) return;
1661 AliVCaloCells* cells = 0;
1662 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1663 else cells = GetPHOSCells();
1665 // First recalculate energy in case non linearity was applied
1668 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1670 Int_t id = clus->GetCellsAbsId()[ipos];
1672 //Recalibrate cell energy if needed
1673 Float_t amp = cells->GetCellAmplitude(id);
1674 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
1684 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1688 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1689 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1690 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1692 //Get the ratio and log ratio to all cells in cluster
1693 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1694 Int_t id = clus->GetCellsAbsId()[ipos];
1696 //Recalibrate cell energy if needed
1697 Float_t amp = cells->GetCellAmplitude(id);
1698 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
1700 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1701 fhECellClusterRatio ->Fill(energy,amp/energy);
1702 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1705 //Recalculate shower shape for different W0
1706 if(fCalorimeter=="EMCAL"){
1708 Float_t l0org = clus->GetM02();
1709 Float_t l1org = clus->GetM20();
1710 Float_t dorg = clus->GetDispersion();
1712 for(Int_t iw = 0; iw < 14; iw++){
1714 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1715 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1717 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1718 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1720 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1724 // Set the original values back
1725 clus->SetM02(l0org);
1726 clus->SetM20(l1org);
1727 clus->SetDispersion(dorg);