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, const Int_t mcTag, const 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" && GetModuleNumber(cluster) > 5)
267 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
268 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
269 fhDispETRD[pidIndex]->Fill(energy,disp);
272 if(!fFillOnlySimpleSSHisto)
276 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
277 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
278 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
282 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
283 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
284 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
287 if(fCalorimeter == "EMCAL")
289 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
290 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
291 fhDispEtaE [pidIndex]-> Fill(energy,dEta);
292 fhDispPhiE [pidIndex]-> Fill(energy,dPhi);
293 fhSumEtaE [pidIndex]-> Fill(energy,sEta);
294 fhSumPhiE [pidIndex]-> Fill(energy,sPhi);
295 fhSumEtaPhiE [pidIndex]-> Fill(energy,sEtaPhi);
296 fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
297 if(dEta+dPhi>0)fhSphericityE [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
299 if (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
300 else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
301 else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
302 else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
303 else fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
310 AliVCaloCells* cells = 0;
311 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
312 else cells = GetPHOSCells();
314 //Fill histograms to check shape of embedded clusters
315 Float_t fraction = 0;
316 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
318 Float_t clusterE = 0; // recalculate in case corrections applied.
320 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
321 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
323 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
326 //Fraction of total energy due to the embedded signal
329 if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
331 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
333 } // embedded fraction
335 // Get the fraction of the cluster energy that carries the cell with highest energy
337 Float_t maxCellFraction = 0.;
339 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
341 // Check the origin and fill histograms
342 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
343 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
344 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
345 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
349 }//photon no conversion
350 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
351 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
353 index = kmcssElectron;
355 if(!GetReader()->IsEmbeddedClusterSelectionOn())
357 //Check particle overlaps in cluster
359 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
360 Int_t ancPDG = 0, ancStatus = -1;
361 TLorentzVector momentum; TVector3 prodVertex;
364 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
365 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
366 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
370 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
372 else if(noverlaps == 2){
373 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
375 else if(noverlaps > 2){
376 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
379 printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
383 //Fill histograms to check shape of embedded clusters
384 if(GetReader()->IsEmbeddedClusterSelectionOn())
388 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
390 else if(fraction > 0.5)
392 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
394 else if(fraction > 0.1)
396 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
400 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
404 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
405 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
407 index = kmcssConversion;
409 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
413 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
422 fhMCELambda0[pidIndex][index] ->Fill(energy, lambda0);
424 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
426 fhMCEDispEta [pidIndex][index]-> Fill(energy,dEta);
427 fhMCEDispPhi [pidIndex][index]-> Fill(energy,dPhi);
428 fhMCESumEtaPhi [pidIndex][index]-> Fill(energy,sEtaPhi);
429 fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
430 if(dEta+dPhi>0)fhMCESphericity [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
439 //_____________________________________________
440 TObjString * AliAnaElectron::GetAnalysisCuts()
442 //Save parameters used for analysis
443 TString parList ; //this will be list of parameters used for this analysis.
444 const Int_t buffersize = 255;
445 char onePar[buffersize] ;
447 snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
449 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
451 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
453 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f \n",fEOverPMin, fEOverPMax) ;
455 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
457 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
459 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
462 //Get parameters set in base class.
463 parList += GetBaseParametersList() ;
465 //Get parameters set in PID class.
466 parList += GetCaloPID()->GetPIDParametersList() ;
468 //Get parameters set in FiducialCut class (not available yet)
469 //parlist += GetFidCut()->GetFidCutParametersList()
471 return new TObjString(parList) ;
474 //_______________________________________________
475 TList * AliAnaElectron::GetCreateOutputObjects()
477 // Create histograms to be saved in output file and
478 // store them in outputContainer
479 TList * outputContainer = new TList() ;
480 outputContainer->SetName("ElectronHistos") ;
482 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
483 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
484 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
485 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
486 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
487 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
488 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
489 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
492 // MC labels, titles, for originator particles
493 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
494 TString pnamess[] = { "Photon","Hadron" ,"Pi0" ,"Eta" ,"Conversion" ,"Electron"} ;
495 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
496 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
498 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
499 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
502 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
503 fhdEdxvsE->SetXTitle("E (GeV)");
504 fhdEdxvsE->SetYTitle("<dE/dx>");
505 outputContainer->Add(fhdEdxvsE);
507 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
508 fhdEdxvsP->SetXTitle("P (GeV/c)");
509 fhdEdxvsP->SetYTitle("<dE/dx>");
510 outputContainer->Add(fhdEdxvsP);
512 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
513 fhEOverPvsE->SetXTitle("E (GeV)");
514 fhEOverPvsE->SetYTitle("E/p");
515 outputContainer->Add(fhEOverPvsE);
517 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
518 fhEOverPvsP->SetXTitle("P (GeV/c)");
519 fhEOverPvsP->SetYTitle("E/p");
520 outputContainer->Add(fhEOverPvsP);
523 fhdEdxvsECutM02 = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
524 fhdEdxvsECutM02->SetXTitle("E (GeV)");
525 fhdEdxvsECutM02->SetYTitle("<dE/dx>");
526 outputContainer->Add(fhdEdxvsECutM02);
528 fhdEdxvsPCutM02 = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
529 fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
530 fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
531 outputContainer->Add(fhdEdxvsPCutM02);
533 fhEOverPvsECutM02 = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
534 fhEOverPvsECutM02->SetXTitle("E (GeV)");
535 fhEOverPvsECutM02->SetYTitle("E/p");
536 outputContainer->Add(fhEOverPvsECutM02);
538 fhEOverPvsPCutM02 = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
539 fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
540 fhEOverPvsPCutM02->SetYTitle("E/p");
541 outputContainer->Add(fhEOverPvsPCutM02);
544 fhdEdxvsECutEOverP = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
545 fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
546 fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
547 outputContainer->Add(fhdEdxvsECutEOverP);
549 fhdEdxvsPCutEOverP = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
550 fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
551 fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
552 outputContainer->Add(fhdEdxvsPCutEOverP);
554 fhEOverPvsECutM02CutdEdx = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
555 fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
556 fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
557 outputContainer->Add(fhEOverPvsECutM02CutdEdx);
559 fhEOverPvsPCutM02CutdEdx = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
560 fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
561 fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
562 outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
566 for(Int_t i = 0; i < fNOriginHistograms; i++)
568 fhMCdEdxvsE[i] = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
569 Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
570 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
571 fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
572 fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
573 outputContainer->Add(fhMCdEdxvsE[i]) ;
575 fhMCdEdxvsP[i] = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
576 Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
577 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
578 fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
579 fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
580 outputContainer->Add(fhMCdEdxvsP[i]) ;
583 fhMCEOverPvsE[i] = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
584 Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
585 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
586 fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
587 fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
588 outputContainer->Add(fhMCEOverPvsE[i]) ;
590 fhMCEOverPvsP[i] = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
591 Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
592 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
593 fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
594 fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
595 outputContainer->Add(fhMCEOverPvsP[i]) ;
600 TString pidParticle[] = {"Electron","ChargedHadron"} ;
602 if(fFillWeightHistograms)
605 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
606 nptbins,ptmin,ptmax, 100,0,1.);
607 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
608 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
609 outputContainer->Add(fhECellClusterRatio);
611 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
612 nptbins,ptmin,ptmax, 100,-10,0);
613 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
614 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
615 outputContainer->Add(fhECellClusterLogRatio);
617 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
618 nptbins,ptmin,ptmax, 100,0,1.);
619 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
620 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
621 outputContainer->Add(fhEMaxCellClusterRatio);
623 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
624 nptbins,ptmin,ptmax, 100,-10,0);
625 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
626 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
627 outputContainer->Add(fhEMaxCellClusterLogRatio);
629 for(Int_t iw = 0; iw < 14; iw++)
631 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),
632 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
633 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
634 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
635 outputContainer->Add(fhLambda0ForW0[iw]);
637 // 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),
638 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
639 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
640 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
641 // outputContainer->Add(fhLambda1ForW0[iw]);
646 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
649 if(fFillSSHistograms)
651 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
652 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
653 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
654 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
655 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
656 outputContainer->Add(fhLam0E[pidIndex]);
658 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
659 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
660 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
661 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
662 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
663 outputContainer->Add(fhLam1E[pidIndex]);
665 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
666 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
667 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
668 fhDispE[pidIndex]->SetYTitle("D^{2}");
669 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
670 outputContainer->Add(fhDispE[pidIndex]);
672 if(fCalorimeter == "EMCAL")
674 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
675 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
676 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
677 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
678 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
679 outputContainer->Add(fhLam0ETRD[pidIndex]);
681 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
682 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
683 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
684 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
685 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
686 outputContainer->Add(fhLam1ETRD[pidIndex]);
688 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
689 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
690 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
691 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
692 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
693 outputContainer->Add(fhDispETRD[pidIndex]);
696 if(!fFillOnlySimpleSSHisto)
699 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
700 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
701 nbins,nmin, nmax, ssbins,ssmin,ssmax);
702 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
703 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
704 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
706 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
707 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
708 nbins,nmin, nmax, ssbins,ssmin,ssmax);
709 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
710 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
711 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
714 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
715 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
716 netabins,etamin,etamax, ssbins,ssmin,ssmax);
717 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
718 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
719 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
721 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
722 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
723 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
724 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
725 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
726 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
728 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
729 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
730 netabins,etamin,etamax, ssbins,ssmin,ssmax);
731 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
732 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
733 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
735 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
736 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
737 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
738 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
739 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
740 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
742 if(fCalorimeter == "EMCAL")
744 fhDispEtaE[pidIndex] = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
745 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
746 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
747 fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
748 fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
749 outputContainer->Add(fhDispEtaE[pidIndex]);
751 fhDispPhiE[pidIndex] = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
752 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
753 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
754 fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
755 fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
756 outputContainer->Add(fhDispPhiE[pidIndex]);
758 fhSumEtaE[pidIndex] = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
759 Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
760 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
761 fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
762 fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
763 outputContainer->Add(fhSumEtaE[pidIndex]);
765 fhSumPhiE[pidIndex] = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
766 Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),
767 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
768 fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
769 fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
770 outputContainer->Add(fhSumPhiE[pidIndex]);
772 fhSumEtaPhiE[pidIndex] = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
773 Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),
774 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
775 fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
776 fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
777 outputContainer->Add(fhSumEtaPhiE[pidIndex]);
779 fhDispEtaPhiDiffE[pidIndex] = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
780 Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
781 nptbins,ptmin,ptmax,200, -10,10);
782 fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
783 fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
784 outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
786 fhSphericityE[pidIndex] = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
787 Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),
788 nptbins,ptmin,ptmax, 200, -1,1);
789 fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
790 fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
791 outputContainer->Add(fhSphericityE[pidIndex]);
793 Int_t bin[] = {0,2,4,6,10,1000};
794 for(Int_t i = 0; i < 5; i++)
796 fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
797 Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
798 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
799 fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
800 fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
801 outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
809 if(fFillSSHistograms)
811 for(Int_t i = 0; i < 6; i++)
813 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
814 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
815 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
816 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
817 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
818 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
820 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
822 fhMCEDispEta[pidIndex][i] = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
823 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()),
824 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
825 fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
826 fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
827 outputContainer->Add(fhMCEDispEta[pidIndex][i]);
829 fhMCEDispPhi[pidIndex][i] = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
830 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()),
831 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
832 fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
833 fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
834 outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
836 fhMCESumEtaPhi[pidIndex][i] = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
837 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()),
838 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
839 fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
840 fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
841 outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
843 fhMCEDispEtaPhiDiff[pidIndex][i] = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
844 Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
845 nptbins,ptmin,ptmax,200,-10,10);
846 fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
847 fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
848 outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
850 fhMCESphericity[pidIndex][i] = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
851 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()),
852 nptbins,ptmin,ptmax, 200,-1,1);
853 fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
854 fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
855 outputContainer->Add(fhMCESphericity[pidIndex][i]);
862 //if(IsCaloPIDOn() && pidIndex > 0) continue;
864 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
865 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
866 nptbins,ptmin,ptmax, nbins,nmin,nmax);
867 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
868 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
869 outputContainer->Add(fhNCellsE[pidIndex]);
871 fhNLME[pidIndex] = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
872 Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
873 nptbins,ptmin,ptmax, 10,0,10);
874 fhNLME[pidIndex]->SetXTitle("E (GeV)");
875 fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
876 outputContainer->Add(fhNLME[pidIndex]);
878 fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
879 Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
880 ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
881 fhTimeE[pidIndex]->SetXTitle("E (GeV)");
882 fhTimeE[pidIndex]->SetYTitle(" t (ns)");
883 outputContainer->Add(fhTimeE[pidIndex]);
885 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
886 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
887 nptbins,ptmin,ptmax, 500,0,1.);
888 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
889 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
890 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
892 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
893 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
894 nptbins,ptmin,ptmax);
895 fhE[pidIndex]->SetYTitle("N");
896 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
897 outputContainer->Add(fhE[pidIndex]) ;
899 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
900 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
901 nptbins,ptmin,ptmax);
902 fhPt[pidIndex]->SetYTitle("N");
903 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
904 outputContainer->Add(fhPt[pidIndex]) ;
906 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
907 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
908 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
909 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
910 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
911 outputContainer->Add(fhPhi[pidIndex]) ;
913 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
914 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
915 nptbins,ptmin,ptmax,netabins,etamin,etamax);
916 fhEta[pidIndex]->SetYTitle("#eta");
917 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
918 outputContainer->Add(fhEta[pidIndex]) ;
920 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
921 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
922 netabins,etamin,etamax,nphibins,phimin,phimax);
923 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
924 fhEtaPhi[pidIndex]->SetXTitle("#eta");
925 outputContainer->Add(fhEtaPhi[pidIndex]) ;
928 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
929 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
930 netabins,etamin,etamax,nphibins,phimin,phimax);
931 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
932 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
933 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
939 for(Int_t i = 0; i < fNOriginHistograms; i++)
941 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
942 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
943 nptbins,ptmin,ptmax);
944 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
945 outputContainer->Add(fhMCE[pidIndex][i]) ;
947 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
948 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
949 nptbins,ptmin,ptmax);
950 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
951 outputContainer->Add(fhMCPt[pidIndex][i]) ;
953 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
954 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
955 nptbins,ptmin,ptmax,netabins,etamin,etamax);
956 fhMCEta[pidIndex][i]->SetYTitle("#eta");
957 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
958 outputContainer->Add(fhMCEta[pidIndex][i]) ;
960 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
961 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
962 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
963 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
964 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
965 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
968 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
969 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
970 nptbins,ptmin,ptmax, 200,-50,50);
971 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
972 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
974 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
975 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
976 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
977 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
978 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
979 outputContainer->Add(fhMC2E[pidIndex][i]);
987 if(fFillSSHistograms)
991 if(!GetReader()->IsEmbeddedClusterSelectionOn())
993 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
994 "cluster from Electron : E vs #lambda_{0}^{2}",
995 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
996 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
997 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
998 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
1000 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
1001 "cluster from Electron : E vs #lambda_{0}^{2}",
1002 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1003 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1004 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
1005 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
1007 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
1008 "cluster from Electron : E vs #lambda_{0}^{2}",
1009 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1010 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1011 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
1012 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
1016 //Fill histograms to check shape of embedded clusters
1017 if(GetReader()->IsEmbeddedClusterSelectionOn())
1020 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1021 "Energy Fraction of embedded signal versus cluster energy",
1022 nptbins,ptmin,ptmax,100,0.,1.);
1023 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1024 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1025 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1027 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
1028 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1029 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1030 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1031 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1032 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1034 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1035 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1036 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1037 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1038 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1039 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1041 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1042 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1043 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1044 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1045 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1046 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1048 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
1049 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1050 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1051 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1052 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1053 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1056 }// embedded histograms
1060 }// Fill SS MC histograms
1062 return outputContainer ;
1066 //_________________________
1067 void AliAnaElectron::Init()
1072 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1073 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1076 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1077 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1083 //___________________________________
1084 void AliAnaElectron::InitParameters()
1087 //Initialize the parameters of the analysis.
1088 AddToHistogramsName("AnaElectron_");
1090 fCalorimeter = "EMCAL" ;
1096 fTimeCutMax = 9999999;
1099 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1100 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1102 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1103 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1107 //_________________________________________
1108 void AliAnaElectron::MakeAnalysisFillAOD()
1110 //Do photon analysis and fill aods
1113 Double_t v[3] = {0,0,0}; //vertex ;
1114 GetReader()->GetVertex(v);
1116 //Select the Calorimeter of the photon
1117 TObjArray * pl = 0x0;
1118 if(fCalorimeter == "PHOS")
1119 pl = GetPHOSClusters();
1120 else if (fCalorimeter == "EMCAL")
1121 pl = GetEMCALClusters();
1124 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1128 //Init arrays, variables, get number of clusters
1129 TLorentzVector mom, mom2 ;
1130 Int_t nCaloClusters = pl->GetEntriesFast();
1131 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1133 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
1135 //----------------------------------------------------
1136 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1137 //----------------------------------------------------
1139 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
1141 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1142 //printf("calo %d, %f\n",icalo,calo->E());
1144 //Get the index where the cluster comes, to retrieve the corresponding vertex
1145 Int_t evtIndex = 0 ;
1146 if (GetMixedEvent()) {
1147 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1148 //Get the vertex and check it is not too large in z
1149 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1152 //Cluster selection, not charged, with photon id and in fiducial cut
1153 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1154 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1156 Double_t vertex[]={0,0,0};
1157 calo->GetMomentum(mom,vertex) ;
1160 //--------------------------------------
1161 // Cluster selection
1162 //--------------------------------------
1163 AliVCaloCells* cells = 0;
1164 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1165 else cells = GetPHOSCells();
1167 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1168 if(!ClusterSelected(calo,mom,nMaxima)) continue;
1170 //----------------------------
1171 //Create AOD for analysis
1172 //----------------------------
1173 AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
1175 //...............................................
1176 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1177 Int_t label = calo->GetLabel();
1178 aodpart.SetLabel(label);
1179 aodpart.SetCaloLabel(calo->GetID(),-1);
1180 aodpart.SetDetector(fCalorimeter);
1181 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1183 //...............................................
1184 //Set bad channel distance bit
1185 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1186 if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1187 else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1188 else aodpart.SetDistToBad(0) ;
1189 //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1191 //-------------------------------------
1192 //PID selection via dEdx
1193 //-------------------------------------
1195 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1199 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
1203 //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1204 Float_t dEdx = track->GetTPCsignal();
1205 Float_t eOverp = calo->E()/track->P();
1207 fhdEdxvsE->Fill(calo->E(), dEdx);
1208 fhdEdxvsP->Fill(track->P(),dEdx);
1210 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1212 fhdEdxvsECutEOverP ->Fill(calo->E(), dEdx);
1213 fhdEdxvsPCutEOverP ->Fill(track->P(),dEdx);
1216 //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1217 Float_t m02 = calo->GetM02();
1218 if(m02 > 0.1 && m02 < 0.4)
1220 fhdEdxvsECutM02 ->Fill(calo->E(), dEdx);
1221 fhdEdxvsPCutM02 ->Fill(track->P(),dEdx);
1222 fhEOverPvsECutM02->Fill(calo->E(), eOverp);
1223 fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1226 Int_t pid = AliCaloPID::kChargedHadron;
1228 if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1230 fhEOverPvsE->Fill(calo->E(), eOverp);
1231 fhEOverPvsP->Fill(track->P(), eOverp);
1233 if(m02 > 0.1 && m02 < 0.4)
1235 fhEOverPvsECutM02CutdEdx->Fill(calo->E(), eOverp);
1236 fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1239 if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1241 pid = AliCaloPID::kElectron;
1246 aodpart.SetIdentifiedParticleType(pid);
1248 Int_t pidIndex = 0;// Electron
1249 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1251 //--------------------------------------------------------------------------------------
1252 //Play with the MC stack if available
1253 //--------------------------------------------------------------------------------------
1255 //Check origin of the candidates
1258 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
1259 aodpart.SetTag(tag);
1262 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodpart.GetTag());
1264 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1266 fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx);
1267 fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx);
1268 fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1269 fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1271 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1273 fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx);
1274 fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx);
1275 fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1276 fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1278 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1279 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1281 fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1282 fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx);
1283 fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1284 fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1286 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1287 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1289 fhMCdEdxvsE [kmcOtherDecay]->Fill(calo ->E(), dEdx);
1290 fhMCdEdxvsP [kmcOtherDecay]->Fill(track->P(), dEdx);
1291 fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
1292 fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
1294 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1296 fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx);
1297 fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx);
1298 fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1299 fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1301 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1303 fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx);
1304 fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx);
1305 fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1306 fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1309 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1311 fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1312 fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx);
1313 fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1314 fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1316 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1318 fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx);
1319 fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx);
1320 fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1321 fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1323 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1325 fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx);
1326 fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx);
1327 fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1328 fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1330 else if( fhMCE[pidIndex][kmcOther])
1332 fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx);
1333 fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx);
1334 fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1335 fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1337 }// set MC tag and fill Histograms with MC
1339 //---------------------------------
1340 //Fill some shower shape histograms
1341 //---------------------------------
1343 FillShowerShapeHistograms(calo,aodpart.GetTag(),pid);
1345 if(pid == AliCaloPID::kElectron)
1346 WeightHistograms(calo);
1348 //-----------------------------------------
1349 //PID Shower Shape selection or bit setting
1350 //-----------------------------------------
1352 // Data, PID check on
1355 // Get most probable PID, 2 options check bayesian PID weights or redo PID
1356 // By default, redo PID
1358 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1362 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodpart.GetIdentifiedParticleType());
1366 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
1367 aodpart.Pt(), aodpart.GetIdentifiedParticleType());
1369 //FIXME, this to MakeAnalysisFillHistograms ...
1371 Float_t maxCellFraction = 0;
1373 absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1374 fhMaxCellDiffClusterE[pidIndex]->Fill(aodpart.E(),maxCellFraction);
1375 fhNCellsE[pidIndex] ->Fill(aodpart.E(),calo->GetNCells());
1376 fhNLME[pidIndex] ->Fill(aodpart.E(),nMaxima);
1377 fhTimeE[pidIndex] ->Fill(aodpart.E(),calo->GetTOF()*1.e9);
1379 //Add AOD with electron/hadron object to aod branch
1380 if ( pid == fAODParticle || fAODParticle == 0 )
1382 AddAODParticle(aodpart);
1387 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1391 //________________________________________________
1392 void AliAnaElectron::MakeAnalysisFillHistograms()
1396 //-------------------------------------------------------------------
1397 // Access MC information in stack if requested, check that it exists.
1398 AliStack * stack = 0x0;
1399 TParticle * primary = 0x0;
1400 TClonesArray * mcparticles = 0x0;
1401 AliAODMCParticle * aodprimary = 0x0;
1405 if(GetReader()->ReadStack()){
1406 stack = GetMCStack() ;
1408 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1413 else if(GetReader()->ReadAODMCParticles()){
1415 //Get the list of MC particles
1416 mcparticles = GetReader()->GetAODMCParticles();
1417 if(!mcparticles && GetDebug() > 0) {
1418 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1425 Double_t v[3] = {0,0,0}; //vertex ;
1426 GetReader()->GetVertex(v);
1427 //fhVertex->Fill(v[0],v[1],v[2]);
1428 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1430 //----------------------------------
1431 //Loop on stored AOD photons
1432 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1433 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1435 for(Int_t iaod = 0; iaod < naod ; iaod++)
1437 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1438 Int_t pdg = ph->GetIdentifiedParticleType();
1440 Int_t pidIndex = 0;// Electron
1441 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1442 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1445 if(ph->GetDetector() != fCalorimeter) continue;
1448 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1450 //................................
1451 //Fill photon histograms
1452 Float_t ptcluster = ph->Pt();
1453 Float_t phicluster = ph->Phi();
1454 Float_t etacluster = ph->Eta();
1455 Float_t ecluster = ph->E();
1457 fhE[pidIndex] ->Fill(ecluster);
1458 fhPt[pidIndex] ->Fill(ptcluster);
1459 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1460 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1461 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1462 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1464 //.......................................
1465 //Play with the MC data if available
1468 //....................................................................
1469 // Access MC information in stack if requested, check that it exists.
1470 Int_t label =ph->GetLabel();
1472 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1478 if(GetReader()->ReadStack()){
1480 if(label >= stack->GetNtrack()) {
1481 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1485 primary = stack->Particle(label);
1487 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1491 eprim = primary->Energy();
1492 ptprim = primary->Pt();
1495 else if(GetReader()->ReadAODMCParticles()){
1496 //Check which is the input
1497 if(ph->GetInputFileIndex() == 0){
1498 if(!mcparticles) continue;
1499 if(label >= mcparticles->GetEntriesFast()) {
1500 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1501 label, mcparticles->GetEntriesFast());
1505 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1510 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1514 eprim = aodprimary->E();
1515 ptprim = aodprimary->Pt();
1519 Int_t tag =ph->GetTag();
1521 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1523 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1524 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1525 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1526 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1528 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1529 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1531 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1533 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1534 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1535 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1536 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1538 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1539 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1542 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1543 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1545 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1546 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1547 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1548 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1550 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1551 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1553 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1554 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1556 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1557 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1558 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1559 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1561 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1562 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1564 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1566 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1567 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1568 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1569 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1571 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1572 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1575 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1577 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1578 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1579 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1580 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1582 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1583 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1587 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1589 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1590 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1591 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1592 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1594 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1595 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1598 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1600 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1601 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1602 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1603 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1605 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1606 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1609 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1611 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1612 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1613 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1614 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1616 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1617 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1620 else if( fhMCE[pidIndex][kmcOther]){
1621 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1622 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1623 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1624 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1626 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1627 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1631 }//Histograms with MC
1637 //____________________________________________________
1638 void AliAnaElectron::Print(const Option_t * opt) const
1640 //Print some relevant parameters set for the analysis
1645 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1646 AliAnaCaloTrackCorrBaseClass::Print(" ");
1648 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1649 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1650 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1651 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1652 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1653 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1654 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1655 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1660 //______________________________________________________
1661 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1663 // Calculate weights and fill histograms
1665 if(!fFillWeightHistograms || GetMixedEvent()) return;
1667 AliVCaloCells* cells = 0;
1668 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1669 else cells = GetPHOSCells();
1671 // First recalculate energy in case non linearity was applied
1674 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1676 Int_t id = clus->GetCellsAbsId()[ipos];
1678 //Recalibrate cell energy if needed
1679 Float_t amp = cells->GetCellAmplitude(id);
1680 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
1690 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1694 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1695 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1696 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1698 //Get the ratio and log ratio to all cells in cluster
1699 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1700 Int_t id = clus->GetCellsAbsId()[ipos];
1702 //Recalibrate cell energy if needed
1703 Float_t amp = cells->GetCellAmplitude(id);
1704 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
1706 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1707 fhECellClusterRatio ->Fill(energy,amp/energy);
1708 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1711 //Recalculate shower shape for different W0
1712 if(fCalorimeter=="EMCAL"){
1714 Float_t l0org = clus->GetM02();
1715 Float_t l1org = clus->GetM20();
1716 Float_t dorg = clus->GetDispersion();
1718 for(Int_t iw = 0; iw < 14; iw++){
1720 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1721 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1723 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1724 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1726 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1730 // Set the original values back
1731 clus->SetM02(l0org);
1732 clus->SetM20(l1org);
1733 clus->SetDispersion(dorg);