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),
64 fhPtEta(0), fhPtPhi(0), fhEtaPhi(0),
65 fhEtaPhiEMCALBC0(0), fhEtaPhiEMCALBC1(0), fhEtaPhiEMCALBCN(0),
66 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
67 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
68 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
69 fhPtCentrality(), fhPtEventPlane(0),
70 fhPtReject(0), fhEReject(0),
71 fhEEtaReject(0), fhEPhiReject(0), fhEtaPhiReject(0),
72 fhMass(0), fhMassPt(0), fhMassSplitPt(0),
73 fhSelectedMass(0), fhSelectedMassPt(0), fhSelectedMassSplitPt(0),
74 fhMassNoOverlap(0), fhMassPtNoOverlap(0), fhMassSplitPtNoOverlap(0),
75 fhSelectedMassNoOverlap(0), fhSelectedMassPtNoOverlap(0), fhSelectedMassSplitPtNoOverlap(0),
76 fhMCPi0PtRecoPtPrim(0), fhMCEtaPtRecoPtPrim(0),
77 fhMCPi0PtRecoPtPrimNoOverlap(0), fhMCEtaPtRecoPtPrimNoOverlap(0),
78 fhMCPi0SplitPtRecoPtPrim(0), fhMCEtaSplitPtRecoPtPrim(0),
79 fhMCPi0SplitPtRecoPtPrimNoOverlap(0), fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
80 fhMCPi0SelectedPtRecoPtPrim(0), fhMCEtaSelectedPtRecoPtPrim(0),
81 fhMCPi0SelectedPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
82 fhMCPi0SelectedSplitPtRecoPtPrim(0), fhMCEtaSelectedSplitPtRecoPtPrim(0),
83 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
84 fhAsymmetry(0), fhSelectedAsymmetry(0),
85 fhSplitE(0), fhSplitPt(0),
86 fhSplitPtEta(0), fhSplitPtPhi(0),
88 fhPtDecay(0), fhEDecay(0),
89 // Shower shape histos
90 fhEDispersion(0), fhELambda0(0), fhELambda1(0),
91 fhELambda0NoTRD(0), fhELambda0FracMaxCellCut(0),
92 fhEFracMaxCell(0), fhEFracMaxCellNoTRD(0),
93 fhENCells(0), fhETime(0), fhEPairDiffTime(0),
94 fhDispEtaE(0), fhDispPhiE(0),
95 fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
96 fhDispEtaPhiDiffE(0), fhSphericityE(0),
100 fhMCPhi(), fhMCEta(),
101 fhMCEReject(), fhMCPtReject(),
103 fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
104 fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
105 fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
107 fhMassPairMCPi0(0), fhMassPairMCEta(0),
108 fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
110 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
111 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
112 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
113 fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
114 fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
115 fhTrackMatchedMCParticleE(0),
116 fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
117 fhdEdx(0), fhEOverP(0), fhEOverPNoTRD(0),
118 // Number of local maxima in cluster
119 fhNLocMaxE(0), fhNLocMaxPt(0), fhNLocMaxPtReject(0),
121 fhTimePtNoCut(0), fhTimePtSPD(0), fhTimePtSPDMulti(0),
122 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
123 fhTimeNPileUpVertContributors(0),
124 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
125 fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
126 fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
127 fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
131 for(Int_t i = 0; i < 6; i++)
137 fhMCPtCentrality [i] = 0;
141 fhMCSplitPtPhi [i] = 0;
142 fhMCSplitPtEta [i] = 0;
144 fhMCNLocMaxPt [i] = 0;
145 fhMCNLocMaxSplitPt [i] = 0;
146 fhMCNLocMaxPtReject[i] = 0;
148 fhEMCLambda0 [i] = 0;
149 fhEMCLambda0NoTRD [i] = 0;
150 fhEMCLambda0FracMaxCellCut[i]= 0;
151 fhEMCFracMaxCell [i] = 0;
152 fhEMCLambda1 [i] = 0;
153 fhEMCDispersion [i] = 0;
155 fhMCEDispEta [i] = 0;
156 fhMCEDispPhi [i] = 0;
157 fhMCESumEtaPhi [i] = 0;
158 fhMCEDispEtaPhiDiff[i] = 0;
159 fhMCESphericity [i] = 0;
160 fhMCEAsymmetry [i] = 0;
163 fhMCMassSplitPt [i]=0;
164 fhMCSelectedMassPt [i]=0;
165 fhMCSelectedMassSplitPt[i]=0;
167 fhMCMassPtNoOverlap [i]=0;
168 fhMCMassSplitPtNoOverlap [i]=0;
169 fhMCSelectedMassPtNoOverlap [i]=0;
170 fhMCSelectedMassSplitPtNoOverlap[i]=0;
172 for(Int_t j = 0; j < 7; j++)
174 fhMCLambda0DispEta [j][i] = 0;
175 fhMCLambda0DispPhi [j][i] = 0;
176 fhMCDispEtaDispPhi [j][i] = 0;
177 fhMCAsymmetryLambda0 [j][i] = 0;
178 fhMCAsymmetryDispEta [j][i] = 0;
179 fhMCAsymmetryDispPhi [j][i] = 0;
183 for(Int_t j = 0; j < 7; j++)
185 fhLambda0DispEta [j] = 0;
186 fhLambda0DispPhi [j] = 0;
187 fhDispEtaDispPhi [j] = 0;
188 fhAsymmetryLambda0 [j] = 0;
189 fhAsymmetryDispEta [j] = 0;
190 fhAsymmetryDispPhi [j] = 0;
195 for(Int_t i = 0; i < 3; i++)
197 fhELambda0LocMax [i] = 0;
198 fhELambda1LocMax [i] = 0;
199 fhEDispersionLocMax [i] = 0;
200 fhEDispEtaLocMax [i] = 0;
201 fhEDispPhiLocMax [i] = 0;
202 fhESumEtaPhiLocMax [i] = 0;
203 fhEDispEtaPhiDiffLocMax[i] = 0;
204 fhESphericityLocMax [i] = 0;
205 fhEAsymmetryLocMax [i] = 0;
209 for(Int_t i =0; i < 14; i++){
210 fhLambda0ForW0[i] = 0;
211 //fhLambda1ForW0[i] = 0;
212 if(i<8)fhMassPairLocMax[i] = 0;
215 for(Int_t i = 0; i < 11; i++)
217 fhEtaPhiTriggerEMCALBC [i] = 0 ;
218 fhTimeTriggerEMCALBC [i] = 0 ;
219 fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
221 fhEtaPhiTriggerEMCALBCUM [i] = 0 ;
222 fhTimeTriggerEMCALBCUM [i] = 0 ;
226 //Initialize parameters
231 //___________________________________________________________________________________
232 void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
234 // Fill some histograms to understand pile-up
235 if(!fFillPileUpHistograms) return;
237 //printf("E %f, time %f\n",energy,time);
238 AliVEvent * event = GetReader()->GetInputEvent();
240 fhTimePtNoCut->Fill(pt,time);
241 if(GetReader()->IsPileUpFromSPD())
243 if(GetReader()->IsPileUpFromSPD()) { fhPtPileUp[0]->Fill(pt); fhTimePtSPD ->Fill(pt,time); }
244 if(GetReader()->IsPileUpFromEMCal()) fhPtPileUp[1]->Fill(pt);
245 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPileUp[2]->Fill(pt);
246 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPileUp[3]->Fill(pt);
247 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPileUp[4]->Fill(pt);
248 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPileUp[5]->Fill(pt);
249 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
251 if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
255 AliVCaloCells* cells = 0;
256 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
257 else cells = GetPHOSCells();
259 Float_t maxCellFraction = 0.;
260 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
262 Double_t tmax = cells->GetCellTime(absIdMax);
263 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
266 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
267 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
269 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
271 Int_t absId = calo->GetCellsAbsId()[ipos];
273 if( absId == absIdMax ) continue ;
275 Double_t timecell = cells->GetCellTime(absId);
276 Float_t amp = cells->GetCellAmplitude(absId);
277 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
278 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
281 Float_t diff = (tmax-timecell);
283 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
285 if(GetReader()->IsPileUpFromSPD())
287 fhPtCellTimePileUp[0]->Fill(pt, timecell);
288 fhPtTimeDiffPileUp[0]->Fill(pt, diff);
291 if(GetReader()->IsPileUpFromEMCal())
293 fhPtCellTimePileUp[1]->Fill(pt, timecell);
294 fhPtTimeDiffPileUp[1]->Fill(pt, diff);
297 if(GetReader()->IsPileUpFromSPDOrEMCal())
299 fhPtCellTimePileUp[2]->Fill(pt, timecell);
300 fhPtTimeDiffPileUp[2]->Fill(pt, diff);
303 if(GetReader()->IsPileUpFromSPDAndEMCal())
305 fhPtCellTimePileUp[3]->Fill(pt, timecell);
306 fhPtTimeDiffPileUp[3]->Fill(pt, diff);
309 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
311 fhPtCellTimePileUp[4]->Fill(pt, timecell);
312 fhPtTimeDiffPileUp[4]->Fill(pt, diff);
315 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
317 fhPtCellTimePileUp[5]->Fill(pt, timecell);
318 fhPtTimeDiffPileUp[5]->Fill(pt, diff);
321 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
323 fhPtCellTimePileUp[6]->Fill(pt, timecell);
324 fhPtTimeDiffPileUp[6]->Fill(pt, diff);
329 if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
331 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
332 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
334 // N pile up vertices
340 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
341 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
346 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
347 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
350 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
351 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
353 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
354 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
356 if(TMath::Abs(time) < 25)
358 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
359 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
362 if(time < 75 && time > -25)
364 fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
365 fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
368 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
369 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
372 Float_t z1 = -1, z2 = -1;
374 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
378 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
379 ncont=pv->GetNContributors();
380 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
382 diamZ = esdEv->GetDiamondZ();
386 AliAODVertex *pv=aodEv->GetVertex(iVert);
387 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
388 ncont=pv->GetNContributors();
389 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
391 diamZ = aodEv->GetDiamondZ();
394 Double_t distZ = TMath::Abs(z2-z1);
395 diamZ = TMath::Abs(z2-diamZ);
397 fhTimeNPileUpVertContributors ->Fill(time,ncont);
398 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
399 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
405 //______________________________________________________________________________________________
406 void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
408 // Fill histograms that do not pass the identification (SS case only)
410 Float_t ener = mom.E();
411 Float_t pt = mom.Pt();
412 Float_t phi = mom.Phi();
413 if(phi < 0) phi+=TMath::TwoPi();
414 Float_t eta = mom.Eta();
416 fhPtReject ->Fill(pt);
417 fhEReject ->Fill(ener);
419 fhEEtaReject ->Fill(ener,eta);
420 fhEPhiReject ->Fill(ener,phi);
421 fhEtaPhiReject ->Fill(eta,phi);
423 fhNLocMaxPtReject->Fill(pt,nMaxima);
427 Int_t mcIndex = GetMCIndex(mctag);
428 fhMCEReject [mcIndex] ->Fill(ener);
429 fhMCPtReject [mcIndex] ->Fill(pt);
430 fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
434 //___________________________________________________________________________________
435 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Int_t nMaxima,
436 Int_t tag, Float_t asy)
438 // Fill shower shape, timing and other histograms for selected clusters from decay
440 Float_t e = cluster->E();
441 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
442 Float_t l0 = cluster->GetM02();
443 Float_t l1 = cluster->GetM20();
444 Int_t nSM = GetModuleNumber(cluster);
447 if (e < 2 ) ebin = 0;
448 else if (e < 4 ) ebin = 1;
449 else if (e < 6 ) ebin = 2;
450 else if (e < 10) ebin = 3;
451 else if (e < 15) ebin = 4;
452 else if (e < 20) ebin = 5;
456 if (nMaxima==1) indexMax = 0 ;
457 else if(nMaxima==2) indexMax = 1 ;
461 AliVCaloCells * cell = 0x0;
462 if(fCalorimeter == "PHOS")
463 cell = GetPHOSCells();
465 cell = GetEMCALCells();
467 Float_t maxCellFraction = 0;
468 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
469 fhEFracMaxCell->Fill(e,maxCellFraction);
471 FillWeightHistograms(cluster);
473 fhEDispersion->Fill(e, disp);
474 fhELambda0 ->Fill(e, l0 );
475 fhELambda1 ->Fill(e, l1 );
477 Float_t ll0 = 0., ll1 = 0.;
478 Float_t dispp= 0., dEta = 0., dPhi = 0.;
479 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
480 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
482 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
483 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
485 fhDispEtaE -> Fill(e,dEta);
486 fhDispPhiE -> Fill(e,dPhi);
487 fhSumEtaE -> Fill(e,sEta);
488 fhSumPhiE -> Fill(e,sPhi);
489 fhSumEtaPhiE -> Fill(e,sEtaPhi);
490 fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
491 if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
493 fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
494 fhLambda0DispEta[ebin]->Fill(l0 ,dEta);
495 fhLambda0DispPhi[ebin]->Fill(l0 ,dPhi);
497 if (fAnaType==kSSCalo)
499 // Asymmetry histograms
500 fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
501 fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
502 fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
506 fhNLocMaxE ->Fill(e ,nMaxima);
508 fhELambda0LocMax [indexMax]->Fill(e,l0);
509 fhELambda1LocMax [indexMax]->Fill(e,l1);
510 fhEDispersionLocMax[indexMax]->Fill(e,disp);
512 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
514 fhEDispEtaLocMax [indexMax]-> Fill(e,dEta);
515 fhEDispPhiLocMax [indexMax]-> Fill(e,dPhi);
516 fhESumEtaPhiLocMax [indexMax]-> Fill(e,sEtaPhi);
517 fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
518 if(dEta+dPhi>0) fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
519 if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e ,asy);
523 if(fCalorimeter=="EMCAL" && nSM < 6)
525 fhELambda0NoTRD->Fill(e, l0 );
526 fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
529 if(maxCellFraction < 0.5)
530 fhELambda0FracMaxCellCut->Fill(e, l0 );
532 fhETime ->Fill(e, cluster->GetTOF()*1.e9);
533 fhENCells->Fill(e, cluster->GetNCells());
535 // Fill Track matching control histograms
538 Float_t dZ = cluster->GetTrackDz();
539 Float_t dR = cluster->GetTrackDx();
541 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
543 dR = 2000., dZ = 2000.;
544 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
546 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
548 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
550 Bool_t positive = kFALSE;
551 if(track) positive = (track->Charge()>0);
553 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
555 fhTrackMatchedDEta->Fill(e,dZ);
556 fhTrackMatchedDPhi->Fill(e,dR);
557 if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
563 fhTrackMatchedDEtaPos->Fill(cluster->E(),dZ);
564 fhTrackMatchedDPhiPos->Fill(cluster->E(),dR);
565 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
569 fhTrackMatchedDEtaNeg->Fill(cluster->E(),dZ);
570 fhTrackMatchedDPhiNeg->Fill(cluster->E(),dR);
571 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
575 // Check dEdx and E/p of matched clusters
577 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
581 Float_t dEdx = track->GetTPCsignal();
582 fhdEdx->Fill(e, dEdx);
584 Float_t eOverp = e/track->P();
585 fhEOverP->Fill(e, eOverp);
587 if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e, eOverp);
591 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
596 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
598 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
599 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 2.5 ;
600 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 0.5 ;
601 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 1.5 ;
607 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
608 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 6.5 ;
609 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 4.5 ;
610 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 5.5 ;
614 fhTrackMatchedMCParticleE ->Fill(e , mctag);
615 fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
616 fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
620 }// Track matching histograms
624 Int_t mcIndex = GetMCIndex(tag);
626 fhEMCLambda0[mcIndex] ->Fill(e, l0);
627 fhEMCLambda1[mcIndex] ->Fill(e, l1);
628 fhEMCDispersion[mcIndex] ->Fill(e, disp);
629 fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
631 if(fCalorimeter=="EMCAL" && nSM < 6)
632 fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0 );
634 if(maxCellFraction < 0.5)
635 fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0 );
637 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
639 fhMCEDispEta [mcIndex]-> Fill(e,dEta);
640 fhMCEDispPhi [mcIndex]-> Fill(e,dPhi);
641 fhMCESumEtaPhi [mcIndex]-> Fill(e,sEtaPhi);
642 fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
643 if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
645 if (fAnaType==kSSCalo)
647 fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
648 fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
649 fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
652 fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
653 fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0 ,dEta);
654 fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0 ,dPhi);
662 //________________________________________________________
663 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
665 // Calculate weights and fill histograms
667 if(!fFillWeightHistograms || GetMixedEvent()) return;
669 AliVCaloCells* cells = 0;
670 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
671 else cells = GetPHOSCells();
673 // First recalculate energy in case non linearity was applied
676 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
679 Int_t id = clus->GetCellsAbsId()[ipos];
681 //Recalibrate cell energy if needed
682 Float_t amp = cells->GetCellAmplitude(id);
683 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
694 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
698 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
699 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
701 //Get the ratio and log ratio to all cells in cluster
702 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
704 Int_t id = clus->GetCellsAbsId()[ipos];
706 //Recalibrate cell energy if needed
707 Float_t amp = cells->GetCellAmplitude(id);
708 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
710 fhECellClusterRatio ->Fill(energy,amp/energy);
711 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
714 //Recalculate shower shape for different W0
715 if(fCalorimeter=="EMCAL"){
717 Float_t l0org = clus->GetM02();
718 Float_t l1org = clus->GetM20();
719 Float_t dorg = clus->GetDispersion();
721 for(Int_t iw = 0; iw < 14; iw++)
723 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
724 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
726 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
727 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
731 // Set the original values back
734 clus->SetDispersion(dorg);
739 //__________________________________________
740 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
742 //Save parameters used for analysis
743 TString parList ; //this will be list of parameters used for this analysis.
744 const Int_t buffersize = 255;
745 char onePar[buffersize] ;
747 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
749 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
752 if(fAnaType == kSSCalo)
754 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
756 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
758 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
760 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
764 //Get parameters set in base class.
765 parList += GetBaseParametersList() ;
767 //Get parameters set in PID class.
768 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
770 return new TObjString(parList) ;
773 //_____________________________________________
774 TList * AliAnaPi0EbE::GetCreateOutputObjects()
776 // Create histograms to be saved in output file and
777 // store them in outputContainer
778 TList * outputContainer = new TList() ;
779 outputContainer->SetName("Pi0EbEHistos") ;
781 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
782 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
783 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
784 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
785 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
786 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
787 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
789 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
790 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
791 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
793 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
794 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
795 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
796 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
797 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
798 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
800 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
801 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
802 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
803 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
804 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
805 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
807 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();
808 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
809 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
811 TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
812 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
813 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
814 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
816 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
817 fhPt->SetYTitle("N");
818 fhPt->SetXTitle("p_{T} (GeV/c)");
819 outputContainer->Add(fhPt) ;
821 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
823 fhE->SetXTitle("E (GeV)");
824 outputContainer->Add(fhE) ;
827 ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
828 fhEPhi->SetYTitle("#phi (rad)");
829 fhEPhi->SetXTitle("E (GeV)");
830 outputContainer->Add(fhEPhi) ;
833 ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
834 fhEEta->SetYTitle("#eta");
835 fhEEta->SetXTitle("E (GeV)");
836 outputContainer->Add(fhEEta) ;
839 ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
840 fhPtPhi->SetYTitle("#phi (rad)");
841 fhPtPhi->SetXTitle("p_{T} (GeV/c)");
842 outputContainer->Add(fhPtPhi) ;
845 ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
846 fhPtEta->SetYTitle("#eta");
847 fhPtEta->SetXTitle("p_{T} (GeV/c)");
848 outputContainer->Add(fhPtEta) ;
851 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
852 fhEtaPhi->SetYTitle("#phi (rad)");
853 fhEtaPhi->SetXTitle("#eta");
854 outputContainer->Add(fhEtaPhi) ;
856 if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
858 fhEtaPhiEMCALBC0 = new TH2F
859 ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
860 fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
861 fhEtaPhiEMCALBC0->SetXTitle("#eta");
862 outputContainer->Add(fhEtaPhiEMCALBC0) ;
864 fhEtaPhiEMCALBC1 = new TH2F
865 ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
866 fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
867 fhEtaPhiEMCALBC1->SetXTitle("#eta");
868 outputContainer->Add(fhEtaPhiEMCALBC1) ;
870 fhEtaPhiEMCALBCN = new TH2F
871 ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
872 fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
873 fhEtaPhiEMCALBCN->SetXTitle("#eta");
874 outputContainer->Add(fhEtaPhiEMCALBCN) ;
876 for(Int_t i = 0; i < 11; i++)
878 fhEtaPhiTriggerEMCALBC[i] = new TH2F
879 (Form("hEtaPhiTriggerEMCALBC%d",i-5),
880 Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
881 netabins,etamin,etamax,nphibins,phimin,phimax);
882 fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
883 fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
884 outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
886 fhTimeTriggerEMCALBC[i] = new TH2F
887 (Form("hTimeTriggerEMCALBC%d",i-5),
888 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
889 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
890 fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
891 fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
892 outputContainer->Add(fhTimeTriggerEMCALBC[i]);
894 fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
895 (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
896 Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
897 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
898 fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
899 fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
900 outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
902 fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
903 (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
904 Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
905 netabins,etamin,etamax,nphibins,phimin,phimax);
906 fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
907 fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
908 outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
910 fhTimeTriggerEMCALBCUM[i] = new TH2F
911 (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
912 Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
913 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
914 fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
915 fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
916 outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
920 fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
921 "cluster time vs E of clusters, no match, rematch open time",
922 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
923 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("E (GeV)");
924 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("time (ns)");
925 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
928 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
929 "cluster time vs E of clusters, no match, rematch with neigbour parches",
930 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
931 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("E (GeV)");
932 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("time (ns)");
933 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
935 fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
936 "cluster time vs E of clusters, no match, rematch open time and neigbour",
937 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
938 fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("E (GeV)");
939 fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("time (ns)");
940 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
944 fhPtCentrality = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
945 fhPtCentrality->SetYTitle("centrality");
946 fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
947 outputContainer->Add(fhPtCentrality) ;
949 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
950 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
951 fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
952 outputContainer->Add(fhPtEventPlane) ;
954 if(fAnaType == kSSCalo)
956 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
957 fhPtReject->SetYTitle("N");
958 fhPtReject->SetXTitle("p_{T} (GeV/c)");
959 outputContainer->Add(fhPtReject) ;
961 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
962 fhEReject->SetYTitle("N");
963 fhEReject->SetXTitle("E (GeV)");
964 outputContainer->Add(fhEReject) ;
966 fhEPhiReject = new TH2F
967 ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
968 fhEPhiReject->SetYTitle("#phi (rad)");
969 fhEPhiReject->SetXTitle("E (GeV)");
970 outputContainer->Add(fhEPhiReject) ;
972 fhEEtaReject = new TH2F
973 ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
974 fhEEtaReject->SetYTitle("#eta");
975 fhEEtaReject->SetXTitle("E (GeV)");
976 outputContainer->Add(fhEEtaReject) ;
978 fhEtaPhiReject = new TH2F
979 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
980 fhEtaPhiReject->SetYTitle("#phi (rad)");
981 fhEtaPhiReject->SetXTitle("#eta");
982 outputContainer->Add(fhEtaPhiReject) ;
986 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
987 fhMass->SetYTitle("mass (GeV/c^{2})");
988 fhMass->SetXTitle("E (GeV)");
989 outputContainer->Add(fhMass) ;
991 fhSelectedMass = new TH2F
992 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
993 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
994 fhSelectedMass->SetXTitle("E (GeV)");
995 outputContainer->Add(fhSelectedMass) ;
998 ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
999 fhMassPt->SetYTitle("mass (GeV/c^{2})");
1000 fhMassPt->SetXTitle("p_{T} (GeV/c)");
1001 outputContainer->Add(fhMassPt) ;
1003 fhSelectedMassPt = new TH2F
1004 ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1005 fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
1006 fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
1007 outputContainer->Add(fhSelectedMassPt) ;
1009 if(IsDataMC() && fAnaType == kSSCalo)
1011 fhMassNoOverlap = new TH2F
1012 ("hMassNoOverlap","all pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1013 fhMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1014 fhMassNoOverlap->SetXTitle("E (GeV)");
1015 outputContainer->Add(fhMassNoOverlap) ;
1017 fhSelectedMassNoOverlap = new TH2F
1018 ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs mass: E vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1019 fhSelectedMassNoOverlap->SetYTitle("mass (GeV/c^{2})");
1020 fhSelectedMassNoOverlap->SetXTitle("E (GeV)");
1021 outputContainer->Add(fhSelectedMassNoOverlap) ;
1023 fhMassPtNoOverlap = new TH2F
1024 ("hMassPtNoOverlap","all pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1025 fhMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1026 fhMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1027 outputContainer->Add(fhMassPtNoOverlap) ;
1029 fhSelectedMassPtNoOverlap = new TH2F
1030 ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1031 fhSelectedMassPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1032 fhSelectedMassPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1033 outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1036 if(fAnaType != kSSCalo)
1038 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1039 fhPtDecay->SetYTitle("N");
1040 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
1041 outputContainer->Add(fhPtDecay) ;
1043 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1044 fhEDecay->SetYTitle("N");
1045 fhEDecay->SetXTitle("E (GeV)");
1046 outputContainer->Add(fhEDecay) ;
1051 if( fFillSelectClHisto )
1054 fhEDispersion = new TH2F
1055 ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1056 fhEDispersion->SetYTitle("D^{2}");
1057 fhEDispersion->SetXTitle("E (GeV)");
1058 outputContainer->Add(fhEDispersion) ;
1060 fhELambda0 = new TH2F
1061 ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1062 fhELambda0->SetYTitle("#lambda_{0}^{2}");
1063 fhELambda0->SetXTitle("E (GeV)");
1064 outputContainer->Add(fhELambda0) ;
1066 fhELambda1 = new TH2F
1067 ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1068 fhELambda1->SetYTitle("#lambda_{1}^{2}");
1069 fhELambda1->SetXTitle("E (GeV)");
1070 outputContainer->Add(fhELambda1) ;
1072 fhELambda0FracMaxCellCut = new TH2F
1073 ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1074 fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1075 fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
1076 outputContainer->Add(fhELambda0FracMaxCellCut) ;
1078 fhEFracMaxCell = new TH2F
1079 ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1080 fhEFracMaxCell->SetYTitle("Fraction");
1081 fhEFracMaxCell->SetXTitle("E (GeV)");
1082 outputContainer->Add(fhEFracMaxCell) ;
1084 if(fCalorimeter=="EMCAL")
1086 fhELambda0NoTRD = new TH2F
1087 ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1088 fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1089 fhELambda0NoTRD->SetXTitle("E (GeV)");
1090 outputContainer->Add(fhELambda0NoTRD) ;
1092 fhEFracMaxCellNoTRD = new TH2F
1093 ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1094 fhEFracMaxCellNoTRD->SetYTitle("Fraction");
1095 fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
1096 outputContainer->Add(fhEFracMaxCellNoTRD) ;
1098 if(!fFillOnlySimpleSSHisto)
1100 fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1101 fhDispEtaE->SetXTitle("E (GeV)");
1102 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
1103 outputContainer->Add(fhDispEtaE);
1105 fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1106 fhDispPhiE->SetXTitle("E (GeV)");
1107 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
1108 outputContainer->Add(fhDispPhiE);
1110 fhSumEtaE = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1111 fhSumEtaE->SetXTitle("E (GeV)");
1112 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
1113 outputContainer->Add(fhSumEtaE);
1115 fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
1116 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1117 fhSumPhiE->SetXTitle("E (GeV)");
1118 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
1119 outputContainer->Add(fhSumPhiE);
1121 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
1122 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1123 fhSumEtaPhiE->SetXTitle("E (GeV)");
1124 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
1125 outputContainer->Add(fhSumEtaPhiE);
1127 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
1128 nptbins,ptmin,ptmax,200, -10,10);
1129 fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
1130 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1131 outputContainer->Add(fhDispEtaPhiDiffE);
1133 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
1134 nptbins,ptmin,ptmax, 200, -1,1);
1135 fhSphericityE->SetXTitle("E (GeV)");
1136 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1137 outputContainer->Add(fhSphericityE);
1139 for(Int_t i = 0; i < 7; i++)
1141 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]),
1142 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1143 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1144 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1145 outputContainer->Add(fhDispEtaDispPhi[i]);
1147 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]),
1148 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1149 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1150 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1151 outputContainer->Add(fhLambda0DispEta[i]);
1153 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]),
1154 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1155 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1156 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1157 outputContainer->Add(fhLambda0DispPhi[i]);
1163 fhNLocMaxE = new TH2F("hNLocMaxE","Number of local maxima in cluster",
1164 nptbins,ptmin,ptmax,10,0,10);
1165 fhNLocMaxE ->SetYTitle("N maxima");
1166 fhNLocMaxE ->SetXTitle("E (GeV)");
1167 outputContainer->Add(fhNLocMaxE) ;
1169 if(fAnaType == kSSCalo)
1171 fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1172 nptbins,ptmin,ptmax,10,0,10);
1173 fhNLocMaxPt ->SetYTitle("N maxima");
1174 fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
1175 outputContainer->Add(fhNLocMaxPt) ;
1177 fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1178 nptbins,ptmin,ptmax,10,0,10);
1179 fhNLocMaxPtReject ->SetYTitle("N maxima");
1180 fhNLocMaxPtReject ->SetXTitle("p_{T} (GeV/c)");
1181 outputContainer->Add(fhNLocMaxPtReject) ;
1184 for (Int_t i = 0; i < 3; i++)
1186 fhELambda0LocMax[i] = new TH2F(Form("hELambda0LocMax%d",i+1),
1187 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
1188 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1189 fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1190 fhELambda0LocMax[i]->SetXTitle("E (GeV)");
1191 outputContainer->Add(fhELambda0LocMax[i]) ;
1193 fhELambda1LocMax[i] = new TH2F(Form("hELambda1LocMax%d",i+1),
1194 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
1195 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1196 fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1197 fhELambda1LocMax[i]->SetXTitle("E (GeV)");
1198 outputContainer->Add(fhELambda1LocMax[i]) ;
1200 fhEDispersionLocMax[i] = new TH2F(Form("hEDispersionLocMax%d",i+1),
1201 Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
1202 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1203 fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1204 fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
1205 outputContainer->Add(fhEDispersionLocMax[i]) ;
1207 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1209 fhEDispEtaLocMax[i] = new TH2F(Form("hEDispEtaLocMax%d",i+1),
1210 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1211 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1212 fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1213 fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
1214 outputContainer->Add(fhEDispEtaLocMax[i]) ;
1216 fhEDispPhiLocMax[i] = new TH2F(Form("hEDispPhiLocMax%d",i+1),
1217 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1218 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1219 fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1220 fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
1221 outputContainer->Add(fhEDispPhiLocMax[i]) ;
1223 fhESumEtaPhiLocMax[i] = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
1224 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1225 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1226 fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1227 fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
1228 outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
1230 fhEDispEtaPhiDiffLocMax[i] = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
1231 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1232 nptbins,ptmin,ptmax,200, -10,10);
1233 fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1234 fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
1235 outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
1237 fhESphericityLocMax[i] = new TH2F(Form("hESphericityLocMax%d",i+1),
1238 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1239 nptbins,ptmin,ptmax,200, -1,1);
1240 fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1241 fhESphericityLocMax[i]->SetXTitle("E (GeV)");
1242 outputContainer->Add(fhESphericityLocMax[i]) ;
1247 fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1248 fhENCells->SetXTitle("E (GeV)");
1249 fhENCells->SetYTitle("# of cells in cluster");
1250 outputContainer->Add(fhENCells);
1252 fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1253 fhETime->SetXTitle("E (GeV)");
1254 fhETime->SetYTitle("t (ns)");
1255 outputContainer->Add(fhETime);
1260 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1261 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
1262 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1263 outputContainer->Add(fhEPairDiffTime);
1265 if(fAnaType == kIMCalo)
1267 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1268 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1269 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1270 "2 Local Maxima paired with more than 2 Local Maxima",
1271 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1273 for (Int_t i = 0; i < 8 ; i++)
1276 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1278 fhMassPairLocMax[i] = new TH2F
1279 (Form("MassPairLocMax%s",combiName[i].Data()),
1280 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1281 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1282 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
1283 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
1284 outputContainer->Add(fhMassPairLocMax[i]) ;
1290 fhTrackMatchedDEta = new TH2F
1291 ("hTrackMatchedDEta",
1292 "d#eta of cluster-track vs cluster energy",
1293 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1294 fhTrackMatchedDEta->SetYTitle("d#eta");
1295 fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
1297 fhTrackMatchedDPhi = new TH2F
1298 ("hTrackMatchedDPhi",
1299 "d#phi of cluster-track vs cluster energy",
1300 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1301 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1302 fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
1304 fhTrackMatchedDEtaDPhi = new TH2F
1305 ("hTrackMatchedDEtaDPhi",
1306 "d#eta vs d#phi of cluster-track vs cluster energy",
1307 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1308 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1309 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1311 outputContainer->Add(fhTrackMatchedDEta) ;
1312 outputContainer->Add(fhTrackMatchedDPhi) ;
1313 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1315 fhTrackMatchedDEtaPos = new TH2F
1316 ("hTrackMatchedDEtaPos",
1317 "d#eta of cluster-track vs cluster energy",
1318 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1319 fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1320 fhTrackMatchedDEtaPos->SetXTitle("E_{cluster} (GeV)");
1322 fhTrackMatchedDPhiPos = new TH2F
1323 ("hTrackMatchedDPhiPos",
1324 "d#phi of cluster-track vs cluster energy",
1325 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1326 fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1327 fhTrackMatchedDPhiPos->SetXTitle("E_{cluster} (GeV)");
1329 fhTrackMatchedDEtaDPhiPos = new TH2F
1330 ("hTrackMatchedDEtaDPhiPos",
1331 "d#eta vs d#phi of cluster-track vs cluster energy",
1332 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1333 fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1334 fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1336 outputContainer->Add(fhTrackMatchedDEtaPos) ;
1337 outputContainer->Add(fhTrackMatchedDPhiPos) ;
1338 outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1340 fhTrackMatchedDEtaNeg = new TH2F
1341 ("hTrackMatchedDEtaNeg",
1342 "d#eta of cluster-track vs cluster energy",
1343 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1344 fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1345 fhTrackMatchedDEtaNeg->SetXTitle("E_{cluster} (GeV)");
1347 fhTrackMatchedDPhiNeg = new TH2F
1348 ("hTrackMatchedDPhiNeg",
1349 "d#phi of cluster-track vs cluster energy",
1350 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1351 fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1352 fhTrackMatchedDPhiNeg->SetXTitle("E_{cluster} (GeV)");
1354 fhTrackMatchedDEtaDPhiNeg = new TH2F
1355 ("hTrackMatchedDEtaDPhiNeg",
1356 "d#eta vs d#phi of cluster-track vs cluster energy",
1357 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1358 fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1359 fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1361 outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1362 outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1363 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1365 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1366 fhdEdx->SetXTitle("E (GeV)");
1367 fhdEdx->SetYTitle("<dE/dx>");
1368 outputContainer->Add(fhdEdx);
1370 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1371 fhEOverP->SetXTitle("E (GeV)");
1372 fhEOverP->SetYTitle("E/p");
1373 outputContainer->Add(fhEOverP);
1375 if(fCalorimeter=="EMCAL")
1377 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1378 fhEOverPNoTRD->SetXTitle("E (GeV)");
1379 fhEOverPNoTRD->SetYTitle("E/p");
1380 outputContainer->Add(fhEOverPNoTRD);
1383 if(IsDataMC() && fFillTMHisto)
1385 fhTrackMatchedMCParticleE = new TH2F
1386 ("hTrackMatchedMCParticleE",
1387 "Origin of particle vs energy",
1388 nptbins,ptmin,ptmax,8,0,8);
1389 fhTrackMatchedMCParticleE->SetXTitle("E (GeV)");
1390 //fhTrackMatchedMCParticleE->SetYTitle("Particle type");
1392 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(1 ,"Photon");
1393 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(2 ,"Electron");
1394 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1395 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(4 ,"Rest");
1396 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1397 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1398 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1399 fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1401 outputContainer->Add(fhTrackMatchedMCParticleE);
1403 fhTrackMatchedMCParticleDEta = new TH2F
1404 ("hTrackMatchedMCParticleDEta",
1405 "Origin of particle vs #eta residual",
1406 nresetabins,resetamin,resetamax,8,0,8);
1407 fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1408 //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1410 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1411 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1412 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1413 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1414 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1415 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1416 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1417 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1419 outputContainer->Add(fhTrackMatchedMCParticleDEta);
1421 fhTrackMatchedMCParticleDPhi = new TH2F
1422 ("hTrackMatchedMCParticleDPhi",
1423 "Origin of particle vs #phi residual",
1424 nresphibins,resphimin,resphimax,8,0,8);
1425 fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1426 //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1428 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1429 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1430 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1431 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1432 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1433 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1434 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1435 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1437 outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1443 if(fFillWeightHistograms)
1445 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1446 nptbins,ptmin,ptmax, 100,0,1.);
1447 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1448 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1449 outputContainer->Add(fhECellClusterRatio);
1451 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1452 nptbins,ptmin,ptmax, 100,-10,0);
1453 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1454 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1455 outputContainer->Add(fhECellClusterLogRatio);
1457 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1458 nptbins,ptmin,ptmax, 100,0,1.);
1459 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1460 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1461 outputContainer->Add(fhEMaxCellClusterRatio);
1463 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1464 nptbins,ptmin,ptmax, 100,-10,0);
1465 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1466 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1467 outputContainer->Add(fhEMaxCellClusterLogRatio);
1469 for(Int_t iw = 0; iw < 14; iw++)
1471 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),
1472 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1473 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1474 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1475 outputContainer->Add(fhLambda0ForW0[iw]);
1477 // 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),
1478 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1479 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1480 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1481 // outputContainer->Add(fhLambda1ForW0[iw]);
1488 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1490 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",
1491 nptbins,ptmin,ptmax,200,0,2);
1492 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1493 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1494 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1496 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1497 nptbins,ptmin,ptmax,200,0,2);
1498 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1499 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1500 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1502 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1503 fhMCPi0DecayPt->SetYTitle("N");
1504 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1505 outputContainer->Add(fhMCPi0DecayPt) ;
1507 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}",
1508 nptbins,ptmin,ptmax,100,0,1);
1509 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1510 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1511 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1513 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1514 fhMCEtaDecayPt->SetYTitle("N");
1515 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1516 outputContainer->Add(fhMCEtaDecayPt) ;
1518 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
1519 nptbins,ptmin,ptmax,100,0,1);
1520 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1521 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1522 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1524 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1525 fhMCOtherDecayPt->SetYTitle("N");
1526 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1527 outputContainer->Add(fhMCOtherDecayPt) ;
1531 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1532 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1535 fhAnglePairMCPi0 = new TH2F
1537 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1538 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1539 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1540 outputContainer->Add(fhAnglePairMCPi0) ;
1542 if (fAnaType!= kSSCalo)
1544 fhAnglePairMCEta = new TH2F
1546 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1547 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1548 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1549 outputContainer->Add(fhAnglePairMCEta) ;
1551 fhMassPairMCPi0 = new TH2F
1553 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1554 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1555 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1556 outputContainer->Add(fhMassPairMCPi0) ;
1558 fhMassPairMCEta = new TH2F
1560 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1561 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1562 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1563 outputContainer->Add(fhMassPairMCEta) ;
1566 for(Int_t i = 0; i < 6; i++)
1570 (Form("hE_MC%s",pname[i].Data()),
1571 Form("Identified as #pi^{0} (#eta), cluster from %s",
1573 nptbins,ptmin,ptmax);
1574 fhMCE[i]->SetYTitle("N");
1575 fhMCE[i]->SetXTitle("E (GeV)");
1576 outputContainer->Add(fhMCE[i]) ;
1578 fhMCPt[i] = new TH1F
1579 (Form("hPt_MC%s",pname[i].Data()),
1580 Form("Identified as #pi^{0} (#eta), cluster from %s",
1582 nptbins,ptmin,ptmax);
1583 fhMCPt[i]->SetYTitle("N");
1584 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1585 outputContainer->Add(fhMCPt[i]) ;
1587 fhMCPtCentrality[i] = new TH2F
1588 (Form("hPtCentrality_MC%s",pname[i].Data()),
1589 Form("Identified as #pi^{0} (#eta), cluster from %s",
1591 nptbins,ptmin,ptmax, 100,0,100);
1592 fhMCPtCentrality[i]->SetYTitle("centrality");
1593 fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1594 outputContainer->Add(fhMCPtCentrality[i]) ;
1596 if(fAnaType == kSSCalo)
1598 fhMCNLocMaxPt[i] = new TH2F
1599 (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1600 Form("cluster from %s, pT of cluster vs NLM, accepted",ptype[i].Data()),
1601 nptbins,ptmin,ptmax,10,0,10);
1602 fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
1603 fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
1604 outputContainer->Add(fhMCNLocMaxPt[i]) ;
1606 fhMCNLocMaxPtReject[i] = new TH2F
1607 (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1608 Form("cluster from %s, pT of cluster vs NLM, rejected",ptype[i].Data()),
1609 nptbins,ptmin,ptmax,10,0,10);
1610 fhMCNLocMaxPtReject[i] ->SetYTitle("N maxima");
1611 fhMCNLocMaxPtReject[i] ->SetXTitle("p_{T} (GeV/c)");
1612 outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1614 fhMCEReject[i] = new TH1F
1615 (Form("hEReject_MC%s",pname[i].Data()),
1616 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1618 nptbins,ptmin,ptmax);
1619 fhMCEReject[i]->SetYTitle("N");
1620 fhMCEReject[i]->SetXTitle("E (GeV)");
1621 outputContainer->Add(fhMCEReject[i]) ;
1623 fhMCPtReject[i] = new TH1F
1624 (Form("hPtReject_MC%s",pname[i].Data()),
1625 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1627 nptbins,ptmin,ptmax);
1628 fhMCPtReject[i]->SetYTitle("N");
1629 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1630 outputContainer->Add(fhMCPtReject[i]) ;
1633 fhMCPhi[i] = new TH2F
1634 (Form("hPhi_MC%s",pname[i].Data()),
1635 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1636 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1637 fhMCPhi[i]->SetYTitle("#phi");
1638 fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1639 outputContainer->Add(fhMCPhi[i]) ;
1641 fhMCEta[i] = new TH2F
1642 (Form("hEta_MC%s",pname[i].Data()),
1643 Form("Identified as #pi^{0} (#eta), cluster from %s",
1644 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1645 fhMCEta[i]->SetYTitle("#eta");
1646 fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1647 outputContainer->Add(fhMCEta[i]) ;
1649 fhMCMassPt[i] = new TH2F
1650 (Form("hMassPt_MC%s",pname[i].Data()),
1651 Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1652 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1653 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1654 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1655 outputContainer->Add(fhMCMassPt[i]) ;
1657 fhMCSelectedMassPt[i] = new TH2F
1658 (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1659 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1660 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1661 fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1662 fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1663 outputContainer->Add(fhMCSelectedMassPt[i]) ;
1665 if(fAnaType == kSSCalo)
1667 fhMCMassPtNoOverlap[i] = new TH2F
1668 (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1669 Form("all pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1670 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1671 fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1672 fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1673 outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1675 fhMCSelectedMassPtNoOverlap[i] = new TH2F
1676 (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1677 Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s, no overlap",ptype[i].Data()),
1678 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1679 fhMCSelectedMassPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
1680 fhMCSelectedMassPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
1681 outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1684 if( fFillSelectClHisto )
1686 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1687 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1688 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1689 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1690 fhEMCLambda0[i]->SetXTitle("E (GeV)");
1691 outputContainer->Add(fhEMCLambda0[i]) ;
1693 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1694 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1695 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1696 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1697 fhEMCLambda1[i]->SetXTitle("E (GeV)");
1698 outputContainer->Add(fhEMCLambda1[i]) ;
1700 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1701 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1702 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1703 fhEMCDispersion[i]->SetYTitle("D^{2}");
1704 fhEMCDispersion[i]->SetXTitle("E (GeV)");
1705 outputContainer->Add(fhEMCDispersion[i]) ;
1707 if(fCalorimeter=="EMCAL")
1709 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1710 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1711 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1712 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1713 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1714 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
1716 if(!fFillOnlySimpleSSHisto)
1718 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1719 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1720 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1721 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1722 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1723 outputContainer->Add(fhMCEDispEta[i]);
1725 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1726 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1727 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1728 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1729 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1730 outputContainer->Add(fhMCEDispPhi[i]);
1732 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1733 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
1734 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1735 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1736 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1737 outputContainer->Add(fhMCESumEtaPhi[i]);
1739 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1740 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1741 nptbins,ptmin,ptmax,200,-10,10);
1742 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1743 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1744 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1746 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1747 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()),
1748 nptbins,ptmin,ptmax, 200,-1,1);
1749 fhMCESphericity[i]->SetXTitle("E (GeV)");
1750 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1751 outputContainer->Add(fhMCESphericity[i]);
1753 for(Int_t ie = 0; ie < 7; ie++)
1755 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1756 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]),
1757 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1758 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1759 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1760 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1762 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1763 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]),
1764 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1765 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1766 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1767 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1769 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1770 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]),
1771 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1772 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1773 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1774 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1780 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1781 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1782 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1783 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1784 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1785 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1787 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1788 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1789 nptbins,ptmin,ptmax,100,0,1);
1790 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1791 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1792 outputContainer->Add(fhEMCFracMaxCell[i]) ;
1795 } // shower shape histo
1800 if(fAnaType==kSSCalo)
1802 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1803 nptbins,ptmin,ptmax, 200, -1,1);
1804 fhAsymmetry->SetXTitle("E (GeV)");
1805 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1806 outputContainer->Add(fhAsymmetry);
1808 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1809 nptbins,ptmin,ptmax, 200, -1,1);
1810 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1811 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1812 outputContainer->Add(fhSelectedAsymmetry);
1815 ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1816 fhSplitE->SetYTitle("counts");
1817 fhSplitE->SetXTitle("E (GeV)");
1818 outputContainer->Add(fhSplitE) ;
1820 fhSplitPt = new TH1F
1821 ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
1822 fhSplitPt->SetYTitle("counts");
1823 fhSplitPt->SetXTitle("p_{T} (GeV/c)");
1824 outputContainer->Add(fhSplitPt) ;
1827 fhSplitPtPhi = new TH2F
1828 ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1829 fhSplitPtPhi->SetYTitle("#phi (rad)");
1830 fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
1831 outputContainer->Add(fhSplitPtPhi) ;
1833 fhSplitPtEta = new TH2F
1834 ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1835 fhSplitPtEta->SetYTitle("#eta");
1836 fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
1837 outputContainer->Add(fhSplitPtEta) ;
1840 fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1841 nptbins,ptmin,ptmax,10,0,10);
1842 fhNLocMaxSplitPt ->SetYTitle("N maxima");
1843 fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
1844 outputContainer->Add(fhNLocMaxSplitPt) ;
1847 fhMassSplitPt = new TH2F
1848 ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",
1849 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1850 fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1851 fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1852 outputContainer->Add(fhMassSplitPt) ;
1854 fhSelectedMassSplitPt = new TH2F
1855 ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",
1856 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1857 fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1858 fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1859 outputContainer->Add(fhSelectedMassSplitPt) ;
1863 fhMassSplitPtNoOverlap = new TH2F
1864 ("hMassSplitPtNoOverlap","all pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1865 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1866 fhMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1867 fhMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1868 outputContainer->Add(fhMassSplitPtNoOverlap) ;
1870 fhSelectedMassSplitPtNoOverlap = new TH2F
1871 ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass, no overlap",
1872 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1873 fhSelectedMassSplitPtNoOverlap->SetYTitle("mass (GeV/c^{2})");
1874 fhSelectedMassSplitPtNoOverlap->SetXTitle("p_{T} (GeV/c)");
1875 outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
1878 fhMCPi0PtRecoPtPrim = new TH2F
1879 ("hMCPi0PtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1880 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1881 fhMCPi0PtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1882 fhMCPi0PtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1883 outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
1885 fhMCPi0PtRecoPtPrimNoOverlap = new TH2F
1886 ("hMCPi0PtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1887 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1888 fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1889 fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1890 outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
1892 fhMCPi0SelectedPtRecoPtPrim = new TH2F
1893 ("hMCPi0SelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1894 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1895 fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1896 fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1897 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
1899 fhMCPi0SelectedPtRecoPtPrimNoOverlap = new TH2F
1900 ("hMCPi0SelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1901 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1902 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1903 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1904 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
1907 fhMCPi0SplitPtRecoPtPrim = new TH2F
1908 ("hMCPi0SplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1909 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1910 fhMCPi0SplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1911 fhMCPi0SplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1912 outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
1914 fhMCPi0SplitPtRecoPtPrimNoOverlap = new TH2F
1915 ("hMCPi0SplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1916 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1917 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1918 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1919 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
1921 fhMCPi0SelectedSplitPtRecoPtPrim = new TH2F
1922 ("hMCPi0SelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1923 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1924 fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1925 fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1926 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
1928 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap = new TH2F
1929 ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1930 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1931 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1932 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1933 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
1935 fhMCEtaPtRecoPtPrim = new TH2F
1936 ("hMCEtaPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1937 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1938 fhMCEtaPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1939 fhMCEtaPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1940 outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
1942 fhMCEtaPtRecoPtPrimNoOverlap = new TH2F
1943 ("hMCEtaPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1944 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1945 fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1946 fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1947 outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
1949 fhMCEtaSelectedPtRecoPtPrim = new TH2F
1950 ("hMCEtaSelectedPtRecoPtPrim","p_{T,reco} vs p_{T,gen}",
1951 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1952 fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1953 fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1954 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
1956 fhMCEtaSelectedPtRecoPtPrimNoOverlap = new TH2F
1957 ("hMCEtaSelectedPtRecoPtPrimNoOverlap","p_{T,reco} vs p_{T,gen}, no overlap",
1958 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1959 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1960 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1961 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
1964 fhMCEtaSplitPtRecoPtPrim = new TH2F
1965 ("hMCEtaSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1966 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1967 fhMCEtaSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1968 fhMCEtaSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1969 outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
1971 fhMCEtaSplitPtRecoPtPrimNoOverlap = new TH2F
1972 ("hMCEtaSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1973 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1974 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1975 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1976 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
1978 fhMCEtaSelectedSplitPtRecoPtPrim = new TH2F
1979 ("hMCEtaSelectedSplitPtRecoPtPrim","p_{T,reco} (split sum) vs p_{T,gen}",
1980 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1981 fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("p_{T,gen} (GeV/c)");
1982 fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("p_{T,reco} (GeV/c)");
1983 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
1985 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap = new TH2F
1986 ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","p_{T,reco} (split sum) vs p_{T,gen}, no overlap",
1987 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1988 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("p_{T,gen} (GeV/c)");
1989 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("p_{T,reco} (GeV/c)");
1990 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
1992 for(Int_t i = 0; i< 6; i++)
1994 fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1995 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1996 nptbins,ptmin,ptmax, 200,-1,1);
1997 fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1998 fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1999 outputContainer->Add(fhMCEAsymmetry[i]);
2001 fhMCSplitE[i] = new TH1F
2002 (Form("hSplitE_MC%s",pname[i].Data()),
2003 Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2004 nptbins,ptmin,ptmax);
2005 fhMCSplitE[i]->SetYTitle("counts");
2006 fhMCSplitE[i]->SetXTitle("E (GeV)");
2007 outputContainer->Add(fhMCSplitE[i]) ;
2009 fhMCSplitPt[i] = new TH1F
2010 (Form("hSplitPt_MC%s",pname[i].Data()),
2011 Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
2012 nptbins,ptmin,ptmax);
2013 fhMCSplitPt[i]->SetYTitle("counts");
2014 fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2015 outputContainer->Add(fhMCSplitPt[i]) ;
2018 fhMCSplitPtPhi[i] = new TH2F
2019 (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2020 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2021 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2022 fhMCSplitPtPhi[i]->SetYTitle("#phi");
2023 fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
2024 outputContainer->Add(fhMCSplitPtPhi[i]) ;
2026 fhMCSplitPtEta[i] = new TH2F
2027 (Form("hSplitPtEta_MC%s",pname[i].Data()),
2028 Form("Identified as #pi^{0} (#eta), cluster from %s",
2029 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2030 fhMCSplitPtEta[i]->SetYTitle("#eta");
2031 fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
2032 outputContainer->Add(fhMCSplitPtEta[i]) ;
2035 fhMCNLocMaxSplitPt[i] = new TH2F
2036 (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2037 Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
2038 nptbins,ptmin,ptmax,10,0,10);
2039 fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
2040 fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
2041 outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2043 fhMCMassSplitPt[i] = new TH2F
2044 (Form("hMassSplitPt_MC%s",pname[i].Data()),
2045 Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2046 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2047 fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2048 fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2049 outputContainer->Add(fhMCMassSplitPt[i]) ;
2051 fhMCSelectedMassSplitPt[i] = new TH2F
2052 (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2053 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
2054 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2055 fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
2056 fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
2057 outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2059 fhMCMassSplitPtNoOverlap[i] = new TH2F
2060 (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2061 Form("all pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2062 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2063 fhMCMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2064 fhMCMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2065 outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2067 fhMCSelectedMassSplitPtNoOverlap[i] = new TH2F
2068 (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2069 Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s, no overlap",ptype[i].Data()),
2070 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2071 fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("mass (GeV/c^{2})");
2072 fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("p_{T} (GeV/c)");
2073 outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2078 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2082 for(Int_t i = 0; i< 3; i++)
2084 fhEAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2085 Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
2086 nptbins,ptmin,ptmax,200, -1,1);
2087 fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2088 fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
2089 outputContainer->Add(fhEAsymmetryLocMax[i]) ;
2092 for(Int_t ie = 0; ie< 7; ie++)
2095 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2096 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2097 ssbins,ssmin,ssmax , 200,-1,1);
2098 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2099 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2100 outputContainer->Add(fhAsymmetryLambda0[ie]);
2102 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2103 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2104 ssbins,ssmin,ssmax , 200,-1,1);
2105 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2106 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2107 outputContainer->Add(fhAsymmetryDispEta[ie]);
2109 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2110 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
2111 ssbins,ssmin,ssmax , 200,-1,1);
2112 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2113 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2114 outputContainer->Add(fhAsymmetryDispPhi[ie]);
2120 for(Int_t i = 0; i< 6; i++)
2122 for(Int_t ie = 0; ie < 7; ie++)
2124 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2125 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2126 ssbins,ssmin,ssmax , 200,-1,1);
2127 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2128 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2129 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2131 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2132 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2133 ssbins,ssmin,ssmax , 200,-1,1);
2134 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2135 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2136 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2138 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2139 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
2140 ssbins,ssmin,ssmax , 200,-1,1);
2141 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2142 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2143 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2149 if(fFillPileUpHistograms)
2152 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2154 for(Int_t i = 0 ; i < 7 ; i++)
2156 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2157 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2158 fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2159 outputContainer->Add(fhPtPileUp[i]);
2161 fhPtCellTimePileUp[i] = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2162 Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2163 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2164 fhPtCellTimePileUp[i]->SetXTitle("p_{T} (GeV/c)");
2165 fhPtCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
2166 outputContainer->Add(fhPtCellTimePileUp[i]);
2168 fhPtTimeDiffPileUp[i] = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2169 Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2170 nptbins,ptmin,ptmax,400,-200,200);
2171 fhPtTimeDiffPileUp[i]->SetXTitle("p_{T} (GeV/c");
2172 fhPtTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2173 outputContainer->Add(fhPtTimeDiffPileUp[i]);
2177 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2178 fhTimePtNoCut->SetXTitle("p_{T} (GeV/c)");
2179 fhTimePtNoCut->SetYTitle("time (ns)");
2180 outputContainer->Add(fhTimePtNoCut);
2182 fhTimePtSPD = new TH2F ("hTimePt_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2183 fhTimePtSPD->SetXTitle("p_{T} (GeV/c)");
2184 fhTimePtSPD->SetYTitle("time (ns)");
2185 outputContainer->Add(fhTimePtSPD);
2187 fhTimePtSPDMulti = new TH2F ("hTimePt_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2188 fhTimePtSPDMulti->SetXTitle("p_{T} (GeV/c)");
2189 fhTimePtSPDMulti->SetYTitle("time (ns)");
2190 outputContainer->Add(fhTimePtSPDMulti);
2192 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
2193 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2194 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
2195 outputContainer->Add(fhTimeNPileUpVertSPD);
2197 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
2198 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2199 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
2200 outputContainer->Add(fhTimeNPileUpVertTrack);
2202 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
2203 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2204 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
2205 outputContainer->Add(fhTimeNPileUpVertContributors);
2207 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
2208 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
2209 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
2210 outputContainer->Add(fhTimePileUpMainVertexZDistance);
2212 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
2213 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
2214 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
2215 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2217 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
2218 nptbins,ptmin,ptmax,20,0,20);
2219 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2220 fhPtNPileUpSPDVtx->SetXTitle("p_{T} (GeV/c)");
2221 outputContainer->Add(fhPtNPileUpSPDVtx);
2223 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
2224 nptbins,ptmin,ptmax, 20,0,20 );
2225 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2226 fhPtNPileUpTrkVtx->SetXTitle("p_{T} (GeV/c)");
2227 outputContainer->Add(fhPtNPileUpTrkVtx);
2229 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2230 nptbins,ptmin,ptmax,20,0,20);
2231 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2232 fhPtNPileUpSPDVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2233 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2235 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2236 nptbins,ptmin,ptmax, 20,0,20 );
2237 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2238 fhPtNPileUpTrkVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2239 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2241 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2242 nptbins,ptmin,ptmax,20,0,20);
2243 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2244 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2245 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2247 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2248 nptbins,ptmin,ptmax, 20,0,20 );
2249 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2250 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2251 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2255 //Keep neutral meson selection histograms if requiered
2256 //Setting done in AliNeutralMesonSelection
2258 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2260 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2262 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2263 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2268 return outputContainer ;
2272 //_____________________________________________
2273 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2276 // Assign mc index depending on MC bit set
2278 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
2282 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2286 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2287 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2289 return kmcConversion ;
2290 }//conversion photon
2291 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2294 }//photon no conversion
2295 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2297 return kmcElectron ;
2306 //__________________________________________________________________
2307 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2308 AliAODPWG4Particle * photon2,
2309 Int_t & label, Int_t & tag)
2311 // Check the labels of pare in case mother was same pi0 or eta
2312 // Set the new AOD accordingly
2314 Int_t label1 = photon1->GetLabel();
2315 Int_t label2 = photon2->GetLabel();
2317 if(label1 < 0 || label2 < 0 ) return ;
2319 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2320 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2321 Int_t tag1 = photon1->GetTag();
2322 Int_t tag2 = photon2->GetTag();
2324 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2325 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2326 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
2327 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2328 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
2332 //Check if pi0/eta mother is the same
2333 if(GetReader()->ReadStack())
2337 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2338 label1 = mother1->GetFirstMother();
2339 //mother1 = GetMCStack()->Particle(label1);//pi0
2343 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2344 label2 = mother2->GetFirstMother();
2345 //mother2 = GetMCStack()->Particle(label2);//pi0
2348 else if(GetReader()->ReadAODMCParticles())
2349 {//&& (input > -1)){
2352 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2353 label1 = mother1->GetMother();
2354 //mother1 = GetMCStack()->Particle(label1);//pi0
2358 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2359 label2 = mother2->GetMother();
2360 //mother2 = GetMCStack()->Particle(label2);//pi0
2364 //printf("mother1 %d, mother2 %d\n",label1,label2);
2365 if( label1 == label2 && label1>=0 )
2370 TLorentzVector mom1 = *(photon1->Momentum());
2371 TLorentzVector mom2 = *(photon2->Momentum());
2373 Double_t angle = mom2.Angle(mom1.Vect());
2374 Double_t mass = (mom1+mom2).M();
2375 Double_t epair = (mom1+mom2).E();
2377 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2379 fhMassPairMCPi0 ->Fill(epair,mass);
2380 fhAnglePairMCPi0->Fill(epair,angle);
2381 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2385 fhMassPairMCEta ->Fill(epair,mass);
2386 fhAnglePairMCEta->Fill(epair,angle);
2387 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2391 } // both from eta or pi0 decay
2395 //____________________________________________________________________________
2396 void AliAnaPi0EbE::Init()
2400 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2401 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2404 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2405 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2411 //____________________________________________________________________________
2412 void AliAnaPi0EbE::InitParameters()
2414 //Initialize the parameters of the analysis.
2415 AddToHistogramsName("AnaPi0EbE_");
2417 fInputAODGammaConvName = "PhotonsCTS" ;
2418 fAnaType = kIMCalo ;
2419 fCalorimeter = "EMCAL" ;
2424 fNLMECutMin[0] = 10.;
2425 fNLMECutMin[1] = 6. ;
2426 fNLMECutMin[2] = 6. ;
2430 //__________________________________________________________________
2431 void AliAnaPi0EbE::MakeAnalysisFillAOD()
2433 //Do analysis and fill aods
2438 MakeInvMassInCalorimeter();
2442 MakeShowerShapeIdentification();
2446 MakeInvMassInCalorimeterAndCTS();
2452 //____________________________________________
2453 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
2455 //Do analysis and fill aods
2456 //Search for the photon decay in calorimeters
2457 //Read photon list from AOD, produced in class AliAnaPhoton
2458 //Check if 2 photons have the mass of the pi0.
2460 TLorentzVector mom1;
2461 TLorentzVector mom2;
2462 TLorentzVector mom ;
2467 if(!GetInputAODBranch()){
2468 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2472 //Get shower shape information of clusters
2473 TObjArray *clusters = 0;
2474 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2475 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2477 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2478 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2480 //Vertex cut in case of mixed events
2481 Int_t evtIndex1 = 0 ;
2483 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2484 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
2485 mom1 = *(photon1->Momentum());
2487 //Get original cluster, to recover some information
2489 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2492 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2496 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2498 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2500 Int_t evtIndex2 = 0 ;
2502 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2504 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2507 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
2509 mom2 = *(photon2->Momentum());
2511 //Get original cluster, to recover some information
2513 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2517 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2521 Float_t e1 = photon1->E();
2522 Float_t e2 = photon2->E();
2524 //Select clusters with good time window difference
2525 Float_t tof1 = cluster1->GetTOF()*1e9;;
2526 Float_t tof2 = cluster2->GetTOF()*1e9;;
2527 Double_t t12diff = tof1-tof2;
2528 fhEPairDiffTime->Fill(e1+e2, t12diff);
2529 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2531 //Play with the MC stack if available
2532 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2534 // Check the invariant mass for different selection on the local maxima
2535 // Name of AOD method TO BE FIXED
2536 Int_t nMaxima1 = photon1->GetFiducialArea();
2537 Int_t nMaxima2 = photon2->GetFiducialArea();
2539 Double_t mass = (mom1+mom2).M();
2540 Double_t epair = (mom1+mom2).E();
2542 if(nMaxima1==nMaxima2)
2544 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2545 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2546 else fhMassPairLocMax[2]->Fill(epair,mass);
2548 else if(nMaxima1==1 || nMaxima2==1)
2550 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2551 else fhMassPairLocMax[4]->Fill(epair,mass);
2554 fhMassPairLocMax[5]->Fill(epair,mass);
2556 // combinations with SS axis cut and NLM cut
2557 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2558 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2559 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2560 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2562 //Skip events with too few or too many NLM
2563 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2565 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2568 fhMass->Fill(epair,(mom1+mom2).M());
2570 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2571 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2574 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());
2576 //Fill some histograms about shower shape
2577 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2579 FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
2580 FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
2583 // Tag both photons as decay
2584 photon1->SetTagged(kTRUE);
2585 photon2->SetTagged(kTRUE);
2587 fhPtDecay->Fill(photon1->Pt());
2588 fhEDecay ->Fill(photon1->E() );
2590 fhPtDecay->Fill(photon2->Pt());
2591 fhEDecay ->Fill(photon2->E() );
2593 //Create AOD for analysis
2596 //Mass of selected pairs
2597 fhSelectedMass->Fill(epair,mom.M());
2599 // Fill histograms to undertand pile-up before other cuts applied
2600 // Remember to relax time cuts in the reader
2601 FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2603 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2605 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2606 pi0.SetDetector(photon1->GetDetector());
2609 pi0.SetLabel(label);
2612 //Set the indeces of the original caloclusters
2613 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2614 //pi0.SetInputFileIndex(input);
2616 AddAODParticle(pi0);
2624 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2628 //__________________________________________________
2629 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2631 //Do analysis and fill aods
2632 //Search for the photon decay in calorimeters
2633 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2634 //Check if 2 photons have the mass of the pi0.
2636 TLorentzVector mom1;
2637 TLorentzVector mom2;
2638 TLorentzVector mom ;
2643 // Check calorimeter input
2644 if(!GetInputAODBranch()){
2645 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2649 // Get the array with conversion photons
2650 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2651 if(!inputAODGammaConv) {
2653 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2655 if(!inputAODGammaConv) {
2656 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2662 //Get shower shape information of clusters
2663 TObjArray *clusters = 0;
2664 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2665 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2667 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
2668 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2669 if(nCTS<=0 || nCalo <=0)
2671 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2676 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2678 // Do the loop, first calo, second CTS
2679 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2680 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2681 mom1 = *(photon1->Momentum());
2683 //Get original cluster, to recover some information
2685 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2687 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2688 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2690 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2691 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2693 mom2 = *(photon2->Momentum());
2695 Double_t mass = (mom1+mom2).M();
2696 Double_t epair = (mom1+mom2).E();
2698 Int_t nMaxima = photon1->GetFiducialArea();
2699 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2700 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2701 else fhMassPairLocMax[2]->Fill(epair,mass);
2703 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2704 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2706 //Play with the MC stack if available
2709 Int_t label2 = photon2->GetLabel();
2710 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2712 HasPairSameMCMother(photon1, photon2, label, tag) ;
2715 //Mass of selected pairs
2716 fhMass->Fill(epair,(mom1+mom2).M());
2718 //Select good pair (good phi, pt cuts, aperture and invariant mass)
2719 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2721 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());
2723 //Fill some histograms about shower shape
2724 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2726 FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
2729 // Tag both photons as decay
2730 photon1->SetTagged(kTRUE);
2731 photon2->SetTagged(kTRUE);
2733 fhPtDecay->Fill(photon1->Pt());
2734 fhEDecay ->Fill(photon1->E() );
2736 //Create AOD for analysis
2740 //Mass of selected pairs
2741 fhSelectedMass->Fill(epair,mom.M());
2743 // Fill histograms to undertand pile-up before other cuts applied
2744 // Remember to relax time cuts in the reader
2745 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2747 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2749 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2750 pi0.SetDetector(photon1->GetDetector());
2753 pi0.SetLabel(label);
2756 //Set the indeces of the original tracks or caloclusters
2757 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2758 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2759 //pi0.SetInputFileIndex(input);
2761 AddAODParticle(pi0);
2768 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2773 //_________________________________________________
2774 void AliAnaPi0EbE::MakeShowerShapeIdentification()
2776 //Search for pi0 in fCalorimeter with shower shape analysis
2778 TObjArray * pl = 0x0;
2779 AliVCaloCells * cells = 0x0;
2780 //Select the Calorimeter of the photon
2781 if (fCalorimeter == "PHOS" )
2783 pl = GetPHOSClusters();
2784 cells = GetPHOSCells();
2786 else if (fCalorimeter == "EMCAL")
2788 pl = GetEMCALClusters();
2789 cells = GetEMCALCells();
2794 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2798 TLorentzVector mom ;
2799 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2801 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2803 Int_t evtIndex = 0 ;
2804 if (GetMixedEvent())
2806 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2809 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2811 //Get Momentum vector,
2812 Double_t vertex[]={0,0,0};
2813 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2815 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2816 }//Assume that come from vertex in straight line
2819 calo->GetMomentum(mom,vertex) ;
2822 //If too small or big pt, skip it
2823 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2825 //Check acceptance selection
2826 if(IsFiducialCutOn())
2828 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2829 if(! in ) continue ;
2833 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());
2835 //Play with the MC stack if available
2836 //Check origin of the candidates
2840 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
2841 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2842 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2845 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
2847 //Check Distance to Bad channel, set bit.
2848 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2849 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2850 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
2851 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2854 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2856 //If too low number of cells, skip it
2857 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
2859 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2864 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
2865 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
2867 //.......................................
2868 // TOF cut, BE CAREFUL WITH THIS CUT
2869 Double_t tof = calo->GetTOF()*1e9;
2870 if(tof < fTimeCutMin || tof > fTimeCutMax)
2872 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2877 //PID selection or bit setting
2879 Double_t mass = 0, angle = 0;
2880 Int_t absId1 =-1, absId2 =-1;
2881 Float_t distbad1 =-1, distbad2 =-1;
2882 Bool_t fidcut1 = 0, fidcut2 = 0;
2883 TLorentzVector l1, l2;
2885 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2886 GetVertex(evtIndex),nMaxima,
2887 mass,angle,l1,l2,absId1,absId2,
2888 distbad1,distbad2,fidcut1,fidcut2) ;
2891 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2894 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
2895 if( (fCheckSplitDistToBad) &&
2896 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
2899 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
2900 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
2902 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2906 //Skip events with too few or too many NLM
2907 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2909 //FillRejectedClusterHistograms(mom,tag,nMaxima);
2914 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2916 //Skip matched clusters with tracks
2917 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
2919 FillRejectedClusterHistograms(mom,tag,nMaxima);
2923 Float_t e1 = l1.Energy();
2924 Float_t e2 = l2.Energy();
2925 TLorentzVector l12 = l1+l2;
2926 Float_t ptSplit = l12.Pt();
2927 Float_t eSplit = e1+e2;
2930 Int_t noverlaps = 0;
2934 mcIndex = GetMCIndex(tag);
2937 Int_t mcLabel = calo->GetLabel();
2939 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
2941 Int_t mesonLabel = -1;
2943 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
2945 if(mcIndex == kmcPi0)
2947 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
2948 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
2952 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
2953 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
2957 const UInt_t nlabels = calo->GetNLabels();
2958 Int_t overpdg[nlabels];
2959 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
2963 if ( nMaxima == 1 ) inlm = 0;
2964 else if( nMaxima == 2 ) inlm = 1;
2965 else if( nMaxima < 1 )
2967 Info("MakeShowerShapeIdentification","Wrong number of maxima %d\n",nMaxima);
2971 //mass of all clusters
2972 fhMass ->Fill(mom.E() ,mass);
2973 fhMassPt ->Fill(mom.Pt(),mass);
2974 fhMassSplitPt->Fill(ptSplit ,mass);
2978 fhMCMassPt[mcIndex] ->Fill(mom.Pt(),mass);
2979 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
2982 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
2983 fhMCPi0SplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
2985 else if(mcIndex==kmcEta)
2987 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
2988 fhMCEtaSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
2995 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
2996 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
2998 else if(mcIndex==kmcEta)
3000 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3001 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3004 fhMassNoOverlap ->Fill(mom.E() ,mass);
3005 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3006 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3008 fhMCMassPtNoOverlap[mcIndex] ->Fill(mom.Pt(),mass);
3009 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3013 // Asymmetry of all clusters
3016 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3017 fhAsymmetry->Fill(mom.E(),asy);
3021 fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
3024 // If cluster does not pass pid, not pi0/eta, skip it.
3025 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3027 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3028 FillRejectedClusterHistograms(mom,tag,nMaxima);
3032 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3034 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3035 FillRejectedClusterHistograms(mom,tag,nMaxima);
3040 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3041 mom.Pt(), idPartType);
3043 //Mass and asymmetry of selected pairs
3044 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3045 fhSelectedMass ->Fill(mom.E() ,mass);
3046 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3047 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3053 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3054 fhMCPi0SelectedSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3056 else if(mcIndex==kmcEta)
3058 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3059 fhMCEtaSelectedSplitPtRecoPtPrim->Fill(ptSplit ,ptprim);
3064 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3065 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3066 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3070 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3071 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3073 else if(mcIndex==kmcEta)
3075 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3076 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3081 fhSplitE ->Fill( eSplit);
3082 fhSplitPt ->Fill(ptSplit);
3083 Float_t phi = mom.Phi();
3084 if(phi<0) phi+=TMath::TwoPi();
3085 fhSplitPtPhi ->Fill(ptSplit,phi);
3086 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3087 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3088 fhNLocMaxPt ->Fill(mom.Pt(),nMaxima);
3090 //Check split-clusters with good time window difference
3091 Double_t tof1 = cells->GetCellTime(absId1);
3092 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3095 Double_t tof2 = cells->GetCellTime(absId2);
3096 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3099 Double_t t12diff = tof1-tof2;
3100 fhEPairDiffTime->Fill(e1+e2, t12diff);
3104 fhMCSplitE [mcIndex]->Fill( eSplit);
3105 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3106 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3107 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3108 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3109 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3111 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3112 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3115 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3116 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3120 //-----------------------
3121 //Create AOD for analysis
3123 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3124 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3125 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3127 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3128 aodpi0.SetLabel(calo->GetLabel());
3130 //Set the indeces of the original caloclusters
3131 aodpi0.SetCaloLabel(calo->GetID(),-1);
3132 aodpi0.SetDetector(fCalorimeter);
3134 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3135 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3136 else aodpi0.SetDistToBad(0) ;
3138 // Check if cluster is pi0 via cluster splitting
3139 aodpi0.SetIdentifiedParticleType(idPartType);
3141 // Add number of local maxima to AOD, method name in AOD to be FIXED
3142 aodpi0.SetFiducialArea(nMaxima);
3146 //Fill some histograms about shower shape
3147 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3149 FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
3152 // Fill histograms to undertand pile-up before other cuts applied
3153 // Remember to relax time cuts in the reader
3154 Double_t tofcluster = calo->GetTOF()*1e9;
3155 Double_t tofclusterUS = TMath::Abs(tofcluster);
3157 FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
3159 Int_t id = GetReader()->GetTriggerClusterId();
3160 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
3162 Float_t phicluster = aodpi0.Phi();
3163 if(phicluster < 0) phicluster+=TMath::TwoPi();
3167 if (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
3168 else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
3169 else fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
3172 Int_t bc = GetReader()->GetTriggerClusterBC();
3173 if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
3175 if(GetReader()->IsTriggerMatched())
3177 if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
3178 fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
3179 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
3183 if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
3184 fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
3188 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(calo->E(), tofcluster);
3189 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
3190 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(calo->E(), tofcluster);
3194 else if(TMath::Abs(bc) >= 6)
3195 Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
3198 //Add AOD with pi0 object to aod branch
3199 AddAODParticle(aodpi0);
3203 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3206 //______________________________________________
3207 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3209 //Do analysis and fill histograms
3211 if(!GetOutputAODBranch())
3213 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3215 //Loop on stored AOD pi0
3216 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3217 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3219 Float_t cen = GetEventCentrality();
3220 Float_t ep = GetEventPlaneAngle();
3222 for(Int_t iaod = 0; iaod < naod ; iaod++)
3225 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3226 Int_t pdg = pi0->GetIdentifiedParticleType();
3228 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
3230 //Fill pi0 histograms
3231 Float_t ener = pi0->E();
3232 Float_t pt = pi0->Pt();
3233 Float_t phi = pi0->Phi();
3234 if(phi < 0) phi+=TMath::TwoPi();
3235 Float_t eta = pi0->Eta();
3240 fhEEta ->Fill(ener,eta);
3241 fhEPhi ->Fill(ener,phi);
3242 fhPtEta ->Fill(pt ,eta);
3243 fhPtPhi ->Fill(pt ,phi);
3244 fhEtaPhi ->Fill(eta ,phi);
3246 fhPtCentrality ->Fill(pt,cen) ;
3247 fhPtEventPlane ->Fill(pt,ep ) ;
3251 Int_t tag = pi0->GetTag();
3252 Int_t mcIndex = GetMCIndex(tag);
3254 fhMCE [mcIndex] ->Fill(ener);
3255 fhMCPt [mcIndex] ->Fill(pt);
3256 fhMCPhi[mcIndex] ->Fill(pt,phi);
3257 fhMCEta[mcIndex] ->Fill(pt,eta);
3259 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3261 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
3263 Float_t efracMC = 0;
3264 Int_t label = pi0->GetLabel();
3265 Int_t momlabel = -1;
3268 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3271 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3273 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3274 if(grandmom.E() > 0 && ok)
3276 efracMC = grandmom.E()/ener;
3277 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3280 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3282 fhMCPi0DecayPt->Fill(pt);
3283 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3284 if(grandmom.E() > 0 && ok)
3286 efracMC = mom.E()/grandmom.E();
3287 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3290 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3292 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3293 if(grandmom.E() > 0 && ok)
3295 efracMC = grandmom.E()/ener;
3296 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3299 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3301 fhMCEtaDecayPt->Fill(pt);
3302 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3303 if(grandmom.E() > 0 && ok)
3305 efracMC = mom.E()/grandmom.E();
3306 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3309 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3311 fhMCOtherDecayPt->Fill(pt);
3316 }//Histograms with MC
3322 //__________________________________________________________________
3323 void AliAnaPi0EbE::Print(const Option_t * opt) const
3325 //Print some relevant parameters set for the analysis
3329 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3330 AliAnaCaloTrackCorrBaseClass::Print("");
3331 printf("Analysis Type = %d \n", fAnaType) ;
3332 if(fAnaType == kSSCalo)
3334 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3335 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3336 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3337 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);