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),
109 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
110 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
111 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
112 fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
113 fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
114 fhTrackMatchedMCParticlePt(0),
115 fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
116 fhdEdx(0), fhEOverP(0), fhEOverPNoTRD(0),
117 // Number of local maxima in cluster
118 fhNLocMaxPt(0), fhNLocMaxPtReject(0),
120 fhTimePtNoCut(0), fhTimePtSPD(0), fhTimePtSPDMulti(0),
121 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
122 fhTimeNPileUpVertContributors(0),
123 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
124 fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
125 fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
126 fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
130 for(Int_t i = 0; i < 6; i++)
136 fhMCPtCentrality [i] = 0;
140 fhMCSplitPtPhi [i] = 0;
141 fhMCSplitPtEta [i] = 0;
143 fhMCNLocMaxPt [i] = 0;
144 fhMCNLocMaxSplitPt [i] = 0;
145 fhMCNLocMaxPtReject[i] = 0;
147 fhMCPtLambda0 [i] = 0;
148 fhMCPtLambda0NoTRD [i] = 0;
149 fhMCPtLambda0FracMaxCellCut[i]= 0;
150 fhMCPtFracMaxCell [i] = 0;
151 fhMCPtLambda1 [i] = 0;
152 fhMCPtDispersion [i] = 0;
154 fhMCPtDispEta [i] = 0;
155 fhMCPtDispPhi [i] = 0;
156 fhMCPtSumEtaPhi [i] = 0;
157 fhMCPtDispEtaPhiDiff[i] = 0;
158 fhMCPtSphericity [i] = 0;
159 fhMCPtAsymmetry [i] = 0;
162 fhMCMassSplitPt [i]=0;
163 fhMCSelectedMassPt [i]=0;
164 fhMCSelectedMassSplitPt[i]=0;
166 fhMCMassPtNoOverlap [i]=0;
167 fhMCMassSplitPtNoOverlap [i]=0;
168 fhMCSelectedMassPtNoOverlap [i]=0;
169 fhMCSelectedMassSplitPtNoOverlap[i]=0;
171 for(Int_t j = 0; j < 7; j++)
173 fhMCLambda0DispEta [j][i] = 0;
174 fhMCLambda0DispPhi [j][i] = 0;
175 fhMCDispEtaDispPhi [j][i] = 0;
176 fhMCAsymmetryLambda0 [j][i] = 0;
177 fhMCAsymmetryDispEta [j][i] = 0;
178 fhMCAsymmetryDispPhi [j][i] = 0;
182 for(Int_t j = 0; j < 7; j++)
184 fhLambda0DispEta [j] = 0;
185 fhLambda0DispPhi [j] = 0;
186 fhDispEtaDispPhi [j] = 0;
187 fhAsymmetryLambda0 [j] = 0;
188 fhAsymmetryDispEta [j] = 0;
189 fhAsymmetryDispPhi [j] = 0;
194 for(Int_t i = 0; i < 3; i++)
196 fhPtLambda0LocMax [i] = 0;
197 fhPtLambda1LocMax [i] = 0;
198 fhPtDispersionLocMax [i] = 0;
199 fhPtDispEtaLocMax [i] = 0;
200 fhPtDispPhiLocMax [i] = 0;
201 fhPtSumEtaPhiLocMax [i] = 0;
202 fhPtDispEtaPhiDiffLocMax[i] = 0;
203 fhPtSphericityLocMax [i] = 0;
204 fhPtAsymmetryLocMax [i] = 0;
205 fhMassPtLocMax [i] = 0;
206 fhSelectedMassPtLocMax [i] = 0;
207 for(Int_t ipart = 0; ipart<6; ipart++)
209 fhMCPtLambda0LocMax [ipart][i] = 0;
210 fhMCSelectedMassPtLocMax[ipart][i] = 0;
215 for(Int_t i =0; i < 14; i++){
216 fhLambda0ForW0[i] = 0;
217 //fhLambda1ForW0[i] = 0;
218 if(i<8)fhMassPairLocMax[i] = 0;
221 for(Int_t i = 0; i < 11; i++)
223 fhEtaPhiTriggerEMCALBC [i] = 0 ;
224 fhTimeTriggerEMCALBC [i] = 0 ;
225 fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
227 fhEtaPhiTriggerEMCALBCUM [i] = 0 ;
228 fhTimeTriggerEMCALBCUM [i] = 0 ;
232 //Initialize parameters
237 //___________________________________________________________________________________
238 void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
240 // Fill some histograms to understand pile-up
241 if(!fFillPileUpHistograms) return;
243 //printf("E %f, time %f\n",energy,time);
244 AliVEvent * event = GetReader()->GetInputEvent();
246 fhTimePtNoCut->Fill(pt,time);
247 if(GetReader()->IsPileUpFromSPD())
249 if(GetReader()->IsPileUpFromSPD()) { fhPtPileUp[0]->Fill(pt); fhTimePtSPD ->Fill(pt,time); }
250 if(GetReader()->IsPileUpFromEMCal()) fhPtPileUp[1]->Fill(pt);
251 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPileUp[2]->Fill(pt);
252 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPileUp[3]->Fill(pt);
253 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPileUp[4]->Fill(pt);
254 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPileUp[5]->Fill(pt);
255 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
257 if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
261 AliVCaloCells* cells = 0;
262 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
263 else cells = GetPHOSCells();
265 Float_t maxCellFraction = 0.;
266 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
268 Double_t tmax = cells->GetCellTime(absIdMax);
269 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
272 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
273 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
275 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
277 Int_t absId = calo->GetCellsAbsId()[ipos];
279 if( absId == absIdMax ) continue ;
281 Double_t timecell = cells->GetCellTime(absId);
282 Float_t amp = cells->GetCellAmplitude(absId);
283 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
284 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
287 Float_t diff = (tmax-timecell);
289 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
291 if(GetReader()->IsPileUpFromSPD())
293 fhPtCellTimePileUp[0]->Fill(pt, timecell);
294 fhPtTimeDiffPileUp[0]->Fill(pt, diff);
297 if(GetReader()->IsPileUpFromEMCal())
299 fhPtCellTimePileUp[1]->Fill(pt, timecell);
300 fhPtTimeDiffPileUp[1]->Fill(pt, diff);
303 if(GetReader()->IsPileUpFromSPDOrEMCal())
305 fhPtCellTimePileUp[2]->Fill(pt, timecell);
306 fhPtTimeDiffPileUp[2]->Fill(pt, diff);
309 if(GetReader()->IsPileUpFromSPDAndEMCal())
311 fhPtCellTimePileUp[3]->Fill(pt, timecell);
312 fhPtTimeDiffPileUp[3]->Fill(pt, diff);
315 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
317 fhPtCellTimePileUp[4]->Fill(pt, timecell);
318 fhPtTimeDiffPileUp[4]->Fill(pt, diff);
321 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
323 fhPtCellTimePileUp[5]->Fill(pt, timecell);
324 fhPtTimeDiffPileUp[5]->Fill(pt, diff);
327 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
329 fhPtCellTimePileUp[6]->Fill(pt, timecell);
330 fhPtTimeDiffPileUp[6]->Fill(pt, diff);
335 if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
337 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
338 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
340 // N pile up vertices
346 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
347 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
352 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
353 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
356 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
357 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
359 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
360 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
362 if(TMath::Abs(time) < 25)
364 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
365 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
368 if(time < 75 && time > -25)
370 fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
371 fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
374 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
375 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
378 Float_t z1 = -1, z2 = -1;
380 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
384 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
385 ncont=pv->GetNContributors();
386 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
388 diamZ = esdEv->GetDiamondZ();
392 AliAODVertex *pv=aodEv->GetVertex(iVert);
393 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
394 ncont=pv->GetNContributors();
395 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
397 diamZ = aodEv->GetDiamondZ();
400 Double_t distZ = TMath::Abs(z2-z1);
401 diamZ = TMath::Abs(z2-diamZ);
403 fhTimeNPileUpVertContributors ->Fill(time,ncont);
404 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
405 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
411 //______________________________________________________________________________________________
412 void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
414 // Fill histograms that do not pass the identification (SS case only)
416 Float_t ener = mom.E();
417 Float_t pt = mom.Pt();
418 Float_t phi = mom.Phi();
419 if(phi < 0) phi+=TMath::TwoPi();
420 Float_t eta = mom.Eta();
422 fhPtReject ->Fill(pt);
423 fhEReject ->Fill(ener);
425 fhPtEtaReject ->Fill(ener,eta);
426 fhPtPhiReject ->Fill(ener,phi);
427 fhEtaPhiReject ->Fill(eta,phi);
429 fhNLocMaxPtReject->Fill(pt,nMaxima);
433 Int_t mcIndex = GetMCIndex(mctag);
434 fhMCEReject [mcIndex] ->Fill(ener);
435 fhMCPtReject [mcIndex] ->Fill(pt);
436 fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
440 //___________________________________________________________________________________
441 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
442 Int_t tag, Float_t asy)
444 // Fill shower shape, timing and other histograms for selected clusters from decay
446 Float_t ener = cluster->E();
447 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
448 Float_t l0 = cluster->GetM02();
449 Float_t l1 = cluster->GetM20();
450 Int_t nSM = GetModuleNumber(cluster);
453 if (pt < 2 ) ptbin = 0;
454 else if (pt < 4 ) ptbin = 1;
455 else if (pt < 6 ) ptbin = 2;
456 else if (pt < 10) ptbin = 3;
457 else if (pt < 15) ptbin = 4;
458 else if (pt < 20) ptbin = 5;
462 if (nMaxima==1) indexMax = 0 ;
463 else if(nMaxima==2) indexMax = 1 ;
467 AliVCaloCells * cell = 0x0;
468 if(fCalorimeter == "PHOS")
469 cell = GetPHOSCells();
471 cell = GetEMCALCells();
473 Float_t maxCellFraction = 0;
474 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
475 fhPtFracMaxCell->Fill(pt,maxCellFraction);
477 FillWeightHistograms(cluster);
479 fhPtDispersion->Fill(pt, disp);
480 fhPtLambda0 ->Fill(pt, l0 );
481 fhPtLambda1 ->Fill(pt, l1 );
483 Float_t ll0 = 0., ll1 = 0.;
484 Float_t dispp= 0., dEta = 0., dPhi = 0.;
485 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
486 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
488 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
489 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
491 fhPtDispEta -> Fill(pt,dEta);
492 fhPtDispPhi -> Fill(pt,dPhi);
493 fhPtSumEta -> Fill(pt,sEta);
494 fhPtSumPhi -> Fill(pt,sPhi);
495 fhPtSumEtaPhi -> Fill(pt,sEtaPhi);
496 fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
497 if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
499 fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
500 fhLambda0DispEta[ptbin]->Fill(l0 ,dEta);
501 fhLambda0DispPhi[ptbin]->Fill(l0 ,dPhi);
503 if (fAnaType==kSSCalo)
505 // Asymmetry histograms
506 fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
507 fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
508 fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
512 fhNLocMaxPt->Fill(pt,nMaxima);
514 fhPtLambda0LocMax [indexMax]->Fill(pt,l0);
515 fhPtLambda1LocMax [indexMax]->Fill(pt,l1);
516 fhPtDispersionLocMax[indexMax]->Fill(pt,disp);
518 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
520 fhPtDispEtaLocMax [indexMax]-> Fill(pt,dEta);
521 fhPtDispPhiLocMax [indexMax]-> Fill(pt,dPhi);
522 fhPtSumEtaPhiLocMax [indexMax]-> Fill(pt,sEtaPhi);
523 fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
524 if(dEta+dPhi>0) fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
525 if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt ,asy);
529 if(fCalorimeter=="EMCAL" && nSM < 6)
531 fhPtLambda0NoTRD ->Fill(pt, l0 );
532 fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
535 if(maxCellFraction < 0.5)
536 fhPtLambda0FracMaxCellCut->Fill(pt, l0 );
538 fhPtTime ->Fill(pt, cluster->GetTOF()*1.e9);
539 fhPtNCells->Fill(pt, cluster->GetNCells());
541 // Fill Track matching control histograms
544 Float_t dZ = cluster->GetTrackDz();
545 Float_t dR = cluster->GetTrackDx();
547 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
549 dR = 2000., dZ = 2000.;
550 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
552 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
554 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
556 Bool_t positive = kFALSE;
557 if(track) positive = (track->Charge()>0);
559 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
561 fhTrackMatchedDEta->Fill(pt,dZ);
562 fhTrackMatchedDPhi->Fill(pt,dR);
563 if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
569 fhTrackMatchedDEtaPos->Fill(pt,dZ);
570 fhTrackMatchedDPhiPos->Fill(pt,dR);
571 if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
575 fhTrackMatchedDEtaNeg->Fill(pt,dZ);
576 fhTrackMatchedDPhiNeg->Fill(pt,dR);
577 if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
581 // Check dEdx and E/p of matched clusters
583 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
587 Float_t dEdx = track->GetTPCsignal();
588 fhdEdx->Fill(pt, dEdx);
590 Float_t eOverp = cluster->E()/track->P();
591 fhEOverP->Fill(pt, eOverp);
593 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
594 if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(pt, eOverp);
598 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
603 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
605 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
606 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 2.5 ;
607 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 0.5 ;
608 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 1.5 ;
614 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
615 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 6.5 ;
616 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 4.5 ;
617 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 5.5 ;
621 fhTrackMatchedMCParticlePt ->Fill(pt, mctag);
622 fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
623 fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
627 }// Track matching histograms
631 Int_t mcIndex = GetMCIndex(tag);
633 fhMCPtLambda0[mcIndex] ->Fill(pt, l0);
634 fhMCPtLambda1[mcIndex] ->Fill(pt, l1);
635 fhMCPtDispersion[mcIndex] ->Fill(pt, disp);
636 fhMCPtFracMaxCell[mcIndex]->Fill(pt,maxCellFraction);
638 fhMCPtLambda0LocMax [mcIndex][indexMax]->Fill(pt,l0);
640 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
641 if(fCalorimeter=="EMCAL" && nSM < 6)
642 fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0 );
644 if(maxCellFraction < 0.5)
645 fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0 );
647 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
649 fhMCPtDispEta [mcIndex]-> Fill(pt,dEta);
650 fhMCPtDispPhi [mcIndex]-> Fill(pt,dPhi);
651 fhMCPtSumEtaPhi [mcIndex]-> Fill(pt,sEtaPhi);
652 fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
653 if(dEta+dPhi>0)fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
655 if (fAnaType==kSSCalo)
657 fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
658 fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
659 fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
662 fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
663 fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0 ,dEta);
664 fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0 ,dPhi);
672 //________________________________________________________
673 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
675 // Calculate weights and fill histograms
677 if(!fFillWeightHistograms || GetMixedEvent()) return;
679 AliVCaloCells* cells = 0;
680 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
681 else cells = GetPHOSCells();
683 // First recalculate energy in case non linearity was applied
686 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
689 Int_t id = clus->GetCellsAbsId()[ipos];
691 //Recalibrate cell energy if needed
692 Float_t amp = cells->GetCellAmplitude(id);
693 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
704 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
708 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
709 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
711 //Get the ratio and log ratio to all cells in cluster
712 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
714 Int_t id = clus->GetCellsAbsId()[ipos];
716 //Recalibrate cell energy if needed
717 Float_t amp = cells->GetCellAmplitude(id);
718 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
720 fhECellClusterRatio ->Fill(energy,amp/energy);
721 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
724 //Recalculate shower shape for different W0
725 if(fCalorimeter=="EMCAL"){
727 Float_t l0org = clus->GetM02();
728 Float_t l1org = clus->GetM20();
729 Float_t dorg = clus->GetDispersion();
731 for(Int_t iw = 0; iw < 14; iw++)
733 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
734 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
736 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
737 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
741 // Set the original values back
744 clus->SetDispersion(dorg);
749 //__________________________________________
750 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
752 //Save parameters used for analysis
753 TString parList ; //this will be list of parameters used for this analysis.
754 const Int_t buffersize = 255;
755 char onePar[buffersize] ;
757 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
759 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
762 if(fAnaType == kSSCalo)
764 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
766 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
768 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
770 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
774 //Get parameters set in base class.
775 parList += GetBaseParametersList() ;
777 //Get parameters set in PID class.
778 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
780 return new TObjString(parList) ;
783 //_____________________________________________
784 TList * AliAnaPi0EbE::GetCreateOutputObjects()
786 // Create histograms to be saved in output file and
787 // store them in outputContainer
788 TList * outputContainer = new TList() ;
789 outputContainer->SetName("Pi0EbEHistos") ;
791 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
792 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
793 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
794 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
795 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
796 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
797 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
799 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
800 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
801 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
803 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
804 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
805 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
806 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
807 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
808 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
810 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
811 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
812 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
813 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
814 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
815 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
817 Int_t ntimptbins= GetHistogramRanges()->GetHistoTimeBins();
818 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
819 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
821 TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
822 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
823 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
824 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
826 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
827 fhPt->SetYTitle("N");
828 fhPt->SetXTitle("p_{T} (GeV/c)");
829 outputContainer->Add(fhPt) ;
831 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
833 fhE->SetXTitle("E (GeV)");
834 outputContainer->Add(fhE) ;
837 ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
838 fhPtPhi->SetYTitle("#phi (rad)");
839 fhPtPhi->SetXTitle("E (GeV)");
840 outputContainer->Add(fhPtPhi) ;
843 ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
844 fhPtEta->SetYTitle("#eta");
845 fhPtEta->SetXTitle("E (GeV)");
846 outputContainer->Add(fhPtEta) ;
849 ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
850 fhPtPhi->SetYTitle("#phi (rad)");
851 fhPtPhi->SetXTitle("p_{T} (GeV/c)");
852 outputContainer->Add(fhPtPhi) ;
855 ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
856 fhPtEta->SetYTitle("#eta");
857 fhPtEta->SetXTitle("p_{T} (GeV/c)");
858 outputContainer->Add(fhPtEta) ;
861 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
862 fhEtaPhi->SetYTitle("#phi (rad)");
863 fhEtaPhi->SetXTitle("#eta");
864 outputContainer->Add(fhEtaPhi) ;
866 if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
868 fhEtaPhiEMCALBC0 = new TH2F
869 ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
870 fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
871 fhEtaPhiEMCALBC0->SetXTitle("#eta");
872 outputContainer->Add(fhEtaPhiEMCALBC0) ;
874 fhEtaPhiEMCALBC1 = new TH2F
875 ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
876 fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
877 fhEtaPhiEMCALBC1->SetXTitle("#eta");
878 outputContainer->Add(fhEtaPhiEMCALBC1) ;
880 fhEtaPhiEMCALBCN = new TH2F
881 ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
882 fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
883 fhEtaPhiEMCALBCN->SetXTitle("#eta");
884 outputContainer->Add(fhEtaPhiEMCALBCN) ;
886 for(Int_t i = 0; i < 11; i++)
888 fhEtaPhiTriggerEMCALBC[i] = new TH2F
889 (Form("hEtaPhiTriggerEMCALBC%d",i-5),
890 Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
891 netabins,etamin,etamax,nphibins,phimin,phimax);
892 fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
893 fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
894 outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
896 fhTimeTriggerEMCALBC[i] = new TH2F
897 (Form("hTimeTriggerEMCALBC%d",i-5),
898 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
899 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
900 fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
901 fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
902 outputContainer->Add(fhTimeTriggerEMCALBC[i]);
904 fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
905 (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
906 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
907 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
908 fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
909 fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
910 outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
912 fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
913 (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
914 Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
915 netabins,etamin,etamax,nphibins,phimin,phimax);
916 fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
917 fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
918 outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
920 fhTimeTriggerEMCALBCUM[i] = new TH2F
921 (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
922 Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
923 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
924 fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
925 fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
926 outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
930 fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
931 "cluster time vs E of clusters, no match, rematch open time",
932 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
933 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("E (GeV)");
934 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("time (ns)");
935 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
938 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
939 "cluster time vs E of clusters, no match, rematch with neigbour parches",
940 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
941 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("E (GeV)");
942 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("time (ns)");
943 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
945 fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
946 "cluster time vs E of clusters, no match, rematch open time and neigbour",
947 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
948 fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("E (GeV)");
949 fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("time (ns)");
950 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
954 fhPtCentrality = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
955 fhPtCentrality->SetYTitle("centrality");
956 fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
957 outputContainer->Add(fhPtCentrality) ;
959 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
960 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
961 fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
962 outputContainer->Add(fhPtEventPlane) ;
964 if(fAnaType == kSSCalo)
966 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
967 fhPtReject->SetYTitle("N");
968 fhPtReject->SetXTitle("p_{T} (GeV/c)");
969 outputContainer->Add(fhPtReject) ;
971 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
972 fhEReject->SetYTitle("N");
973 fhEReject->SetXTitle("E (GeV)");
974 outputContainer->Add(fhEReject) ;
976 fhPtPhiReject = new TH2F
977 ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
978 fhPtPhiReject->SetYTitle("#phi (rad)");
979 fhPtPhiReject->SetXTitle("p_{T} (GeV/c)");
980 outputContainer->Add(fhPtPhiReject) ;
982 fhPtEtaReject = new TH2F
983 ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
984 fhPtEtaReject->SetYTitle("#eta");
985 fhPtEtaReject->SetXTitle("p_{T} (GeV/c)");
986 outputContainer->Add(fhPtEtaReject) ;
988 fhEtaPhiReject = new TH2F
989 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
990 fhEtaPhiReject->SetYTitle("#phi (rad)");
991 fhEtaPhiReject->SetXTitle("#eta");
992 outputContainer->Add(fhEtaPhiReject) ;
996 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
997 fhMass->SetYTitle("mass (GeV/c^{2})");
998 fhMass->SetXTitle("E (GeV)");
999 outputContainer->Add(fhMass) ;
1001 fhSelectedMass = new TH2F
1002 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1003 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
1004 fhSelectedMass->SetXTitle("E (GeV)");
1005 outputContainer->Add(fhSelectedMass) ;
1007 if(fAnaType == kSSCalo)
1011 ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1012 fhMassPt->SetYTitle("mass (GeV/c^{2})");
1013 fhMassPt->SetXTitle("p_{T} (GeV/c)");
1014 outputContainer->Add(fhMassPt) ;
1016 fhSelectedMassPt = new TH2F
1017 ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1018 fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
1019 fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
1020 outputContainer->Add(fhSelectedMassPt) ;
1022 for(Int_t inlm = 0; inlm < 3; inlm++)
1024 fhMassPtLocMax[inlm] = new TH2F
1025 (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs mass: p_{T} vs mass and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1026 fhMassPtLocMax[inlm]->SetYTitle("mass (GeV/c^{2})");
1027 fhMassPtLocMax[inlm]->SetXTitle("p_{T} (GeV/c)");
1028 outputContainer->Add(fhMassPtLocMax[inlm]) ;
1030 fhSelectedMassPtLocMax[inlm] = new TH2F
1031 (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1032 fhSelectedMassPtLocMax[inlm]->SetYTitle("mass (GeV/c^{2})");
1033 fhSelectedMassPtLocMax[inlm]->SetXTitle("p_{T} (GeV/c)");
1034 outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
1038 for(Int_t ipart = 0; ipart < 6; ipart++)
1040 fhMCSelectedMassPtLocMax[ipart][inlm] = new TH2F
1041 (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
1042 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, NLM=%s, %s",nlm[inlm].Data(),pname[ipart].Data()),
1043 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1044 fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("mass (GeV/c^{2})");
1045 fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("p_{T} (GeV/c)");
1046 outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
1053 fhMassNoOverlap = new TH2F
1054 ("hMassNoOverlap","all pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1055 fhMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1056 fhMassNoOverlap->SetXTitle("E (GeV)");
1057 outputContainer->Add(fhMassNoOverlap) ;
1059 fhSelectedMassNoOverlap = new TH2F
1060 ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1061 fhSelectedMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1062 fhSelectedMassNoOverlap->SetXTitle("E (GeV)");
1063 outputContainer->Add(fhSelectedMassNoOverlap) ;
1065 fhMassPtNoOverlap = new TH2F
1066 ("hMassPtNoOverlap","all pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1067 fhMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1068 fhMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1069 outputContainer->Add(fhMassPtNoOverlap) ;
1071 fhSelectedMassPtNoOverlap = new TH2F
1072 ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1073 fhSelectedMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1074 fhSelectedMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1075 outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1079 if(fAnaType != kSSCalo)
1081 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1082 fhPtDecay->SetYTitle("N");
1083 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
1084 outputContainer->Add(fhPtDecay) ;
1086 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1087 fhEDecay->SetYTitle("N");
1088 fhEDecay->SetXTitle("E (GeV)");
1089 outputContainer->Add(fhEDecay) ;
1094 if( fFillSelectClHisto )
1096 fhPtDispersion = new TH2F
1097 ("hPtDispersion","Selected #pi^{0} (#eta) pairs: p_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1098 fhPtDispersion->SetYTitle("D^{2}");
1099 fhPtDispersion->SetXTitle("p_{T} (GeV/c)");
1100 outputContainer->Add(fhPtDispersion) ;
1102 fhPtLambda0 = new TH2F
1103 ("hPtLambda0","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1104 fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1105 fhPtLambda0->SetXTitle("p_{T} (GeV/c)");
1106 outputContainer->Add(fhPtLambda0) ;
1108 fhPtLambda1 = new TH2F
1109 ("hPtLambda1","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1110 fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1111 fhPtLambda1->SetXTitle("p_{T} (GeV/c)");
1112 outputContainer->Add(fhPtLambda1) ;
1114 fhPtLambda0FracMaxCellCut = new TH2F
1115 ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1116 fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1117 fhPtLambda0FracMaxCellCut->SetXTitle("p_{T} (GeV/c)");
1118 outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1120 fhPtFracMaxCell = new TH2F
1121 ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1122 fhPtFracMaxCell->SetYTitle("Fraction");
1123 fhPtFracMaxCell->SetXTitle("p_{T} (GeV/c)");
1124 outputContainer->Add(fhPtFracMaxCell) ;
1126 if(fCalorimeter=="EMCAL")
1128 fhPtLambda0NoTRD = new TH2F
1129 ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1130 fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1131 fhPtLambda0NoTRD->SetXTitle("p_{T} (GeV/c)");
1132 outputContainer->Add(fhPtLambda0NoTRD) ;
1134 fhPtFracMaxCellNoTRD = new TH2F
1135 ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1136 fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1137 fhPtFracMaxCellNoTRD->SetXTitle("p_{T} (GeV/c)");
1138 outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1140 if(!fFillOnlySimpleSSHisto)
1142 fhPtDispEta = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1143 fhPtDispEta->SetXTitle("p_{T} (GeV/c)");
1144 fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1145 outputContainer->Add(fhPtDispEta);
1147 fhPtDispPhi = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1148 fhPtDispPhi->SetXTitle("p_{T} (GeV/c)");
1149 fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1150 outputContainer->Add(fhPtDispPhi);
1152 fhPtSumEta = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs p_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1153 fhPtSumEta->SetXTitle("p_{T} (GeV/c)");
1154 fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1155 outputContainer->Add(fhPtSumEta);
1157 fhPtSumPhi = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs p_{T}",
1158 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1159 fhPtSumPhi->SetXTitle("p_{T} (GeV/c)");
1160 fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1161 outputContainer->Add(fhPtSumPhi);
1163 fhPtSumEtaPhi = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs p_{T}",
1164 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1165 fhPtSumEtaPhi->SetXTitle("p_{T} (GeV/c)");
1166 fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1167 outputContainer->Add(fhPtSumEtaPhi);
1169 fhPtDispEtaPhiDiff = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs p_{T}",
1170 nptbins,ptmin,ptmax,200, -10,10);
1171 fhPtDispEtaPhiDiff->SetXTitle("p_{T} (GeV/c)");
1172 fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1173 outputContainer->Add(fhPtDispEtaPhiDiff);
1175 fhPtSphericity = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs p_{T} (GeV/c)",
1176 nptbins,ptmin,ptmax, 200, -1,1);
1177 fhPtSphericity->SetXTitle("p_{T} (GeV/c)");
1178 fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1179 outputContainer->Add(fhPtSphericity);
1181 for(Int_t i = 0; i < 7; i++)
1183 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]),
1184 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1185 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1186 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1187 outputContainer->Add(fhDispEtaDispPhi[i]);
1189 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]),
1190 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1191 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1192 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1193 outputContainer->Add(fhLambda0DispEta[i]);
1195 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]),
1196 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1197 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1198 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1199 outputContainer->Add(fhLambda0DispPhi[i]);
1205 fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1206 nptbins,ptmin,ptmax,20,0,20);
1207 fhNLocMaxPt ->SetYTitle("N maxima");
1208 fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
1209 outputContainer->Add(fhNLocMaxPt) ;
1211 if(fAnaType == kSSCalo)
1214 fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1215 nptbins,ptmin,ptmax,20,0,20);
1216 fhNLocMaxPtReject ->SetYTitle("N maxima");
1217 fhNLocMaxPtReject ->SetXTitle("p_{T} (GeV/c)");
1218 outputContainer->Add(fhNLocMaxPtReject) ;
1221 for (Int_t i = 0; i < 3; i++)
1223 fhPtLambda0LocMax[i] = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1224 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1225 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1226 fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1227 fhPtLambda0LocMax[i]->SetXTitle("p_{T} (GeV/c)");
1228 outputContainer->Add(fhPtLambda0LocMax[i]) ;
1232 for(Int_t ipart = 0; ipart < 6; ipart++)
1234 fhMCPtLambda0LocMax[ipart][i] = new TH2F
1235 (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1236 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),pname[ipart].Data()),
1237 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1238 fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1239 fhMCPtLambda0LocMax[ipart][i]->SetXTitle("p_{T} (GeV/c)");
1240 outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1244 fhPtLambda1LocMax[i] = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1245 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1246 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1247 fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1248 fhPtLambda1LocMax[i]->SetXTitle("p_{T} (GeV/c)");
1249 outputContainer->Add(fhPtLambda1LocMax[i]) ;
1251 fhPtDispersionLocMax[i] = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1252 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1253 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1254 fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1255 fhPtDispersionLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1256 outputContainer->Add(fhPtDispersionLocMax[i]) ;
1258 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1260 fhPtDispEtaLocMax[i] = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1261 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1262 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1263 fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1264 fhPtDispEtaLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1265 outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1267 fhPtDispPhiLocMax[i] = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1268 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1269 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1270 fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1271 fhPtDispPhiLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1272 outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1274 fhPtSumEtaPhiLocMax[i] = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1275 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1276 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1277 fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1278 fhPtSumEtaPhiLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1279 outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1281 fhPtDispEtaPhiDiffLocMax[i] = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1282 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1283 nptbins,ptmin,ptmax,200, -10,10);
1284 fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1285 fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1286 outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1288 fhPtSphericityLocMax[i] = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1289 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1290 nptbins,ptmin,ptmax,200, -1,1);
1291 fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1292 fhPtSphericityLocMax[i]->SetXTitle("p_{T} (GeV/c)");
1293 outputContainer->Add(fhPtSphericityLocMax[i]) ;
1298 fhPtNCells = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1299 fhPtNCells->SetXTitle("p_{T} (GeV/c)");
1300 fhPtNCells->SetYTitle("# of cells in cluster");
1301 outputContainer->Add(fhPtNCells);
1303 fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1304 fhPtTime->SetXTitle("p_{T} (GeV/c)");
1305 fhPtTime->SetYTitle("t (ns)");
1306 outputContainer->Add(fhPtTime);
1311 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1312 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
1313 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1314 outputContainer->Add(fhEPairDiffTime);
1316 if(fAnaType == kIMCalo)
1318 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1319 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1320 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1321 "2 Local Maxima paired with more than 2 Local Maxima",
1322 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1324 for (Int_t i = 0; i < 8 ; i++)
1327 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1329 fhMassPairLocMax[i] = new TH2F
1330 (Form("MassPairLocMax%s",combiName[i].Data()),
1331 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1332 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1333 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
1334 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
1335 outputContainer->Add(fhMassPairLocMax[i]) ;
1341 fhTrackMatchedDEta = new TH2F
1342 ("hTrackMatchedDEta",
1343 "d#eta of cluster-track vs cluster p_{T}",
1344 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1345 fhTrackMatchedDEta->SetYTitle("d#eta");
1346 fhTrackMatchedDEta->SetXTitle("p_{T} (GeV/c)");
1348 fhTrackMatchedDPhi = new TH2F
1349 ("hTrackMatchedDPhi",
1350 "d#phi of cluster-track vs cluster p_{T}",
1351 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1352 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1353 fhTrackMatchedDPhi->SetXTitle("p_{T} (GeV/c)");
1355 fhTrackMatchedDEtaDPhi = new TH2F
1356 ("hTrackMatchedDEtaDPhi",
1357 "d#eta vs d#phi of cluster-track",
1358 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1359 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1360 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1362 outputContainer->Add(fhTrackMatchedDEta) ;
1363 outputContainer->Add(fhTrackMatchedDPhi) ;
1364 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1366 fhTrackMatchedDEtaPos = new TH2F
1367 ("hTrackMatchedDEtaPos",
1368 "d#eta of cluster-track vs cluster p_{T}",
1369 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1370 fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1371 fhTrackMatchedDEtaPos->SetXTitle("p_{T} (GeV/c)");
1373 fhTrackMatchedDPhiPos = new TH2F
1374 ("hTrackMatchedDPhiPos",
1375 "d#phi of cluster-track vs cluster p_{T}",
1376 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1377 fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1378 fhTrackMatchedDPhiPos->SetXTitle("p_{T} (GeV/c)");
1380 fhTrackMatchedDEtaDPhiPos = new TH2F
1381 ("hTrackMatchedDEtaDPhiPos",
1382 "d#eta vs d#phi of cluster-track",
1383 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1384 fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1385 fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1387 outputContainer->Add(fhTrackMatchedDEtaPos) ;
1388 outputContainer->Add(fhTrackMatchedDPhiPos) ;
1389 outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1391 fhTrackMatchedDEtaNeg = new TH2F
1392 ("hTrackMatchedDEtaNeg",
1393 "d#eta of cluster-track vs cluster p_{T}",
1394 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1395 fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1396 fhTrackMatchedDEtaNeg->SetXTitle("p_{T} (GeV/c)");
1398 fhTrackMatchedDPhiNeg = new TH2F
1399 ("hTrackMatchedDPhiNeg",
1400 "d#phi of cluster-track vs cluster p_{T}",
1401 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1402 fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1403 fhTrackMatchedDPhiNeg->SetXTitle("p_{T} (GeV/c)");
1405 fhTrackMatchedDEtaDPhiNeg = new TH2F
1406 ("hTrackMatchedDEtaDPhiNeg",
1407 "d#eta vs d#phi of cluster-track",
1408 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1409 fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1410 fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1412 outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1413 outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1414 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1416 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster p_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1417 fhdEdx->SetXTitle("p_{T} (GeV/c)");
1418 fhdEdx->SetYTitle("<dE/dx>");
1419 outputContainer->Add(fhdEdx);
1421 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster p_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1422 fhEOverP->SetXTitle("p_{T} (GeV/c)");
1423 fhEOverP->SetYTitle("E/p");
1424 outputContainer->Add(fhEOverP);
1426 if(fCalorimeter=="EMCAL")
1428 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1429 fhEOverPNoTRD->SetXTitle("E (GeV)");
1430 fhEOverPNoTRD->SetYTitle("E/p");
1431 outputContainer->Add(fhEOverPNoTRD);
1434 if(IsDataMC() && fFillTMHisto)
1436 fhTrackMatchedMCParticlePt = new TH2F
1437 ("hTrackMatchedMCParticlePt",
1438 "Origin of particle vs energy",
1439 nptbins,ptmin,ptmax,8,0,8);
1440 fhTrackMatchedMCParticlePt->SetXTitle("p_{T} (GeV/c)");
1441 //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1443 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1444 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1445 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1446 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1447 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1448 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1449 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1450 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1452 outputContainer->Add(fhTrackMatchedMCParticlePt);
1454 fhTrackMatchedMCParticleDEta = new TH2F
1455 ("hTrackMatchedMCParticleDEta",
1456 "Origin of particle vs #eta residual",
1457 nresetabins,resetamin,resetamax,8,0,8);
1458 fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1459 //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1461 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1462 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1463 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1464 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1465 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1466 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1467 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1468 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1470 outputContainer->Add(fhTrackMatchedMCParticleDEta);
1472 fhTrackMatchedMCParticleDPhi = new TH2F
1473 ("hTrackMatchedMCParticleDPhi",
1474 "Origin of particle vs #phi residual",
1475 nresphibins,resphimin,resphimax,8,0,8);
1476 fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1477 //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1479 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1480 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1481 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1482 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1483 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1484 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1485 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1486 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1488 outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1494 if(fFillWeightHistograms)
1496 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1497 nptbins,ptmin,ptmax, 100,0,1.);
1498 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1499 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1500 outputContainer->Add(fhECellClusterRatio);
1502 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1503 nptbins,ptmin,ptmax, 100,-10,0);
1504 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1505 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1506 outputContainer->Add(fhECellClusterLogRatio);
1508 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1509 nptbins,ptmin,ptmax, 100,0,1.);
1510 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1511 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1512 outputContainer->Add(fhEMaxCellClusterRatio);
1514 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1515 nptbins,ptmin,ptmax, 100,-10,0);
1516 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1517 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1518 outputContainer->Add(fhEMaxCellClusterLogRatio);
1520 for(Int_t iw = 0; iw < 14; iw++)
1522 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),
1523 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1524 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1525 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1526 outputContainer->Add(fhLambda0ForW0[iw]);
1528 // 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),
1529 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1530 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1531 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1532 // outputContainer->Add(fhLambda1ForW0[iw]);
1539 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1541 fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
1542 nptbins,ptmin,ptmax,200,0,2);
1543 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1544 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1545 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1547 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1548 nptbins,ptmin,ptmax,200,0,2);
1549 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1550 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1551 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1553 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1554 fhMCPi0DecayPt->SetYTitle("N");
1555 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1556 outputContainer->Add(fhMCPi0DecayPt) ;
1558 fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #pi^{0}",
1559 nptbins,ptmin,ptmax,100,0,1);
1560 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1561 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1562 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1564 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1565 fhMCEtaDecayPt->SetYTitle("N");
1566 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1567 outputContainer->Add(fhMCEtaDecayPt) ;
1569 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
1570 nptbins,ptmin,ptmax,100,0,1);
1571 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1572 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1573 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1575 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1576 fhMCOtherDecayPt->SetYTitle("N");
1577 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1578 outputContainer->Add(fhMCOtherDecayPt) ;
1582 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1583 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1586 fhAnglePairMCPi0 = new TH2F
1588 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1589 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1590 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1591 outputContainer->Add(fhAnglePairMCPi0) ;
1593 if (fAnaType!= kSSCalo)
1595 fhAnglePairMCEta = new TH2F
1597 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1598 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1599 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1600 outputContainer->Add(fhAnglePairMCEta) ;
1602 fhMassPairMCPi0 = new TH2F
1604 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1605 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1606 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1607 outputContainer->Add(fhMassPairMCPi0) ;
1609 fhMassPairMCEta = new TH2F
1611 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1612 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1613 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1614 outputContainer->Add(fhMassPairMCEta) ;
1617 for(Int_t i = 0; i < 6; i++)
1621 (Form("hE_MC%s",pname[i].Data()),
1622 Form("Identified as #pi^{0} (#eta), cluster from %s",
1624 nptbins,ptmin,ptmax);
1625 fhMCE[i]->SetYTitle("N");
1626 fhMCE[i]->SetXTitle("E (GeV)");
1627 outputContainer->Add(fhMCE[i]) ;
1629 fhMCPt[i] = new TH1F
1630 (Form("hPt_MC%s",pname[i].Data()),
1631 Form("Identified as #pi^{0} (#eta), cluster from %s",
1633 nptbins,ptmin,ptmax);
1634 fhMCPt[i]->SetYTitle("N");
1635 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1636 outputContainer->Add(fhMCPt[i]) ;
1638 fhMCPtCentrality[i] = new TH2F
1639 (Form("hPtCentrality_MC%s",pname[i].Data()),
1640 Form("Identified as #pi^{0} (#eta), cluster from %s",
1642 nptbins,ptmin,ptmax, 100,0,100);
1643 fhMCPtCentrality[i]->SetYTitle("centrality");
1644 fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1645 outputContainer->Add(fhMCPtCentrality[i]) ;
1647 if(fAnaType == kSSCalo)
1649 fhMCNLocMaxPt[i] = new TH2F
1650 (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1651 Form("cluster from %s, pT of cluster vs NLM, accepted",ptype[i].Data()),
1652 nptbins,ptmin,ptmax,20,0,20);
1653 fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
1654 fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
1655 outputContainer->Add(fhMCNLocMaxPt[i]) ;
1657 fhMCNLocMaxPtReject[i] = new TH2F
1658 (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1659 Form("cluster from %s, pT of cluster vs NLM, rejected",ptype[i].Data()),
1660 nptbins,ptmin,ptmax,20,0,20);
1661 fhMCNLocMaxPtReject[i] ->SetYTitle("N maxima");
1662 fhMCNLocMaxPtReject[i] ->SetXTitle("p_{T} (GeV/c)");
1663 outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1665 fhMCEReject[i] = new TH1F
1666 (Form("hEReject_MC%s",pname[i].Data()),
1667 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1669 nptbins,ptmin,ptmax);
1670 fhMCEReject[i]->SetYTitle("N");
1671 fhMCEReject[i]->SetXTitle("E (GeV)");
1672 outputContainer->Add(fhMCEReject[i]) ;
1674 fhMCPtReject[i] = new TH1F
1675 (Form("hPtReject_MC%s",pname[i].Data()),
1676 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1678 nptbins,ptmin,ptmax);
1679 fhMCPtReject[i]->SetYTitle("N");
1680 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1681 outputContainer->Add(fhMCPtReject[i]) ;
1684 fhMCPtPhi[i] = new TH2F
1685 (Form("hPtPhi_MC%s",pname[i].Data()),
1686 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1687 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1688 fhMCPtPhi[i]->SetYTitle("#phi");
1689 fhMCPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
1690 outputContainer->Add(fhMCPtPhi[i]) ;
1692 fhMCPtEta[i] = new TH2F
1693 (Form("hPtEta_MC%s",pname[i].Data()),
1694 Form("Identified as #pi^{0} (#eta), cluster from %s",
1695 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1696 fhMCPtEta[i]->SetYTitle("#eta");
1697 fhMCPtEta[i]->SetXTitle("p_{T} (GeV/c)");
1698 outputContainer->Add(fhMCPtEta[i]) ;
1700 fhMCMassPt[i] = new TH2F
1701 (Form("hMassPt_MC%s",pname[i].Data()),
1702 Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1703 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1704 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1705 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1706 outputContainer->Add(fhMCMassPt[i]) ;
1708 fhMCSelectedMassPt[i] = new TH2F
1709 (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1710 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1711 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1712 fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1713 fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1714 outputContainer->Add(fhMCSelectedMassPt[i]) ;
1716 if(fAnaType == kSSCalo)
1718 fhMCMassPtNoOverlap[i] = new TH2F
1719 (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1720 Form("all pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1721 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1722 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1723 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1724 outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1726 fhMCSelectedMassPtNoOverlap[i] = new TH2F
1727 (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1728 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1729 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1730 fhMCSelectedMassPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
1731 fhMCSelectedMassPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
1732 outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1735 if( fFillSelectClHisto )
1737 fhMCPtLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1738 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1739 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1740 fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1741 fhMCPtLambda0[i]->SetXTitle("p_{T} (GeV/c)");
1742 outputContainer->Add(fhMCPtLambda0[i]) ;
1744 fhMCPtLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1745 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1746 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1747 fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1748 fhMCPtLambda1[i]->SetXTitle("p_{T} (GeV/c)");
1749 outputContainer->Add(fhMCPtLambda1[i]) ;
1751 fhMCPtDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1752 Form("Selected pair, cluster from %s : p_{T} vs dispersion^{2}",ptype[i].Data()),
1753 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1754 fhMCPtDispersion[i]->SetYTitle("D^{2}");
1755 fhMCPtDispersion[i]->SetXTitle("p_{T} (GeV/c)");
1756 outputContainer->Add(fhMCPtDispersion[i]) ;
1758 if(fCalorimeter=="EMCAL")
1760 fhMCPtLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1761 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1762 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1763 fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1764 fhMCPtLambda0NoTRD[i]->SetXTitle("p_{T} (GeV/c)");
1765 outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1767 if(!fFillOnlySimpleSSHisto)
1769 fhMCPtDispEta[i] = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1770 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs p_{T}",ptype[i].Data()),
1771 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1772 fhMCPtDispEta[i]->SetXTitle("p_{T} (GeV/c)");
1773 fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1774 outputContainer->Add(fhMCPtDispEta[i]);
1776 fhMCPtDispPhi[i] = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1777 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs p_{T}",ptype[i].Data()),
1778 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1779 fhMCPtDispPhi[i]->SetXTitle("p_{T} (GeV/c)");
1780 fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1781 outputContainer->Add(fhMCPtDispPhi[i]);
1783 fhMCPtSumEtaPhi[i] = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1784 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs p_{T}",ptype[i].Data()),
1785 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1786 fhMCPtSumEtaPhi[i]->SetXTitle("p_{T} (GeV/c)");
1787 fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1788 outputContainer->Add(fhMCPtSumEtaPhi[i]);
1790 fhMCPtDispEtaPhiDiff[i] = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1791 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs p_{T}",ptype[i].Data()),
1792 nptbins,ptmin,ptmax,200,-10,10);
1793 fhMCPtDispEtaPhiDiff[i]->SetXTitle("p_{T} (GeV/c)");
1794 fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1795 outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1797 fhMCPtSphericity[i] = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1798 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()),
1799 nptbins,ptmin,ptmax, 200,-1,1);
1800 fhMCPtSphericity[i]->SetXTitle("p_{T} (GeV/c)");
1801 fhMCPtSphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1802 outputContainer->Add(fhMCPtSphericity[i]);
1804 for(Int_t ie = 0; ie < 7; ie++)
1806 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1807 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]),
1808 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1809 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1810 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1811 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1813 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1814 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]),
1815 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1816 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1817 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1818 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1820 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1821 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]),
1822 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1823 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1824 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1825 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1831 fhMCPtLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1832 Form("Selected pair, cluster from %s : p_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1833 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1834 fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1835 fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1836 outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
1838 fhMCPtFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1839 Form("Selected pair, cluster from %s : p_{T} vs Max cell fraction of energy",ptype[i].Data()),
1840 nptbins,ptmin,ptmax,100,0,1);
1841 fhMCPtFracMaxCell[i]->SetYTitle("Fraction");
1842 fhMCPtFracMaxCell[i]->SetXTitle("E (GeV)");
1843 outputContainer->Add(fhMCPtFracMaxCell[i]) ;
1846 } // shower shape histo
1851 if(fAnaType==kSSCalo)
1853 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1854 nptbins,ptmin,ptmax, 200, -1,1);
1855 fhAsymmetry->SetXTitle("E (GeV)");
1856 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1857 outputContainer->Add(fhAsymmetry);
1859 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1860 nptbins,ptmin,ptmax, 200, -1,1);
1861 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1862 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1863 outputContainer->Add(fhSelectedAsymmetry);
1866 ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1867 fhSplitE->SetYTitle("counts");
1868 fhSplitE->SetXTitle("E (GeV)");
1869 outputContainer->Add(fhSplitE) ;
1871 fhSplitPt = new TH1F
1872 ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
1873 fhSplitPt->SetYTitle("counts");
1874 fhSplitPt->SetXTitle("p_{T} (GeV/c)");
1875 outputContainer->Add(fhSplitPt) ;
1878 fhSplitPtPhi = new TH2F
1879 ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1880 fhSplitPtPhi->SetYTitle("#phi (rad)");
1881 fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
1882 outputContainer->Add(fhSplitPtPhi) ;
1884 fhSplitPtEta = new TH2F
1885 ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1886 fhSplitPtEta->SetYTitle("#eta");
1887 fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
1888 outputContainer->Add(fhSplitPtEta) ;
1891 fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1892 nptbins,ptmin,ptmax,20,0,20);
1893 fhNLocMaxSplitPt ->SetYTitle("N maxima");
1894 fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
1895 outputContainer->Add(fhNLocMaxSplitPt) ;
1898 fhMassSplitPt = new TH2F
1899 ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",
1900 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1901 fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1902 fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1903 outputContainer->Add(fhMassSplitPt) ;
1905 fhSelectedMassSplitPt = new TH2F
1906 ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",
1907 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1908 fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1909 fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1910 outputContainer->Add(fhSelectedMassSplitPt) ;
1914 fhMassSplitPtNoOverlap = new TH2F
1915 ("hMassSplitPtNoOverlap","all pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1916 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1917 fhMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1918 fhMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1919 outputContainer->Add(fhMassSplitPtNoOverlap) ;
1921 fhSelectedMassSplitPtNoOverlap = new TH2F
1922 ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1923 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1924 fhSelectedMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1925 fhSelectedMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1926 outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
1929 fhMCPi0PtRecoPtPrim = new TH2F
1930 ("hMCPi0PtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1931 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1932 fhMCPi0PtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1933 fhMCPi0PtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1934 outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
1936 fhMCPi0PtRecoPtPrimNoOverlap = new TH2F
1937 ("hMCPi0PtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1938 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1939 fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1940 fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1941 outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
1943 fhMCPi0SelectedPtRecoPtPrim = new TH2F
1944 ("hMCPi0SelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1945 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1946 fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1947 fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1948 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
1950 fhMCPi0SelectedPtRecoPtPrimNoOverlap = new TH2F
1951 ("hMCPi0SelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1952 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1953 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1954 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1955 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
1958 fhMCPi0SplitPtRecoPtPrim = new TH2F
1959 ("hMCPi0SplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1960 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1961 fhMCPi0SplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1962 fhMCPi0SplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1963 outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
1965 fhMCPi0SplitPtRecoPtPrimNoOverlap = new TH2F
1966 ("hMCPi0SplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1967 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1968 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1969 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1970 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
1972 fhMCPi0SelectedSplitPtRecoPtPrim = new TH2F
1973 ("hMCPi0SelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1974 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1975 fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1976 fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1977 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
1979 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap = new TH2F
1980 ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1981 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1982 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1983 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1984 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
1986 fhMCEtaPtRecoPtPrim = new TH2F
1987 ("hMCEtaPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1988 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1989 fhMCEtaPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1990 fhMCEtaPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1991 outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
1993 fhMCEtaPtRecoPtPrimNoOverlap = new TH2F
1994 ("hMCEtaPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1995 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1996 fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1997 fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1998 outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2000 fhMCEtaSelectedPtRecoPtPrim = new TH2F
2001 ("hMCEtaSelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
2002 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2003 fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2004 fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2005 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2007 fhMCEtaSelectedPtRecoPtPrimNoOverlap = new TH2F
2008 ("hMCEtaSelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
2009 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2010 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2011 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2012 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2015 fhMCEtaSplitPtRecoPtPrim = new TH2F
2016 ("hMCEtaSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
2017 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2018 fhMCEtaSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2019 fhMCEtaSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2020 outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2022 fhMCEtaSplitPtRecoPtPrimNoOverlap = new TH2F
2023 ("hMCEtaSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2024 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2025 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2026 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2027 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2029 fhMCEtaSelectedSplitPtRecoPtPrim = new TH2F
2030 ("hMCEtaSelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
2031 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2032 fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
2033 fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
2034 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2036 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2037 ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
2038 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2039 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
2040 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
2041 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2043 for(Int_t i = 0; i< 6; i++)
2045 fhMCPtAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2046 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
2047 nptbins,ptmin,ptmax, 200,-1,1);
2048 fhMCPtAsymmetry[i]->SetXTitle("E (GeV)");
2049 fhMCPtAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2050 outputContainer->Add(fhMCPtAsymmetry[i]);
2052 fhMCSplitE[i] = new TH1F
2053 (Form("hSplitE_MC%s",pname[i].Data()),
2054 Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2055 nptbins,ptmin,ptmax);
2056 fhMCSplitE[i]->SetYTitle("counts");
2057 fhMCSplitE[i]->SetXTitle("E (GeV)");
2058 outputContainer->Add(fhMCSplitE[i]) ;
2060 fhMCSplitPt[i] = new TH1F
2061 (Form("hSplitPt_MC%s",pname[i].Data()),
2062 Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
2063 nptbins,ptmin,ptmax);
2064 fhMCSplitPt[i]->SetYTitle("counts");
2065 fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2066 outputContainer->Add(fhMCSplitPt[i]) ;
2069 fhMCSplitPtPhi[i] = new TH2F
2070 (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2071 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2072 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2073 fhMCSplitPtPhi[i]->SetYTitle("#phi");
2074 fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
2075 outputContainer->Add(fhMCSplitPtPhi[i]) ;
2077 fhMCSplitPtEta[i] = new TH2F
2078 (Form("hSplitPtEta_MC%s",pname[i].Data()),
2079 Form("Identified as #pi^{0} (#eta), cluster from %s",
2080 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2081 fhMCSplitPtEta[i]->SetYTitle("#eta");
2082 fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
2083 outputContainer->Add(fhMCSplitPtEta[i]) ;
2086 fhMCNLocMaxSplitPt[i] = new TH2F
2087 (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2088 Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
2089 nptbins,ptmin,ptmax,20,0,20);
2090 fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
2091 fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
2092 outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2094 fhMCMassSplitPt[i] = new TH2F
2095 (Form("hMassSplitPt_MC%s",pname[i].Data()),
2096 Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2097 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2098 fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2099 fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2100 outputContainer->Add(fhMCMassSplitPt[i]) ;
2102 fhMCSelectedMassSplitPt[i] = new TH2F
2103 (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2104 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2105 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2106 fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2107 fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2108 outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2110 fhMCMassSplitPtNoOverlap[i] = new TH2F
2111 (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2112 Form("all pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2113 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2114 fhMCMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2115 fhMCMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2116 outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2118 fhMCSelectedMassSplitPtNoOverlap[i] = new TH2F
2119 (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2120 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2121 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2122 fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2123 fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2124 outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2129 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2133 for(Int_t i = 0; i< 3; i++)
2135 fhPtAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2136 Form("Selected #pi^{0} (#eta) pairs: p_{T} vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
2137 nptbins,ptmin,ptmax,200, -1,1);
2138 fhPtAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2139 fhPtAsymmetryLocMax[i]->SetXTitle("p_{T} (GeV/c)");
2140 outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2143 for(Int_t ie = 0; ie< 7; ie++)
2146 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2147 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2148 ssbins,ssmin,ssmax , 200,-1,1);
2149 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2150 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2151 outputContainer->Add(fhAsymmetryLambda0[ie]);
2153 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2154 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2155 ssbins,ssmin,ssmax , 200,-1,1);
2156 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2157 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2158 outputContainer->Add(fhAsymmetryDispEta[ie]);
2160 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2161 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2162 ssbins,ssmin,ssmax , 200,-1,1);
2163 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2164 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2165 outputContainer->Add(fhAsymmetryDispPhi[ie]);
2171 for(Int_t i = 0; i< 6; i++)
2173 for(Int_t ie = 0; ie < 7; ie++)
2175 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2176 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2177 ssbins,ssmin,ssmax , 200,-1,1);
2178 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2179 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2180 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2182 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2183 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2184 ssbins,ssmin,ssmax , 200,-1,1);
2185 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2186 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2187 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2189 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2190 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2191 ssbins,ssmin,ssmax , 200,-1,1);
2192 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2193 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2194 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2200 if(fFillPileUpHistograms)
2203 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2205 for(Int_t i = 0 ; i < 7 ; i++)
2207 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2208 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2209 fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2210 outputContainer->Add(fhPtPileUp[i]);
2212 fhPtCellTimePileUp[i] = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2213 Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2214 nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2215 fhPtCellTimePileUp[i]->SetXTitle("p_{T} (GeV/c)");
2216 fhPtCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
2217 outputContainer->Add(fhPtCellTimePileUp[i]);
2219 fhPtTimeDiffPileUp[i] = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2220 Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2221 nptbins,ptmin,ptmax,400,-200,200);
2222 fhPtTimeDiffPileUp[i]->SetXTitle("p_{T} (GeV/c");
2223 fhPtTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2224 outputContainer->Add(fhPtTimeDiffPileUp[i]);
2228 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2229 fhTimePtNoCut->SetXTitle("p_{T} (GeV/c)");
2230 fhTimePtNoCut->SetYTitle("time (ns)");
2231 outputContainer->Add(fhTimePtNoCut);
2233 fhTimePtSPD = new TH2F ("hTimePt_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2234 fhTimePtSPD->SetXTitle("p_{T} (GeV/c)");
2235 fhTimePtSPD->SetYTitle("time (ns)");
2236 outputContainer->Add(fhTimePtSPD);
2238 fhTimePtSPDMulti = new TH2F ("hTimePt_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2239 fhTimePtSPDMulti->SetXTitle("p_{T} (GeV/c)");
2240 fhTimePtSPDMulti->SetYTitle("time (ns)");
2241 outputContainer->Add(fhTimePtSPDMulti);
2243 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2244 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2245 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
2246 outputContainer->Add(fhTimeNPileUpVertSPD);
2248 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2249 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2250 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
2251 outputContainer->Add(fhTimeNPileUpVertTrack);
2253 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2254 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2255 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
2256 outputContainer->Add(fhTimeNPileUpVertContributors);
2258 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
2259 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
2260 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
2261 outputContainer->Add(fhTimePileUpMainVertexZDistance);
2263 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
2264 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
2265 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
2266 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2268 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
2269 nptbins,ptmin,ptmax,20,0,20);
2270 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2271 fhPtNPileUpSPDVtx->SetXTitle("p_{T} (GeV/c)");
2272 outputContainer->Add(fhPtNPileUpSPDVtx);
2274 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
2275 nptbins,ptmin,ptmax, 20,0,20 );
2276 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2277 fhPtNPileUpTrkVtx->SetXTitle("p_{T} (GeV/c)");
2278 outputContainer->Add(fhPtNPileUpTrkVtx);
2280 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2281 nptbins,ptmin,ptmax,20,0,20);
2282 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2283 fhPtNPileUpSPDVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2284 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2286 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2287 nptbins,ptmin,ptmax, 20,0,20 );
2288 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2289 fhPtNPileUpTrkVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2290 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2292 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2293 nptbins,ptmin,ptmax,20,0,20);
2294 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2295 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2296 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2298 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2299 nptbins,ptmin,ptmax, 20,0,20 );
2300 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2301 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2302 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2306 //Keep neutral meson selection histograms if requiered
2307 //Setting done in AliNeutralMesonSelection
2309 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2311 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2313 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2314 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2319 return outputContainer ;
2323 //_____________________________________________
2324 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2327 // Assign mc index depending on MC bit set
2329 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
2333 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2337 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2338 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2340 return kmcConversion ;
2341 }//conversion photon
2342 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2345 }//photon no conversion
2346 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2348 return kmcElectron ;
2357 //__________________________________________________________________
2358 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2359 AliAODPWG4Particle * photon2,
2360 Int_t & label, Int_t & tag)
2362 // Check the labels of pare in case mother was same pi0 or eta
2363 // Set the new AOD accordingly
2365 Int_t label1 = photon1->GetLabel();
2366 Int_t label2 = photon2->GetLabel();
2368 if(label1 < 0 || label2 < 0 ) return ;
2370 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2371 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2372 Int_t tag1 = photon1->GetTag();
2373 Int_t tag2 = photon2->GetTag();
2375 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2376 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2377 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
2378 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2379 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
2383 //Check if pi0/eta mother is the same
2384 if(GetReader()->ReadStack())
2388 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2389 label1 = mother1->GetFirstMother();
2390 //mother1 = GetMCStack()->Particle(label1);//pi0
2394 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2395 label2 = mother2->GetFirstMother();
2396 //mother2 = GetMCStack()->Particle(label2);//pi0
2399 else if(GetReader()->ReadAODMCParticles())
2400 {//&& (input > -1)){
2403 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2404 label1 = mother1->GetMother();
2405 //mother1 = GetMCStack()->Particle(label1);//pi0
2409 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2410 label2 = mother2->GetMother();
2411 //mother2 = GetMCStack()->Particle(label2);//pi0
2415 //printf("mother1 %d, mother2 %d\n",label1,label2);
2416 if( label1 == label2 && label1>=0 )
2421 TLorentzVector mom1 = *(photon1->Momentum());
2422 TLorentzVector mom2 = *(photon2->Momentum());
2424 Double_t angle = mom2.Angle(mom1.Vect());
2425 Double_t mass = (mom1+mom2).M();
2426 Double_t epair = (mom1+mom2).E();
2428 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2430 fhMassPairMCPi0 ->Fill(epair,mass);
2431 fhAnglePairMCPi0->Fill(epair,angle);
2432 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2436 fhMassPairMCEta ->Fill(epair,mass);
2437 fhAnglePairMCEta->Fill(epair,angle);
2438 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2442 } // both from eta or pi0 decay
2446 //____________________________________________________________________________
2447 void AliAnaPi0EbE::Init()
2451 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2452 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2455 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2456 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2462 //____________________________________________________________________________
2463 void AliAnaPi0EbE::InitParameters()
2465 //Initialize the parameters of the analysis.
2466 AddToHistogramsName("AnaPi0EbE_");
2468 fInputAODGammaConvName = "PhotonsCTS" ;
2469 fAnaType = kIMCalo ;
2470 fCalorimeter = "EMCAL" ;
2475 fNLMECutMin[0] = 10.;
2476 fNLMECutMin[1] = 6. ;
2477 fNLMECutMin[2] = 6. ;
2481 //__________________________________________________________________
2482 void AliAnaPi0EbE::MakeAnalysisFillAOD()
2484 //Do analysis and fill aods
2489 MakeInvMassInCalorimeter();
2493 MakeShowerShapeIdentification();
2497 MakeInvMassInCalorimeterAndCTS();
2503 //____________________________________________
2504 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
2506 //Do analysis and fill aods
2507 //Search for the photon decay in calorimeters
2508 //Read photon list from AOD, produced in class AliAnaPhoton
2509 //Check if 2 photons have the mass of the pi0.
2511 TLorentzVector mom1;
2512 TLorentzVector mom2;
2513 TLorentzVector mom ;
2518 if(!GetInputAODBranch()){
2519 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2523 //Get shower shape information of clusters
2524 TObjArray *clusters = 0;
2525 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2526 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2528 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2529 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2531 //Vertex cut in case of mixed events
2532 Int_t evtIndex1 = 0 ;
2534 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2535 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
2536 mom1 = *(photon1->Momentum());
2538 //Get original cluster, to recover some information
2540 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2543 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2547 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2549 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2551 Int_t evtIndex2 = 0 ;
2553 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2555 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2558 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
2560 mom2 = *(photon2->Momentum());
2562 //Get original cluster, to recover some information
2564 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2568 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2572 Float_t e1 = photon1->E();
2573 Float_t e2 = photon2->E();
2575 //Select clusters with good time window difference
2576 Float_t tof1 = cluster1->GetTOF()*1e9;;
2577 Float_t tof2 = cluster2->GetTOF()*1e9;;
2578 Double_t t12diff = tof1-tof2;
2579 fhEPairDiffTime->Fill(e1+e2, t12diff);
2580 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2582 //Play with the MC stack if available
2583 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2585 // Check the invariant mass for different selection on the local maxima
2586 // Name of AOD method TO BE FIXED
2587 Int_t nMaxima1 = photon1->GetFiducialArea();
2588 Int_t nMaxima2 = photon2->GetFiducialArea();
2590 Double_t mass = (mom1+mom2).M();
2591 Double_t epair = (mom1+mom2).E();
2593 if(nMaxima1==nMaxima2)
2595 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2596 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2597 else fhMassPairLocMax[2]->Fill(epair,mass);
2599 else if(nMaxima1==1 || nMaxima2==1)
2601 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2602 else fhMassPairLocMax[4]->Fill(epair,mass);
2605 fhMassPairLocMax[5]->Fill(epair,mass);
2607 // combinations with SS axis cut and NLM cut
2608 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2609 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2610 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2611 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2613 //Skip events with too few or too many NLM
2614 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2616 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2619 fhMass->Fill(epair,(mom1+mom2).M());
2621 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2622 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2625 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());
2627 //Fill some histograms about shower shape
2628 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2630 FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2631 FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2634 // Tag both photons as decay
2635 photon1->SetTagged(kTRUE);
2636 photon2->SetTagged(kTRUE);
2638 fhPtDecay->Fill(photon1->Pt());
2639 fhEDecay ->Fill(photon1->E() );
2641 fhPtDecay->Fill(photon2->Pt());
2642 fhEDecay ->Fill(photon2->E() );
2644 //Create AOD for analysis
2647 //Mass of selected pairs
2648 fhSelectedMass->Fill(epair,mom.M());
2650 // Fill histograms to undertand pile-up before other cuts applied
2651 // Remember to relax time cuts in the reader
2652 FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2654 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2656 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2657 pi0.SetDetector(photon1->GetDetector());
2660 pi0.SetLabel(label);
2663 //Set the indeces of the original caloclusters
2664 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2665 //pi0.SetInputFileIndex(input);
2667 AddAODParticle(pi0);
2675 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2679 //__________________________________________________
2680 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2682 //Do analysis and fill aods
2683 //Search for the photon decay in calorimeters
2684 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2685 //Check if 2 photons have the mass of the pi0.
2687 TLorentzVector mom1;
2688 TLorentzVector mom2;
2689 TLorentzVector mom ;
2694 // Check calorimeter input
2695 if(!GetInputAODBranch()){
2696 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2700 // Get the array with conversion photons
2701 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2702 if(!inputAODGammaConv) {
2704 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2706 if(!inputAODGammaConv) {
2707 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2713 //Get shower shape information of clusters
2714 TObjArray *clusters = 0;
2715 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2716 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2718 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
2719 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2720 if(nCTS<=0 || nCalo <=0)
2722 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2727 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2729 // Do the loop, first calo, second CTS
2730 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2731 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2732 mom1 = *(photon1->Momentum());
2734 //Get original cluster, to recover some information
2736 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2738 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2739 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2741 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2742 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2744 mom2 = *(photon2->Momentum());
2746 Double_t mass = (mom1+mom2).M();
2747 Double_t epair = (mom1+mom2).E();
2749 Int_t nMaxima = photon1->GetFiducialArea();
2750 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2751 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2752 else fhMassPairLocMax[2]->Fill(epair,mass);
2754 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2755 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2757 //Play with the MC stack if available
2760 Int_t label2 = photon2->GetLabel();
2761 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2763 HasPairSameMCMother(photon1, photon2, label, tag) ;
2766 //Mass of selected pairs
2767 fhMass->Fill(epair,(mom1+mom2).M());
2769 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2770 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2772 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());
2774 //Fill some histograms about shower shape
2775 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2777 FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
2780 // Tag both photons as decay
2781 photon1->SetTagged(kTRUE);
2782 photon2->SetTagged(kTRUE);
2784 fhPtDecay->Fill(photon1->Pt());
2785 fhEDecay ->Fill(photon1->E() );
2787 //Create AOD for analysis
2791 //Mass of selected pairs
2792 fhSelectedMass->Fill(epair,mom.M());
2794 // Fill histograms to undertand pile-up before other cuts applied
2795 // Remember to relax time cuts in the reader
2796 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2798 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2800 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2801 pi0.SetDetector(photon1->GetDetector());
2804 pi0.SetLabel(label);
2807 //Set the indeces of the original tracks or caloclusters
2808 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2809 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2810 //pi0.SetInputFileIndex(input);
2812 AddAODParticle(pi0);
2819 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2824 //_________________________________________________
2825 void AliAnaPi0EbE::MakeShowerShapeIdentification()
2827 //Search for pi0 in fCalorimeter with shower shape analysis
2829 TObjArray * pl = 0x0;
2830 AliVCaloCells * cells = 0x0;
2831 //Select the Calorimeter of the photon
2832 if (fCalorimeter == "PHOS" )
2834 pl = GetPHOSClusters();
2835 cells = GetPHOSCells();
2837 else if (fCalorimeter == "EMCAL")
2839 pl = GetEMCALClusters();
2840 cells = GetEMCALCells();
2845 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2849 TLorentzVector mom ;
2850 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2852 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2854 Int_t evtIndex = 0 ;
2855 if (GetMixedEvent())
2857 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2860 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2862 //Get Momentum vector,
2863 Double_t vertex[]={0,0,0};
2864 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2866 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2867 }//Assume that come from vertex in straight line
2870 calo->GetMomentum(mom,vertex) ;
2873 //If too small or big pt, skip it
2874 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2876 //Check acceptance selection
2877 if(IsFiducialCutOn())
2879 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2880 if(! in ) continue ;
2884 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());
2886 //Play with the MC stack if available
2887 //Check origin of the candidates
2891 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
2892 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2893 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2896 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
2898 //Check Distance to Bad channel, set bit.
2899 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2900 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2901 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
2902 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2906 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2908 //If too low number of cells, skip it
2909 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
2911 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2916 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
2917 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
2919 //.......................................
2920 // TOF cut, BE CAREFUL WITH THIS CUT
2921 Double_t tof = calo->GetTOF()*1e9;
2922 if(tof < fTimeCutMin || tof > fTimeCutMax)
2924 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2929 //PID selection or bit setting
2931 Double_t mass = 0, angle = 0;
2932 Int_t absId1 =-1, absId2 =-1;
2933 Float_t distbad1 =-1, distbad2 =-1;
2934 Bool_t fidcut1 = 0, fidcut2 = 0;
2935 TLorentzVector l1, l2;
2937 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2938 GetVertex(evtIndex),nMaxima,
2939 mass,angle,l1,l2,absId1,absId2,
2940 distbad1,distbad2,fidcut1,fidcut2) ;
2943 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2946 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
2947 if( (fCheckSplitDistToBad) &&
2948 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
2951 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
2952 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
2954 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2958 //Skip events with too few or too many NLM
2959 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2961 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2966 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2968 //Skip matched clusters with tracks
2969 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
2971 FillRejectedClusterHistograms(mom,tag,nMaxima);
2975 Float_t e1 = l1.Energy();
2976 Float_t e2 = l2.Energy();
2977 TLorentzVector l12 = l1+l2;
2978 Float_t ptSplit = l12.Pt();
2979 Float_t eSplit = e1+e2;
2982 Int_t noverlaps = 0;
2986 mcIndex = GetMCIndex(tag);
2989 Int_t mcLabel = calo->GetLabel();
2991 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
2993 Int_t mesonLabel = -1;
2995 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
2997 if(mcIndex == kmcPi0)
2999 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3000 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3004 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3005 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3009 const UInt_t nlabels = calo->GetNLabels();
3010 Int_t overpdg[nlabels];
3011 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3014 //mass of all clusters
3015 fhMass ->Fill(mom.E() ,mass);
3016 fhMassPt ->Fill(mom.Pt(),mass);
3017 fhMassSplitPt->Fill(ptSplit ,mass);
3019 Int_t indexMax = -1;
3020 if (nMaxima==1) indexMax = 0 ;
3021 else if(nMaxima==2) indexMax = 1 ;
3023 fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3027 fhMCMassPt[mcIndex] ->Fill(mom.Pt(),mass);
3028 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3031 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3032 fhMCPi0SplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3034 else if(mcIndex==kmcEta)
3036 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3037 fhMCEtaSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3044 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3045 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3047 else if(mcIndex==kmcEta)
3049 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3050 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3053 fhMassNoOverlap ->Fill(mom.E() ,mass);
3054 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3055 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3057 fhMCMassPtNoOverlap[mcIndex] ->Fill(mom.Pt(),mass);
3058 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3062 // Asymmetry of all clusters
3065 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3066 fhAsymmetry->Fill(mom.E(),asy);
3070 fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3073 // If cluster does not pass pid, not pi0/eta, skip it.
3074 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3076 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3077 FillRejectedClusterHistograms(mom,tag,nMaxima);
3081 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3083 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3084 FillRejectedClusterHistograms(mom,tag,nMaxima);
3089 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3090 mom.Pt(), idPartType);
3092 //Mass and asymmetry of selected pairs
3093 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3094 fhSelectedMass ->Fill(mom.E() ,mass);
3095 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3096 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3097 fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3103 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3104 fhMCPi0SelectedSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3106 else if(mcIndex==kmcEta)
3108 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3109 fhMCEtaSelectedSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3114 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3115 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3116 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3120 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3121 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3123 else if(mcIndex==kmcEta)
3125 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3126 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3131 fhSplitE ->Fill( eSplit);
3132 fhSplitPt ->Fill(ptSplit);
3133 Float_t phi = mom.Phi();
3134 if(phi<0) phi+=TMath::TwoPi();
3135 fhSplitPtPhi ->Fill(ptSplit,phi);
3136 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3137 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3139 //Check split-clusters with good time window difference
3140 Double_t tof1 = cells->GetCellTime(absId1);
3141 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3144 Double_t tof2 = cells->GetCellTime(absId2);
3145 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3148 Double_t t12diff = tof1-tof2;
3149 fhEPairDiffTime->Fill(e1+e2, t12diff);
3153 fhMCSplitE [mcIndex]->Fill( eSplit);
3154 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3155 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3156 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3157 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3158 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3160 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3161 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3162 fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3166 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3167 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3171 //-----------------------
3172 //Create AOD for analysis
3174 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3175 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3176 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3178 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3179 aodpi0.SetLabel(calo->GetLabel());
3181 //Set the indeces of the original caloclusters
3182 aodpi0.SetCaloLabel(calo->GetID(),-1);
3183 aodpi0.SetDetector(fCalorimeter);
3185 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3186 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3187 else aodpi0.SetDistToBad(0) ;
3189 // Check if cluster is pi0 via cluster splitting
3190 aodpi0.SetIdentifiedParticleType(idPartType);
3192 // Add number of local maxima to AOD, method name in AOD to be FIXED
3193 aodpi0.SetFiducialArea(nMaxima);
3197 //Fill some histograms about shower shape
3198 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3200 FillSelectedClusterHistograms(calo, aodpi0.Pt(), nMaxima, tag, asy);
3203 // Fill histograms to undertand pile-up before other cuts applied
3204 // Remember to relax time cuts in the reader
3205 Double_t tofcluster = calo->GetTOF()*1e9;
3206 Double_t tofclusterUS = TMath::Abs(tofcluster);
3208 FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
3210 Int_t id = GetReader()->GetTriggerClusterId();
3211 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
3213 Float_t phicluster = aodpi0.Phi();
3214 if(phicluster < 0) phicluster+=TMath::TwoPi();
3218 if (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
3219 else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
3220 else fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
3223 Int_t bc = GetReader()->GetTriggerClusterBC();
3224 if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
3226 if(GetReader()->IsTriggerMatched())
3228 if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
3229 fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
3230 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
3234 if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
3235 fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
3239 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(calo->E(), tofcluster);
3240 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
3241 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(calo->E(), tofcluster);
3245 else if(TMath::Abs(bc) >= 6)
3246 Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
3249 //Add AOD with pi0 object to aod branch
3250 AddAODParticle(aodpi0);
3254 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3257 //______________________________________________
3258 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3260 //Do analysis and fill histograms
3262 if(!GetOutputAODBranch())
3264 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3266 //Loop on stored AOD pi0
3267 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3268 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3270 Float_t cen = GetEventCentrality();
3271 Float_t ep = GetEventPlaneAngle();
3273 for(Int_t iaod = 0; iaod < naod ; iaod++)
3276 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3277 Int_t pdg = pi0->GetIdentifiedParticleType();
3279 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
3281 //Fill pi0 histograms
3282 Float_t ener = pi0->E();
3283 Float_t pt = pi0->Pt();
3284 Float_t phi = pi0->Phi();
3285 if(phi < 0) phi+=TMath::TwoPi();
3286 Float_t eta = pi0->Eta();
3291 fhPtEta ->Fill(pt ,eta);
3292 fhPtPhi ->Fill(pt ,phi);
3293 fhEtaPhi ->Fill(eta ,phi);
3295 fhPtCentrality ->Fill(pt,cen) ;
3296 fhPtEventPlane ->Fill(pt,ep ) ;
3300 Int_t tag = pi0->GetTag();
3301 Int_t mcIndex = GetMCIndex(tag);
3303 fhMCE [mcIndex] ->Fill(ener);
3304 fhMCPt [mcIndex] ->Fill(pt);
3305 fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3306 fhMCPtEta[mcIndex] ->Fill(pt,eta);
3308 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3310 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
3312 Float_t efracMC = 0;
3313 Int_t label = pi0->GetLabel();
3314 Int_t momlabel = -1;
3317 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3320 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3322 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3323 if(grandmom.E() > 0 && ok)
3325 efracMC = grandmom.E()/ener;
3326 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3329 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3331 fhMCPi0DecayPt->Fill(pt);
3332 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3333 if(grandmom.E() > 0 && ok)
3335 efracMC = mom.E()/grandmom.E();
3336 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3339 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3341 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3342 if(grandmom.E() > 0 && ok)
3344 efracMC = grandmom.E()/ener;
3345 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3348 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3350 fhMCEtaDecayPt->Fill(pt);
3351 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3352 if(grandmom.E() > 0 && ok)
3354 efracMC = mom.E()/grandmom.E();
3355 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3358 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3360 fhMCOtherDecayPt->Fill(pt);
3365 }//Histograms with MC
3371 //__________________________________________________________________
3372 void AliAnaPi0EbE::Print(const Option_t * opt) const
3374 //Print some relevant parameters set for the analysis
3378 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3379 AliAnaCaloTrackCorrBaseClass::Print("");
3380 printf("Analysis Type = %d \n", fAnaType) ;
3381 if(fAnaType == kSSCalo)
3383 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3384 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3385 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3386 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);