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 is 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 //_________________________________________________________________________
17 // Class for the analysis of high pT pi0 event by event
18 // Pi0/Eta identified by one of the following:
19 // -Invariant mass of 2 cluster in calorimeter
20 // -Shower shape analysis in calorimeter
21 // -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
23 // -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
29 #include <TClonesArray.h>
30 #include <TObjString.h>
32 // --- Analysis system ---
33 #include "AliAnaPi0EbE.h"
34 #include "AliCaloTrackReader.h"
35 #include "AliIsolationCut.h"
36 #include "AliNeutralMesonSelection.h"
37 #include "AliCaloPID.h"
38 #include "AliMCAnalysisUtils.h"
40 #include "AliFiducialCut.h"
41 #include "TParticle.h"
42 #include "AliVCluster.h"
43 #include "AliESDEvent.h"
44 #include "AliAODEvent.h"
45 #include "AliAODMCParticle.h"
47 ClassImp(AliAnaPi0EbE)
49 //____________________________
50 AliAnaPi0EbE::AliAnaPi0EbE() :
51 AliAnaCaloTrackCorrBaseClass(), fAnaType(kIMCalo), fCalorimeter(""),
52 fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
53 fNLMCutMin(-1), fNLMCutMax(10),
54 fTimeCutMin(-10000), fTimeCutMax(10000),
55 fRejectTrackMatch(kTRUE),
56 fFillPileUpHistograms(0),
57 fFillWeightHistograms(kFALSE), fFillTMHisto(0),
58 fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1), fFillEMCALBCHistograms(0),
59 fInputAODGammaConvName(""),
60 fCheckSplitDistToBad(0), fNSuperModules(0),
63 fhPtEta(0), fhPtPhi(0), fhEtaPhi(0),
64 fhEtaPhiEMCALBC0(0), fhEtaPhiEMCALBC1(0), fhEtaPhiEMCALBCN(0),
65 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
66 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
67 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
68 fhPtCentrality(), fhPtEventPlane(0),
69 fhPtReject(0), fhEReject(0),
70 fhPtEtaReject(0), fhPtPhiReject(0), fhEtaPhiReject(0),
71 fhMass(0), fhMassPt(0), fhMassSplitPt(0),
72 fhSelectedMass(), fhSelectedMassPt(0), fhSelectedMassSplitPt(0),
73 fhMassNoOverlap(0), fhMassPtNoOverlap(0), fhMassSplitPtNoOverlap(0),
74 fhSelectedMassNoOverlap(0), fhSelectedMassPtNoOverlap(0), fhSelectedMassSplitPtNoOverlap(0),
75 fhMCPi0PtRecoPtPrim(0), fhMCEtaPtRecoPtPrim(0),
76 fhMCPi0PtRecoPtPrimNoOverlap(0), fhMCEtaPtRecoPtPrimNoOverlap(0),
77 fhMCPi0SplitPtRecoPtPrim(0), fhMCEtaSplitPtRecoPtPrim(0),
78 fhMCPi0SplitPtRecoPtPrimNoOverlap(0), fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
79 fhMCPi0SelectedPtRecoPtPrim(0), fhMCEtaSelectedPtRecoPtPrim(0),
80 fhMCPi0SelectedPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
81 fhMCPi0SelectedSplitPtRecoPtPrim(0), fhMCEtaSelectedSplitPtRecoPtPrim(0),
82 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
83 fhAsymmetry(0), fhSelectedAsymmetry(0),
84 fhSplitE(0), fhSplitPt(0),
85 fhSplitPtEta(0), fhSplitPtPhi(0),
87 fhPtDecay(0), fhEDecay(0),
88 // Shower shape histos
89 fhPtDispersion(0), fhPtLambda0(0), fhPtLambda1(0),
90 fhPtLambda0NoTRD(0), fhPtLambda0FracMaxCellCut(0),
91 fhPtFracMaxCell(0), fhPtFracMaxCellNoTRD(0),
92 fhPtNCells(0), fhPtTime(0), fhEPairDiffTime(0),
93 fhPtDispEta(0), fhPtDispPhi(0),
94 fhPtSumEta(0), fhPtSumPhi(0), fhPtSumEtaPhi(0),
95 fhPtDispEtaPhiDiff(0), fhPtSphericity(0),
99 fhMCPtPhi(), fhMCPtEta(),
100 fhMCEReject(), fhMCPtReject(),
102 fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
103 fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
104 fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
106 fhMassPairMCPi0(0), fhMassPairMCEta(0),
107 fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
109 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
110 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
111 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
112 fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
113 fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
114 fhTrackMatchedMCParticlePt(0),
115 fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
116 fhdEdx(0), fhEOverP(0), fhEOverPNoTRD(0),
117 // Number of local maxima in cluster
118 fhNLocMaxPt(0), fhNLocMaxPtReject(0),
120 fhTimePtNoCut(0), fhTimePtSPD(0), fhTimePtSPDMulti(0),
121 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
122 fhTimeNPileUpVertContributors(0),
123 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
124 fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
125 fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
126 fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
130 for(Int_t i = 0; i < 6; i++)
136 fhMCPtCentrality [i] = 0;
140 fhMCSplitPtPhi [i] = 0;
141 fhMCSplitPtEta [i] = 0;
143 fhMCNLocMaxPt [i] = 0;
144 fhMCNLocMaxSplitPt [i] = 0;
145 fhMCNLocMaxPtReject[i] = 0;
147 fhMCPtLambda0 [i] = 0;
148 fhMCPtLambda0NoTRD [i] = 0;
149 fhMCPtLambda0FracMaxCellCut[i]= 0;
150 fhMCPtFracMaxCell [i] = 0;
151 fhMCPtLambda1 [i] = 0;
152 fhMCPtDispersion [i] = 0;
154 fhMCPtDispEta [i] = 0;
155 fhMCPtDispPhi [i] = 0;
156 fhMCPtSumEtaPhi [i] = 0;
157 fhMCPtDispEtaPhiDiff[i] = 0;
158 fhMCPtSphericity [i] = 0;
159 fhMCPtAsymmetry [i] = 0;
162 fhMCMassSplitPt [i]=0;
163 fhMCSelectedMassPt [i]=0;
164 fhMCSelectedMassSplitPt[i]=0;
166 fhMCMassPtNoOverlap [i]=0;
167 fhMCMassSplitPtNoOverlap [i]=0;
168 fhMCSelectedMassPtNoOverlap [i]=0;
169 fhMCSelectedMassSplitPtNoOverlap[i]=0;
171 for(Int_t j = 0; j < 7; j++)
173 fhMCLambda0DispEta [j][i] = 0;
174 fhMCLambda0DispPhi [j][i] = 0;
175 fhMCDispEtaDispPhi [j][i] = 0;
176 fhMCAsymmetryLambda0 [j][i] = 0;
177 fhMCAsymmetryDispEta [j][i] = 0;
178 fhMCAsymmetryDispPhi [j][i] = 0;
182 for(Int_t j = 0; j < 7; j++)
184 fhLambda0DispEta [j] = 0;
185 fhLambda0DispPhi [j] = 0;
186 fhDispEtaDispPhi [j] = 0;
187 fhAsymmetryLambda0 [j] = 0;
188 fhAsymmetryDispEta [j] = 0;
189 fhAsymmetryDispPhi [j] = 0;
194 for(Int_t i = 0; i < 3; i++)
196 fhPtLambda0LocMax [i] = 0;
197 fhPtLambda1LocMax [i] = 0;
198 fhPtDispersionLocMax [i] = 0;
199 fhPtDispEtaLocMax [i] = 0;
200 fhPtDispPhiLocMax [i] = 0;
201 fhPtSumEtaPhiLocMax [i] = 0;
202 fhPtDispEtaPhiDiffLocMax[i] = 0;
203 fhPtSphericityLocMax [i] = 0;
204 fhPtAsymmetryLocMax [i] = 0;
205 fhMassPtLocMax [i] = 0;
206 fhSelectedMassPtLocMax [i] = 0;
207 for(Int_t ipart = 0; ipart<6; ipart++)
209 fhMCPtLambda0LocMax [ipart][i] = 0;
210 fhMCSelectedMassPtLocMax[ipart][i] = 0;
213 fhMCPi0PtRecoPtPrimLocMax [i] = 0;
214 fhMCEtaPtRecoPtPrimLocMax [i] = 0;
215 fhMCPi0SplitPtRecoPtPrimLocMax [i] = 0;
216 fhMCEtaSplitPtRecoPtPrimLocMax [i] = 0;
218 fhMCPi0SelectedPtRecoPtPrimLocMax [i] = 0;
219 fhMCEtaSelectedPtRecoPtPrimLocMax [i] = 0;
220 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[i] = 0;
221 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[i] = 0;
226 for(Int_t i =0; i < 14; i++){
227 fhLambda0ForW0[i] = 0;
228 //fhLambda1ForW0[i] = 0;
229 if(i<8)fhMassPairLocMax[i] = 0;
232 for(Int_t i = 0; i < 11; i++)
234 fhEtaPhiTriggerEMCALBC [i] = 0 ;
235 fhTimeTriggerEMCALBC [i] = 0 ;
236 fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
238 fhEtaPhiTriggerEMCALBCUM [i] = 0 ;
239 fhTimeTriggerEMCALBCUM [i] = 0 ;
243 for(Int_t iSM = 0; iSM < 22; iSM++)
245 fhNLocMaxPtSM[iSM] = 0;
246 for(Int_t inlm = 0; inlm < 3; inlm++)
248 fhSelectedMassPtLocMaxSM [inlm][iSM] = 0;
249 fhSelectedLambda0PtLocMaxSM [inlm][iSM] = 0;
252 //Initialize parameters
257 //___________________________________________________________________________________
258 void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
260 // Fill some histograms to understand pile-up
261 if(!fFillPileUpHistograms) return;
263 //printf("E %f, time %f\n",energy,time);
264 AliVEvent * event = GetReader()->GetInputEvent();
266 fhTimePtNoCut->Fill(pt,time);
267 if(GetReader()->IsPileUpFromSPD())
269 if(GetReader()->IsPileUpFromSPD()) { fhPtPileUp[0]->Fill(pt); fhTimePtSPD ->Fill(pt,time); }
270 if(GetReader()->IsPileUpFromEMCal()) fhPtPileUp[1]->Fill(pt);
271 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPileUp[2]->Fill(pt);
272 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPileUp[3]->Fill(pt);
273 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPileUp[4]->Fill(pt);
274 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPileUp[5]->Fill(pt);
275 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
277 if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
281 AliVCaloCells* cells = 0;
282 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
283 else cells = GetPHOSCells();
285 Float_t maxCellFraction = 0.;
286 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
288 Double_t tmax = cells->GetCellTime(absIdMax);
289 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
292 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
293 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
295 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
297 Int_t absId = calo->GetCellsAbsId()[ipos];
299 if( absId == absIdMax ) continue ;
301 Double_t timecell = cells->GetCellTime(absId);
302 Float_t amp = cells->GetCellAmplitude(absId);
303 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
304 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
307 Float_t diff = (tmax-timecell);
309 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
311 if(GetReader()->IsPileUpFromSPD())
313 fhPtCellTimePileUp[0]->Fill(pt, timecell);
314 fhPtTimeDiffPileUp[0]->Fill(pt, diff);
317 if(GetReader()->IsPileUpFromEMCal())
319 fhPtCellTimePileUp[1]->Fill(pt, timecell);
320 fhPtTimeDiffPileUp[1]->Fill(pt, diff);
323 if(GetReader()->IsPileUpFromSPDOrEMCal())
325 fhPtCellTimePileUp[2]->Fill(pt, timecell);
326 fhPtTimeDiffPileUp[2]->Fill(pt, diff);
329 if(GetReader()->IsPileUpFromSPDAndEMCal())
331 fhPtCellTimePileUp[3]->Fill(pt, timecell);
332 fhPtTimeDiffPileUp[3]->Fill(pt, diff);
335 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
337 fhPtCellTimePileUp[4]->Fill(pt, timecell);
338 fhPtTimeDiffPileUp[4]->Fill(pt, diff);
341 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
343 fhPtCellTimePileUp[5]->Fill(pt, timecell);
344 fhPtTimeDiffPileUp[5]->Fill(pt, diff);
347 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
349 fhPtCellTimePileUp[6]->Fill(pt, timecell);
350 fhPtTimeDiffPileUp[6]->Fill(pt, diff);
355 if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
357 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
358 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
360 // N pile up vertices
366 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
367 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
372 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
373 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
376 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
377 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
379 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
380 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
382 if(TMath::Abs(time) < 25)
384 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
385 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
388 if(time < 75 && time > -25)
390 fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
391 fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
394 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
395 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
398 Float_t z1 = -1, z2 = -1;
400 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
404 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
405 ncont=pv->GetNContributors();
406 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
408 diamZ = esdEv->GetDiamondZ();
412 AliAODVertex *pv=aodEv->GetVertex(iVert);
413 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
414 ncont=pv->GetNContributors();
415 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
417 diamZ = aodEv->GetDiamondZ();
420 Double_t distZ = TMath::Abs(z2-z1);
421 diamZ = TMath::Abs(z2-diamZ);
423 fhTimeNPileUpVertContributors ->Fill(time,ncont);
424 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
425 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
431 //______________________________________________________________________________________________
432 void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
434 // Fill histograms that do not pass the identification (SS case only)
436 Float_t ener = mom.E();
437 Float_t pt = mom.Pt();
438 Float_t phi = mom.Phi();
439 if(phi < 0) phi+=TMath::TwoPi();
440 Float_t eta = mom.Eta();
442 fhPtReject ->Fill(pt);
443 fhEReject ->Fill(ener);
445 fhPtEtaReject ->Fill(ener,eta);
446 fhPtPhiReject ->Fill(ener,phi);
447 fhEtaPhiReject ->Fill(eta,phi);
449 fhNLocMaxPtReject->Fill(pt,nMaxima);
453 Int_t mcIndex = GetMCIndex(mctag);
454 fhMCEReject [mcIndex] ->Fill(ener);
455 fhMCPtReject [mcIndex] ->Fill(pt);
456 fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
460 //___________________________________________________________________________________
461 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
462 Int_t tag, Float_t asy)
464 // Fill shower shape, timing and other histograms for selected clusters from decay
466 Float_t ener = cluster->E();
467 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
468 Float_t l0 = cluster->GetM02();
469 Float_t l1 = cluster->GetM20();
470 Int_t nSM = GetModuleNumber(cluster);
473 if (pt < 2 ) ptbin = 0;
474 else if (pt < 4 ) ptbin = 1;
475 else if (pt < 6 ) ptbin = 2;
476 else if (pt < 10) ptbin = 3;
477 else if (pt < 15) ptbin = 4;
478 else if (pt < 20) ptbin = 5;
482 if (nMaxima==1) indexMax = 0 ;
483 else if(nMaxima==2) indexMax = 1 ;
487 AliVCaloCells * cell = 0x0;
488 if(fCalorimeter == "PHOS")
489 cell = GetPHOSCells();
491 cell = GetEMCALCells();
493 Float_t maxCellFraction = 0;
494 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
495 fhPtFracMaxCell->Fill(pt,maxCellFraction);
497 FillWeightHistograms(cluster);
499 fhPtDispersion->Fill(pt, disp);
500 fhPtLambda0 ->Fill(pt, l0 );
501 fhPtLambda1 ->Fill(pt, l1 );
503 Float_t ll0 = 0., ll1 = 0.;
504 Float_t dispp= 0., dEta = 0., dPhi = 0.;
505 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
506 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
508 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
509 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
511 fhPtDispEta -> Fill(pt,dEta);
512 fhPtDispPhi -> Fill(pt,dPhi);
513 fhPtSumEta -> Fill(pt,sEta);
514 fhPtSumPhi -> Fill(pt,sPhi);
515 fhPtSumEtaPhi -> Fill(pt,sEtaPhi);
516 fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
517 if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
519 fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
520 fhLambda0DispEta[ptbin]->Fill(l0 ,dEta);
521 fhLambda0DispPhi[ptbin]->Fill(l0 ,dPhi);
523 if (fAnaType==kSSCalo)
525 // Asymmetry histograms
526 fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
527 fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
528 fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
532 fhNLocMaxPt->Fill(pt,nMaxima);
534 if(nSM < fNSuperModules && nSM >=0)
535 fhNLocMaxPtSM[nSM]->Fill(pt,nMaxima);
537 fhPtLambda0LocMax [indexMax]->Fill(pt,l0);
538 fhPtLambda1LocMax [indexMax]->Fill(pt,l1);
539 fhPtDispersionLocMax[indexMax]->Fill(pt,disp);
541 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
543 fhPtDispEtaLocMax [indexMax]-> Fill(pt,dEta);
544 fhPtDispPhiLocMax [indexMax]-> Fill(pt,dPhi);
545 fhPtSumEtaPhiLocMax [indexMax]-> Fill(pt,sEtaPhi);
546 fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
547 if(dEta+dPhi>0) fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
548 if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt ,asy);
552 if(fCalorimeter=="EMCAL" && nSM < 6) // CAREFUL FOR 2012-13 runs change 6 to 4, -1 for 2015 ...
554 fhPtLambda0NoTRD ->Fill(pt, l0 );
555 fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
558 if(maxCellFraction < 0.5)
559 fhPtLambda0FracMaxCellCut->Fill(pt, l0 );
561 fhPtTime ->Fill(pt, cluster->GetTOF()*1.e9);
562 fhPtNCells->Fill(pt, cluster->GetNCells());
564 // Fill Track matching control histograms
567 Float_t dZ = cluster->GetTrackDz();
568 Float_t dR = cluster->GetTrackDx();
570 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
572 dR = 2000., dZ = 2000.;
573 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
575 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
577 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
579 Bool_t positive = kFALSE;
580 if(track) positive = (track->Charge()>0);
582 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
584 fhTrackMatchedDEta->Fill(pt,dZ);
585 fhTrackMatchedDPhi->Fill(pt,dR);
586 if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
592 fhTrackMatchedDEtaPos->Fill(pt,dZ);
593 fhTrackMatchedDPhiPos->Fill(pt,dR);
594 if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
598 fhTrackMatchedDEtaNeg->Fill(pt,dZ);
599 fhTrackMatchedDPhiNeg->Fill(pt,dR);
600 if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
604 // Check dEdx and E/p of matched clusters
606 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
610 Float_t dEdx = track->GetTPCsignal();
611 fhdEdx->Fill(pt, dEdx);
613 Float_t eOverp = cluster->E()/track->P();
614 fhEOverP->Fill(pt, eOverp);
616 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
617 if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(pt, eOverp);
621 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
626 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
628 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
629 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 2.5 ;
630 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 0.5 ;
631 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 1.5 ;
637 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
638 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 6.5 ;
639 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 4.5 ;
640 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 5.5 ;
644 fhTrackMatchedMCParticlePt ->Fill(pt, mctag);
645 fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
646 fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
650 }// Track matching histograms
654 Int_t mcIndex = GetMCIndex(tag);
656 fhMCPtLambda0[mcIndex] ->Fill(pt, l0);
657 fhMCPtLambda1[mcIndex] ->Fill(pt, l1);
658 fhMCPtDispersion[mcIndex] ->Fill(pt, disp);
659 fhMCPtFracMaxCell[mcIndex]->Fill(pt,maxCellFraction);
661 fhMCPtLambda0LocMax [mcIndex][indexMax]->Fill(pt,l0);
663 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
664 if(fCalorimeter=="EMCAL" && nSM < 6)
665 fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0 );
667 if(maxCellFraction < 0.5)
668 fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0 );
670 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
672 fhMCPtDispEta [mcIndex]-> Fill(pt,dEta);
673 fhMCPtDispPhi [mcIndex]-> Fill(pt,dPhi);
674 fhMCPtSumEtaPhi [mcIndex]-> Fill(pt,sEtaPhi);
675 fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
676 if(dEta+dPhi>0)fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
678 if (fAnaType==kSSCalo)
680 fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
681 fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
682 fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
685 fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
686 fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0 ,dEta);
687 fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0 ,dPhi);
695 //________________________________________________________
696 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
698 // Calculate weights and fill histograms
700 if(!fFillWeightHistograms || GetMixedEvent()) return;
702 AliVCaloCells* cells = 0;
703 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
704 else cells = GetPHOSCells();
706 // First recalculate energy in case non linearity was applied
709 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
712 Int_t id = clus->GetCellsAbsId()[ipos];
714 //Recalibrate cell energy if needed
715 Float_t amp = cells->GetCellAmplitude(id);
716 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
727 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
731 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
732 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
734 //Get the ratio and log ratio to all cells in cluster
735 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
737 Int_t id = clus->GetCellsAbsId()[ipos];
739 //Recalibrate cell energy if needed
740 Float_t amp = cells->GetCellAmplitude(id);
741 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
743 fhECellClusterRatio ->Fill(energy,amp/energy);
744 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
747 //Recalculate shower shape for different W0
748 if(fCalorimeter=="EMCAL"){
750 Float_t l0org = clus->GetM02();
751 Float_t l1org = clus->GetM20();
752 Float_t dorg = clus->GetDispersion();
754 for(Int_t iw = 0; iw < 14; iw++)
756 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
757 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
759 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
760 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
764 // Set the original values back
767 clus->SetDispersion(dorg);
772 //__________________________________________
773 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
775 //Save parameters used for analysis
776 TString parList ; //this will be list of parameters used for this analysis.
777 const Int_t buffersize = 255;
778 char onePar[buffersize] ;
780 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
782 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
785 if(fAnaType == kSSCalo)
787 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
789 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
791 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
793 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
797 //Get parameters set in base class.
798 parList += GetBaseParametersList() ;
800 //Get parameters set in PID class.
801 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
803 return new TObjString(parList) ;
806 //_____________________________________________
807 TList * AliAnaPi0EbE::GetCreateOutputObjects()
809 // Create histograms to be saved in output file and
810 // store them in outputContainer
811 TList * outputContainer = new TList() ;
812 outputContainer->SetName("Pi0EbEHistos") ;
814 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
815 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
816 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
817 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
818 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
819 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
820 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
822 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
823 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
824 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
826 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
827 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
828 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
829 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
830 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
831 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
833 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
834 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
835 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
836 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
837 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
838 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
840 Int_t ntimptbins= GetHistogramRanges()->GetHistoTimeBins();
841 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
842 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
844 TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
845 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
846 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
847 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
849 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
850 fhPt->SetYTitle("N");
851 fhPt->SetXTitle("p_{T} (GeV/c)");
852 outputContainer->Add(fhPt) ;
854 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
856 fhE->SetXTitle("E (GeV)");
857 outputContainer->Add(fhE) ;
860 ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
861 fhPtPhi->SetYTitle("#phi (rad)");
862 fhPtPhi->SetXTitle("E (GeV)");
863 outputContainer->Add(fhPtPhi) ;
866 ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
867 fhPtEta->SetYTitle("#eta");
868 fhPtEta->SetXTitle("E (GeV)");
869 outputContainer->Add(fhPtEta) ;
872 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
873 fhEtaPhi->SetYTitle("#phi (rad)");
874 fhEtaPhi->SetXTitle("#eta");
875 outputContainer->Add(fhEtaPhi) ;
877 if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
879 fhEtaPhiEMCALBC0 = new TH2F
880 ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
881 fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
882 fhEtaPhiEMCALBC0->SetXTitle("#eta");
883 outputContainer->Add(fhEtaPhiEMCALBC0) ;
885 fhEtaPhiEMCALBC1 = new TH2F
886 ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
887 fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
888 fhEtaPhiEMCALBC1->SetXTitle("#eta");
889 outputContainer->Add(fhEtaPhiEMCALBC1) ;
891 fhEtaPhiEMCALBCN = new TH2F
892 ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
893 fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
894 fhEtaPhiEMCALBCN->SetXTitle("#eta");
895 outputContainer->Add(fhEtaPhiEMCALBCN) ;
897 for(Int_t i = 0; i < 11; i++)
899 fhEtaPhiTriggerEMCALBC[i] = new TH2F
900 (Form("hEtaPhiTriggerEMCALBC%d",i-5),
901 Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
902 netabins,etamin,etamax,nphibins,phimin,phimax);
903 fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
904 fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
905 outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
907 fhTimeTriggerEMCALBC[i] = new TH2F
908 (Form("hTimeTriggerEMCALBC%d",i-5),
909 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
910 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
911 fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
912 fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
913 outputContainer->Add(fhTimeTriggerEMCALBC[i]);
915 fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
916 (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
917 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
918 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
919 fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
920 fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
921 outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
923 fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
924 (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
925 Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
926 netabins,etamin,etamax,nphibins,phimin,phimax);
927 fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
928 fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
929 outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
931 fhTimeTriggerEMCALBCUM[i] = new TH2F
932 (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
933 Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
934 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
935 fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
936 fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
937 outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
941 fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
942 "cluster time vs E of clusters, no match, rematch open time",
943 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
944 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("E (GeV)");
945 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("time (ns)");
946 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
949 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
950 "cluster time vs E of clusters, no match, rematch with neigbour parches",
951 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
952 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("E (GeV)");
953 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("time (ns)");
954 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
956 fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
957 "cluster time vs E of clusters, no match, rematch open time and neigbour",
958 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
959 fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("E (GeV)");
960 fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("time (ns)");
961 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
965 fhPtCentrality = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
966 fhPtCentrality->SetYTitle("centrality");
967 fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
968 outputContainer->Add(fhPtCentrality) ;
970 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
971 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
972 fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
973 outputContainer->Add(fhPtEventPlane) ;
975 if(fAnaType == kSSCalo)
977 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
978 fhPtReject->SetYTitle("N");
979 fhPtReject->SetXTitle("p_{T} (GeV/c)");
980 outputContainer->Add(fhPtReject) ;
982 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
983 fhEReject->SetYTitle("N");
984 fhEReject->SetXTitle("E (GeV)");
985 outputContainer->Add(fhEReject) ;
987 fhPtPhiReject = new TH2F
988 ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
989 fhPtPhiReject->SetYTitle("#phi (rad)");
990 fhPtPhiReject->SetXTitle("p_{T} (GeV/c)");
991 outputContainer->Add(fhPtPhiReject) ;
993 fhPtEtaReject = new TH2F
994 ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
995 fhPtEtaReject->SetYTitle("#eta");
996 fhPtEtaReject->SetXTitle("p_{T} (GeV/c)");
997 outputContainer->Add(fhPtEtaReject) ;
999 fhEtaPhiReject = new TH2F
1000 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1001 fhEtaPhiReject->SetYTitle("#phi (rad)");
1002 fhEtaPhiReject->SetXTitle("#eta");
1003 outputContainer->Add(fhEtaPhiReject) ;
1007 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1008 fhMass->SetYTitle("mass (GeV/c^{2})");
1009 fhMass->SetXTitle("E (GeV)");
1010 outputContainer->Add(fhMass) ;
1012 fhSelectedMass = new TH2F
1013 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1014 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
1015 fhSelectedMass->SetXTitle("E (GeV)");
1016 outputContainer->Add(fhSelectedMass) ;
1018 if(fAnaType == kSSCalo)
1022 ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1023 fhMassPt->SetYTitle("mass (GeV/c^{2})");
1024 fhMassPt->SetXTitle("p_{T} (GeV/c)");
1025 outputContainer->Add(fhMassPt) ;
1027 fhSelectedMassPt = new TH2F
1028 ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1029 fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
1030 fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
1031 outputContainer->Add(fhSelectedMassPt) ;
1033 for(Int_t inlm = 0; inlm < 3; inlm++)
1035 fhMassPtLocMax[inlm] = new TH2F
1036 (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs mass: p_{T} vs mass and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1037 fhMassPtLocMax[inlm]->SetYTitle("mass (GeV/c^{2})");
1038 fhMassPtLocMax[inlm]->SetXTitle("p_{T} (GeV/c)");
1039 outputContainer->Add(fhMassPtLocMax[inlm]) ;
1041 fhSelectedMassPtLocMax[inlm] = new TH2F
1042 (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1043 fhSelectedMassPtLocMax[inlm]->SetYTitle("mass (GeV/c^{2})");
1044 fhSelectedMassPtLocMax[inlm]->SetXTitle("p_{T} (GeV/c)");
1045 outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
1047 for(Int_t iSM = 0; iSM < fNSuperModules; iSM++)
1049 fhSelectedMassPtLocMaxSM[inlm][iSM] = new TH2F
1050 (Form("hSelectedMassPtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1051 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetYTitle("mass (GeV/c^{2})");
1052 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetXTitle("p_{T} (GeV/c)");
1053 outputContainer->Add(fhSelectedMassPtLocMaxSM[inlm][iSM]) ;
1055 fhSelectedLambda0PtLocMaxSM[inlm][iSM] = new TH2F
1056 (Form("hSelectedLambda0PtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #lambda_{0}^{2}: p_{T} vs mass, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1057 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetYTitle("#lambda_{0}^{2}");
1058 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetXTitle("p_{T} (GeV/c)");
1059 outputContainer->Add(fhSelectedLambda0PtLocMaxSM[inlm][iSM]) ;
1064 for(Int_t ipart = 0; ipart < 6; ipart++)
1066 fhMCSelectedMassPtLocMax[ipart][inlm] = new TH2F
1067 (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
1068 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, NLM=%s, %s",nlm[inlm].Data(),pname[ipart].Data()),
1069 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1070 fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("mass (GeV/c^{2})");
1071 fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("p_{T} (GeV/c)");
1072 outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
1079 fhMassNoOverlap = new TH2F
1080 ("hMassNoOverlap","all pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1081 fhMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1082 fhMassNoOverlap->SetXTitle("E (GeV)");
1083 outputContainer->Add(fhMassNoOverlap) ;
1085 fhSelectedMassNoOverlap = new TH2F
1086 ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1087 fhSelectedMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1088 fhSelectedMassNoOverlap->SetXTitle("E (GeV)");
1089 outputContainer->Add(fhSelectedMassNoOverlap) ;
1091 fhMassPtNoOverlap = new TH2F
1092 ("hMassPtNoOverlap","all pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1093 fhMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1094 fhMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1095 outputContainer->Add(fhMassPtNoOverlap) ;
1097 fhSelectedMassPtNoOverlap = new TH2F
1098 ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1099 fhSelectedMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1100 fhSelectedMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1101 outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1105 if(fAnaType != kSSCalo)
1107 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1108 fhPtDecay->SetYTitle("N");
1109 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
1110 outputContainer->Add(fhPtDecay) ;
1112 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1113 fhEDecay->SetYTitle("N");
1114 fhEDecay->SetXTitle("E (GeV)");
1115 outputContainer->Add(fhEDecay) ;
1120 if( fFillSelectClHisto )
1122 fhPtDispersion = new TH2F
1123 ("hPtDispersion","Selected #pi^{0} (#eta) pairs: p_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1124 fhPtDispersion->SetYTitle("D^{2}");
1125 fhPtDispersion->SetXTitle("p_{T} (GeV/c)");
1126 outputContainer->Add(fhPtDispersion) ;
1128 fhPtLambda0 = new TH2F
1129 ("hPtLambda0","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1130 fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1131 fhPtLambda0->SetXTitle("p_{T} (GeV/c)");
1132 outputContainer->Add(fhPtLambda0) ;
1134 fhPtLambda1 = new TH2F
1135 ("hPtLambda1","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1136 fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1137 fhPtLambda1->SetXTitle("p_{T} (GeV/c)");
1138 outputContainer->Add(fhPtLambda1) ;
1140 fhPtLambda0FracMaxCellCut = new TH2F
1141 ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1142 fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1143 fhPtLambda0FracMaxCellCut->SetXTitle("p_{T} (GeV/c)");
1144 outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1146 fhPtFracMaxCell = new TH2F
1147 ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1148 fhPtFracMaxCell->SetYTitle("Fraction");
1149 fhPtFracMaxCell->SetXTitle("p_{T} (GeV/c)");
1150 outputContainer->Add(fhPtFracMaxCell) ;
1152 if(fCalorimeter=="EMCAL")
1154 fhPtLambda0NoTRD = new TH2F
1155 ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1156 fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1157 fhPtLambda0NoTRD->SetXTitle("p_{T} (GeV/c)");
1158 outputContainer->Add(fhPtLambda0NoTRD) ;
1160 fhPtFracMaxCellNoTRD = new TH2F
1161 ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1162 fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1163 fhPtFracMaxCellNoTRD->SetXTitle("p_{T} (GeV/c)");
1164 outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1166 if(!fFillOnlySimpleSSHisto)
1168 fhPtDispEta = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1169 fhPtDispEta->SetXTitle("p_{T} (GeV/c)");
1170 fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1171 outputContainer->Add(fhPtDispEta);
1173 fhPtDispPhi = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1174 fhPtDispPhi->SetXTitle("p_{T} (GeV/c)");
1175 fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1176 outputContainer->Add(fhPtDispPhi);
1178 fhPtSumEta = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1179 fhPtSumEta->SetXTitle("p_{T} (GeV/c)");
1180 fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1181 outputContainer->Add(fhPtSumEta);
1183 fhPtSumPhi = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs p_{T}",
1184 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1185 fhPtSumPhi->SetXTitle("p_{T} (GeV/c)");
1186 fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1187 outputContainer->Add(fhPtSumPhi);
1189 fhPtSumEtaPhi = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs p_{T}",
1190 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1191 fhPtSumEtaPhi->SetXTitle("p_{T} (GeV/c)");
1192 fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1193 outputContainer->Add(fhPtSumEtaPhi);
1195 fhPtDispEtaPhiDiff = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs p_{T}",
1196 nptbins,ptmin,ptmax,200, -10,10);
1197 fhPtDispEtaPhiDiff->SetXTitle("p_{T} (GeV/c)");
1198 fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1199 outputContainer->Add(fhPtDispEtaPhiDiff);
1201 fhPtSphericity = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs p_{T} (GeV/c)",
1202 nptbins,ptmin,ptmax, 200, -1,1);
1203 fhPtSphericity->SetXTitle("p_{T} (GeV/c)");
1204 fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1205 outputContainer->Add(fhPtSphericity);
1207 for(Int_t i = 0; i < 7; i++)
1209 fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1210 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1211 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1212 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1213 outputContainer->Add(fhDispEtaDispPhi[i]);
1215 fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1216 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1217 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1218 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1219 outputContainer->Add(fhLambda0DispEta[i]);
1221 fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
1222 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1223 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1224 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1225 outputContainer->Add(fhLambda0DispPhi[i]);
1231 fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1232 nptbins,ptmin,ptmax,20,0,20);
1233 fhNLocMaxPt ->SetYTitle("N maxima");
1234 fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
1235 outputContainer->Add(fhNLocMaxPt) ;
1237 for(Int_t iSM = 0; iSM < fNSuperModules; iSM++)
1239 fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
1240 nptbins,ptmin,ptmax,20,0,20);
1241 fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
1242 fhNLocMaxPtSM[iSM] ->SetXTitle("p_{T} (GeV/c)");
1243 outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
1246 if(fAnaType == kSSCalo)
1249 fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1250 nptbins,ptmin,ptmax,20,0,20);
1251 fhNLocMaxPtReject ->SetYTitle("N maxima");
1252 fhNLocMaxPtReject ->SetXTitle("p_{T} (GeV/c)");
1253 outputContainer->Add(fhNLocMaxPtReject) ;
1256 for (Int_t i = 0; i < 3; i++)
1258 fhPtLambda0LocMax[i] = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1259 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1260 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1261 fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1262 fhPtLambda0LocMax[i]->SetXTitle("p_{T} (GeV/c)");
1263 outputContainer->Add(fhPtLambda0LocMax[i]) ;
1267 for(Int_t ipart = 0; ipart < 6; ipart++)
1269 fhMCPtLambda0LocMax[ipart][i] = new TH2F
1270 (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1271 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),pname[ipart].Data()),
1272 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1273 fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1274 fhMCPtLambda0LocMax[ipart][i]->SetXTitle("p_{T} (GeV/c)");
1275 outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1279 fhPtLambda1LocMax[i] = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1280 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1281 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1282 fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1283 fhPtLambda1LocMax[i]->SetXTitle("p_{T} (GeV/c)");
1284 outputContainer->Add(fhPtLambda1LocMax[i]) ;
1286 fhPtDispersionLocMax[i] = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1287 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1288 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1289 fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1290 fhPtDispersionLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1291 outputContainer->Add(fhPtDispersionLocMax[i]) ;
1293 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1295 fhPtDispEtaLocMax[i] = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1296 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1297 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1298 fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1299 fhPtDispEtaLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1300 outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1302 fhPtDispPhiLocMax[i] = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1303 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1304 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1305 fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1306 fhPtDispPhiLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1307 outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1309 fhPtSumEtaPhiLocMax[i] = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1310 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1311 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1312 fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1313 fhPtSumEtaPhiLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1314 outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1316 fhPtDispEtaPhiDiffLocMax[i] = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1317 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1318 nptbins,ptmin,ptmax,200, -10,10);
1319 fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1320 fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1321 outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1323 fhPtSphericityLocMax[i] = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1324 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1325 nptbins,ptmin,ptmax,200, -1,1);
1326 fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1327 fhPtSphericityLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1328 outputContainer->Add(fhPtSphericityLocMax[i]) ;
1333 fhPtNCells = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1334 fhPtNCells->SetXTitle("p_{T} (GeV/c)");
1335 fhPtNCells->SetYTitle("# of cells in cluster");
1336 outputContainer->Add(fhPtNCells);
1338 fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1339 fhPtTime->SetXTitle("p_{T} (GeV/c)");
1340 fhPtTime->SetYTitle("t (ns)");
1341 outputContainer->Add(fhPtTime);
1346 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1347 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
1348 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1349 outputContainer->Add(fhEPairDiffTime);
1351 if(fAnaType == kIMCalo)
1353 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1354 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1355 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1356 "2 Local Maxima paired with more than 2 Local Maxima",
1357 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1359 for (Int_t i = 0; i < 8 ; i++)
1362 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1364 fhMassPairLocMax[i] = new TH2F
1365 (Form("MassPairLocMax%s",combiName[i].Data()),
1366 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1367 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1368 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
1369 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
1370 outputContainer->Add(fhMassPairLocMax[i]) ;
1376 fhTrackMatchedDEta = new TH2F
1377 ("hTrackMatchedDEta",
1378 "d#eta of cluster-track vs cluster p_{T}",
1379 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1380 fhTrackMatchedDEta->SetYTitle("d#eta");
1381 fhTrackMatchedDEta->SetXTitle("p_{T} (GeV/c)");
1383 fhTrackMatchedDPhi = new TH2F
1384 ("hTrackMatchedDPhi",
1385 "d#phi of cluster-track vs cluster p_{T}",
1386 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1387 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1388 fhTrackMatchedDPhi->SetXTitle("p_{T} (GeV/c)");
1390 fhTrackMatchedDEtaDPhi = new TH2F
1391 ("hTrackMatchedDEtaDPhi",
1392 "d#eta vs d#phi of cluster-track",
1393 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1394 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1395 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1397 outputContainer->Add(fhTrackMatchedDEta) ;
1398 outputContainer->Add(fhTrackMatchedDPhi) ;
1399 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1401 fhTrackMatchedDEtaPos = new TH2F
1402 ("hTrackMatchedDEtaPos",
1403 "d#eta of cluster-track vs cluster p_{T}",
1404 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1405 fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1406 fhTrackMatchedDEtaPos->SetXTitle("p_{T} (GeV/c)");
1408 fhTrackMatchedDPhiPos = new TH2F
1409 ("hTrackMatchedDPhiPos",
1410 "d#phi of cluster-track vs cluster p_{T}",
1411 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1412 fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1413 fhTrackMatchedDPhiPos->SetXTitle("p_{T} (GeV/c)");
1415 fhTrackMatchedDEtaDPhiPos = new TH2F
1416 ("hTrackMatchedDEtaDPhiPos",
1417 "d#eta vs d#phi of cluster-track",
1418 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1419 fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1420 fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1422 outputContainer->Add(fhTrackMatchedDEtaPos) ;
1423 outputContainer->Add(fhTrackMatchedDPhiPos) ;
1424 outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1426 fhTrackMatchedDEtaNeg = new TH2F
1427 ("hTrackMatchedDEtaNeg",
1428 "d#eta of cluster-track vs cluster p_{T}",
1429 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1430 fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1431 fhTrackMatchedDEtaNeg->SetXTitle("p_{T} (GeV/c)");
1433 fhTrackMatchedDPhiNeg = new TH2F
1434 ("hTrackMatchedDPhiNeg",
1435 "d#phi of cluster-track vs cluster p_{T}",
1436 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1437 fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1438 fhTrackMatchedDPhiNeg->SetXTitle("p_{T} (GeV/c)");
1440 fhTrackMatchedDEtaDPhiNeg = new TH2F
1441 ("hTrackMatchedDEtaDPhiNeg",
1442 "d#eta vs d#phi of cluster-track",
1443 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1444 fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1445 fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1447 outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1448 outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1449 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1451 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster p_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1452 fhdEdx->SetXTitle("p_{T} (GeV/c)");
1453 fhdEdx->SetYTitle("<dE/dx>");
1454 outputContainer->Add(fhdEdx);
1456 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster p_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1457 fhEOverP->SetXTitle("p_{T} (GeV/c)");
1458 fhEOverP->SetYTitle("E/p");
1459 outputContainer->Add(fhEOverP);
1461 if(fCalorimeter=="EMCAL")
1463 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1464 fhEOverPNoTRD->SetXTitle("E (GeV)");
1465 fhEOverPNoTRD->SetYTitle("E/p");
1466 outputContainer->Add(fhEOverPNoTRD);
1469 if(IsDataMC() && fFillTMHisto)
1471 fhTrackMatchedMCParticlePt = new TH2F
1472 ("hTrackMatchedMCParticlePt",
1473 "Origin of particle vs energy",
1474 nptbins,ptmin,ptmax,8,0,8);
1475 fhTrackMatchedMCParticlePt->SetXTitle("p_{T} (GeV/c)");
1476 //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1478 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1479 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1480 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1481 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1482 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1483 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1484 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1485 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1487 outputContainer->Add(fhTrackMatchedMCParticlePt);
1489 fhTrackMatchedMCParticleDEta = new TH2F
1490 ("hTrackMatchedMCParticleDEta",
1491 "Origin of particle vs #eta residual",
1492 nresetabins,resetamin,resetamax,8,0,8);
1493 fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1494 //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1496 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1497 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1498 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1499 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1500 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1501 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1502 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1503 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1505 outputContainer->Add(fhTrackMatchedMCParticleDEta);
1507 fhTrackMatchedMCParticleDPhi = new TH2F
1508 ("hTrackMatchedMCParticleDPhi",
1509 "Origin of particle vs #phi residual",
1510 nresphibins,resphimin,resphimax,8,0,8);
1511 fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1512 //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1514 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1515 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1516 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1517 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1518 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1519 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1520 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1521 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1523 outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1529 if(fFillWeightHistograms)
1531 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1532 nptbins,ptmin,ptmax, 100,0,1.);
1533 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1534 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1535 outputContainer->Add(fhECellClusterRatio);
1537 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1538 nptbins,ptmin,ptmax, 100,-10,0);
1539 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1540 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1541 outputContainer->Add(fhECellClusterLogRatio);
1543 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1544 nptbins,ptmin,ptmax, 100,0,1.);
1545 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1546 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1547 outputContainer->Add(fhEMaxCellClusterRatio);
1549 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1550 nptbins,ptmin,ptmax, 100,-10,0);
1551 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1552 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1553 outputContainer->Add(fhEMaxCellClusterLogRatio);
1555 for(Int_t iw = 0; iw < 14; iw++)
1557 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
1558 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1559 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1560 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1561 outputContainer->Add(fhLambda0ForW0[iw]);
1563 // fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
1564 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1565 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1566 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1567 // outputContainer->Add(fhLambda1ForW0[iw]);
1574 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1576 fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
1577 nptbins,ptmin,ptmax,200,0,2);
1578 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1579 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1580 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1582 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1583 nptbins,ptmin,ptmax,200,0,2);
1584 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1585 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1586 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1588 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1589 fhMCPi0DecayPt->SetYTitle("N");
1590 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1591 outputContainer->Add(fhMCPi0DecayPt) ;
1593 fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #pi^{0}",
1594 nptbins,ptmin,ptmax,100,0,1);
1595 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1596 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1597 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1599 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1600 fhMCEtaDecayPt->SetYTitle("N");
1601 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1602 outputContainer->Add(fhMCEtaDecayPt) ;
1604 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
1605 nptbins,ptmin,ptmax,100,0,1);
1606 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1607 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1608 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1610 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1611 fhMCOtherDecayPt->SetYTitle("N");
1612 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1613 outputContainer->Add(fhMCOtherDecayPt) ;
1617 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1618 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1621 fhAnglePairMCPi0 = new TH2F
1623 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1624 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1625 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1626 outputContainer->Add(fhAnglePairMCPi0) ;
1628 if (fAnaType!= kSSCalo)
1630 fhAnglePairMCEta = new TH2F
1632 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1633 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1634 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1635 outputContainer->Add(fhAnglePairMCEta) ;
1637 fhMassPairMCPi0 = new TH2F
1639 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1640 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1641 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1642 outputContainer->Add(fhMassPairMCPi0) ;
1644 fhMassPairMCEta = new TH2F
1646 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1647 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1648 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1649 outputContainer->Add(fhMassPairMCEta) ;
1652 for(Int_t i = 0; i < 6; i++)
1656 (Form("hE_MC%s",pname[i].Data()),
1657 Form("Identified as #pi^{0} (#eta), cluster from %s",
1659 nptbins,ptmin,ptmax);
1660 fhMCE[i]->SetYTitle("N");
1661 fhMCE[i]->SetXTitle("E (GeV)");
1662 outputContainer->Add(fhMCE[i]) ;
1664 fhMCPt[i] = new TH1F
1665 (Form("hPt_MC%s",pname[i].Data()),
1666 Form("Identified as #pi^{0} (#eta), cluster from %s",
1668 nptbins,ptmin,ptmax);
1669 fhMCPt[i]->SetYTitle("N");
1670 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1671 outputContainer->Add(fhMCPt[i]) ;
1673 fhMCPtCentrality[i] = new TH2F
1674 (Form("hPtCentrality_MC%s",pname[i].Data()),
1675 Form("Identified as #pi^{0} (#eta), cluster from %s",
1677 nptbins,ptmin,ptmax, 100,0,100);
1678 fhMCPtCentrality[i]->SetYTitle("centrality");
1679 fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1680 outputContainer->Add(fhMCPtCentrality[i]) ;
1682 if(fAnaType == kSSCalo)
1684 fhMCNLocMaxPt[i] = new TH2F
1685 (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1686 Form("cluster from %s, pT of cluster vs NLM, accepted",ptype[i].Data()),
1687 nptbins,ptmin,ptmax,20,0,20);
1688 fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
1689 fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
1690 outputContainer->Add(fhMCNLocMaxPt[i]) ;
1692 fhMCNLocMaxPtReject[i] = new TH2F
1693 (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1694 Form("cluster from %s, pT of cluster vs NLM, rejected",ptype[i].Data()),
1695 nptbins,ptmin,ptmax,20,0,20);
1696 fhMCNLocMaxPtReject[i] ->SetYTitle("N maxima");
1697 fhMCNLocMaxPtReject[i] ->SetXTitle("p_{T} (GeV/c)");
1698 outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1700 fhMCEReject[i] = new TH1F
1701 (Form("hEReject_MC%s",pname[i].Data()),
1702 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1704 nptbins,ptmin,ptmax);
1705 fhMCEReject[i]->SetYTitle("N");
1706 fhMCEReject[i]->SetXTitle("E (GeV)");
1707 outputContainer->Add(fhMCEReject[i]) ;
1709 fhMCPtReject[i] = new TH1F
1710 (Form("hPtReject_MC%s",pname[i].Data()),
1711 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1713 nptbins,ptmin,ptmax);
1714 fhMCPtReject[i]->SetYTitle("N");
1715 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1716 outputContainer->Add(fhMCPtReject[i]) ;
1719 fhMCPtPhi[i] = new TH2F
1720 (Form("hPtPhi_MC%s",pname[i].Data()),
1721 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1722 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1723 fhMCPtPhi[i]->SetYTitle("#phi");
1724 fhMCPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
1725 outputContainer->Add(fhMCPtPhi[i]) ;
1727 fhMCPtEta[i] = new TH2F
1728 (Form("hPtEta_MC%s",pname[i].Data()),
1729 Form("Identified as #pi^{0} (#eta), cluster from %s",
1730 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1731 fhMCPtEta[i]->SetYTitle("#eta");
1732 fhMCPtEta[i]->SetXTitle("p_{T} (GeV/c)");
1733 outputContainer->Add(fhMCPtEta[i]) ;
1735 fhMCMassPt[i] = new TH2F
1736 (Form("hMassPt_MC%s",pname[i].Data()),
1737 Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1738 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1739 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1740 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1741 outputContainer->Add(fhMCMassPt[i]) ;
1743 fhMCSelectedMassPt[i] = new TH2F
1744 (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1745 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1746 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1747 fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1748 fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1749 outputContainer->Add(fhMCSelectedMassPt[i]) ;
1751 if(fAnaType == kSSCalo)
1753 fhMCMassPtNoOverlap[i] = new TH2F
1754 (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1755 Form("all pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1756 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1757 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1758 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1759 outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1761 fhMCSelectedMassPtNoOverlap[i] = new TH2F
1762 (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1763 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1764 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1765 fhMCSelectedMassPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
1766 fhMCSelectedMassPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
1767 outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1770 if( fFillSelectClHisto )
1772 fhMCPtLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1773 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1774 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1775 fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1776 fhMCPtLambda0[i]->SetXTitle("p_{T} (GeV/c)");
1777 outputContainer->Add(fhMCPtLambda0[i]) ;
1779 fhMCPtLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1780 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1781 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1782 fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1783 fhMCPtLambda1[i]->SetXTitle("p_{T} (GeV/c)");
1784 outputContainer->Add(fhMCPtLambda1[i]) ;
1786 fhMCPtDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1787 Form("Selected pair, cluster from %s : p_{T} vs dispersion^{2}",ptype[i].Data()),
1788 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1789 fhMCPtDispersion[i]->SetYTitle("D^{2}");
1790 fhMCPtDispersion[i]->SetXTitle("p_{T} (GeV/c)");
1791 outputContainer->Add(fhMCPtDispersion[i]) ;
1793 if(fCalorimeter=="EMCAL")
1795 fhMCPtLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1796 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1797 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1798 fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1799 fhMCPtLambda0NoTRD[i]->SetXTitle("p_{T} (GeV/c)");
1800 outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1802 if(!fFillOnlySimpleSSHisto)
1804 fhMCPtDispEta[i] = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1805 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs p_{T}",ptype[i].Data()),
1806 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1807 fhMCPtDispEta[i]->SetXTitle("p_{T} (GeV/c)");
1808 fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1809 outputContainer->Add(fhMCPtDispEta[i]);
1811 fhMCPtDispPhi[i] = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1812 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs p_{T}",ptype[i].Data()),
1813 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1814 fhMCPtDispPhi[i]->SetXTitle("p_{T} (GeV/c)");
1815 fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1816 outputContainer->Add(fhMCPtDispPhi[i]);
1818 fhMCPtSumEtaPhi[i] = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1819 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs p_{T}",ptype[i].Data()),
1820 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1821 fhMCPtSumEtaPhi[i]->SetXTitle("p_{T} (GeV/c)");
1822 fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1823 outputContainer->Add(fhMCPtSumEtaPhi[i]);
1825 fhMCPtDispEtaPhiDiff[i] = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1826 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs p_{T}",ptype[i].Data()),
1827 nptbins,ptmin,ptmax,200,-10,10);
1828 fhMCPtDispEtaPhiDiff[i]->SetXTitle("p_{T} (GeV/c)");
1829 fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1830 outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1832 fhMCPtSphericity[i] = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1833 Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
1834 nptbins,ptmin,ptmax, 200,-1,1);
1835 fhMCPtSphericity[i]->SetXTitle("p_{T} (GeV/c)");
1836 fhMCPtSphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1837 outputContainer->Add(fhMCPtSphericity[i]);
1839 for(Int_t ie = 0; ie < 7; ie++)
1841 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1842 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1843 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1844 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1845 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1846 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1848 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1849 Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1850 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1851 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1852 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1853 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1855 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1856 Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1857 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1858 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1859 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1860 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1866 fhMCPtLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1867 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1868 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1869 fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1870 fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1871 outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
1873 fhMCPtFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1874 Form("Selected pair, cluster from %s : p_{T} vs Max cell fraction of energy",ptype[i].Data()),
1875 nptbins,ptmin,ptmax,100,0,1);
1876 fhMCPtFracMaxCell[i]->SetYTitle("Fraction");
1877 fhMCPtFracMaxCell[i]->SetXTitle("E (GeV)");
1878 outputContainer->Add(fhMCPtFracMaxCell[i]) ;
1881 } // shower shape histo
1886 if(fAnaType==kSSCalo)
1888 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1889 nptbins,ptmin,ptmax, 200, -1,1);
1890 fhAsymmetry->SetXTitle("E (GeV)");
1891 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1892 outputContainer->Add(fhAsymmetry);
1894 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1895 nptbins,ptmin,ptmax, 200, -1,1);
1896 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1897 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1898 outputContainer->Add(fhSelectedAsymmetry);
1901 ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1902 fhSplitE->SetYTitle("counts");
1903 fhSplitE->SetXTitle("E (GeV)");
1904 outputContainer->Add(fhSplitE) ;
1906 fhSplitPt = new TH1F
1907 ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
1908 fhSplitPt->SetYTitle("counts");
1909 fhSplitPt->SetXTitle("p_{T} (GeV/c)");
1910 outputContainer->Add(fhSplitPt) ;
1913 fhSplitPtPhi = new TH2F
1914 ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1915 fhSplitPtPhi->SetYTitle("#phi (rad)");
1916 fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
1917 outputContainer->Add(fhSplitPtPhi) ;
1919 fhSplitPtEta = new TH2F
1920 ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1921 fhSplitPtEta->SetYTitle("#eta");
1922 fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
1923 outputContainer->Add(fhSplitPtEta) ;
1926 fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1927 nptbins,ptmin,ptmax,20,0,20);
1928 fhNLocMaxSplitPt ->SetYTitle("N maxima");
1929 fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
1930 outputContainer->Add(fhNLocMaxSplitPt) ;
1933 fhMassSplitPt = new TH2F
1934 ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",
1935 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1936 fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1937 fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1938 outputContainer->Add(fhMassSplitPt) ;
1940 fhSelectedMassSplitPt = new TH2F
1941 ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",
1942 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1943 fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1944 fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1945 outputContainer->Add(fhSelectedMassSplitPt) ;
1949 fhMassSplitPtNoOverlap = new TH2F
1950 ("hMassSplitPtNoOverlap","all pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1951 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1952 fhMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1953 fhMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1954 outputContainer->Add(fhMassSplitPtNoOverlap) ;
1956 fhSelectedMassSplitPtNoOverlap = new TH2F
1957 ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1958 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1959 fhSelectedMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1960 fhSelectedMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1961 outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
1964 fhMCPi0PtRecoPtPrim = new TH2F
1965 ("hMCPi0PtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1966 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1967 fhMCPi0PtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1968 fhMCPi0PtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1969 outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
1971 fhMCPi0PtRecoPtPrimNoOverlap = new TH2F
1972 ("hMCPi0PtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1973 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1974 fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1975 fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1976 outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
1978 fhMCPi0SelectedPtRecoPtPrim = new TH2F
1979 ("hMCPi0SelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1980 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1981 fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1982 fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1983 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
1985 fhMCPi0SelectedPtRecoPtPrimNoOverlap = new TH2F
1986 ("hMCPi0SelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1987 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1988 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1989 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1990 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
1993 fhMCPi0SplitPtRecoPtPrim = new TH2F
1994 ("hMCPi0SplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1995 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1996 fhMCPi0SplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1997 fhMCPi0SplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1998 outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
2000 fhMCPi0SplitPtRecoPtPrimNoOverlap = new TH2F
2001 ("hMCPi0SplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2002 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2003 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2004 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2005 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
2007 fhMCPi0SelectedSplitPtRecoPtPrim = new TH2F
2008 ("hMCPi0SelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
2009 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2010 fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2011 fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2012 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
2014 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2015 ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2016 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2017 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2018 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2019 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
2021 fhMCEtaPtRecoPtPrim = new TH2F
2022 ("hMCEtaPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
2023 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2024 fhMCEtaPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2025 fhMCEtaPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2026 outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
2028 fhMCEtaPtRecoPtPrimNoOverlap = new TH2F
2029 ("hMCEtaPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
2030 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2031 fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2032 fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2033 outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2035 fhMCEtaSelectedPtRecoPtPrim = new TH2F
2036 ("hMCEtaSelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
2037 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2038 fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2039 fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2040 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2042 fhMCEtaSelectedPtRecoPtPrimNoOverlap = new TH2F
2043 ("hMCEtaSelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
2044 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2045 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2046 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2047 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2050 fhMCEtaSplitPtRecoPtPrim = new TH2F
2051 ("hMCEtaSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
2052 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2053 fhMCEtaSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2054 fhMCEtaSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2055 outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2057 fhMCEtaSplitPtRecoPtPrimNoOverlap = new TH2F
2058 ("hMCEtaSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2059 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2060 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2061 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2062 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2064 fhMCEtaSelectedSplitPtRecoPtPrim = new TH2F
2065 ("hMCEtaSelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
2066 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2067 fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2068 fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2069 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2071 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2072 ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2073 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2074 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2075 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2076 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2079 for(Int_t inlm = 0; inlm < 3; inlm++)
2081 fhMCPi0PtRecoPtPrimLocMax[inlm] = new TH2F
2082 (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} vs p_{T,gen}, %s",nlm[inlm].Data()),
2083 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2084 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2085 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2086 outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
2088 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2089 (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} vs p_{T,gen}, %s",nlm[inlm].Data()),
2090 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2091 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2092 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2093 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
2095 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] = new TH2F
2096 (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} (split sum) vs p_{T,gen}, %s",nlm[inlm].Data()),
2097 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2098 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2099 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2100 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
2102 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2103 (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} (split sum) vs p_{T,gen}, %s",nlm[inlm].Data()),
2104 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2105 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2106 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2107 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2109 fhMCEtaPtRecoPtPrimLocMax[inlm] = new TH2F
2110 (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} vs p_{T,gen}, %s",nlm[inlm].Data()),
2111 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2112 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2113 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2114 outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
2116 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2117 (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} vs p_{T,gen}, %s",nlm[inlm].Data()),
2118 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2119 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2120 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2121 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
2123 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2124 (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} (split sum) vs p_{T,gen}, %s",nlm[inlm].Data()),
2125 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2126 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2127 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2128 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
2130 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2131 (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("p_{T,reco} (split sum) vs p_{T,gen}, %s",nlm[inlm].Data()),
2132 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2133 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("p_{T,gen} (GeV/c)");
2134 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("p_{T,reco} (GeV/c)");
2135 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2139 for(Int_t i = 0; i< 6; i++)
2141 fhMCPtAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2142 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
2143 nptbins,ptmin,ptmax, 200,-1,1);
2144 fhMCPtAsymmetry[i]->SetXTitle("E (GeV)");
2145 fhMCPtAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2146 outputContainer->Add(fhMCPtAsymmetry[i]);
2148 fhMCSplitE[i] = new TH1F
2149 (Form("hSplitE_MC%s",pname[i].Data()),
2150 Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2151 nptbins,ptmin,ptmax);
2152 fhMCSplitE[i]->SetYTitle("counts");
2153 fhMCSplitE[i]->SetXTitle("E (GeV)");
2154 outputContainer->Add(fhMCSplitE[i]) ;
2156 fhMCSplitPt[i] = new TH1F
2157 (Form("hSplitPt_MC%s",pname[i].Data()),
2158 Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
2159 nptbins,ptmin,ptmax);
2160 fhMCSplitPt[i]->SetYTitle("counts");
2161 fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2162 outputContainer->Add(fhMCSplitPt[i]) ;
2165 fhMCSplitPtPhi[i] = new TH2F
2166 (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2167 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2168 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2169 fhMCSplitPtPhi[i]->SetYTitle("#phi");
2170 fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
2171 outputContainer->Add(fhMCSplitPtPhi[i]) ;
2173 fhMCSplitPtEta[i] = new TH2F
2174 (Form("hSplitPtEta_MC%s",pname[i].Data()),
2175 Form("Identified as #pi^{0} (#eta), cluster from %s",
2176 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2177 fhMCSplitPtEta[i]->SetYTitle("#eta");
2178 fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
2179 outputContainer->Add(fhMCSplitPtEta[i]) ;
2182 fhMCNLocMaxSplitPt[i] = new TH2F
2183 (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2184 Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
2185 nptbins,ptmin,ptmax,20,0,20);
2186 fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
2187 fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
2188 outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2190 fhMCMassSplitPt[i] = new TH2F
2191 (Form("hMassSplitPt_MC%s",pname[i].Data()),
2192 Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2193 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2194 fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2195 fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2196 outputContainer->Add(fhMCMassSplitPt[i]) ;
2198 fhMCSelectedMassSplitPt[i] = new TH2F
2199 (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2200 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2201 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2202 fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2203 fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2204 outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2206 fhMCMassSplitPtNoOverlap[i] = new TH2F
2207 (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2208 Form("all pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2209 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2210 fhMCMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2211 fhMCMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2212 outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2214 fhMCSelectedMassSplitPtNoOverlap[i] = new TH2F
2215 (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2216 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2217 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2218 fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2219 fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2220 outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2225 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2229 for(Int_t i = 0; i< 3; i++)
2231 fhPtAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2232 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
2233 nptbins,ptmin,ptmax,200, -1,1);
2234 fhPtAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2235 fhPtAsymmetryLocMax[i]->SetXTitle("p_{T} (GeV/c)");
2236 outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2239 for(Int_t ie = 0; ie< 7; ie++)
2242 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2243 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2244 ssbins,ssmin,ssmax , 200,-1,1);
2245 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2246 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2247 outputContainer->Add(fhAsymmetryLambda0[ie]);
2249 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2250 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2251 ssbins,ssmin,ssmax , 200,-1,1);
2252 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2253 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2254 outputContainer->Add(fhAsymmetryDispEta[ie]);
2256 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2257 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2258 ssbins,ssmin,ssmax , 200,-1,1);
2259 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2260 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2261 outputContainer->Add(fhAsymmetryDispPhi[ie]);
2267 for(Int_t i = 0; i< 6; i++)
2269 for(Int_t ie = 0; ie < 7; ie++)
2271 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2272 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2273 ssbins,ssmin,ssmax , 200,-1,1);
2274 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2275 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2276 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2278 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2279 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2280 ssbins,ssmin,ssmax , 200,-1,1);
2281 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2282 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2283 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2285 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2286 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2287 ssbins,ssmin,ssmax , 200,-1,1);
2288 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2289 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2290 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2296 if(fFillPileUpHistograms)
2299 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2301 for(Int_t i = 0 ; i < 7 ; i++)
2303 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2304 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2305 fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2306 outputContainer->Add(fhPtPileUp[i]);
2308 fhPtCellTimePileUp[i] = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2309 Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2310 nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2311 fhPtCellTimePileUp[i]->SetXTitle("p_{T} (GeV/c)");
2312 fhPtCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
2313 outputContainer->Add(fhPtCellTimePileUp[i]);
2315 fhPtTimeDiffPileUp[i] = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2316 Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2317 nptbins,ptmin,ptmax,400,-200,200);
2318 fhPtTimeDiffPileUp[i]->SetXTitle("p_{T} (GeV/c");
2319 fhPtTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2320 outputContainer->Add(fhPtTimeDiffPileUp[i]);
2324 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2325 fhTimePtNoCut->SetXTitle("p_{T} (GeV/c)");
2326 fhTimePtNoCut->SetYTitle("time (ns)");
2327 outputContainer->Add(fhTimePtNoCut);
2329 fhTimePtSPD = new TH2F ("hTimePt_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2330 fhTimePtSPD->SetXTitle("p_{T} (GeV/c)");
2331 fhTimePtSPD->SetYTitle("time (ns)");
2332 outputContainer->Add(fhTimePtSPD);
2334 fhTimePtSPDMulti = new TH2F ("hTimePt_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2335 fhTimePtSPDMulti->SetXTitle("p_{T} (GeV/c)");
2336 fhTimePtSPDMulti->SetYTitle("time (ns)");
2337 outputContainer->Add(fhTimePtSPDMulti);
2339 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2340 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2341 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
2342 outputContainer->Add(fhTimeNPileUpVertSPD);
2344 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2345 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2346 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
2347 outputContainer->Add(fhTimeNPileUpVertTrack);
2349 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2350 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2351 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
2352 outputContainer->Add(fhTimeNPileUpVertContributors);
2354 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
2355 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
2356 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
2357 outputContainer->Add(fhTimePileUpMainVertexZDistance);
2359 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
2360 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
2361 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
2362 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2364 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
2365 nptbins,ptmin,ptmax,20,0,20);
2366 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2367 fhPtNPileUpSPDVtx->SetXTitle("p_{T} (GeV/c)");
2368 outputContainer->Add(fhPtNPileUpSPDVtx);
2370 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
2371 nptbins,ptmin,ptmax, 20,0,20 );
2372 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2373 fhPtNPileUpTrkVtx->SetXTitle("p_{T} (GeV/c)");
2374 outputContainer->Add(fhPtNPileUpTrkVtx);
2376 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2377 nptbins,ptmin,ptmax,20,0,20);
2378 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2379 fhPtNPileUpSPDVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2380 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2382 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2383 nptbins,ptmin,ptmax, 20,0,20 );
2384 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2385 fhPtNPileUpTrkVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2386 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2388 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2389 nptbins,ptmin,ptmax,20,0,20);
2390 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2391 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2392 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2394 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2395 nptbins,ptmin,ptmax, 20,0,20 );
2396 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2397 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2398 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2402 //Keep neutral meson selection histograms if requiered
2403 //Setting done in AliNeutralMesonSelection
2405 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2407 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2409 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2410 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2415 return outputContainer ;
2419 //_____________________________________________
2420 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2423 // Assign mc index depending on MC bit set
2425 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
2429 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2433 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2434 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2436 return kmcConversion ;
2437 }//conversion photon
2438 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2441 }//photon no conversion
2442 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2444 return kmcElectron ;
2453 //__________________________________________________________________
2454 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2455 AliAODPWG4Particle * photon2,
2456 Int_t & label, Int_t & tag)
2458 // Check the labels of pare in case mother was same pi0 or eta
2459 // Set the new AOD accordingly
2461 Int_t label1 = photon1->GetLabel();
2462 Int_t label2 = photon2->GetLabel();
2464 if(label1 < 0 || label2 < 0 ) return ;
2466 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2467 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2468 Int_t tag1 = photon1->GetTag();
2469 Int_t tag2 = photon2->GetTag();
2471 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2472 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2473 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
2474 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2475 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
2479 //Check if pi0/eta mother is the same
2480 if(GetReader()->ReadStack())
2484 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2485 label1 = mother1->GetFirstMother();
2486 //mother1 = GetMCStack()->Particle(label1);//pi0
2490 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2491 label2 = mother2->GetFirstMother();
2492 //mother2 = GetMCStack()->Particle(label2);//pi0
2495 else if(GetReader()->ReadAODMCParticles())
2496 {//&& (input > -1)){
2499 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2500 label1 = mother1->GetMother();
2501 //mother1 = GetMCStack()->Particle(label1);//pi0
2505 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2506 label2 = mother2->GetMother();
2507 //mother2 = GetMCStack()->Particle(label2);//pi0
2511 //printf("mother1 %d, mother2 %d\n",label1,label2);
2512 if( label1 == label2 && label1>=0 )
2517 TLorentzVector mom1 = *(photon1->Momentum());
2518 TLorentzVector mom2 = *(photon2->Momentum());
2520 Double_t angle = mom2.Angle(mom1.Vect());
2521 Double_t mass = (mom1+mom2).M();
2522 Double_t epair = (mom1+mom2).E();
2524 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2526 fhMassPairMCPi0 ->Fill(epair,mass);
2527 fhAnglePairMCPi0->Fill(epair,angle);
2528 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2532 fhMassPairMCEta ->Fill(epair,mass);
2533 fhAnglePairMCEta->Fill(epair,angle);
2534 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2538 } // both from eta or pi0 decay
2542 //____________________________________________________________________________
2543 void AliAnaPi0EbE::Init()
2547 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2548 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2551 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2552 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2558 //____________________________________________________________________________
2559 void AliAnaPi0EbE::InitParameters()
2561 //Initialize the parameters of the analysis.
2562 AddToHistogramsName("AnaPi0EbE_");
2564 fInputAODGammaConvName = "PhotonsCTS" ;
2565 fAnaType = kIMCalo ;
2566 fCalorimeter = "EMCAL" ;
2571 fNLMECutMin[0] = 10.;
2572 fNLMECutMin[1] = 6. ;
2573 fNLMECutMin[2] = 6. ;
2575 fNSuperModules = 10;
2579 //__________________________________________________________________
2580 void AliAnaPi0EbE::MakeAnalysisFillAOD()
2582 //Do analysis and fill aods
2587 MakeInvMassInCalorimeter();
2591 MakeShowerShapeIdentification();
2595 MakeInvMassInCalorimeterAndCTS();
2601 //____________________________________________
2602 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
2604 //Do analysis and fill aods
2605 //Search for the photon decay in calorimeters
2606 //Read photon list from AOD, produced in class AliAnaPhoton
2607 //Check if 2 photons have the mass of the pi0.
2609 TLorentzVector mom1;
2610 TLorentzVector mom2;
2611 TLorentzVector mom ;
2616 if(!GetInputAODBranch()){
2617 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2621 //Get shower shape information of clusters
2622 TObjArray *clusters = 0;
2623 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2624 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2626 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2627 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2629 //Vertex cut in case of mixed events
2630 Int_t evtIndex1 = 0 ;
2632 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2633 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
2634 mom1 = *(photon1->Momentum());
2636 //Get original cluster, to recover some information
2638 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2641 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2645 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2647 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2649 Int_t evtIndex2 = 0 ;
2651 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2653 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2656 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
2658 mom2 = *(photon2->Momentum());
2660 //Get original cluster, to recover some information
2662 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2666 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2670 Float_t e1 = photon1->E();
2671 Float_t e2 = photon2->E();
2673 //Select clusters with good time window difference
2674 Float_t tof1 = cluster1->GetTOF()*1e9;;
2675 Float_t tof2 = cluster2->GetTOF()*1e9;;
2676 Double_t t12diff = tof1-tof2;
2677 fhEPairDiffTime->Fill(e1+e2, t12diff);
2678 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2680 //Play with the MC stack if available
2681 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2683 // Check the invariant mass for different selection on the local maxima
2684 // Name of AOD method TO BE FIXED
2685 Int_t nMaxima1 = photon1->GetFiducialArea();
2686 Int_t nMaxima2 = photon2->GetFiducialArea();
2688 Double_t mass = (mom1+mom2).M();
2689 Double_t epair = (mom1+mom2).E();
2691 if(nMaxima1==nMaxima2)
2693 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2694 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2695 else fhMassPairLocMax[2]->Fill(epair,mass);
2697 else if(nMaxima1==1 || nMaxima2==1)
2699 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2700 else fhMassPairLocMax[4]->Fill(epair,mass);
2703 fhMassPairLocMax[5]->Fill(epair,mass);
2705 // combinations with SS axis cut and NLM cut
2706 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2707 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2708 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2709 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2711 //Skip events with too few or too many NLM
2712 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2714 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2717 fhMass->Fill(epair,(mom1+mom2).M());
2719 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2720 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2723 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2725 //Fill some histograms about shower shape
2726 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2728 FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2729 FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2732 // Tag both photons as decay
2733 photon1->SetTagged(kTRUE);
2734 photon2->SetTagged(kTRUE);
2736 fhPtDecay->Fill(photon1->Pt());
2737 fhEDecay ->Fill(photon1->E() );
2739 fhPtDecay->Fill(photon2->Pt());
2740 fhEDecay ->Fill(photon2->E() );
2742 //Create AOD for analysis
2745 //Mass of selected pairs
2746 fhSelectedMass->Fill(epair,mom.M());
2748 // Fill histograms to undertand pile-up before other cuts applied
2749 // Remember to relax time cuts in the reader
2750 FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2752 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2754 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2755 pi0.SetDetector(photon1->GetDetector());
2758 pi0.SetLabel(label);
2761 //Set the indeces of the original caloclusters
2762 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2763 //pi0.SetInputFileIndex(input);
2765 AddAODParticle(pi0);
2773 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2777 //__________________________________________________
2778 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2780 //Do analysis and fill aods
2781 //Search for the photon decay in calorimeters
2782 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2783 //Check if 2 photons have the mass of the pi0.
2785 TLorentzVector mom1;
2786 TLorentzVector mom2;
2787 TLorentzVector mom ;
2792 // Check calorimeter input
2793 if(!GetInputAODBranch()){
2794 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2798 // Get the array with conversion photons
2799 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2800 if(!inputAODGammaConv) {
2802 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2804 if(!inputAODGammaConv) {
2805 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2811 //Get shower shape information of clusters
2812 TObjArray *clusters = 0;
2813 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2814 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2816 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
2817 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2818 if(nCTS<=0 || nCalo <=0)
2820 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2825 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2827 // Do the loop, first calo, second CTS
2828 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2829 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2830 mom1 = *(photon1->Momentum());
2832 //Get original cluster, to recover some information
2834 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2836 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2837 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2839 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2840 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2842 mom2 = *(photon2->Momentum());
2844 Double_t mass = (mom1+mom2).M();
2845 Double_t epair = (mom1+mom2).E();
2847 Int_t nMaxima = photon1->GetFiducialArea();
2848 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2849 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2850 else fhMassPairLocMax[2]->Fill(epair,mass);
2852 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2853 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2855 //Play with the MC stack if available
2858 Int_t label2 = photon2->GetLabel();
2859 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2861 HasPairSameMCMother(photon1, photon2, label, tag) ;
2864 //Mass of selected pairs
2865 fhMass->Fill(epair,(mom1+mom2).M());
2867 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2868 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2870 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2872 //Fill some histograms about shower shape
2873 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2875 FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
2878 // Tag both photons as decay
2879 photon1->SetTagged(kTRUE);
2880 photon2->SetTagged(kTRUE);
2882 fhPtDecay->Fill(photon1->Pt());
2883 fhEDecay ->Fill(photon1->E() );
2885 //Create AOD for analysis
2889 //Mass of selected pairs
2890 fhSelectedMass->Fill(epair,mom.M());
2892 // Fill histograms to undertand pile-up before other cuts applied
2893 // Remember to relax time cuts in the reader
2894 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2896 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2898 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2899 pi0.SetDetector(photon1->GetDetector());
2902 pi0.SetLabel(label);
2905 //Set the indeces of the original tracks or caloclusters
2906 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2907 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2908 //pi0.SetInputFileIndex(input);
2910 AddAODParticle(pi0);
2917 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2922 //_________________________________________________
2923 void AliAnaPi0EbE::MakeShowerShapeIdentification()
2925 //Search for pi0 in fCalorimeter with shower shape analysis
2927 TObjArray * pl = 0x0;
2928 AliVCaloCells * cells = 0x0;
2929 //Select the Calorimeter of the photon
2930 if (fCalorimeter == "PHOS" )
2932 pl = GetPHOSClusters();
2933 cells = GetPHOSCells();
2935 else if (fCalorimeter == "EMCAL")
2937 pl = GetEMCALClusters();
2938 cells = GetEMCALCells();
2943 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2947 TLorentzVector mom ;
2948 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2950 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2952 Int_t evtIndex = 0 ;
2953 if (GetMixedEvent())
2955 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2958 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2960 //Get Momentum vector,
2961 Double_t vertex[]={0,0,0};
2962 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2964 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2965 }//Assume that come from vertex in straight line
2968 calo->GetMomentum(mom,vertex) ;
2971 //If too small or big pt, skip it
2972 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2974 //Check acceptance selection
2975 if(IsFiducialCutOn())
2977 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2978 if(! in ) continue ;
2982 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
2984 //Play with the MC stack if available
2985 //Check origin of the candidates
2989 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
2990 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2991 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2994 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
2996 //Check Distance to Bad channel, set bit.
2997 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2998 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2999 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3000 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3004 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
3006 //If too low number of cells, skip it
3007 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3009 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3014 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
3015 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
3017 //.......................................
3018 // TOF cut, BE CAREFUL WITH THIS CUT
3019 Double_t tof = calo->GetTOF()*1e9;
3020 if(tof < fTimeCutMin || tof > fTimeCutMax)
3022 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3027 //PID selection or bit setting
3029 Double_t mass = 0, angle = 0;
3030 Int_t absId1 =-1, absId2 =-1;
3031 Float_t distbad1 =-1, distbad2 =-1;
3032 Bool_t fidcut1 = 0, fidcut2 = 0;
3033 TLorentzVector l1, l2;
3035 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3036 GetVertex(evtIndex),nMaxima,
3037 mass,angle,l1,l2,absId1,absId2,
3038 distbad1,distbad2,fidcut1,fidcut2) ;
3041 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
3044 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3045 if( (fCheckSplitDistToBad) &&
3046 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3049 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
3050 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
3052 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3056 //Skip events with too few or too many NLM
3057 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3059 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3064 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
3066 //Skip matched clusters with tracks
3067 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3069 FillRejectedClusterHistograms(mom,tag,nMaxima);
3073 Float_t e1 = l1.Energy();
3074 Float_t e2 = l2.Energy();
3075 TLorentzVector l12 = l1+l2;
3076 Float_t ptSplit = l12.Pt();
3077 Float_t eSplit = e1+e2;
3080 Int_t noverlaps = 0;
3084 mcIndex = GetMCIndex(tag);
3087 Int_t mcLabel = calo->GetLabel();
3089 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3091 Int_t mesonLabel = -1;
3093 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3095 if(mcIndex == kmcPi0)
3097 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3098 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3102 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3103 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3107 const UInt_t nlabels = calo->GetNLabels();
3108 Int_t overpdg[nlabels];
3109 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3112 //mass of all clusters
3113 fhMass ->Fill(mom.E() ,mass);
3114 fhMassPt ->Fill(mom.Pt(),mass);
3115 fhMassSplitPt->Fill(ptSplit ,mass);
3117 Int_t indexMax = -1;
3118 if (nMaxima==1) indexMax = 0 ;
3119 else if(nMaxima==2) indexMax = 1 ;
3121 fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3125 fhMCMassPt[mcIndex] ->Fill(mom.Pt(),mass);
3126 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3129 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3130 fhMCPi0SplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3131 fhMCPi0PtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3132 fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3135 else if(mcIndex==kmcEta)
3137 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3138 fhMCEtaSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3139 fhMCEtaPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3140 fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3147 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3148 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3150 else if(mcIndex==kmcEta)
3152 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3153 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3156 fhMassNoOverlap ->Fill(mom.E() ,mass);
3157 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3158 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3160 fhMCMassPtNoOverlap[mcIndex] ->Fill(mom.Pt(),mass);
3161 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3165 // Asymmetry of all clusters
3168 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3169 fhAsymmetry->Fill(mom.E(),asy);
3173 fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3176 // If cluster does not pass pid, not pi0/eta, skip it.
3177 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3179 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3180 FillRejectedClusterHistograms(mom,tag,nMaxima);
3184 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3186 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3187 FillRejectedClusterHistograms(mom,tag,nMaxima);
3192 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3193 mom.Pt(), idPartType);
3195 //Mass and asymmetry of selected pairs
3196 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3197 fhSelectedMass ->Fill(mom.E() ,mass);
3198 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3199 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3200 fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3202 Int_t nSM = GetModuleNumber(calo);
3203 if(nSM < fNSuperModules && nSM >=0)
3205 fhSelectedMassPtLocMaxSM [indexMax][nSM]->Fill(mom.Pt(),mass);
3206 fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),calo->GetM02());
3213 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3214 fhMCPi0SelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3215 fhMCPi0SelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3216 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3218 else if(mcIndex==kmcEta)
3220 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3221 fhMCEtaSelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3222 fhMCEtaSelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3223 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3228 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3229 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3230 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3234 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3235 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3237 else if(mcIndex==kmcEta)
3239 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3240 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3245 fhSplitE ->Fill( eSplit);
3246 fhSplitPt ->Fill(ptSplit);
3247 Float_t phi = mom.Phi();
3248 if(phi<0) phi+=TMath::TwoPi();
3249 fhSplitPtPhi ->Fill(ptSplit,phi);
3250 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3251 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3253 //Check split-clusters with good time window difference
3254 Double_t tof1 = cells->GetCellTime(absId1);
3255 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3258 Double_t tof2 = cells->GetCellTime(absId2);
3259 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3262 Double_t t12diff = tof1-tof2;
3263 fhEPairDiffTime->Fill(e1+e2, t12diff);
3267 fhMCSplitE [mcIndex]->Fill( eSplit);
3268 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3269 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3270 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3271 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3272 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3274 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3275 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3276 fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3280 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3281 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3285 //-----------------------
3286 //Create AOD for analysis
3288 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3289 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3290 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3292 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3293 aodpi0.SetLabel(calo->GetLabel());
3295 //Set the indeces of the original caloclusters
3296 aodpi0.SetCaloLabel(calo->GetID(),-1);
3297 aodpi0.SetDetector(fCalorimeter);
3299 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3300 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3301 else aodpi0.SetDistToBad(0) ;
3303 // Check if cluster is pi0 via cluster splitting
3304 aodpi0.SetIdentifiedParticleType(idPartType);
3306 // Add number of local maxima to AOD, method name in AOD to be FIXED
3307 aodpi0.SetFiducialArea(nMaxima);
3311 //Fill some histograms about shower shape
3312 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3314 FillSelectedClusterHistograms(calo, aodpi0.Pt(), nMaxima, tag, asy);
3317 // Fill histograms to undertand pile-up before other cuts applied
3318 // Remember to relax time cuts in the reader
3319 Double_t tofcluster = calo->GetTOF()*1e9;
3320 Double_t tofclusterUS = TMath::Abs(tofcluster);
3322 FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
3324 Int_t id = GetReader()->GetTriggerClusterId();
3325 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
3327 Float_t phicluster = aodpi0.Phi();
3328 if(phicluster < 0) phicluster+=TMath::TwoPi();
3332 if (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
3333 else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
3334 else fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
3337 Int_t bc = GetReader()->GetTriggerClusterBC();
3338 if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
3340 if(GetReader()->IsTriggerMatched())
3342 if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
3343 fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
3344 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
3348 if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
3349 fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
3353 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(calo->E(), tofcluster);
3354 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
3355 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(calo->E(), tofcluster);
3359 else if(TMath::Abs(bc) >= 6)
3360 Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
3363 //Add AOD with pi0 object to aod branch
3364 AddAODParticle(aodpi0);
3368 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3371 //______________________________________________
3372 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3374 //Do analysis and fill histograms
3376 if(!GetOutputAODBranch())
3378 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3380 //Loop on stored AOD pi0
3381 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3382 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3384 Float_t cen = GetEventCentrality();
3385 Float_t ep = GetEventPlaneAngle();
3387 for(Int_t iaod = 0; iaod < naod ; iaod++)
3390 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3391 Int_t pdg = pi0->GetIdentifiedParticleType();
3393 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
3395 //Fill pi0 histograms
3396 Float_t ener = pi0->E();
3397 Float_t pt = pi0->Pt();
3398 Float_t phi = pi0->Phi();
3399 if(phi < 0) phi+=TMath::TwoPi();
3400 Float_t eta = pi0->Eta();
3405 fhPtEta ->Fill(pt ,eta);
3406 fhPtPhi ->Fill(pt ,phi);
3407 fhEtaPhi ->Fill(eta ,phi);
3409 fhPtCentrality ->Fill(pt,cen) ;
3410 fhPtEventPlane ->Fill(pt,ep ) ;
3414 Int_t tag = pi0->GetTag();
3415 Int_t mcIndex = GetMCIndex(tag);
3417 fhMCE [mcIndex] ->Fill(ener);
3418 fhMCPt [mcIndex] ->Fill(pt);
3419 fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3420 fhMCPtEta[mcIndex] ->Fill(pt,eta);
3422 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3424 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
3426 Float_t efracMC = 0;
3427 Int_t label = pi0->GetLabel();
3428 Int_t momlabel = -1;
3431 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3434 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3436 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3437 if(grandmom.E() > 0 && ok)
3439 efracMC = grandmom.E()/ener;
3440 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3443 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3445 fhMCPi0DecayPt->Fill(pt);
3446 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3447 if(grandmom.E() > 0 && ok)
3449 efracMC = mom.E()/grandmom.E();
3450 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3453 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3455 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3456 if(grandmom.E() > 0 && ok)
3458 efracMC = grandmom.E()/ener;
3459 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3462 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3464 fhMCEtaDecayPt->Fill(pt);
3465 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3466 if(grandmom.E() > 0 && ok)
3468 efracMC = mom.E()/grandmom.E();
3469 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3472 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3474 fhMCOtherDecayPt->Fill(pt);
3479 }//Histograms with MC
3485 //__________________________________________________________________
3486 void AliAnaPi0EbE::Print(const Option_t * opt) const
3488 //Print some relevant parameters set for the analysis
3492 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3493 AliAnaCaloTrackCorrBaseClass::Print("");
3494 printf("Analysis Type = %d \n", fAnaType) ;
3495 if(fAnaType == kSSCalo)
3497 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3498 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3499 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3500 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);