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),
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),
108 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
109 fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0),
112 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
113 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
114 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
115 fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
116 fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
117 fhTrackMatchedMCParticlePt(0),
118 fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
119 fhdEdx(0), fhEOverP(0), fhEOverPNoTRD(0),
120 // Number of local maxima in cluster
121 fhNLocMaxPt(0), fhNLocMaxPtReject(0),
123 fhTimePtNoCut(0), fhTimePtSPD(0), fhTimePtSPDMulti(0),
124 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
125 fhTimeNPileUpVertContributors(0),
126 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
127 fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
128 fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
129 fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
133 for(Int_t i = 0; i < 6; i++)
139 fhMCPtCentrality [i] = 0;
143 fhMCSplitPtPhi [i] = 0;
144 fhMCSplitPtEta [i] = 0;
146 fhMCNLocMaxPt [i] = 0;
147 fhMCNLocMaxSplitPt [i] = 0;
148 fhMCNLocMaxPtReject[i] = 0;
150 fhMCPtLambda0 [i] = 0;
151 fhMCPtLambda0NoTRD [i] = 0;
152 fhMCPtLambda0FracMaxCellCut[i]= 0;
153 fhMCPtFracMaxCell [i] = 0;
154 fhMCPtLambda1 [i] = 0;
155 fhMCPtDispersion [i] = 0;
157 fhMCPtDispEta [i] = 0;
158 fhMCPtDispPhi [i] = 0;
159 fhMCPtSumEtaPhi [i] = 0;
160 fhMCPtDispEtaPhiDiff[i] = 0;
161 fhMCPtSphericity [i] = 0;
162 fhMCPtAsymmetry [i] = 0;
165 fhMCMassSplitPt [i]=0;
166 fhMCSelectedMassPt [i]=0;
167 fhMCSelectedMassSplitPt[i]=0;
169 fhMCMassPtNoOverlap [i]=0;
170 fhMCMassSplitPtNoOverlap [i]=0;
171 fhMCSelectedMassPtNoOverlap [i]=0;
172 fhMCSelectedMassSplitPtNoOverlap[i]=0;
174 for(Int_t j = 0; j < 7; j++)
176 fhMCLambda0DispEta [j][i] = 0;
177 fhMCLambda0DispPhi [j][i] = 0;
178 fhMCDispEtaDispPhi [j][i] = 0;
179 fhMCAsymmetryLambda0 [j][i] = 0;
180 fhMCAsymmetryDispEta [j][i] = 0;
181 fhMCAsymmetryDispPhi [j][i] = 0;
185 for(Int_t j = 0; j < 7; j++)
187 fhLambda0DispEta [j] = 0;
188 fhLambda0DispPhi [j] = 0;
189 fhDispEtaDispPhi [j] = 0;
190 fhAsymmetryLambda0 [j] = 0;
191 fhAsymmetryDispEta [j] = 0;
192 fhAsymmetryDispPhi [j] = 0;
197 for(Int_t i = 0; i < 3; i++)
199 fhPtLambda0LocMax [i] = 0;
200 fhPtLambda1LocMax [i] = 0;
201 fhPtDispersionLocMax [i] = 0;
202 fhPtDispEtaLocMax [i] = 0;
203 fhPtDispPhiLocMax [i] = 0;
204 fhPtSumEtaPhiLocMax [i] = 0;
205 fhPtDispEtaPhiDiffLocMax[i] = 0;
206 fhPtSphericityLocMax [i] = 0;
207 fhPtAsymmetryLocMax [i] = 0;
208 fhMassPtLocMax [i] = 0;
209 fhSelectedMassPtLocMax [i] = 0;
210 for(Int_t ipart = 0; ipart<6; ipart++)
212 fhMCPtLambda0LocMax [ipart][i] = 0;
213 fhMCSelectedMassPtLocMax[ipart][i] = 0;
216 fhMCPi0PtRecoPtPrimLocMax [i] = 0;
217 fhMCEtaPtRecoPtPrimLocMax [i] = 0;
218 fhMCPi0SplitPtRecoPtPrimLocMax [i] = 0;
219 fhMCEtaSplitPtRecoPtPrimLocMax [i] = 0;
221 fhMCPi0SelectedPtRecoPtPrimLocMax [i] = 0;
222 fhMCEtaSelectedPtRecoPtPrimLocMax [i] = 0;
223 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[i] = 0;
224 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[i] = 0;
229 for(Int_t i =0; i < 14; i++){
230 fhLambda0ForW0[i] = 0;
231 //fhLambda1ForW0[i] = 0;
232 if(i<8)fhMassPairLocMax[i] = 0;
235 for(Int_t i = 0; i < 11; i++)
237 fhEtaPhiTriggerEMCALBC [i] = 0 ;
238 fhTimeTriggerEMCALBC [i] = 0 ;
239 fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
241 fhEtaPhiTriggerEMCALBCUM [i] = 0 ;
242 fhTimeTriggerEMCALBCUM [i] = 0 ;
246 for(Int_t iSM = 0; iSM < 22; iSM++)
248 fhNLocMaxPtSM[iSM] = 0;
249 for(Int_t inlm = 0; inlm < 3; inlm++)
251 fhSelectedMassPtLocMaxSM [inlm][iSM] = 0;
252 fhSelectedLambda0PtLocMaxSM [inlm][iSM] = 0;
255 //Initialize parameters
260 //___________________________________________________________________________________
261 void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
263 // Fill some histograms to understand pile-up
264 if(!fFillPileUpHistograms) return;
266 //printf("E %f, time %f\n",energy,time);
267 AliVEvent * event = GetReader()->GetInputEvent();
269 fhTimePtNoCut->Fill(pt,time);
270 if(GetReader()->IsPileUpFromSPD())
272 if(GetReader()->IsPileUpFromSPD()) { fhPtPileUp[0]->Fill(pt); fhTimePtSPD ->Fill(pt,time); }
273 if(GetReader()->IsPileUpFromEMCal()) fhPtPileUp[1]->Fill(pt);
274 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPileUp[2]->Fill(pt);
275 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPileUp[3]->Fill(pt);
276 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPileUp[4]->Fill(pt);
277 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPileUp[5]->Fill(pt);
278 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
280 if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
284 AliVCaloCells* cells = 0;
285 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
286 else cells = GetPHOSCells();
288 Float_t maxCellFraction = 0.;
289 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
291 Double_t tmax = cells->GetCellTime(absIdMax);
292 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
295 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
296 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
298 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
300 Int_t absId = calo->GetCellsAbsId()[ipos];
302 if( absId == absIdMax ) continue ;
304 Double_t timecell = cells->GetCellTime(absId);
305 Float_t amp = cells->GetCellAmplitude(absId);
306 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
307 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
310 Float_t diff = (tmax-timecell);
312 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
314 if(GetReader()->IsPileUpFromSPD())
316 fhPtCellTimePileUp[0]->Fill(pt, timecell);
317 fhPtTimeDiffPileUp[0]->Fill(pt, diff);
320 if(GetReader()->IsPileUpFromEMCal())
322 fhPtCellTimePileUp[1]->Fill(pt, timecell);
323 fhPtTimeDiffPileUp[1]->Fill(pt, diff);
326 if(GetReader()->IsPileUpFromSPDOrEMCal())
328 fhPtCellTimePileUp[2]->Fill(pt, timecell);
329 fhPtTimeDiffPileUp[2]->Fill(pt, diff);
332 if(GetReader()->IsPileUpFromSPDAndEMCal())
334 fhPtCellTimePileUp[3]->Fill(pt, timecell);
335 fhPtTimeDiffPileUp[3]->Fill(pt, diff);
338 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
340 fhPtCellTimePileUp[4]->Fill(pt, timecell);
341 fhPtTimeDiffPileUp[4]->Fill(pt, diff);
344 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
346 fhPtCellTimePileUp[5]->Fill(pt, timecell);
347 fhPtTimeDiffPileUp[5]->Fill(pt, diff);
350 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
352 fhPtCellTimePileUp[6]->Fill(pt, timecell);
353 fhPtTimeDiffPileUp[6]->Fill(pt, diff);
358 if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
360 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
361 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
363 // N pile up vertices
369 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
370 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
375 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
376 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
379 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
380 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
382 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
383 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
385 if(TMath::Abs(time) < 25)
387 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
388 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
391 if(time < 75 && time > -25)
393 fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
394 fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
397 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
398 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
401 Float_t z1 = -1, z2 = -1;
403 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
407 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
408 ncont=pv->GetNContributors();
409 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
411 diamZ = esdEv->GetDiamondZ();
415 AliAODVertex *pv=aodEv->GetVertex(iVert);
416 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
417 ncont=pv->GetNContributors();
418 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
420 diamZ = aodEv->GetDiamondZ();
423 Double_t distZ = TMath::Abs(z2-z1);
424 diamZ = TMath::Abs(z2-diamZ);
426 fhTimeNPileUpVertContributors ->Fill(time,ncont);
427 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
428 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
434 //______________________________________________________________________________________________
435 void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
437 // Fill histograms that do not pass the identification (SS case only)
439 Float_t ener = mom.E();
440 Float_t pt = mom.Pt();
441 Float_t phi = mom.Phi();
442 if(phi < 0) phi+=TMath::TwoPi();
443 Float_t eta = mom.Eta();
445 fhPtReject ->Fill(pt);
446 fhEReject ->Fill(ener);
448 fhPtEtaReject ->Fill(ener,eta);
449 fhPtPhiReject ->Fill(ener,phi);
450 fhEtaPhiReject ->Fill(eta,phi);
452 fhNLocMaxPtReject->Fill(pt,nMaxima);
456 Int_t mcIndex = GetMCIndex(mctag);
457 fhMCEReject [mcIndex] ->Fill(ener);
458 fhMCPtReject [mcIndex] ->Fill(pt);
459 fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
463 //___________________________________________________________________________________
464 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
465 Int_t tag, Float_t asy)
467 // Fill shower shape, timing and other histograms for selected clusters from decay
469 Float_t ener = cluster->E();
470 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
471 Float_t l0 = cluster->GetM02();
472 Float_t l1 = cluster->GetM20();
473 Int_t nSM = GetModuleNumber(cluster);
476 if (pt < 2 ) ptbin = 0;
477 else if (pt < 4 ) ptbin = 1;
478 else if (pt < 6 ) ptbin = 2;
479 else if (pt < 10) ptbin = 3;
480 else if (pt < 15) ptbin = 4;
481 else if (pt < 20) ptbin = 5;
485 if (nMaxima==1) indexMax = 0 ;
486 else if(nMaxima==2) indexMax = 1 ;
490 AliVCaloCells * cell = 0x0;
491 if(fCalorimeter == "PHOS")
492 cell = GetPHOSCells();
494 cell = GetEMCALCells();
496 Float_t maxCellFraction = 0;
497 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
498 fhPtFracMaxCell->Fill(pt,maxCellFraction);
500 FillWeightHistograms(cluster);
502 fhPtDispersion->Fill(pt, disp);
503 fhPtLambda0 ->Fill(pt, l0 );
504 fhPtLambda1 ->Fill(pt, l1 );
506 Float_t ll0 = 0., ll1 = 0.;
507 Float_t dispp= 0., dEta = 0., dPhi = 0.;
508 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
509 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
511 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
512 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
514 fhPtDispEta -> Fill(pt,dEta);
515 fhPtDispPhi -> Fill(pt,dPhi);
516 fhPtSumEta -> Fill(pt,sEta);
517 fhPtSumPhi -> Fill(pt,sPhi);
518 fhPtSumEtaPhi -> Fill(pt,sEtaPhi);
519 fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
520 if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
522 fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
523 fhLambda0DispEta[ptbin]->Fill(l0 ,dEta);
524 fhLambda0DispPhi[ptbin]->Fill(l0 ,dPhi);
526 if (fAnaType==kSSCalo)
528 // Asymmetry histograms
529 fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
530 fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
531 fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
535 fhNLocMaxPt->Fill(pt,nMaxima);
537 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
538 fhNLocMaxPtSM[nSM]->Fill(pt,nMaxima);
540 fhPtLambda0LocMax [indexMax]->Fill(pt,l0);
541 fhPtLambda1LocMax [indexMax]->Fill(pt,l1);
542 fhPtDispersionLocMax[indexMax]->Fill(pt,disp);
544 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
546 fhPtDispEtaLocMax [indexMax]-> Fill(pt,dEta);
547 fhPtDispPhiLocMax [indexMax]-> Fill(pt,dPhi);
548 fhPtSumEtaPhiLocMax [indexMax]-> Fill(pt,sEtaPhi);
549 fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
550 if(dEta+dPhi>0) fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
551 if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt ,asy);
555 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
556 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
558 fhPtLambda0NoTRD ->Fill(pt, l0 );
559 fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
562 if(maxCellFraction < 0.5)
563 fhPtLambda0FracMaxCellCut->Fill(pt, l0 );
565 fhPtTime ->Fill(pt, cluster->GetTOF()*1.e9);
566 fhPtNCells->Fill(pt, cluster->GetNCells());
568 // Fill Track matching control histograms
571 Float_t dZ = cluster->GetTrackDz();
572 Float_t dR = cluster->GetTrackDx();
574 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
576 dR = 2000., dZ = 2000.;
577 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
579 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
581 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
583 Bool_t positive = kFALSE;
584 if(track) positive = (track->Charge()>0);
586 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
588 fhTrackMatchedDEta->Fill(pt,dZ);
589 fhTrackMatchedDPhi->Fill(pt,dR);
590 if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
596 fhTrackMatchedDEtaPos->Fill(pt,dZ);
597 fhTrackMatchedDPhiPos->Fill(pt,dR);
598 if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
602 fhTrackMatchedDEtaNeg->Fill(pt,dZ);
603 fhTrackMatchedDPhiNeg->Fill(pt,dR);
604 if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
608 // Check dEdx and E/p of matched clusters
610 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
614 Float_t dEdx = track->GetTPCsignal();
615 fhdEdx->Fill(pt, dEdx);
617 Float_t eOverp = cluster->E()/track->P();
618 fhEOverP->Fill(pt, eOverp);
620 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
621 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
622 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
623 fhEOverPNoTRD->Fill(pt, eOverp);
627 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
632 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
634 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
635 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 2.5 ;
636 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 0.5 ;
637 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 1.5 ;
643 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
644 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 6.5 ;
645 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 4.5 ;
646 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 5.5 ;
650 fhTrackMatchedMCParticlePt ->Fill(pt, mctag);
651 fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
652 fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
656 }// Track matching histograms
660 Int_t mcIndex = GetMCIndex(tag);
662 fhMCPtLambda0[mcIndex] ->Fill(pt, l0);
663 fhMCPtLambda1[mcIndex] ->Fill(pt, l1);
664 fhMCPtDispersion[mcIndex] ->Fill(pt, disp);
665 fhMCPtFracMaxCell[mcIndex]->Fill(pt,maxCellFraction);
667 fhMCPtLambda0LocMax [mcIndex][indexMax]->Fill(pt,l0);
669 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
670 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
671 fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0 );
673 if(maxCellFraction < 0.5)
674 fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0 );
676 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
678 fhMCPtDispEta [mcIndex]-> Fill(pt,dEta);
679 fhMCPtDispPhi [mcIndex]-> Fill(pt,dPhi);
680 fhMCPtSumEtaPhi [mcIndex]-> Fill(pt,sEtaPhi);
681 fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
682 if(dEta+dPhi>0)fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
684 if (fAnaType==kSSCalo)
686 fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
687 fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
688 fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
691 fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
692 fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0 ,dEta);
693 fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0 ,dPhi);
701 //________________________________________________________
702 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
704 // Calculate weights and fill histograms
706 if(!fFillWeightHistograms || GetMixedEvent()) return;
708 AliVCaloCells* cells = 0;
709 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
710 else cells = GetPHOSCells();
712 // First recalculate energy in case non linearity was applied
715 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
718 Int_t id = clus->GetCellsAbsId()[ipos];
720 //Recalibrate cell energy if needed
721 Float_t amp = cells->GetCellAmplitude(id);
722 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
733 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
737 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
738 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
740 //Get the ratio and log ratio to all cells in cluster
741 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
743 Int_t id = clus->GetCellsAbsId()[ipos];
745 //Recalibrate cell energy if needed
746 Float_t amp = cells->GetCellAmplitude(id);
747 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
749 fhECellClusterRatio ->Fill(energy,amp/energy);
750 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
753 //Recalculate shower shape for different W0
754 if(fCalorimeter=="EMCAL"){
756 Float_t l0org = clus->GetM02();
757 Float_t l1org = clus->GetM20();
758 Float_t dorg = clus->GetDispersion();
760 for(Int_t iw = 0; iw < 14; iw++)
762 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
763 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
765 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
766 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
770 // Set the original values back
773 clus->SetDispersion(dorg);
778 //__________________________________________
779 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
781 //Save parameters used for analysis
782 TString parList ; //this will be list of parameters used for this analysis.
783 const Int_t buffersize = 255;
784 char onePar[buffersize] ;
786 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
788 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
791 if(fAnaType == kSSCalo)
793 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
795 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
797 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
799 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
803 //Get parameters set in base class.
804 parList += GetBaseParametersList() ;
806 //Get parameters set in PID class.
807 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
809 return new TObjString(parList) ;
812 //_____________________________________________
813 TList * AliAnaPi0EbE::GetCreateOutputObjects()
815 // Create histograms to be saved in output file and
816 // store them in outputContainer
817 TList * outputContainer = new TList() ;
818 outputContainer->SetName("Pi0EbEHistos") ;
820 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
821 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
822 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
823 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
824 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
825 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
826 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
828 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
829 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
830 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
832 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
833 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
834 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
835 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
836 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
837 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
839 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
840 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
841 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
842 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
843 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
844 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
846 Int_t ntimptbins = GetHistogramRanges()->GetHistoTimeBins();
847 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
848 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
850 TString nlm[] = {"1 Local Maxima","2 Local Maxima", "NLM > 2"};
851 TString ptype[] = {"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
852 TString pname[] = {"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
853 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
855 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
856 fhPt->SetYTitle("#it{N}");
857 fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
858 outputContainer->Add(fhPt) ;
860 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
861 fhE->SetYTitle("#it{N}");
862 fhE->SetXTitle("#it{E} (GeV)");
863 outputContainer->Add(fhE) ;
866 ("hPtPhi","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
867 fhPtPhi->SetYTitle("#phi (rad)");
868 fhPtPhi->SetXTitle("#it{E} (GeV)");
869 outputContainer->Add(fhPtPhi) ;
872 ("hPtEta","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
873 fhPtEta->SetYTitle("#eta");
874 fhPtEta->SetXTitle("#it{E} (GeV)");
875 outputContainer->Add(fhPtEta) ;
878 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
879 fhEtaPhi->SetYTitle("#phi (rad)");
880 fhEtaPhi->SetXTitle("#eta");
881 outputContainer->Add(fhEtaPhi) ;
883 if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
885 fhEtaPhiEMCALBC0 = new TH2F
886 ("hEtaPhiEMCALBC0","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
887 fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
888 fhEtaPhiEMCALBC0->SetXTitle("#eta");
889 outputContainer->Add(fhEtaPhiEMCALBC0) ;
891 fhEtaPhiEMCALBC1 = new TH2F
892 ("hEtaPhiEMCALBC1","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |#it{t}| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
893 fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
894 fhEtaPhiEMCALBC1->SetXTitle("#eta");
895 outputContainer->Add(fhEtaPhiEMCALBC1) ;
897 fhEtaPhiEMCALBCN = new TH2F
898 ("hEtaPhiEMCALBCN","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
899 fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
900 fhEtaPhiEMCALBCN->SetXTitle("#eta");
901 outputContainer->Add(fhEtaPhiEMCALBCN) ;
903 for(Int_t i = 0; i < 11; i++)
905 fhEtaPhiTriggerEMCALBC[i] = new TH2F
906 (Form("hEtaPhiTriggerEMCALBC%d",i-5),
907 Form("meson #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
908 netabins,etamin,etamax,nphibins,phimin,phimax);
909 fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
910 fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
911 outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
913 fhTimeTriggerEMCALBC[i] = new TH2F
914 (Form("hTimeTriggerEMCALBC%d",i-5),
915 Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
916 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
917 fhTimeTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
918 fhTimeTriggerEMCALBC[i]->SetYTitle("#it{t} (ns)");
919 outputContainer->Add(fhTimeTriggerEMCALBC[i]);
921 fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
922 (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
923 Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
924 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
925 fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
926 fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{t} (ns)");
927 outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
929 fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
930 (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
931 Form("meson #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
932 netabins,etamin,etamax,nphibins,phimin,phimax);
933 fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
934 fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
935 outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
937 fhTimeTriggerEMCALBCUM[i] = new TH2F
938 (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
939 Form("meson #it{t} vs #it{E}, unmatched trigger EMCAL-BC=%d",i-5),
940 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
941 fhTimeTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
942 fhTimeTriggerEMCALBCUM[i]->SetYTitle("#it{t} (ns)");
943 outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
947 fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
948 "cluster #it{t} vs #it{E} of clusters, no match, rematch open time",
949 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
950 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
951 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{t} (ns)");
952 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
955 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
956 "cluster #it{t} vs #it{E} of clusters, no match, rematch with neigbour parches",
957 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
958 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
959 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{t} (ns)");
960 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
962 fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
963 "cluster #it{t} vs #it{E} of clusters, no match, rematch open time and neigbour",
964 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
965 fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
966 fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{t} (ns)");
967 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
971 fhPtCentrality = new TH2F("hPtCentrality","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
972 fhPtCentrality->SetYTitle("centrality");
973 fhPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
974 outputContainer->Add(fhPtCentrality) ;
976 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
977 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
978 fhPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
979 outputContainer->Add(fhPtEventPlane) ;
981 if(fAnaType == kSSCalo)
983 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
984 fhPtReject->SetYTitle("#it{N}");
985 fhPtReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
986 outputContainer->Add(fhPtReject) ;
988 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
989 fhEReject->SetYTitle("#it{N}");
990 fhEReject->SetXTitle("#it{E} (GeV)");
991 outputContainer->Add(fhEReject) ;
993 fhPtPhiReject = new TH2F
994 ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
995 fhPtPhiReject->SetYTitle("#phi (rad)");
996 fhPtPhiReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
997 outputContainer->Add(fhPtPhiReject) ;
999 fhPtEtaReject = new TH2F
1000 ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1001 fhPtEtaReject->SetYTitle("#eta");
1002 fhPtEtaReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1003 outputContainer->Add(fhPtEtaReject) ;
1005 fhEtaPhiReject = new TH2F
1006 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1007 fhEtaPhiReject->SetYTitle("#phi (rad)");
1008 fhEtaPhiReject->SetXTitle("#eta");
1009 outputContainer->Add(fhEtaPhiReject) ;
1013 ("hMass","all pairs #it{M}: #it{E} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1014 fhMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1015 fhMass->SetXTitle("#it{E} (GeV)");
1016 outputContainer->Add(fhMass) ;
1018 fhSelectedMass = new TH2F
1019 ("hSelectedMass","Selected #pi^{0} (#eta) pairs #it{M}: E vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1020 fhSelectedMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1021 fhSelectedMass->SetXTitle("#it{E} (GeV)");
1022 outputContainer->Add(fhSelectedMass) ;
1024 if(fAnaType == kSSCalo)
1028 ("hMassPt","all pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1029 fhMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1030 fhMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1031 outputContainer->Add(fhMassPt) ;
1033 fhSelectedMassPt = new TH2F
1034 ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1035 fhSelectedMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1036 fhSelectedMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1037 outputContainer->Add(fhSelectedMassPt) ;
1039 for(Int_t inlm = 0; inlm < 3; inlm++)
1041 fhMassPtLocMax[inlm] = new TH2F
1042 (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs #it{M}: #it{p}_{T} vs #it{M} and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1043 fhMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1044 fhMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1045 outputContainer->Add(fhMassPtLocMax[inlm]) ;
1047 fhSelectedMassPtLocMax[inlm] = new TH2F
1048 (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1049 fhSelectedMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1050 fhSelectedMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1051 outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
1053 for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1055 fhSelectedMassPtLocMaxSM[inlm][iSM] = new TH2F
1056 (Form("hSelectedMassPtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1057 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1058 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1059 outputContainer->Add(fhSelectedMassPtLocMaxSM[inlm][iSM]) ;
1061 fhSelectedLambda0PtLocMaxSM[inlm][iSM] = new TH2F
1062 (Form("hSelectedLambda0PtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #lambda_{0}^{2}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1063 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetYTitle("#lambda_{0}^{2}");
1064 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1065 outputContainer->Add(fhSelectedLambda0PtLocMaxSM[inlm][iSM]) ;
1070 for(Int_t ipart = 0; ipart < 6; ipart++)
1072 fhMCSelectedMassPtLocMax[ipart][inlm] = new TH2F
1073 (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
1074 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s, %s",nlm[inlm].Data(),pname[ipart].Data()),
1075 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1076 fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1077 fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1078 outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
1085 fhMassNoOverlap = new TH2F
1086 ("hMassNoOverlap","all pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1087 fhMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1088 fhMassNoOverlap->SetXTitle("#it{E} (GeV)");
1089 outputContainer->Add(fhMassNoOverlap) ;
1091 fhSelectedMassNoOverlap = new TH2F
1092 ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1093 fhSelectedMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1094 fhSelectedMassNoOverlap->SetXTitle("#it{E} (GeV)");
1095 outputContainer->Add(fhSelectedMassNoOverlap) ;
1097 fhMassPtNoOverlap = new TH2F
1098 ("hMassPtNoOverlap","all pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1099 fhMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1100 fhMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1101 outputContainer->Add(fhMassPtNoOverlap) ;
1103 fhSelectedMassPtNoOverlap = new TH2F
1104 ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1105 fhSelectedMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1106 fhSelectedMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1107 outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1111 if(fAnaType != kSSCalo)
1113 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1114 fhPtDecay->SetYTitle("#it{N}");
1115 fhPtDecay->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1116 outputContainer->Add(fhPtDecay) ;
1118 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1119 fhEDecay->SetYTitle("#it{N}");
1120 fhEDecay->SetXTitle("#it{E} (GeV)");
1121 outputContainer->Add(fhEDecay) ;
1126 if( fFillSelectClHisto )
1128 fhPtDispersion = new TH2F
1129 ("hPtDispersion","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1130 fhPtDispersion->SetYTitle("D^{2}");
1131 fhPtDispersion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1132 outputContainer->Add(fhPtDispersion) ;
1134 fhPtLambda0 = new TH2F
1135 ("hPtLambda0","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1136 fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1137 fhPtLambda0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1138 outputContainer->Add(fhPtLambda0) ;
1140 fhPtLambda1 = new TH2F
1141 ("hPtLambda1","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1142 fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1143 fhPtLambda1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1144 outputContainer->Add(fhPtLambda1) ;
1146 fhPtLambda0FracMaxCellCut = new TH2F
1147 ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1148 fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1149 fhPtLambda0FracMaxCellCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1150 outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1152 fhPtFracMaxCell = new TH2F
1153 ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1154 fhPtFracMaxCell->SetYTitle("Fraction");
1155 fhPtFracMaxCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1156 outputContainer->Add(fhPtFracMaxCell) ;
1158 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >=0 )
1160 fhPtLambda0NoTRD = new TH2F
1161 ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1162 fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1163 fhPtLambda0NoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1164 outputContainer->Add(fhPtLambda0NoTRD) ;
1166 fhPtFracMaxCellNoTRD = new TH2F
1167 ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1168 fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1169 fhPtFracMaxCellNoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1170 outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1172 if(!fFillOnlySimpleSSHisto)
1174 fhPtDispEta = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1175 fhPtDispEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1176 fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1177 outputContainer->Add(fhPtDispEta);
1179 fhPtDispPhi = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1180 fhPtDispPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1181 fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1182 outputContainer->Add(fhPtDispPhi);
1184 fhPtSumEta = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1185 fhPtSumEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1186 fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1187 outputContainer->Add(fhPtSumEta);
1189 fhPtSumPhi = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs #it{p}_{T}",
1190 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1191 fhPtSumPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1192 fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1193 outputContainer->Add(fhPtSumPhi);
1195 fhPtSumEtaPhi = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",
1196 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1197 fhPtSumEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1198 fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1199 outputContainer->Add(fhPtSumEtaPhi);
1201 fhPtDispEtaPhiDiff = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",
1202 nptbins,ptmin,ptmax,200, -10,10);
1203 fhPtDispEtaPhiDiff->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1204 fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1205 outputContainer->Add(fhPtDispEtaPhiDiff);
1207 fhPtSphericity = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs #it{p}_{T} (GeV/#it{c})",
1208 nptbins,ptmin,ptmax, 200, -1,1);
1209 fhPtSphericity->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1210 fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1211 outputContainer->Add(fhPtSphericity);
1213 for(Int_t i = 0; i < 7; i++)
1215 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]),
1216 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1217 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1218 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1219 outputContainer->Add(fhDispEtaDispPhi[i]);
1221 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]),
1222 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1223 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1224 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1225 outputContainer->Add(fhLambda0DispEta[i]);
1227 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]),
1228 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1229 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1230 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1231 outputContainer->Add(fhLambda0DispPhi[i]);
1237 fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1238 nptbins,ptmin,ptmax,20,0,20);
1239 fhNLocMaxPt ->SetYTitle("N maxima");
1240 fhNLocMaxPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1241 outputContainer->Add(fhNLocMaxPt) ;
1243 for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1245 fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
1246 nptbins,ptmin,ptmax,20,0,20);
1247 fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
1248 fhNLocMaxPtSM[iSM] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1249 outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
1252 if(fAnaType == kSSCalo)
1255 fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1256 nptbins,ptmin,ptmax,20,0,20);
1257 fhNLocMaxPtReject ->SetYTitle("N maxima");
1258 fhNLocMaxPtReject ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1259 outputContainer->Add(fhNLocMaxPtReject) ;
1262 for (Int_t i = 0; i < 3; i++)
1264 fhPtLambda0LocMax[i] = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1265 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1266 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1267 fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1268 fhPtLambda0LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1269 outputContainer->Add(fhPtLambda0LocMax[i]) ;
1273 for(Int_t ipart = 0; ipart < 6; ipart++)
1275 fhMCPtLambda0LocMax[ipart][i] = new TH2F
1276 (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1277 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),pname[ipart].Data()),
1278 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1279 fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1280 fhMCPtLambda0LocMax[ipart][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1281 outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1285 fhPtLambda1LocMax[i] = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1286 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1287 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1288 fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1289 fhPtLambda1LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1290 outputContainer->Add(fhPtLambda1LocMax[i]) ;
1292 fhPtDispersionLocMax[i] = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1293 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1294 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1295 fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1296 fhPtDispersionLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1297 outputContainer->Add(fhPtDispersionLocMax[i]) ;
1299 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1301 fhPtDispEtaLocMax[i] = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1302 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1303 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1304 fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1305 fhPtDispEtaLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1306 outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1308 fhPtDispPhiLocMax[i] = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1309 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1310 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1311 fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1312 fhPtDispPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1313 outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1315 fhPtSumEtaPhiLocMax[i] = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1316 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1317 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1318 fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1319 fhPtSumEtaPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1320 outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1322 fhPtDispEtaPhiDiffLocMax[i] = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1323 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1324 nptbins,ptmin,ptmax,200, -10,10);
1325 fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1326 fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1327 outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1329 fhPtSphericityLocMax[i] = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1330 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1331 nptbins,ptmin,ptmax,200, -1,1);
1332 fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1333 fhPtSphericityLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1334 outputContainer->Add(fhPtSphericityLocMax[i]) ;
1339 fhPtNCells = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1340 fhPtNCells->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1341 fhPtNCells->SetYTitle("# of cells in cluster");
1342 outputContainer->Add(fhPtNCells);
1344 fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1345 fhPtTime->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1346 fhPtTime->SetYTitle("t (ns)");
1347 outputContainer->Add(fhPtTime);
1352 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1353 fhEPairDiffTime->SetXTitle("#it{E}_{pair} (GeV)");
1354 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1355 outputContainer->Add(fhEPairDiffTime);
1357 if(fAnaType == kIMCalo)
1359 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1360 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1361 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1362 "2 Local Maxima paired with more than 2 Local Maxima",
1363 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1365 for (Int_t i = 0; i < 8 ; i++)
1368 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1370 fhMassPairLocMax[i] = new TH2F
1371 (Form("MassPairLocMax%s",combiName[i].Data()),
1372 Form("#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1373 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1374 fhMassPairLocMax[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1375 fhMassPairLocMax[i]->SetXTitle("#it{E}_{pair} (GeV)");
1376 outputContainer->Add(fhMassPairLocMax[i]) ;
1382 fhTrackMatchedDEta = new TH2F
1383 ("hTrackMatchedDEta",
1384 "d#eta of cluster-track vs cluster #it{p}_{T}",
1385 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1386 fhTrackMatchedDEta->SetYTitle("d#eta");
1387 fhTrackMatchedDEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1389 fhTrackMatchedDPhi = new TH2F
1390 ("hTrackMatchedDPhi",
1391 "d#phi of cluster-track vs cluster #it{p}_{T}",
1392 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1393 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1394 fhTrackMatchedDPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1396 fhTrackMatchedDEtaDPhi = new TH2F
1397 ("hTrackMatchedDEtaDPhi",
1398 "d#eta vs d#phi of cluster-track",
1399 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1400 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1401 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1403 outputContainer->Add(fhTrackMatchedDEta) ;
1404 outputContainer->Add(fhTrackMatchedDPhi) ;
1405 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1407 fhTrackMatchedDEtaPos = new TH2F
1408 ("hTrackMatchedDEtaPos",
1409 "d#eta of cluster-track vs cluster #it{p}_{T}",
1410 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1411 fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1412 fhTrackMatchedDEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1414 fhTrackMatchedDPhiPos = new TH2F
1415 ("hTrackMatchedDPhiPos",
1416 "d#phi of cluster-track vs cluster #it{p}_{T}",
1417 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1418 fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1419 fhTrackMatchedDPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1421 fhTrackMatchedDEtaDPhiPos = new TH2F
1422 ("hTrackMatchedDEtaDPhiPos",
1423 "d#eta vs d#phi of cluster-track",
1424 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1425 fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1426 fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1428 outputContainer->Add(fhTrackMatchedDEtaPos) ;
1429 outputContainer->Add(fhTrackMatchedDPhiPos) ;
1430 outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1432 fhTrackMatchedDEtaNeg = new TH2F
1433 ("hTrackMatchedDEtaNeg",
1434 "d#eta of cluster-track vs cluster #it{p}_{T}",
1435 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1436 fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1437 fhTrackMatchedDEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1439 fhTrackMatchedDPhiNeg = new TH2F
1440 ("hTrackMatchedDPhiNeg",
1441 "d#phi of cluster-track vs cluster #it{p}_{T}",
1442 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1443 fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1444 fhTrackMatchedDPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1446 fhTrackMatchedDEtaDPhiNeg = new TH2F
1447 ("hTrackMatchedDEtaDPhiNeg",
1448 "d#eta vs d#phi of cluster-track",
1449 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1450 fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1451 fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1453 outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1454 outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1455 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1457 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1458 fhdEdx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1459 fhdEdx->SetYTitle("<#it{dE}/#it{dx}>");
1460 outputContainer->Add(fhdEdx);
1462 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1463 fhEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1464 fhEOverP->SetYTitle("#it{E}/#it{p}");
1465 outputContainer->Add(fhEOverP);
1467 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >=0)
1469 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1470 fhEOverPNoTRD->SetXTitle("#it{E} (GeV)");
1471 fhEOverPNoTRD->SetYTitle("#it{E}/#it{p}");
1472 outputContainer->Add(fhEOverPNoTRD);
1475 if(IsDataMC() && fFillTMHisto)
1477 fhTrackMatchedMCParticlePt = new TH2F
1478 ("hTrackMatchedMCParticlePt",
1479 "Origin of particle vs energy",
1480 nptbins,ptmin,ptmax,8,0,8);
1481 fhTrackMatchedMCParticlePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1482 //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1484 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1485 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1486 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1487 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1488 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1489 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1490 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1491 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1493 outputContainer->Add(fhTrackMatchedMCParticlePt);
1495 fhTrackMatchedMCParticleDEta = new TH2F
1496 ("hTrackMatchedMCParticleDEta",
1497 "Origin of particle vs #eta residual",
1498 nresetabins,resetamin,resetamax,8,0,8);
1499 fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1500 //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1502 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1503 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1504 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1505 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1506 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1507 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1508 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1509 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1511 outputContainer->Add(fhTrackMatchedMCParticleDEta);
1513 fhTrackMatchedMCParticleDPhi = new TH2F
1514 ("hTrackMatchedMCParticleDPhi",
1515 "Origin of particle vs #phi residual",
1516 nresphibins,resphimin,resphimax,8,0,8);
1517 fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1518 //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1520 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1521 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1522 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1523 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1524 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1525 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1526 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1527 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1529 outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1535 if(fFillWeightHistograms)
1537 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1538 nptbins,ptmin,ptmax, 100,0,1.);
1539 fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1540 fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
1541 outputContainer->Add(fhECellClusterRatio);
1543 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1544 nptbins,ptmin,ptmax, 100,-10,0);
1545 fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1546 fhECellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1547 outputContainer->Add(fhECellClusterLogRatio);
1549 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1550 nptbins,ptmin,ptmax, 100,0,1.);
1551 fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1552 fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
1553 outputContainer->Add(fhEMaxCellClusterRatio);
1555 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1556 nptbins,ptmin,ptmax, 100,-10,0);
1557 fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1558 fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1559 outputContainer->Add(fhEMaxCellClusterLogRatio);
1561 for(Int_t iw = 0; iw < 14; iw++)
1563 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),
1564 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1565 fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1566 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1567 outputContainer->Add(fhLambda0ForW0[iw]);
1569 // 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),
1570 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1571 // fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1572 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1573 // outputContainer->Add(fhLambda1ForW0[iw]);
1582 fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1583 fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1584 fhMCPi0PtOrigin->SetYTitle("Origin");
1585 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1586 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1587 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1588 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1589 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1590 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1591 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1592 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1593 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1594 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1595 outputContainer->Add(fhMCPi0PtOrigin) ;
1597 fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1598 fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1599 fhMCEtaPtOrigin->SetYTitle("Origin");
1600 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1601 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1602 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1603 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1604 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1605 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1606 outputContainer->Add(fhMCEtaPtOrigin) ;
1608 fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1609 fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1610 fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1611 outputContainer->Add(fhMCPi0ProdVertex) ;
1613 fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1614 fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1615 fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1616 outputContainer->Add(fhMCEtaProdVertex) ;
1618 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1620 fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #pi^{0} / E reco",
1621 nptbins,ptmin,ptmax,200,0,2);
1622 fhMCPi0PtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1623 fhMCPi0PtGenRecoFraction->SetYTitle("#it{E}^{#pi^{0} mother} / #it{E}^{rec}");
1624 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1626 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),#it{p}_{T} versus E primary #eta / E reco",
1627 nptbins,ptmin,ptmax,200,0,2);
1628 fhMCEtaPtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1629 fhMCEtaPtGenRecoFraction->SetYTitle("#it{E}^{ #eta mother} / #it{E}^{rec}");
1630 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1632 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1633 fhMCPi0DecayPt->SetYTitle("#it{N}");
1634 fhMCPi0DecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1635 outputContainer->Add(fhMCPi0DecayPt) ;
1637 fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #gamma / #it{E} primary #pi^{0}",
1638 nptbins,ptmin,ptmax,100,0,1);
1639 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/#it{c})");
1640 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1641 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1643 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1644 fhMCEtaDecayPt->SetYTitle("#it{N}");
1645 fhMCEtaDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1646 outputContainer->Add(fhMCEtaDecayPt) ;
1648 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #gamma / E primary #eta",
1649 nptbins,ptmin,ptmax,100,0,1);
1650 fhMCEtaDecayPtFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1651 fhMCEtaDecayPtFraction->SetYTitle("#it{E}^{gen} / #it{E}^{gen-mother}");
1652 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1654 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1655 fhMCOtherDecayPt->SetYTitle("#it{N}");
1656 fhMCOtherDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1657 outputContainer->Add(fhMCOtherDecayPt) ;
1661 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1662 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1665 fhAnglePairMCPi0 = new TH2F
1667 "Angle between decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1668 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1669 fhAnglePairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1670 outputContainer->Add(fhAnglePairMCPi0) ;
1672 if (fAnaType!= kSSCalo)
1674 fhAnglePairMCEta = new TH2F
1676 "Angle between decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1677 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1678 fhAnglePairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1679 outputContainer->Add(fhAnglePairMCEta) ;
1681 fhMassPairMCPi0 = new TH2F
1683 "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1684 fhMassPairMCPi0->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1685 fhMassPairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1686 outputContainer->Add(fhMassPairMCPi0) ;
1688 fhMassPairMCEta = new TH2F
1690 "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1691 fhMassPairMCEta->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1692 fhMassPairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1693 outputContainer->Add(fhMassPairMCEta) ;
1696 for(Int_t i = 0; i < 6; i++)
1700 (Form("hE_MC%s",pname[i].Data()),
1701 Form("Identified as #pi^{0} (#eta), cluster from %s",
1703 nptbins,ptmin,ptmax);
1704 fhMCE[i]->SetYTitle("#it{N}");
1705 fhMCE[i]->SetXTitle("#it{E} (GeV)");
1706 outputContainer->Add(fhMCE[i]) ;
1708 fhMCPt[i] = new TH1F
1709 (Form("hPt_MC%s",pname[i].Data()),
1710 Form("Identified as #pi^{0} (#eta), cluster from %s",
1712 nptbins,ptmin,ptmax);
1713 fhMCPt[i]->SetYTitle("#it{N}");
1714 fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1715 outputContainer->Add(fhMCPt[i]) ;
1717 fhMCPtCentrality[i] = new TH2F
1718 (Form("hPtCentrality_MC%s",pname[i].Data()),
1719 Form("Identified as #pi^{0} (#eta), cluster from %s",
1721 nptbins,ptmin,ptmax, 100,0,100);
1722 fhMCPtCentrality[i]->SetYTitle("centrality");
1723 fhMCPtCentrality[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1724 outputContainer->Add(fhMCPtCentrality[i]) ;
1726 if(fAnaType == kSSCalo)
1728 fhMCNLocMaxPt[i] = new TH2F
1729 (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1730 Form("cluster from %s, #it{p}_{T} of cluster vs NLM, accepted",ptype[i].Data()),
1731 nptbins,ptmin,ptmax,20,0,20);
1732 fhMCNLocMaxPt[i] ->SetYTitle("#it{NLM}");
1733 fhMCNLocMaxPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1734 outputContainer->Add(fhMCNLocMaxPt[i]) ;
1736 fhMCNLocMaxPtReject[i] = new TH2F
1737 (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1738 Form("cluster from %s, #it{p}_{T} of cluster vs NLM, rejected",ptype[i].Data()),
1739 nptbins,ptmin,ptmax,20,0,20);
1740 fhMCNLocMaxPtReject[i] ->SetYTitle("#it{NLM}");
1741 fhMCNLocMaxPtReject[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1742 outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1744 fhMCEReject[i] = new TH1F
1745 (Form("hEReject_MC%s",pname[i].Data()),
1746 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1748 nptbins,ptmin,ptmax);
1749 fhMCEReject[i]->SetYTitle("#it{N}");
1750 fhMCEReject[i]->SetXTitle("#it{E} (GeV)");
1751 outputContainer->Add(fhMCEReject[i]) ;
1753 fhMCPtReject[i] = new TH1F
1754 (Form("hPtReject_MC%s",pname[i].Data()),
1755 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1757 nptbins,ptmin,ptmax);
1758 fhMCPtReject[i]->SetYTitle("#it{N}");
1759 fhMCPtReject[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1760 outputContainer->Add(fhMCPtReject[i]) ;
1763 fhMCPtPhi[i] = new TH2F
1764 (Form("hPtPhi_MC%s",pname[i].Data()),
1765 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1766 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1767 fhMCPtPhi[i]->SetYTitle("#phi");
1768 fhMCPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1769 outputContainer->Add(fhMCPtPhi[i]) ;
1771 fhMCPtEta[i] = new TH2F
1772 (Form("hPtEta_MC%s",pname[i].Data()),
1773 Form("Identified as #pi^{0} (#eta), cluster from %s",
1774 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1775 fhMCPtEta[i]->SetYTitle("#eta");
1776 fhMCPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1777 outputContainer->Add(fhMCPtEta[i]) ;
1779 fhMCMassPt[i] = new TH2F
1780 (Form("hMassPt_MC%s",pname[i].Data()),
1781 Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1782 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1783 fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1784 fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1785 outputContainer->Add(fhMCMassPt[i]) ;
1787 fhMCSelectedMassPt[i] = new TH2F
1788 (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1789 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1790 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1791 fhMCSelectedMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1792 fhMCSelectedMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1793 outputContainer->Add(fhMCSelectedMassPt[i]) ;
1795 if(fAnaType == kSSCalo)
1797 fhMCMassPtNoOverlap[i] = new TH2F
1798 (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1799 Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1800 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1801 fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1802 fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1803 outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1805 fhMCSelectedMassPtNoOverlap[i] = new TH2F
1806 (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1807 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1808 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1809 fhMCSelectedMassPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1810 fhMCSelectedMassPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1811 outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1814 if( fFillSelectClHisto )
1816 fhMCPtLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1817 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1818 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1819 fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1820 fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1821 outputContainer->Add(fhMCPtLambda0[i]) ;
1823 fhMCPtLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1824 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1825 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1826 fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1827 fhMCPtLambda1[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1828 outputContainer->Add(fhMCPtLambda1[i]) ;
1830 fhMCPtDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1831 Form("Selected pair, cluster from %s : #it{p}_{T} vs dispersion^{2}",ptype[i].Data()),
1832 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1833 fhMCPtDispersion[i]->SetYTitle("#it{D}^{2}");
1834 fhMCPtDispersion[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1835 outputContainer->Add(fhMCPtDispersion[i]) ;
1837 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0)
1839 fhMCPtLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1840 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1841 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1842 fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1843 fhMCPtLambda0NoTRD[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1844 outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1846 if(!fFillOnlySimpleSSHisto)
1848 fhMCPtDispEta[i] = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1849 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1850 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1851 fhMCPtDispEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1852 fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1853 outputContainer->Add(fhMCPtDispEta[i]);
1855 fhMCPtDispPhi[i] = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1856 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1857 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1858 fhMCPtDispPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1859 fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1860 outputContainer->Add(fhMCPtDispPhi[i]);
1862 fhMCPtSumEtaPhi[i] = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1863 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",ptype[i].Data()),
1864 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1865 fhMCPtSumEtaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1866 fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1867 outputContainer->Add(fhMCPtSumEtaPhi[i]);
1869 fhMCPtDispEtaPhiDiff[i] = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1870 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",ptype[i].Data()),
1871 nptbins,ptmin,ptmax,200,-10,10);
1872 fhMCPtDispEtaPhiDiff[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1873 fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1874 outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1876 fhMCPtSphericity[i] = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1877 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()),
1878 nptbins,ptmin,ptmax, 200,-1,1);
1879 fhMCPtSphericity[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1880 fhMCPtSphericity[i]->SetYTitle("#it{s} = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1881 outputContainer->Add(fhMCPtSphericity[i]);
1883 for(Int_t ie = 0; ie < 7; ie++)
1885 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1886 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]),
1887 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1888 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1889 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1890 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1892 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1893 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]),
1894 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1895 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1896 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1897 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1899 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1900 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]),
1901 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1902 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1903 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1904 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1910 fhMCPtLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1911 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1912 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1913 fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1914 fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("#it{E} (GeV)");
1915 outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
1917 fhMCPtFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1918 Form("Selected pair, cluster from %s : #it{p}_{T} vs Max cell fraction of energy",ptype[i].Data()),
1919 nptbins,ptmin,ptmax,100,0,1);
1920 fhMCPtFracMaxCell[i]->SetYTitle("#it{Fraction}");
1921 fhMCPtFracMaxCell[i]->SetXTitle("#it{E} (GeV)");
1922 outputContainer->Add(fhMCPtFracMaxCell[i]) ;
1925 } // shower shape histo
1930 if(fAnaType==kSSCalo)
1932 fhAsymmetry = new TH2F ("hAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
1933 nptbins,ptmin,ptmax, 200, -1,1);
1934 fhAsymmetry->SetXTitle("#it{E} (GeV)");
1935 fhAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
1936 outputContainer->Add(fhAsymmetry);
1938 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
1939 nptbins,ptmin,ptmax, 200, -1,1);
1940 fhSelectedAsymmetry->SetXTitle("#it{E} (GeV)");
1941 fhSelectedAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
1942 outputContainer->Add(fhSelectedAsymmetry);
1945 ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1946 fhSplitE->SetYTitle("counts");
1947 fhSplitE->SetXTitle("#it{E} (GeV)");
1948 outputContainer->Add(fhSplitE) ;
1950 fhSplitPt = new TH1F
1951 ("hSplitPt","Selected #pi^{0} (#eta) pairs #it{p}_{T} sum of split sub-clusters",nptbins,ptmin,ptmax);
1952 fhSplitPt->SetYTitle("counts");
1953 fhSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1954 outputContainer->Add(fhSplitPt) ;
1957 fhSplitPtPhi = new TH2F
1958 ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1959 fhSplitPtPhi->SetYTitle("#phi (rad)");
1960 fhSplitPtPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1961 outputContainer->Add(fhSplitPtPhi) ;
1963 fhSplitPtEta = new TH2F
1964 ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1965 fhSplitPtEta->SetYTitle("#eta");
1966 fhSplitPtEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1967 outputContainer->Add(fhSplitPtEta) ;
1970 fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1971 nptbins,ptmin,ptmax,20,0,20);
1972 fhNLocMaxSplitPt ->SetYTitle("#it{NLM}");
1973 fhNLocMaxSplitPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1974 outputContainer->Add(fhNLocMaxSplitPt) ;
1977 fhMassSplitPt = new TH2F
1978 ("hMassSplitPt","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
1979 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1980 fhMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1981 fhMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1982 outputContainer->Add(fhMassSplitPt) ;
1984 fhSelectedMassSplitPt = new TH2F
1985 ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
1986 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1987 fhSelectedMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1988 fhSelectedMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1989 outputContainer->Add(fhSelectedMassSplitPt) ;
1993 fhMassSplitPtNoOverlap = new TH2F
1994 ("hMassSplitPtNoOverlap","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
1995 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1996 fhMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1997 fhMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1998 outputContainer->Add(fhMassSplitPtNoOverlap) ;
2000 fhSelectedMassSplitPtNoOverlap = new TH2F
2001 ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2002 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2003 fhSelectedMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2004 fhSelectedMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2005 outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
2008 fhMCPi0PtRecoPtPrim = new TH2F
2009 ("hMCPi0PtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2010 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2011 fhMCPi0PtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2012 fhMCPi0PtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2013 outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
2015 fhMCPi0PtRecoPtPrimNoOverlap = new TH2F
2016 ("hMCPi0PtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2017 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2018 fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2019 fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2020 outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
2022 fhMCPi0SelectedPtRecoPtPrim = new TH2F
2023 ("hMCPi0SelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2024 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2025 fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2026 fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2027 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
2029 fhMCPi0SelectedPtRecoPtPrimNoOverlap = new TH2F
2030 ("hMCPi0SelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2031 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2032 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2033 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2034 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
2037 fhMCPi0SplitPtRecoPtPrim = new TH2F
2038 ("hMCPi0SplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2039 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2040 fhMCPi0SplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2041 fhMCPi0SplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2042 outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
2044 fhMCPi0SplitPtRecoPtPrimNoOverlap = new TH2F
2045 ("hMCPi0SplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2046 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2047 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2048 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2049 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
2051 fhMCPi0SelectedSplitPtRecoPtPrim = new TH2F
2052 ("hMCPi0SelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2053 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2054 fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2055 fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2056 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
2058 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2059 ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2060 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2061 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2062 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2063 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
2065 fhMCEtaPtRecoPtPrim = new TH2F
2066 ("hMCEtaPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2067 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2068 fhMCEtaPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2069 fhMCEtaPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2070 outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
2072 fhMCEtaPtRecoPtPrimNoOverlap = new TH2F
2073 ("hMCEtaPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2074 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2075 fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2076 fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2077 outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2079 fhMCEtaSelectedPtRecoPtPrim = new TH2F
2080 ("hMCEtaSelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2081 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2082 fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2083 fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2084 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2086 fhMCEtaSelectedPtRecoPtPrimNoOverlap = new TH2F
2087 ("hMCEtaSelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2088 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2089 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2090 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2091 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2094 fhMCEtaSplitPtRecoPtPrim = new TH2F
2095 ("hMCEtaSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2096 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2097 fhMCEtaSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2098 fhMCEtaSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2099 outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2101 fhMCEtaSplitPtRecoPtPrimNoOverlap = new TH2F
2102 ("hMCEtaSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2103 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2104 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2105 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2106 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2108 fhMCEtaSelectedSplitPtRecoPtPrim = new TH2F
2109 ("hMCEtaSelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2110 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2111 fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2112 fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2113 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2115 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2116 ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2117 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2118 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2119 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2120 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2123 for(Int_t inlm = 0; inlm < 3; inlm++)
2125 fhMCPi0PtRecoPtPrimLocMax[inlm] = new TH2F
2126 (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2127 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2128 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2129 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2130 outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
2132 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2133 (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2134 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2135 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2136 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2137 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
2139 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] = new TH2F
2140 (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2141 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2142 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2143 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2144 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
2146 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2147 (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2148 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2149 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2150 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2151 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2153 fhMCEtaPtRecoPtPrimLocMax[inlm] = new TH2F
2154 (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2155 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2156 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2157 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2158 outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
2160 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2161 (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2162 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2163 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2164 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2165 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
2167 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2168 (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2169 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2170 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2171 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2172 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
2174 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2175 (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2176 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2177 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2178 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2179 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2183 for(Int_t i = 0; i< 6; i++)
2185 fhMCPtAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2186 Form("cluster from %s : #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",ptype[i].Data()),
2187 nptbins,ptmin,ptmax, 200,-1,1);
2188 fhMCPtAsymmetry[i]->SetXTitle("#it{E} (GeV)");
2189 fhMCPtAsymmetry[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2190 outputContainer->Add(fhMCPtAsymmetry[i]);
2192 fhMCSplitE[i] = new TH1F
2193 (Form("hSplitE_MC%s",pname[i].Data()),
2194 Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2195 nptbins,ptmin,ptmax);
2196 fhMCSplitE[i]->SetYTitle("counts");
2197 fhMCSplitE[i]->SetXTitle("#it{E} (GeV)");
2198 outputContainer->Add(fhMCSplitE[i]) ;
2200 fhMCSplitPt[i] = new TH1F
2201 (Form("hSplitPt_MC%s",pname[i].Data()),
2202 Form("cluster from %s, #it{p}_{T} sum of split sub-clusters",ptype[i].Data()),
2203 nptbins,ptmin,ptmax);
2204 fhMCSplitPt[i]->SetYTitle("counts");
2205 fhMCSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2206 outputContainer->Add(fhMCSplitPt[i]) ;
2209 fhMCSplitPtPhi[i] = new TH2F
2210 (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2211 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2212 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2213 fhMCSplitPtPhi[i]->SetYTitle("#phi");
2214 fhMCSplitPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2215 outputContainer->Add(fhMCSplitPtPhi[i]) ;
2217 fhMCSplitPtEta[i] = new TH2F
2218 (Form("hSplitPtEta_MC%s",pname[i].Data()),
2219 Form("Identified as #pi^{0} (#eta), cluster from %s",
2220 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2221 fhMCSplitPtEta[i]->SetYTitle("#eta");
2222 fhMCSplitPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2223 outputContainer->Add(fhMCSplitPtEta[i]) ;
2226 fhMCNLocMaxSplitPt[i] = new TH2F
2227 (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2228 Form("cluster from %s, #it{p}_{T} sum of split sub-clusters, for NLM",ptype[i].Data()),
2229 nptbins,ptmin,ptmax,20,0,20);
2230 fhMCNLocMaxSplitPt[i] ->SetYTitle("#it{NLM}");
2231 fhMCNLocMaxSplitPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2232 outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2234 fhMCMassSplitPt[i] = new TH2F
2235 (Form("hMassSplitPt_MC%s",pname[i].Data()),
2236 Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2237 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2238 fhMCMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2239 fhMCMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2240 outputContainer->Add(fhMCMassSplitPt[i]) ;
2242 fhMCSelectedMassSplitPt[i] = new TH2F
2243 (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2244 Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2245 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2246 fhMCSelectedMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2247 fhMCSelectedMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2248 outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2250 fhMCMassSplitPtNoOverlap[i] = new TH2F
2251 (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2252 Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2253 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2254 fhMCMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2255 fhMCMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2256 outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2258 fhMCSelectedMassSplitPtNoOverlap[i] = new TH2F
2259 (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2260 Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2261 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2262 fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2263 fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2264 outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2269 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2273 for(Int_t i = 0; i< 3; i++)
2275 fhPtAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2276 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ), %s",nlm[i].Data()),
2277 nptbins,ptmin,ptmax,200, -1,1);
2278 fhPtAsymmetryLocMax[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2279 fhPtAsymmetryLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2280 outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2283 for(Int_t ie = 0; ie< 7; ie++)
2286 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2287 Form("#lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2288 ssbins,ssmin,ssmax , 200,-1,1);
2289 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2290 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2291 outputContainer->Add(fhAsymmetryLambda0[ie]);
2293 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2294 Form("#sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2295 ssbins,ssmin,ssmax , 200,-1,1);
2296 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2297 fhAsymmetryDispEta[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2298 outputContainer->Add(fhAsymmetryDispEta[ie]);
2300 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2301 Form("#sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2302 ssbins,ssmin,ssmax , 200,-1,1);
2303 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2304 fhAsymmetryDispPhi[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2305 outputContainer->Add(fhAsymmetryDispPhi[ie]);
2311 for(Int_t i = 0; i< 6; i++)
2313 for(Int_t ie = 0; ie < 7; ie++)
2315 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2316 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2317 ssbins,ssmin,ssmax , 200,-1,1);
2318 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2319 fhMCAsymmetryLambda0[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2320 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2322 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2323 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2324 ssbins,ssmin,ssmax , 200,-1,1);
2325 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2326 fhMCAsymmetryDispEta[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2327 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2329 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2330 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2331 ssbins,ssmin,ssmax , 200,-1,1);
2332 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2333 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2334 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2340 if(fFillPileUpHistograms)
2343 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2345 for(Int_t i = 0 ; i < 7 ; i++)
2347 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2348 Form("Selected #pi^{0} (#eta) #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2349 fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2350 outputContainer->Add(fhPtPileUp[i]);
2352 fhPtCellTimePileUp[i] = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2353 Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2354 nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2355 fhPtCellTimePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2356 fhPtCellTimePileUp[i]->SetYTitle("#it{t}_{cell} (ns)");
2357 outputContainer->Add(fhPtCellTimePileUp[i]);
2359 fhPtTimeDiffPileUp[i] = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2360 Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2361 nptbins,ptmin,ptmax,400,-200,200);
2362 fhPtTimeDiffPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2363 fhPtTimeDiffPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
2364 outputContainer->Add(fhPtTimeDiffPileUp[i]);
2368 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","#it{t} of cluster vs #it{E} of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2369 fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2370 fhTimePtNoCut->SetYTitle("#it{t} (ns)");
2371 outputContainer->Add(fhTimePtNoCut);
2373 fhTimePtSPD = new TH2F ("hTimePt_SPD","#it{t} of cluster vs #it{E} of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2374 fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2375 fhTimePtSPD->SetYTitle("#it{t} (ns)");
2376 outputContainer->Add(fhTimePtSPD);
2378 fhTimePtSPDMulti = new TH2F ("hTimePt_SPDMulti","time of cluster vs #it{E} of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2379 fhTimePtSPDMulti->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2380 fhTimePtSPDMulti->SetYTitle("#it{t} (ns)");
2381 outputContainer->Add(fhTimePtSPDMulti);
2383 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","#it{t} of cluster vs #it{N} pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2384 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2385 fhTimeNPileUpVertSPD->SetXTitle("#it{t} (ns)");
2386 outputContainer->Add(fhTimeNPileUpVertSPD);
2388 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","#it{t} of cluster vs #it{N} pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2389 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2390 fhTimeNPileUpVertTrack->SetXTitle("#it{t} (ns)");
2391 outputContainer->Add(fhTimeNPileUpVertTrack);
2393 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","#it{t} of cluster vs #it{N} constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2394 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2395 fhTimeNPileUpVertContributors->SetXTitle("#it{t} (ns)");
2396 outputContainer->Add(fhTimeNPileUpVertContributors);
2398 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
2399 fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{Z} (cm) ");
2400 fhTimePileUpMainVertexZDistance->SetXTitle("#it{t} (ns)");
2401 outputContainer->Add(fhTimePileUpMainVertexZDistance);
2403 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
2404 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{Z} (cm) ");
2405 fhTimePileUpMainVertexZDiamond->SetXTitle("#it{t} (ns)");
2406 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2408 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","#it{p}_{T} of cluster vs #it{N} pile-up SPD vertex",
2409 nptbins,ptmin,ptmax,20,0,20);
2410 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2411 fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2412 outputContainer->Add(fhPtNPileUpSPDVtx);
2414 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","#it{p}_{T} of cluster vs #it{N} pile-up Tracks vertex",
2415 nptbins,ptmin,ptmax, 20,0,20 );
2416 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2417 fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2418 outputContainer->Add(fhPtNPileUpTrkVtx);
2420 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","#it{p}_{T} of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2421 nptbins,ptmin,ptmax,20,0,20);
2422 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2423 fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2424 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2426 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","#it{p}_{T} of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2427 nptbins,ptmin,ptmax, 20,0,20 );
2428 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2429 fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2430 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2432 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","#it{p}_{T} of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2433 nptbins,ptmin,ptmax,20,0,20);
2434 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2435 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2436 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2438 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","#it{p}_{T} of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2439 nptbins,ptmin,ptmax, 20,0,20 );
2440 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2441 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2442 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2446 //Keep neutral meson selection histograms if requiered
2447 //Setting done in AliNeutralMesonSelection
2449 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2451 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2453 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2454 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2459 return outputContainer ;
2463 //_____________________________________________
2464 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2467 // Assign mc index depending on MC bit set
2469 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
2473 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2477 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2478 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2480 return kmcConversion ;
2481 }//conversion photon
2482 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2485 }//photon no conversion
2486 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2488 return kmcElectron ;
2497 //__________________________________________________________________
2498 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2499 AliAODPWG4Particle * photon2,
2500 Int_t & label, Int_t & tag)
2502 // Check the labels of pare in case mother was same pi0 or eta
2503 // Set the new AOD accordingly
2505 Int_t label1 = photon1->GetLabel();
2506 Int_t label2 = photon2->GetLabel();
2508 if(label1 < 0 || label2 < 0 ) return ;
2510 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2511 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2512 Int_t tag1 = photon1->GetTag();
2513 Int_t tag2 = photon2->GetTag();
2515 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2516 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2517 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
2518 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2519 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
2523 //Check if pi0/eta mother is the same
2524 if(GetReader()->ReadStack())
2528 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2529 label1 = mother1->GetFirstMother();
2530 //mother1 = GetMCStack()->Particle(label1);//pi0
2534 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2535 label2 = mother2->GetFirstMother();
2536 //mother2 = GetMCStack()->Particle(label2);//pi0
2539 else if(GetReader()->ReadAODMCParticles())
2540 {//&& (input > -1)){
2543 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2544 label1 = mother1->GetMother();
2545 //mother1 = GetMCStack()->Particle(label1);//pi0
2549 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2550 label2 = mother2->GetMother();
2551 //mother2 = GetMCStack()->Particle(label2);//pi0
2555 //printf("mother1 %d, mother2 %d\n",label1,label2);
2556 if( label1 == label2 && label1>=0 )
2561 TLorentzVector mom1 = *(photon1->Momentum());
2562 TLorentzVector mom2 = *(photon2->Momentum());
2564 Double_t angle = mom2.Angle(mom1.Vect());
2565 Double_t mass = (mom1+mom2).M();
2566 Double_t epair = (mom1+mom2).E();
2568 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2570 fhMassPairMCPi0 ->Fill(epair,mass);
2571 fhAnglePairMCPi0->Fill(epair,angle);
2572 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2576 fhMassPairMCEta ->Fill(epair,mass);
2577 fhAnglePairMCEta->Fill(epair,angle);
2578 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2582 } // both from eta or pi0 decay
2586 //____________________________________________________________________________
2587 void AliAnaPi0EbE::Init()
2591 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2592 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2595 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2596 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2602 //____________________________________________________________________________
2603 void AliAnaPi0EbE::InitParameters()
2605 //Initialize the parameters of the analysis.
2606 AddToHistogramsName("AnaPi0EbE_");
2608 fInputAODGammaConvName = "PhotonsCTS" ;
2609 fAnaType = kIMCalo ;
2610 fCalorimeter = "EMCAL" ;
2615 fNLMECutMin[0] = 10.;
2616 fNLMECutMin[1] = 6. ;
2617 fNLMECutMin[2] = 6. ;
2620 //__________________________________________________________________
2621 void AliAnaPi0EbE::MakeAnalysisFillAOD()
2623 //Do analysis and fill aods
2628 MakeInvMassInCalorimeter();
2632 MakeShowerShapeIdentification();
2636 MakeInvMassInCalorimeterAndCTS();
2642 //____________________________________________
2643 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
2645 //Do analysis and fill aods
2646 //Search for the photon decay in calorimeters
2647 //Read photon list from AOD, produced in class AliAnaPhoton
2648 //Check if 2 photons have the mass of the pi0.
2650 TLorentzVector mom1;
2651 TLorentzVector mom2;
2652 TLorentzVector mom ;
2657 if(!GetInputAODBranch()){
2658 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2662 //Get shower shape information of clusters
2663 TObjArray *clusters = 0;
2664 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2665 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2667 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2668 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2670 //Vertex cut in case of mixed events
2671 Int_t evtIndex1 = 0 ;
2673 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2674 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
2675 mom1 = *(photon1->Momentum());
2677 //Get original cluster, to recover some information
2679 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2682 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2686 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2688 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2690 Int_t evtIndex2 = 0 ;
2692 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2694 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2697 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
2699 mom2 = *(photon2->Momentum());
2701 //Get original cluster, to recover some information
2703 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2704 // start new loop from iclus1+1 to gain some time
2708 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2712 Float_t e1 = photon1->E();
2713 Float_t e2 = photon2->E();
2715 //Select clusters with good time window difference
2716 Float_t tof1 = cluster1->GetTOF()*1e9;;
2717 Float_t tof2 = cluster2->GetTOF()*1e9;;
2718 Double_t t12diff = tof1-tof2;
2719 fhEPairDiffTime->Fill(e1+e2, t12diff);
2720 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2722 //Play with the MC stack if available
2723 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2725 // Check the invariant mass for different selection on the local maxima
2726 // Name of AOD method TO BE FIXED
2727 Int_t nMaxima1 = photon1->GetFiducialArea();
2728 Int_t nMaxima2 = photon2->GetFiducialArea();
2732 Double_t mass = mom.M();
2733 Double_t epair = mom.E();
2735 if(nMaxima1==nMaxima2)
2737 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2738 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2739 else fhMassPairLocMax[2]->Fill(epair,mass);
2741 else if(nMaxima1==1 || nMaxima2==1)
2743 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2744 else fhMassPairLocMax[4]->Fill(epair,mass);
2747 fhMassPairLocMax[5]->Fill(epair,mass);
2749 // combinations with SS axis cut and NLM cut
2750 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2751 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2752 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2753 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2755 //Skip events with too few or too many NLM
2756 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2758 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2761 fhMass->Fill(epair,mass);
2763 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2764 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2767 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",
2768 mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2770 //Fill some histograms about shower shape
2771 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2773 FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2774 FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2777 // Tag both photons as decay, !!careful!! since in case of SideBand analysis, also they will be tagged.
2778 photon1->SetTagged(kTRUE);
2779 photon2->SetTagged(kTRUE);
2781 fhPtDecay->Fill(photon1->Pt());
2782 fhEDecay ->Fill(photon1->E() );
2784 fhPtDecay->Fill(photon2->Pt());
2785 fhEDecay ->Fill(photon2->E() );
2787 //Mass of selected pairs
2788 fhSelectedMass->Fill(epair,mass);
2790 // Fill histograms to undertand pile-up before other cuts applied
2791 // Remember to relax time cuts in the reader
2792 FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2794 //Create AOD for analysis
2796 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2798 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2799 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2802 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2805 pi0.SetDetector(photon1->GetDetector());
2808 pi0.SetLabel(label);
2811 //Set the indeces of the original caloclusters
2812 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2813 //pi0.SetInputFileIndex(input);
2815 AddAODParticle(pi0);
2823 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2827 //__________________________________________________
2828 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2830 //Do analysis and fill aods
2831 //Search for the photon decay in calorimeters
2832 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2833 //Check if 2 photons have the mass of the pi0.
2835 TLorentzVector mom1;
2836 TLorentzVector mom2;
2837 TLorentzVector mom ;
2842 // Check calorimeter input
2843 if(!GetInputAODBranch())
2845 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2849 // Get the array with conversion photons
2850 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2851 if(!inputAODGammaConv)
2853 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2855 if(!inputAODGammaConv)
2857 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2863 //Get shower shape information of clusters
2864 TObjArray *clusters = 0;
2865 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2866 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2868 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
2869 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2870 if(nCTS<=0 || nCalo <=0)
2872 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2877 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2879 // Do the loop, first calo, second CTS
2880 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
2882 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2883 mom1 = *(photon1->Momentum());
2885 //Get original cluster, to recover some information
2887 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2889 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++)
2891 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2894 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2895 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2897 mom2 = *(photon2->Momentum());
2901 Double_t mass = mom.M();
2902 Double_t epair = mom.E();
2904 Int_t nMaxima = photon1->GetFiducialArea();
2905 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2906 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2907 else fhMassPairLocMax[2]->Fill(epair,mass);
2909 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2910 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2912 //Play with the MC stack if available
2915 Int_t label2 = photon2->GetLabel();
2916 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2918 HasPairSameMCMother(photon1, photon2, label, tag) ;
2921 //Mass of selected pairs
2922 fhMass->Fill(epair,mass);
2924 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2925 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2927 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",
2928 mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2930 //Fill some histograms about shower shape
2931 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2933 FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
2936 // Tag both photons as decay, !!careful!! since in case of SideBand analysis, also they will be tagged.
2937 photon1->SetTagged(kTRUE);
2938 photon2->SetTagged(kTRUE);
2940 fhPtDecay->Fill(photon1->Pt());
2941 fhEDecay ->Fill(photon1->E() );
2943 //Mass of selected pairs
2944 fhSelectedMass->Fill(epair,mass);
2946 // Fill histograms to undertand pile-up before other cuts applied
2947 // Remember to relax time cuts in the reader
2948 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2950 //Create AOD for analysis
2952 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2954 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2955 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2958 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2961 pi0.SetDetector(photon1->GetDetector());
2964 pi0.SetLabel(label);
2967 //Set the indeces of the original tracks or caloclusters
2968 pi0.SetCaloLabel (photon1->GetCaloLabel(0) , -1);
2969 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2970 //pi0.SetInputFileIndex(input);
2972 AddAODParticle(pi0);
2979 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2984 //_________________________________________________
2985 void AliAnaPi0EbE::MakeShowerShapeIdentification()
2987 //Search for pi0 in fCalorimeter with shower shape analysis
2989 TObjArray * pl = 0x0;
2990 AliVCaloCells * cells = 0x0;
2991 //Select the Calorimeter of the photon
2992 if (fCalorimeter == "PHOS" )
2994 pl = GetPHOSClusters();
2995 cells = GetPHOSCells();
2997 else if (fCalorimeter == "EMCAL")
2999 pl = GetEMCALClusters();
3000 cells = GetEMCALCells();
3005 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3009 TLorentzVector mom ;
3010 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3012 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3014 Int_t evtIndex = 0 ;
3015 if (GetMixedEvent())
3017 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3020 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
3022 //Get Momentum vector,
3023 Double_t vertex[]={0,0,0};
3024 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3026 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
3027 }//Assume that come from vertex in straight line
3030 calo->GetMomentum(mom,vertex) ;
3033 //If too small or big pt, skip it
3034 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
3036 //Check acceptance selection
3037 if(IsFiducialCutOn())
3039 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
3040 if(! in ) continue ;
3044 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());
3046 //Play with the MC stack if available
3047 //Check origin of the candidates
3051 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
3052 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
3053 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
3056 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3058 //Check Distance to Bad channel, set bit.
3059 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3060 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
3061 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3062 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3066 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
3068 //If too low number of cells, skip it
3069 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3071 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3076 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
3077 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
3079 //.......................................
3080 // TOF cut, BE CAREFUL WITH THIS CUT
3081 Double_t tof = calo->GetTOF()*1e9;
3082 if(tof < fTimeCutMin || tof > fTimeCutMax)
3084 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3089 //PID selection or bit setting
3091 Double_t mass = 0, angle = 0;
3092 Int_t absId1 =-1, absId2 =-1;
3093 Float_t distbad1 =-1, distbad2 =-1;
3094 Bool_t fidcut1 = 0, fidcut2 = 0;
3095 TLorentzVector l1, l2;
3097 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3098 GetVertex(evtIndex),nMaxima,
3099 mass,angle,l1,l2,absId1,absId2,
3100 distbad1,distbad2,fidcut1,fidcut2) ;
3103 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
3106 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3107 if( (fCheckSplitDistToBad) &&
3108 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3111 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
3112 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
3114 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3118 //Skip events with too few or too many NLM
3119 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3121 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3126 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
3128 //Skip matched clusters with tracks
3129 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3131 FillRejectedClusterHistograms(mom,tag,nMaxima);
3135 Float_t e1 = l1.Energy();
3136 Float_t e2 = l2.Energy();
3137 TLorentzVector l12 = l1+l2;
3138 Float_t ptSplit = l12.Pt();
3139 Float_t eSplit = e1+e2;
3142 Int_t noverlaps = 0;
3146 mcIndex = GetMCIndex(tag);
3149 Int_t mcLabel = calo->GetLabel();
3151 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3153 Int_t mesonLabel = -1;
3155 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3157 if(mcIndex == kmcPi0)
3159 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3160 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3164 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3165 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3169 const UInt_t nlabels = calo->GetNLabels();
3170 Int_t overpdg[nlabels];
3171 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3174 //mass of all clusters
3175 fhMass ->Fill(mom.E() ,mass);
3176 fhMassPt ->Fill(mom.Pt(),mass);
3177 fhMassSplitPt->Fill(ptSplit ,mass);
3179 Int_t indexMax = -1;
3180 if (nMaxima==1) indexMax = 0 ;
3181 else if(nMaxima==2) indexMax = 1 ;
3183 fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3187 fhMCMassPt[mcIndex] ->Fill(mom.Pt(),mass);
3188 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3191 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3192 fhMCPi0SplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3193 fhMCPi0PtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3194 fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3197 else if(mcIndex==kmcEta)
3199 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3200 fhMCEtaSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3201 fhMCEtaPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3202 fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3209 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3210 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3212 else if(mcIndex==kmcEta)
3214 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3215 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3218 fhMassNoOverlap ->Fill(mom.E() ,mass);
3219 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3220 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3222 fhMCMassPtNoOverlap[mcIndex] ->Fill(mom.Pt(),mass);
3223 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3227 // Asymmetry of all clusters
3230 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3231 fhAsymmetry->Fill(mom.E(),asy);
3235 fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3238 // If cluster does not pass pid, not pi0/eta, skip it.
3239 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3241 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3242 FillRejectedClusterHistograms(mom,tag,nMaxima);
3246 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3248 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3249 FillRejectedClusterHistograms(mom,tag,nMaxima);
3254 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3255 mom.Pt(), idPartType);
3257 //Mass and asymmetry of selected pairs
3258 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3259 fhSelectedMass ->Fill(mom.E() ,mass);
3260 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3261 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3262 fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3264 Int_t nSM = GetModuleNumber(calo);
3265 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
3267 fhSelectedMassPtLocMaxSM [indexMax][nSM]->Fill(mom.Pt(),mass);
3268 fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),calo->GetM02());
3275 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3276 fhMCPi0SelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3277 fhMCPi0SelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3278 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3280 else if(mcIndex==kmcEta)
3282 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3283 fhMCEtaSelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3284 fhMCEtaSelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3285 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3290 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3291 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3292 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3296 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3297 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3299 else if(mcIndex==kmcEta)
3301 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3302 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3307 fhSplitE ->Fill( eSplit);
3308 fhSplitPt ->Fill(ptSplit);
3309 Float_t phi = mom.Phi();
3310 if(phi<0) phi+=TMath::TwoPi();
3311 fhSplitPtPhi ->Fill(ptSplit,phi);
3312 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3313 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3315 //Check split-clusters with good time window difference
3316 Double_t tof1 = cells->GetCellTime(absId1);
3317 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3320 Double_t tof2 = cells->GetCellTime(absId2);
3321 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3324 Double_t t12diff = tof1-tof2;
3325 fhEPairDiffTime->Fill(e1+e2, t12diff);
3329 fhMCSplitE [mcIndex]->Fill( eSplit);
3330 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3331 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3332 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3333 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3334 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3336 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3337 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3338 fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3342 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3343 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3347 //-----------------------
3348 //Create AOD for analysis
3350 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3351 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3352 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3354 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3355 aodpi0.SetLabel(calo->GetLabel());
3357 //Set the indeces of the original caloclusters
3358 aodpi0.SetCaloLabel(calo->GetID(),-1);
3359 aodpi0.SetDetector(fCalorimeter);
3361 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3362 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3363 else aodpi0.SetDistToBad(0) ;
3365 // Check if cluster is pi0 via cluster splitting
3366 aodpi0.SetIdentifiedParticleType(idPartType);
3368 // Add number of local maxima to AOD, method name in AOD to be FIXED
3369 aodpi0.SetFiducialArea(nMaxima);
3373 //Fill some histograms about shower shape
3374 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3376 FillSelectedClusterHistograms(calo, aodpi0.Pt(), nMaxima, tag, asy);
3379 // Fill histograms to undertand pile-up before other cuts applied
3380 // Remember to relax time cuts in the reader
3381 Double_t tofcluster = calo->GetTOF()*1e9;
3382 Double_t tofclusterUS = TMath::Abs(tofcluster);
3384 FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
3386 Int_t id = GetReader()->GetTriggerClusterId();
3387 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
3389 Float_t phicluster = aodpi0.Phi();
3390 if(phicluster < 0) phicluster+=TMath::TwoPi();
3394 if (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
3395 else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
3396 else fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
3399 Int_t bc = GetReader()->GetTriggerClusterBC();
3400 if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
3402 if(GetReader()->IsTriggerMatched())
3404 if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
3405 fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
3406 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
3410 if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
3411 fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
3415 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(calo->E(), tofcluster);
3416 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
3417 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(calo->E(), tofcluster);
3421 else if(TMath::Abs(bc) >= 6)
3422 Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
3425 //Add AOD with pi0 object to aod branch
3426 AddAODParticle(aodpi0);
3430 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3433 //______________________________________________
3434 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3436 //Do analysis and fill histograms
3438 if(!GetOutputAODBranch())
3440 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3443 //Loop on stored AOD pi0
3444 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3445 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3447 Float_t cen = GetEventCentrality();
3448 Float_t ep = GetEventPlaneAngle();
3450 for(Int_t iaod = 0; iaod < naod ; iaod++)
3452 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3453 Int_t pdg = pi0->GetIdentifiedParticleType();
3455 if( ( pdg != AliCaloPID::kPi0 && pdg != AliCaloPID::kEta ) ) continue;
3457 //Fill pi0 histograms
3458 Float_t ener = pi0->E();
3459 Float_t pt = pi0->Pt();
3460 Float_t phi = pi0->Phi();
3461 if(phi < 0) phi+=TMath::TwoPi();
3462 Float_t eta = pi0->Eta();
3467 fhPtEta ->Fill(pt ,eta);
3468 fhPtPhi ->Fill(pt ,phi);
3469 fhEtaPhi ->Fill(eta ,phi);
3471 fhPtCentrality ->Fill(pt,cen) ;
3472 fhPtEventPlane ->Fill(pt,ep ) ;
3476 Int_t tag = pi0->GetTag();
3477 Int_t label = pi0->GetLabel();
3478 Int_t mcIndex = GetMCIndex(tag);
3480 fhMCE [mcIndex] ->Fill(ener);
3481 fhMCPt [mcIndex] ->Fill(pt);
3482 fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3483 fhMCPtEta[mcIndex] ->Fill(pt,eta);
3485 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3487 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
3489 Float_t efracMC = 0;
3490 Int_t momlabel = -1;
3493 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3496 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3498 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3499 if(grandmom.E() > 0 && ok)
3501 efracMC = grandmom.E()/ener;
3502 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3505 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3507 fhMCPi0DecayPt->Fill(pt);
3508 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3509 if(grandmom.E() > 0 && ok)
3511 efracMC = mom.E()/grandmom.E();
3512 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3515 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3517 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3518 if(grandmom.E() > 0 && ok)
3520 efracMC = grandmom.E()/ener;
3521 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3524 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3526 fhMCEtaDecayPt->Fill(pt);
3527 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3528 if(grandmom.E() > 0 && ok)
3530 efracMC = mom.E()/grandmom.E();
3531 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3534 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3536 fhMCOtherDecayPt->Fill(pt);
3541 if( mcIndex==kmcPi0 || mcIndex==kmcEta )
3544 Int_t momindex = -1;
3546 Int_t momstatus = -1;
3548 if(GetReader()->ReadStack())
3550 TParticle* ancestor = GetMCStack()->Particle(label);
3551 momindex = ancestor->GetFirstMother();
3552 if(momindex < 0) return;
3553 TParticle* mother = GetMCStack()->Particle(momindex);
3554 mompdg = TMath::Abs(mother->GetPdgCode());
3555 momstatus = mother->GetStatusCode();
3556 prodR = mother->R();
3560 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3561 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
3562 momindex = ancestor->GetMother();
3563 if(momindex < 0) return;
3564 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3565 mompdg = TMath::Abs(mother->GetPdgCode());
3566 momstatus = mother->GetStatus();
3567 prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3570 if( mcIndex==kmcPi0 )
3572 fhMCPi0ProdVertex->Fill(pt,prodR);
3574 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
3575 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
3576 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
3577 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
3578 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
3579 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
3580 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
3581 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3582 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
3583 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
3584 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
3586 else if (mcIndex==kmcEta )
3588 fhMCEtaProdVertex->Fill(pt,prodR);
3590 if (momstatus == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
3591 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
3592 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
3593 else if(mompdg == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
3594 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
3595 else if(mompdg == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
3596 else if(mompdg == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
3597 else if(mompdg >= 310 && mompdg <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3598 else if(mompdg == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
3599 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
3600 else fhMCEtaPtOrigin->Fill(pt,7.5);//other?
3604 }//Histograms with MC
3610 //__________________________________________________________________
3611 void AliAnaPi0EbE::Print(const Option_t * opt) const
3613 //Print some relevant parameters set for the analysis
3617 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3618 AliAnaCaloTrackCorrBaseClass::Print("");
3619 printf("Analysis Type = %d \n", fAnaType) ;
3620 if(fAnaType == kSSCalo)
3622 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3623 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3624 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3625 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);