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",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2769 //Fill some histograms about shower shape
2770 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2772 FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2773 FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2776 // Tag both photons as decay
2777 photon1->SetTagged(kTRUE);
2778 photon2->SetTagged(kTRUE);
2780 fhPtDecay->Fill(photon1->Pt());
2781 fhEDecay ->Fill(photon1->E() );
2783 fhPtDecay->Fill(photon2->Pt());
2784 fhEDecay ->Fill(photon2->E() );
2786 //Mass of selected pairs
2787 fhSelectedMass->Fill(epair,mass);
2789 // Fill histograms to undertand pile-up before other cuts applied
2790 // Remember to relax time cuts in the reader
2791 FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2793 //Create AOD for analysis
2795 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2797 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2798 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2801 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2804 pi0.SetDetector(photon1->GetDetector());
2807 pi0.SetLabel(label);
2810 //Set the indeces of the original caloclusters
2811 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2812 //pi0.SetInputFileIndex(input);
2814 AddAODParticle(pi0);
2822 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2826 //__________________________________________________
2827 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2829 //Do analysis and fill aods
2830 //Search for the photon decay in calorimeters
2831 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2832 //Check if 2 photons have the mass of the pi0.
2834 TLorentzVector mom1;
2835 TLorentzVector mom2;
2836 TLorentzVector mom ;
2841 // Check calorimeter input
2842 if(!GetInputAODBranch()){
2843 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2847 // Get the array with conversion photons
2848 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2849 if(!inputAODGammaConv) {
2851 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2853 if(!inputAODGammaConv) {
2854 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2860 //Get shower shape information of clusters
2861 TObjArray *clusters = 0;
2862 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2863 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2865 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
2866 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2867 if(nCTS<=0 || nCalo <=0)
2869 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2874 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2876 // Do the loop, first calo, second CTS
2877 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2878 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2879 mom1 = *(photon1->Momentum());
2881 //Get original cluster, to recover some information
2883 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2885 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2886 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2888 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2889 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2891 mom2 = *(photon2->Momentum());
2893 Double_t mass = (mom1+mom2).M();
2894 Double_t epair = (mom1+mom2).E();
2896 Int_t nMaxima = photon1->GetFiducialArea();
2897 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2898 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2899 else fhMassPairLocMax[2]->Fill(epair,mass);
2901 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2902 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2904 //Play with the MC stack if available
2907 Int_t label2 = photon2->GetLabel();
2908 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2910 HasPairSameMCMother(photon1, photon2, label, tag) ;
2913 //Mass of selected pairs
2914 fhMass->Fill(epair,(mom1+mom2).M());
2916 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2917 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2919 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2921 //Fill some histograms about shower shape
2922 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2924 FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
2927 // Tag both photons as decay
2928 photon1->SetTagged(kTRUE);
2929 photon2->SetTagged(kTRUE);
2931 fhPtDecay->Fill(photon1->Pt());
2932 fhEDecay ->Fill(photon1->E() );
2934 //Create AOD for analysis
2938 //Mass of selected pairs
2939 fhSelectedMass->Fill(epair,mom.M());
2941 // Fill histograms to undertand pile-up before other cuts applied
2942 // Remember to relax time cuts in the reader
2943 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2945 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2947 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2948 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2951 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2954 pi0.SetDetector(photon1->GetDetector());
2957 pi0.SetLabel(label);
2960 //Set the indeces of the original tracks or caloclusters
2961 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2962 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2963 //pi0.SetInputFileIndex(input);
2965 AddAODParticle(pi0);
2972 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2977 //_________________________________________________
2978 void AliAnaPi0EbE::MakeShowerShapeIdentification()
2980 //Search for pi0 in fCalorimeter with shower shape analysis
2982 TObjArray * pl = 0x0;
2983 AliVCaloCells * cells = 0x0;
2984 //Select the Calorimeter of the photon
2985 if (fCalorimeter == "PHOS" )
2987 pl = GetPHOSClusters();
2988 cells = GetPHOSCells();
2990 else if (fCalorimeter == "EMCAL")
2992 pl = GetEMCALClusters();
2993 cells = GetEMCALCells();
2998 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3002 TLorentzVector mom ;
3003 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3005 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3007 Int_t evtIndex = 0 ;
3008 if (GetMixedEvent())
3010 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3013 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
3015 //Get Momentum vector,
3016 Double_t vertex[]={0,0,0};
3017 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3019 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
3020 }//Assume that come from vertex in straight line
3023 calo->GetMomentum(mom,vertex) ;
3026 //If too small or big pt, skip it
3027 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
3029 //Check acceptance selection
3030 if(IsFiducialCutOn())
3032 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
3033 if(! in ) continue ;
3037 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());
3039 //Play with the MC stack if available
3040 //Check origin of the candidates
3044 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
3045 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
3046 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
3049 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3051 //Check Distance to Bad channel, set bit.
3052 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3053 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
3054 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3055 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3059 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
3061 //If too low number of cells, skip it
3062 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3064 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3069 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
3070 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
3072 //.......................................
3073 // TOF cut, BE CAREFUL WITH THIS CUT
3074 Double_t tof = calo->GetTOF()*1e9;
3075 if(tof < fTimeCutMin || tof > fTimeCutMax)
3077 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3082 //PID selection or bit setting
3084 Double_t mass = 0, angle = 0;
3085 Int_t absId1 =-1, absId2 =-1;
3086 Float_t distbad1 =-1, distbad2 =-1;
3087 Bool_t fidcut1 = 0, fidcut2 = 0;
3088 TLorentzVector l1, l2;
3090 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3091 GetVertex(evtIndex),nMaxima,
3092 mass,angle,l1,l2,absId1,absId2,
3093 distbad1,distbad2,fidcut1,fidcut2) ;
3096 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
3099 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3100 if( (fCheckSplitDistToBad) &&
3101 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3104 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
3105 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
3107 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3111 //Skip events with too few or too many NLM
3112 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3114 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3119 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
3121 //Skip matched clusters with tracks
3122 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3124 FillRejectedClusterHistograms(mom,tag,nMaxima);
3128 Float_t e1 = l1.Energy();
3129 Float_t e2 = l2.Energy();
3130 TLorentzVector l12 = l1+l2;
3131 Float_t ptSplit = l12.Pt();
3132 Float_t eSplit = e1+e2;
3135 Int_t noverlaps = 0;
3139 mcIndex = GetMCIndex(tag);
3142 Int_t mcLabel = calo->GetLabel();
3144 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3146 Int_t mesonLabel = -1;
3148 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3150 if(mcIndex == kmcPi0)
3152 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3153 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3157 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3158 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3162 const UInt_t nlabels = calo->GetNLabels();
3163 Int_t overpdg[nlabels];
3164 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3167 //mass of all clusters
3168 fhMass ->Fill(mom.E() ,mass);
3169 fhMassPt ->Fill(mom.Pt(),mass);
3170 fhMassSplitPt->Fill(ptSplit ,mass);
3172 Int_t indexMax = -1;
3173 if (nMaxima==1) indexMax = 0 ;
3174 else if(nMaxima==2) indexMax = 1 ;
3176 fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3180 fhMCMassPt[mcIndex] ->Fill(mom.Pt(),mass);
3181 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3184 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3185 fhMCPi0SplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3186 fhMCPi0PtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3187 fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3190 else if(mcIndex==kmcEta)
3192 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3193 fhMCEtaSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3194 fhMCEtaPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3195 fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3202 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3203 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3205 else if(mcIndex==kmcEta)
3207 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3208 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3211 fhMassNoOverlap ->Fill(mom.E() ,mass);
3212 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3213 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3215 fhMCMassPtNoOverlap[mcIndex] ->Fill(mom.Pt(),mass);
3216 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3220 // Asymmetry of all clusters
3223 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3224 fhAsymmetry->Fill(mom.E(),asy);
3228 fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3231 // If cluster does not pass pid, not pi0/eta, skip it.
3232 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3234 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3235 FillRejectedClusterHistograms(mom,tag,nMaxima);
3239 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3241 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3242 FillRejectedClusterHistograms(mom,tag,nMaxima);
3247 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3248 mom.Pt(), idPartType);
3250 //Mass and asymmetry of selected pairs
3251 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3252 fhSelectedMass ->Fill(mom.E() ,mass);
3253 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3254 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3255 fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3257 Int_t nSM = GetModuleNumber(calo);
3258 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
3260 fhSelectedMassPtLocMaxSM [indexMax][nSM]->Fill(mom.Pt(),mass);
3261 fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),calo->GetM02());
3268 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3269 fhMCPi0SelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3270 fhMCPi0SelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3271 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3273 else if(mcIndex==kmcEta)
3275 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3276 fhMCEtaSelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3277 fhMCEtaSelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3278 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3283 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3284 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3285 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3289 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3290 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3292 else if(mcIndex==kmcEta)
3294 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3295 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3300 fhSplitE ->Fill( eSplit);
3301 fhSplitPt ->Fill(ptSplit);
3302 Float_t phi = mom.Phi();
3303 if(phi<0) phi+=TMath::TwoPi();
3304 fhSplitPtPhi ->Fill(ptSplit,phi);
3305 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3306 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3308 //Check split-clusters with good time window difference
3309 Double_t tof1 = cells->GetCellTime(absId1);
3310 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3313 Double_t tof2 = cells->GetCellTime(absId2);
3314 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3317 Double_t t12diff = tof1-tof2;
3318 fhEPairDiffTime->Fill(e1+e2, t12diff);
3322 fhMCSplitE [mcIndex]->Fill( eSplit);
3323 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3324 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3325 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3326 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3327 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3329 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3330 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3331 fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3335 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3336 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3340 //-----------------------
3341 //Create AOD for analysis
3343 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3344 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3345 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3347 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3348 aodpi0.SetLabel(calo->GetLabel());
3350 //Set the indeces of the original caloclusters
3351 aodpi0.SetCaloLabel(calo->GetID(),-1);
3352 aodpi0.SetDetector(fCalorimeter);
3354 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3355 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3356 else aodpi0.SetDistToBad(0) ;
3358 // Check if cluster is pi0 via cluster splitting
3359 aodpi0.SetIdentifiedParticleType(idPartType);
3361 // Add number of local maxima to AOD, method name in AOD to be FIXED
3362 aodpi0.SetFiducialArea(nMaxima);
3366 //Fill some histograms about shower shape
3367 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3369 FillSelectedClusterHistograms(calo, aodpi0.Pt(), nMaxima, tag, asy);
3372 // Fill histograms to undertand pile-up before other cuts applied
3373 // Remember to relax time cuts in the reader
3374 Double_t tofcluster = calo->GetTOF()*1e9;
3375 Double_t tofclusterUS = TMath::Abs(tofcluster);
3377 FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
3379 Int_t id = GetReader()->GetTriggerClusterId();
3380 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
3382 Float_t phicluster = aodpi0.Phi();
3383 if(phicluster < 0) phicluster+=TMath::TwoPi();
3387 if (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
3388 else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
3389 else fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
3392 Int_t bc = GetReader()->GetTriggerClusterBC();
3393 if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
3395 if(GetReader()->IsTriggerMatched())
3397 if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
3398 fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
3399 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
3403 if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
3404 fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
3408 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(calo->E(), tofcluster);
3409 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
3410 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(calo->E(), tofcluster);
3414 else if(TMath::Abs(bc) >= 6)
3415 Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
3418 //Add AOD with pi0 object to aod branch
3419 AddAODParticle(aodpi0);
3423 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3426 //______________________________________________
3427 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3429 //Do analysis and fill histograms
3431 if(!GetOutputAODBranch())
3433 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3435 //Loop on stored AOD pi0
3436 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3437 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3439 Float_t cen = GetEventCentrality();
3440 Float_t ep = GetEventPlaneAngle();
3442 for(Int_t iaod = 0; iaod < naod ; iaod++)
3444 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3445 Int_t pdg = pi0->GetIdentifiedParticleType();
3447 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
3449 //Fill pi0 histograms
3450 Float_t ener = pi0->E();
3451 Float_t pt = pi0->Pt();
3452 Float_t phi = pi0->Phi();
3453 if(phi < 0) phi+=TMath::TwoPi();
3454 Float_t eta = pi0->Eta();
3459 fhPtEta ->Fill(pt ,eta);
3460 fhPtPhi ->Fill(pt ,phi);
3461 fhEtaPhi ->Fill(eta ,phi);
3463 fhPtCentrality ->Fill(pt,cen) ;
3464 fhPtEventPlane ->Fill(pt,ep ) ;
3468 Int_t tag = pi0->GetTag();
3469 Int_t label = pi0->GetLabel();
3470 Int_t mcIndex = GetMCIndex(tag);
3472 fhMCE [mcIndex] ->Fill(ener);
3473 fhMCPt [mcIndex] ->Fill(pt);
3474 fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3475 fhMCPtEta[mcIndex] ->Fill(pt,eta);
3477 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3479 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
3481 Float_t efracMC = 0;
3482 Int_t momlabel = -1;
3485 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3488 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3490 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3491 if(grandmom.E() > 0 && ok)
3493 efracMC = grandmom.E()/ener;
3494 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3497 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3499 fhMCPi0DecayPt->Fill(pt);
3500 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3501 if(grandmom.E() > 0 && ok)
3503 efracMC = mom.E()/grandmom.E();
3504 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3507 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3509 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3510 if(grandmom.E() > 0 && ok)
3512 efracMC = grandmom.E()/ener;
3513 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3516 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3518 fhMCEtaDecayPt->Fill(pt);
3519 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3520 if(grandmom.E() > 0 && ok)
3522 efracMC = mom.E()/grandmom.E();
3523 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3526 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3528 fhMCOtherDecayPt->Fill(pt);
3533 if( mcIndex==kmcPi0 || mcIndex==kmcEta )
3536 Int_t momindex = -1;
3538 Int_t momstatus = -1;
3540 if(GetReader()->ReadStack())
3542 TParticle* ancestor = GetMCStack()->Particle(label);
3543 momindex = ancestor->GetFirstMother();
3544 if(momindex < 0) return;
3545 TParticle* mother = GetMCStack()->Particle(momindex);
3546 mompdg = TMath::Abs(mother->GetPdgCode());
3547 momstatus = mother->GetStatusCode();
3548 prodR = mother->R();
3552 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3553 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
3554 momindex = ancestor->GetMother();
3555 if(momindex < 0) return;
3556 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3557 mompdg = TMath::Abs(mother->GetPdgCode());
3558 momstatus = mother->GetStatus();
3559 prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3562 if( mcIndex==kmcPi0 )
3564 fhMCPi0ProdVertex->Fill(pt,prodR);
3566 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
3567 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
3568 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
3569 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
3570 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
3571 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
3572 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
3573 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3574 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
3575 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
3576 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
3578 else if (mcIndex==kmcEta )
3580 fhMCEtaProdVertex->Fill(pt,prodR);
3582 if (momstatus == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
3583 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
3584 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
3585 else if(mompdg == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
3586 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
3587 else if(mompdg == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
3588 else if(mompdg == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
3589 else if(mompdg >= 310 && mompdg <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3590 else if(mompdg == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
3591 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
3592 else fhMCEtaPtOrigin->Fill(pt,7.5);//other?
3596 }//Histograms with MC
3602 //__________________________________________________________________
3603 void AliAnaPi0EbE::Print(const Option_t * opt) const
3605 //Print some relevant parameters set for the analysis
3609 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3610 AliAnaCaloTrackCorrBaseClass::Print("");
3611 printf("Analysis Type = %d \n", fAnaType) ;
3612 if(fAnaType == kSSCalo)
3614 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3615 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3616 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3617 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);