]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
fix wrong initialization of int with a float
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPi0EbE.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
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
22//
23// -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
24//////////////////////////////////////////////////////////////////////////////
25
26
27// --- ROOT system ---
28#include <TList.h>
29#include <TClonesArray.h>
30#include <TObjString.h>
31
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"
39#include "AliStack.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"
46
47ClassImp(AliAnaPi0EbE)
48
49//____________________________
50AliAnaPi0EbE::AliAnaPi0EbE() :
51AliAnaCaloTrackCorrBaseClass(),
52fAnaType(kIMCalo), fCalorimeter(""),
53fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
54fNLMCutMin(-1), fNLMCutMax(10),
55fTimeCutMin(-10000), fTimeCutMax(10000),
56fRejectTrackMatch(kTRUE), fSelectIsolatedDecay(kFALSE),
57fFillPileUpHistograms(0),
58fFillWeightHistograms(kFALSE), fFillTMHisto(0),
59fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1),
60fFillEMCALBCHistograms(0), fFillHighMultHistograms(0),
61fFillAllNLMHistograms(0),
62fInputAODGammaConvName(""),
63fCheckSplitDistToBad(0),
64// Histograms
65fhPt(0), fhE(0),
66fhPtEta(0), fhPtPhi(0), fhEtaPhi(0),
67fhEtaPhiEMCALBC0(0), fhEtaPhiEMCALBC1(0), fhEtaPhiEMCALBCN(0),
68fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
69fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
70fhTimeTriggerEMCALBC0UMReMatchBoth(0),
71fhPtCentrality(), fhPtEventPlane(0), fhMCPtCentrality(),
72fhPtReject(0), fhEReject(0),
73fhPtEtaReject(0), fhPtPhiReject(0), fhEtaPhiReject(0),
74fhMass(0), fhMassPt(0), fhMassSplitPt(0),
75fhSelectedMass(), fhSelectedMassPt(0), fhSelectedMassSplitPt(0),
76fhMassNoOverlap(0), fhMassPtNoOverlap(0), fhMassSplitPtNoOverlap(0),
77fhSelectedMassNoOverlap(0), fhSelectedMassPtNoOverlap(0), fhSelectedMassSplitPtNoOverlap(0),
78fhMCPi0PtRecoPtPrim(0), fhMCEtaPtRecoPtPrim(0),
79fhMCPi0PtRecoPtPrimNoOverlap(0), fhMCEtaPtRecoPtPrimNoOverlap(0),
80fhMCPi0SplitPtRecoPtPrim(0), fhMCEtaSplitPtRecoPtPrim(0),
81fhMCPi0SplitPtRecoPtPrimNoOverlap(0), fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
82fhMCPi0SelectedPtRecoPtPrim(0), fhMCEtaSelectedPtRecoPtPrim(0),
83fhMCPi0SelectedPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
84fhMCPi0SelectedSplitPtRecoPtPrim(0), fhMCEtaSelectedSplitPtRecoPtPrim(0),
85fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
86fhAsymmetry(0), fhSelectedAsymmetry(0),
87fhSplitE(0), fhSplitPt(0),
88fhSplitPtEta(0), fhSplitPtPhi(0),
89fhNLocMaxSplitPt(0),
90fhPtDecay(0),
91// Shower shape histos
92fhPtDispersion(0), fhPtLambda0(0), fhPtLambda0NoSplitCut(0),
93fhPtLambda1(0), fhPtLambda0NoTRD(0), fhPtLambda0FracMaxCellCut(0),
94fhPtFracMaxCell(0), fhPtFracMaxCellNoTRD(0),
95fhPtNCells(0), fhPtTime(0), fhEPairDiffTime(0),
96fhPtDispEta(0), fhPtDispPhi(0),
97fhPtSumEta(0), fhPtSumPhi(0), fhPtSumEtaPhi(0),
98fhPtDispEtaPhiDiff(0), fhPtSphericity(0),
99
100// MC histos
101fhMCE(), fhMCPt(),
102fhMCPtPhi(), fhMCPtEta(),
103fhMCEReject(), fhMCPtReject(),
104fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
105fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
106fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
107fhMCOtherDecayPt(0),
108fhMassPairMCPi0(0), fhMassPairMCEta(0),
109fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
110fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
111fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0),
112
113// Weight studies
114fhECellClusterRatio(0), fhECellClusterLogRatio(0),
115fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
116fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
117fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), fhTrackMatchedDEtaDPhiPos(0),
118fhTrackMatchedDEtaNeg(0), fhTrackMatchedDPhiNeg(0), fhTrackMatchedDEtaDPhiNeg(0),
119fhTrackMatchedMCParticlePt(0),
120fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
121fhdEdx(0), fhEOverP(0), fhEOverPNoTRD(0),
122// Number of local maxima in cluster
123fhNLocMaxPt(0), fhNLocMaxPtReject(0),
124// PileUp
125fhTimePtNoCut(0), fhTimePtSPD(0), fhTimePtSPDMulti(0),
126fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
127fhTimeNPileUpVertContributors(0),
128fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
129fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
130fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
131fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
132{
133 //default ctor
134
135 for(Int_t i = 0; i < fgkNmcTypes; i++)
136 {
137 fhMCE [i] = 0;
138 fhMCPt [i] = 0;
139 fhMCPtPhi [i] = 0;
140 fhMCPtEta [i] = 0;
141 fhMCPtCentrality [i] = 0;
142
143 fhMCSplitE [i] = 0;
144 fhMCSplitPt [i] = 0;
145 fhMCSplitPtPhi [i] = 0;
146 fhMCSplitPtEta [i] = 0;
147
148 fhMCNLocMaxPt [i] = 0;
149 fhMCNLocMaxSplitPt [i] = 0;
150 fhMCNLocMaxPtReject[i] = 0;
151
152 fhMCPtDecay [i] = 0;
153 fhMCPtLambda0 [i] = 0;
154 fhMCPtLambda0NoTRD [i] = 0;
155 fhMCPtLambda0FracMaxCellCut[i]= 0;
156 fhMCPtFracMaxCell [i] = 0;
157 fhMCPtLambda1 [i] = 0;
158 fhMCPtDispersion [i] = 0;
159
160 fhMCPtDispEta [i] = 0;
161 fhMCPtDispPhi [i] = 0;
162 fhMCPtSumEtaPhi [i] = 0;
163 fhMCPtDispEtaPhiDiff[i] = 0;
164 fhMCPtSphericity [i] = 0;
165 fhMCPtAsymmetry [i] = 0;
166
167 fhMCMassPt [i]=0;
168 fhMCMassSplitPt [i]=0;
169 fhMCSelectedMassPt [i]=0;
170 fhMCSelectedMassSplitPt[i]=0;
171
172 fhMCMassPtNoOverlap [i]=0;
173 fhMCMassSplitPtNoOverlap [i]=0;
174 fhMCSelectedMassPtNoOverlap [i]=0;
175 fhMCSelectedMassSplitPtNoOverlap[i]=0;
176
177 for(Int_t j = 0; j < 7; j++)
178 {
179 fhMCLambda0DispEta [j][i] = 0;
180 fhMCLambda0DispPhi [j][i] = 0;
181 fhMCDispEtaDispPhi [j][i] = 0;
182 fhMCAsymmetryLambda0 [j][i] = 0;
183 fhMCAsymmetryDispEta [j][i] = 0;
184 fhMCAsymmetryDispPhi [j][i] = 0;
185 }
186 }
187
188 for(Int_t j = 0; j < 7; j++)
189 {
190 fhLambda0DispEta [j] = 0;
191 fhLambda0DispPhi [j] = 0;
192 fhDispEtaDispPhi [j] = 0;
193 fhAsymmetryLambda0 [j] = 0;
194 fhAsymmetryDispEta [j] = 0;
195 fhAsymmetryDispPhi [j] = 0;
196
197 fhPtPileUp [j] = 0;
198 }
199
200 for(Int_t i = 0; i < 3; i++)
201 {
202 fhPtLambda0LocMax [i] = 0;
203 fhPtLambda1LocMax [i] = 0;
204 fhPtDispersionLocMax [i] = 0;
205 fhPtDispEtaLocMax [i] = 0;
206 fhPtDispPhiLocMax [i] = 0;
207 fhPtSumEtaPhiLocMax [i] = 0;
208 fhPtDispEtaPhiDiffLocMax[i] = 0;
209 fhPtSphericityLocMax [i] = 0;
210 fhPtAsymmetryLocMax [i] = 0;
211 fhMassPtLocMax [i] = 0;
212 fhSelectedMassPtLocMax [i] = 0;
213 for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
214 {
215 fhMCPtLambda0LocMax [ipart][i] = 0;
216 fhMCSelectedMassPtLocMax[ipart][i] = 0;
217 }
218
219 fhMCPi0PtRecoPtPrimLocMax [i] = 0;
220 fhMCEtaPtRecoPtPrimLocMax [i] = 0;
221 fhMCPi0SplitPtRecoPtPrimLocMax [i] = 0;
222 fhMCEtaSplitPtRecoPtPrimLocMax [i] = 0;
223
224 fhMCPi0SelectedPtRecoPtPrimLocMax [i] = 0;
225 fhMCEtaSelectedPtRecoPtPrimLocMax [i] = 0;
226 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[i] = 0;
227 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[i] = 0;
228
229 }
230
231 //Weight studies
232 for(Int_t i =0; i < 14; i++){
233 fhLambda0ForW0[i] = 0;
234 //fhLambda1ForW0[i] = 0;
235 if(i<8)fhMassPairLocMax[i] = 0;
236 }
237
238 for(Int_t i = 0; i < 11; i++)
239 {
240 fhEtaPhiTriggerEMCALBC [i] = 0 ;
241 fhTimeTriggerEMCALBC [i] = 0 ;
242 fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
243
244 fhEtaPhiTriggerEMCALBCUM [i] = 0 ;
245 fhTimeTriggerEMCALBCUM [i] = 0 ;
246
247 }
248
249 for(Int_t iSM = 0; iSM < 22; iSM++)
250 {
251 fhNLocMaxPtSM[iSM] = 0;
252 for(Int_t inlm = 0; inlm < 3; inlm++)
253 {
254 fhSelectedMassPtLocMaxSM [inlm][iSM] = 0;
255 fhSelectedLambda0PtLocMaxSM [inlm][iSM] = 0;
256 }
257 }
258 //Initialize parameters
259 InitParameters();
260
261}
262
263//______________________________________________________________________________________________
264void AliAnaPi0EbE::FillEMCALBCHistograms(Float_t energy, Float_t eta, Float_t phi, Float_t time)
265{
266 // EMCal trigger cluster BC studies
267
268 Int_t id = GetReader()->GetTriggerClusterId();
269 if( id < 0 ) return;
270
271 Int_t bc = GetReader()->GetTriggerClusterBC();
272 if(TMath::Abs(bc) >= 6)
273 Info("FillEMCALBCHistograms","Trigger BC not expected = %d\n",bc);
274
275 if(phi < 0) phi+=TMath::TwoPi();
276
277 if(energy > 2)
278 {
279 Double_t timeUS = TMath::Abs(time);
280
281 if (timeUS < 25) fhEtaPhiEMCALBC0->Fill(eta, phi);
282 else if (timeUS < 75) fhEtaPhiEMCALBC1->Fill(eta, phi);
283 else fhEtaPhiEMCALBCN->Fill(eta, phi);
284 }
285
286 if(TMath::Abs(bc) >= 6) return ;
287
288 if(GetReader()->IsBadCellTriggerEvent() || GetReader()->IsExoticEvent()) return ;
289
290 if(GetReader()->IsTriggerMatched())
291 {
292 if(energy > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(eta, phi);
293 fhTimeTriggerEMCALBC[bc+5]->Fill(energy, time);
294 if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(energy, time);
295 }
296 else
297 {
298 if(energy > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(eta, phi);
299 fhTimeTriggerEMCALBCUM[bc+5]->Fill(energy, time);
300
301 if(bc==0)
302 {
303 if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime ->Fill(energy, time);
304 if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(energy, time);
305 if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth ->Fill(energy, time);
306 }
307 }
308}
309
310//___________________________________________________________________________________
311void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
312{
313 // Fill some histograms to understand pile-up
314 if(!fFillPileUpHistograms) return;
315
316 //printf("E %f, time %f\n",energy,time);
317 AliVEvent * event = GetReader()->GetInputEvent();
318
319 fhTimePtNoCut->Fill(pt,time);
320 if(GetReader()->IsPileUpFromSPD())
321
322 if(GetReader()->IsPileUpFromSPD()) { fhPtPileUp[0]->Fill(pt); fhTimePtSPD ->Fill(pt,time); }
323 if(GetReader()->IsPileUpFromEMCal()) fhPtPileUp[1]->Fill(pt);
324 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPileUp[2]->Fill(pt);
325 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPileUp[3]->Fill(pt);
326 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPileUp[4]->Fill(pt);
327 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPileUp[5]->Fill(pt);
328 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
329
330 if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
331
332 // cells in cluster
333
334 AliVCaloCells* cells = 0;
335 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
336 else cells = GetPHOSCells();
337
338 Float_t maxCellFraction = 0.;
339 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
340
341 Double_t tmax = cells->GetCellTime(absIdMax);
342 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
343 tmax*=1.e9;
344
345 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
346 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
347 {
348 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
349 {
350 Int_t absId = calo->GetCellsAbsId()[ipos];
351
352 if( absId == absIdMax ) continue ;
353
354 Double_t timecell = cells->GetCellTime(absId);
355 Float_t amp = cells->GetCellAmplitude(absId);
356 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
357 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
358 timecell*=1e9;
359
360 Float_t diff = (tmax-timecell);
361
362 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
363
364 if(GetReader()->IsPileUpFromSPD())
365 {
366 fhPtCellTimePileUp[0]->Fill(pt, timecell);
367 fhPtTimeDiffPileUp[0]->Fill(pt, diff);
368 }
369
370 if(GetReader()->IsPileUpFromEMCal())
371 {
372 fhPtCellTimePileUp[1]->Fill(pt, timecell);
373 fhPtTimeDiffPileUp[1]->Fill(pt, diff);
374 }
375
376 if(GetReader()->IsPileUpFromSPDOrEMCal())
377 {
378 fhPtCellTimePileUp[2]->Fill(pt, timecell);
379 fhPtTimeDiffPileUp[2]->Fill(pt, diff);
380 }
381
382 if(GetReader()->IsPileUpFromSPDAndEMCal())
383 {
384 fhPtCellTimePileUp[3]->Fill(pt, timecell);
385 fhPtTimeDiffPileUp[3]->Fill(pt, diff);
386 }
387
388 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
389 {
390 fhPtCellTimePileUp[4]->Fill(pt, timecell);
391 fhPtTimeDiffPileUp[4]->Fill(pt, diff);
392 }
393
394 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
395 {
396 fhPtCellTimePileUp[5]->Fill(pt, timecell);
397 fhPtTimeDiffPileUp[5]->Fill(pt, diff);
398 }
399
400 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
401 {
402 fhPtCellTimePileUp[6]->Fill(pt, timecell);
403 fhPtTimeDiffPileUp[6]->Fill(pt, diff);
404 }
405 }//loop
406 }
407
408 if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
409
410 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
411 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
412
413 // N pile up vertices
414 Int_t nVtxSPD = -1;
415 Int_t nVtxTrk = -1;
416
417 if (esdEv)
418 {
419 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
420 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
421
422 }//ESD
423 else if (aodEv)
424 {
425 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
426 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
427 }//AOD
428
429 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
430 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
431
432 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
433 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
434
435 if(TMath::Abs(time) < 25)
436 {
437 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
438 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
439 }
440
441 if(time < 75 && time > -25)
442 {
443 fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
444 fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
445 }
446
447 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
448 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
449
450 Int_t ncont = -1;
451 Float_t z1 = -1, z2 = -1;
452 Float_t diamZ = -1;
453 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
454 {
455 if (esdEv)
456 {
457 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
458 ncont=pv->GetNContributors();
459 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
460 z2 = pv->GetZ();
461 diamZ = esdEv->GetDiamondZ();
462 }//ESD
463 else if (aodEv)
464 {
465 AliAODVertex *pv=aodEv->GetVertex(iVert);
466 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
467 ncont=pv->GetNContributors();
468 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
469 z2=pv->GetZ();
470 diamZ = aodEv->GetDiamondZ();
471 }// AOD
472
473 Double_t distZ = TMath::Abs(z2-z1);
474 diamZ = TMath::Abs(z2-diamZ);
475
476 fhTimeNPileUpVertContributors ->Fill(time,ncont);
477 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
478 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
479
480 }// vertex loop
481}
482
483
484//______________________________________________________________________________________________
485void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
486{
487 // Fill histograms that do not pass the identification (SS case only)
488
489 Float_t ener = mom.E();
490 Float_t pt = mom.Pt();
491 Float_t phi = mom.Phi();
492 if(phi < 0) phi+=TMath::TwoPi();
493 Float_t eta = mom.Eta();
494
495 fhPtReject ->Fill(pt);
496 fhEReject ->Fill(ener);
497
498 fhPtEtaReject ->Fill(ener,eta);
499 fhPtPhiReject ->Fill(ener,phi);
500 fhEtaPhiReject ->Fill(eta,phi);
501 fhNLocMaxPtReject->Fill(pt,nMaxima);
502
503 if(IsDataMC())
504 {
505 Int_t mcIndex = GetMCIndex(mctag);
506 fhMCEReject [mcIndex] ->Fill(ener);
507 fhMCPtReject [mcIndex] ->Fill(pt);
508 if(fFillAllNLMHistograms) fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
509 }
510
511}
512
513//___________________________________________________________________________________
514void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
515 Int_t tag, Float_t asy)
516{
517 // Fill shower shape, timing and other histograms for selected clusters from decay
518
519 Float_t ener = cluster->E();
520 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
521 Float_t l0 = cluster->GetM02();
522 Float_t l1 = cluster->GetM20();
523 Int_t nSM = GetModuleNumber(cluster);
524
525 Int_t ptbin = -1;
526 if (pt < 2 ) ptbin = 0;
527 else if (pt < 4 ) ptbin = 1;
528 else if (pt < 6 ) ptbin = 2;
529 else if (pt < 10) ptbin = 3;
530 else if (pt < 15) ptbin = 4;
531 else if (pt < 20) ptbin = 5;
532 else ptbin = 6;
533
534 Int_t indexMax = -1;
535 if (nMaxima==1) indexMax = 0 ;
536 else if(nMaxima==2) indexMax = 1 ;
537 else indexMax = 2 ;
538
539 FillWeightHistograms(cluster);
540
541 fhPtLambda0 ->Fill(pt, l0 );
542 fhPtLambda1 ->Fill(pt, l1 );
543
544 fhNLocMaxPt->Fill(pt,nMaxima);
545
546 if(fFillAllNLMHistograms)
547 {
548 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
549 fhNLocMaxPtSM[nSM]->Fill(pt,nMaxima);
550
551 fhPtLambda0LocMax [indexMax]->Fill(pt,l0);
552 fhPtLambda1LocMax [indexMax]->Fill(pt,l1);
553 }
554
555 Float_t ll0 = 0., ll1 = 0.;
556 Float_t dispp= 0., dEta = 0., dPhi = 0.;
557 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
558 AliVCaloCells * cell = 0x0;
559 Float_t maxCellFraction = 0;
560
561 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
562 {
563 cell = GetEMCALCells();
564
565 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
566 fhPtFracMaxCell->Fill(pt,maxCellFraction);
567
568 if(maxCellFraction < 0.5)
569 fhPtLambda0FracMaxCellCut->Fill(pt, l0 );
570
571 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(),
572 cell, cluster,
573 ll0, ll1, dispp, dEta, dPhi,
574 sEta, sPhi, sEtaPhi);
575 fhPtDispersion -> Fill(pt,disp);
576 fhPtDispEta -> Fill(pt,dEta);
577 fhPtDispPhi -> Fill(pt,dPhi);
578 fhPtSumEta -> Fill(pt,sEta);
579 fhPtSumPhi -> Fill(pt,sPhi);
580 fhPtSumEtaPhi -> Fill(pt,sEtaPhi);
581 fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
582 if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
583
584 fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
585 fhLambda0DispEta[ptbin]->Fill(l0 ,dEta);
586 fhLambda0DispPhi[ptbin]->Fill(l0 ,dPhi);
587
588 if (fAnaType==kSSCalo)
589 {
590 // Asymmetry histograms
591 fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
592 fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
593 fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
594 }
595
596 if(fFillAllNLMHistograms)
597 {
598 fhPtDispersionLocMax [indexMax]->Fill(pt,disp);
599 fhPtDispEtaLocMax [indexMax]-> Fill(pt,dEta);
600 fhPtDispPhiLocMax [indexMax]-> Fill(pt,dPhi);
601 fhPtSumEtaPhiLocMax [indexMax]-> Fill(pt,sEtaPhi);
602 fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
603 if(dEta+dPhi>0) fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
604 if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt ,asy);
605 }
606 }
607
608
609 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
610 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
611 {
612 fhPtLambda0NoTRD ->Fill(pt, l0 );
613 if(!fFillOnlySimpleSSHisto)
614 fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
615 }
616
617 fhPtTime ->Fill(pt, cluster->GetTOF()*1.e9);
618 fhPtNCells->Fill(pt, cluster->GetNCells());
619
620 // Fill Track matching control histograms
621 if(fFillTMHisto)
622 {
623 Float_t dZ = cluster->GetTrackDz();
624 Float_t dR = cluster->GetTrackDx();
625
626 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
627 {
628 dR = 2000., dZ = 2000.;
629 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
630 }
631 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
632
633 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
634
635 Bool_t positive = kFALSE;
636 if(track) positive = (track->Charge()>0);
637
638 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
639 {
640 fhTrackMatchedDEta->Fill(pt,dZ);
641 fhTrackMatchedDPhi->Fill(pt,dR);
642 if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
643
644 if(track)
645 {
646 if(positive)
647 {
648 fhTrackMatchedDEtaPos->Fill(pt,dZ);
649 fhTrackMatchedDPhiPos->Fill(pt,dR);
650 if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
651 }
652 else
653 {
654 fhTrackMatchedDEtaNeg->Fill(pt,dZ);
655 fhTrackMatchedDPhiNeg->Fill(pt,dR);
656 if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
657 }
658 }
659 }
660 // Check dEdx and E/p of matched clusters
661
662 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
663 {
664 if(track)
665 {
666 Float_t dEdx = track->GetTPCsignal();
667 fhdEdx->Fill(pt, dEdx);
668
669 Float_t eOverp = cluster->E()/track->P();
670 fhEOverP->Fill(pt, eOverp);
671
672 // Change nSM for year > 2011 (< 4 in 2012-13, none after)
673 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
674 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
675 fhEOverPNoTRD->Fill(pt, eOverp);
676
677 }
678 //else
679 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
680
681 if(IsDataMC())
682 {
683 Float_t mctag = -1;
684 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
685 {
686 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
687 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 2.5 ;
688 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 0.5 ;
689 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 1.5 ;
690 else mctag = 3.5 ;
691
692 }
693 else
694 {
695 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
696 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mctag = 6.5 ;
697 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) mctag = 4.5 ;
698 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) mctag = 5.5 ;
699 else mctag = 7.5 ;
700 }
701
702 fhTrackMatchedMCParticlePt ->Fill(pt, mctag);
703 fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
704 fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
705
706 } // MC
707 }
708 }// Track matching histograms
709
710 if(IsDataMC())
711 {
712 Int_t mcIndex = GetMCIndex(tag);
713
714 fhMCPtLambda0[mcIndex] ->Fill(pt, l0);
715 fhMCPtLambda1[mcIndex] ->Fill(pt, l1);
716 if(fFillAllNLMHistograms) fhMCPtLambda0LocMax[mcIndex][indexMax]->Fill(pt,l0);
717
718 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
719 GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
720 fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0 );
721
722 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
723 {
724 if(maxCellFraction < 0.5)
725 fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0 );
726
727 fhMCPtDispersion [mcIndex]->Fill(pt, disp);
728 fhMCPtFracMaxCell [mcIndex]->Fill(pt,maxCellFraction);
729
730 fhMCPtDispEta [mcIndex]-> Fill(pt,dEta);
731 fhMCPtDispPhi [mcIndex]-> Fill(pt,dPhi);
732 fhMCPtSumEtaPhi [mcIndex]-> Fill(pt,sEtaPhi);
733 fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
734 if(dEta+dPhi > 0) fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
735
736 if (fAnaType==kSSCalo)
737 {
738 fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
739 fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
740 fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
741 }
742
743 fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
744 fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0 ,dEta);
745 fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0 ,dPhi);
746 }// only SS simple?
747
748 }//MC
749
750}
751
752//________________________________________________________
753void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
754{
755 // Calculate weights and fill histograms
756
757 if(!fFillWeightHistograms || GetMixedEvent()) return;
758
759 AliVCaloCells* cells = 0;
760 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
761 else cells = GetPHOSCells();
762
763 // First recalculate energy in case non linearity was applied
764 Float_t energy = 0;
765 Float_t ampMax = 0;
766 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
767 {
768
769 Int_t id = clus->GetCellsAbsId()[ipos];
770
771 //Recalibrate cell energy if needed
772 Float_t amp = cells->GetCellAmplitude(id);
773 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
774
775 energy += amp;
776
777 if(amp> ampMax)
778 ampMax = amp;
779
780 } // energy loop
781
782 if(energy <=0 )
783 {
784 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
785 return;
786 }
787
788 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
789 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
790
791 //Get the ratio and log ratio to all cells in cluster
792 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
793 {
794 Int_t id = clus->GetCellsAbsId()[ipos];
795
796 //Recalibrate cell energy if needed
797 Float_t amp = cells->GetCellAmplitude(id);
798 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
799
800 fhECellClusterRatio ->Fill(energy,amp/energy);
801 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
802 }
803
804 //Recalculate shower shape for different W0
805 if(fCalorimeter=="EMCAL"){
806
807 Float_t l0org = clus->GetM02();
808 Float_t l1org = clus->GetM20();
809 Float_t dorg = clus->GetDispersion();
810
811 for(Int_t iw = 0; iw < 14; iw++)
812 {
813 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
814 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
815
816 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
817 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
818
819 } // w0 loop
820
821 // Set the original values back
822 clus->SetM02(l0org);
823 clus->SetM20(l1org);
824 clus->SetDispersion(dorg);
825
826 }// EMCAL
827}
828
829//__________________________________________
830TObjString * AliAnaPi0EbE::GetAnalysisCuts()
831{
832 //Save parameters used for analysis
833 TString parList ; //this will be list of parameters used for this analysis.
834 const Int_t buffersize = 255;
835 char onePar[buffersize] ;
836
837 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
838 parList+=onePar ;
839 snprintf(onePar,buffersize,"fAnaType=%d (selection type) \n",fAnaType) ;
840 parList+=onePar ;
841 snprintf(onePar,buffersize,"Calorimeter: %s;",fCalorimeter.Data()) ;
842 parList+=onePar ;
843 snprintf(onePar,buffersize,"Local maxima in cluster: %d < nlm < %d;",fNLMCutMin,fNLMCutMax) ;
844 parList+=onePar ;
845
846 if(fAnaType == kSSCalo)
847 {
848 snprintf(onePar,buffersize,"E cut: %2.2f<E<%2.2f;",GetMinEnergy(),GetMaxEnergy()) ;
849 parList+=onePar ;
850 snprintf(onePar,buffersize,"N cell cut: N > %d;",GetCaloPID()->GetClusterSplittingMinNCells()) ;
851 parList+=onePar ;
852 snprintf(onePar,buffersize,"Min Dist to Bad channel: fMinDist =%2.2f; fMinDist2=%2.2f, fMinDist3=%2.2f;",fMinDist, fMinDist2,fMinDist3) ;
853 parList+=onePar ;
854 snprintf(onePar,buffersize,"Min E cut for NLM cases: 1) %2.2f; 2) %2.2f; 3) %2.2f;",fNLMECutMin[0],fNLMECutMin[1],fNLMECutMin[2]) ;
855 parList+=onePar ;
856 snprintf(onePar,buffersize,"Reject Matched tracks?: %d;",fRejectTrackMatch) ;
857 parList+=onePar ;
858 snprintf(onePar,buffersize,"Reject split cluster close to border or bad?: %d;",fCheckSplitDistToBad) ;
859 parList+=onePar ;
860 snprintf(onePar,buffersize,"Time cut: %2.2f<t<%2.2f;",fTimeCutMin,fTimeCutMax) ;
861 parList+=onePar ;
862 //Get parameters set in PID class.
863 parList += GetCaloPID()->GetPIDParametersList() ;
864 }
865 else if(fAnaType == kIMCalo || fAnaType == kIMCaloTracks)
866 {
867 snprintf(onePar,buffersize,"Select %s;", (GetNeutralMesonSelection()->GetParticle()).Data()) ;
868 parList+=onePar ;
869 snprintf(onePar,buffersize,"Mass cut: %2.2f<M<%2.2f;",GetNeutralMesonSelection()->GetInvMassMinCut() ,GetNeutralMesonSelection()->GetInvMassMaxCut()) ;
870 parList+=onePar ;
871 }
872 else if(fAnaType == kIMCaloTracks)
873 {
874 snprintf(onePar,buffersize,"Photon Conv Array: %s;",fInputAODGammaConvName.Data()) ;
875 parList+=onePar ;
876 }
877 else if(fAnaType == kIMCalo)
878 {
879 snprintf(onePar,buffersize,"Time Diff: %2.2f;",GetPairTimeCut()) ;
880 parList+=onePar ;
881 }
882
883 //Get parameters set in base class.
884 //parList += GetBaseParametersList() ;
885
886 return new TObjString(parList) ;
887}
888
889//_____________________________________________
890TList * AliAnaPi0EbE::GetCreateOutputObjects()
891{
892 // Create histograms to be saved in output file and
893 // store them in outputContainer
894 TList * outputContainer = new TList() ;
895 outputContainer->SetName("Pi0EbEHistos") ;
896
897 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
898 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
899 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
900 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
901 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
902 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
903 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
904
905 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
906 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
907 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
908
909 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
910 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
911 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
912 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
913 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
914 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
915
916 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
917 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
918 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
919 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
920 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
921 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
922
923 Int_t ntimptbins = GetHistogramRanges()->GetHistoTimeBins();
924 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
925 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
926
927 TString nlm[] = {"1 Local Maxima","2 Local Maxima", "NLM > 2"};
928
929 TString ptype [] = {"#gamma (#pi^{0})", "#gamma (#eta)", "#gamma (other)", "#pi^{0}", "#eta", "#gamma (direct)","e^{#pm}" , "hadron/other combinations"};
930 TString pname [] = {"Pi0Decay" , "EtaDecay" , "OtherDecay" , "Pi0" , "Eta" , "Photon" , "Electron", "Hadron"};
931
932 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
933
934 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
935 fhPt->SetYTitle("#it{N}");
936 fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
937 outputContainer->Add(fhPt) ;
938
939 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
940 fhE->SetYTitle("#it{N}");
941 fhE->SetXTitle("#it{E} (GeV)");
942 outputContainer->Add(fhE) ;
943
944 fhPtPhi = new TH2F
945 ("hPtPhi","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
946 fhPtPhi->SetYTitle("#phi (rad)");
947 fhPtPhi->SetXTitle("#it{E} (GeV)");
948 outputContainer->Add(fhPtPhi) ;
949
950 fhPtEta = new TH2F
951 ("hPtEta","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
952 fhPtEta->SetYTitle("#eta");
953 fhPtEta->SetXTitle("#it{E} (GeV)");
954 outputContainer->Add(fhPtEta) ;
955
956 fhEtaPhi = new TH2F
957 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
958 fhEtaPhi->SetYTitle("#phi (rad)");
959 fhEtaPhi->SetXTitle("#eta");
960 outputContainer->Add(fhEtaPhi) ;
961
962 if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
963 {
964 fhEtaPhiEMCALBC0 = new TH2F
965 ("hEtaPhiEMCALBC0","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
966 fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
967 fhEtaPhiEMCALBC0->SetXTitle("#eta");
968 outputContainer->Add(fhEtaPhiEMCALBC0) ;
969
970 fhEtaPhiEMCALBC1 = new TH2F
971 ("hEtaPhiEMCALBC1","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |#it{t}| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
972 fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
973 fhEtaPhiEMCALBC1->SetXTitle("#eta");
974 outputContainer->Add(fhEtaPhiEMCALBC1) ;
975
976 fhEtaPhiEMCALBCN = new TH2F
977 ("hEtaPhiEMCALBCN","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
978 fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
979 fhEtaPhiEMCALBCN->SetXTitle("#eta");
980 outputContainer->Add(fhEtaPhiEMCALBCN) ;
981
982 for(Int_t i = 0; i < 11; i++)
983 {
984 fhEtaPhiTriggerEMCALBC[i] = new TH2F
985 (Form("hEtaPhiTriggerEMCALBC%d",i-5),
986 Form("meson #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
987 netabins,etamin,etamax,nphibins,phimin,phimax);
988 fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
989 fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
990 outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
991
992 fhTimeTriggerEMCALBC[i] = new TH2F
993 (Form("hTimeTriggerEMCALBC%d",i-5),
994 Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
995 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
996 fhTimeTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
997 fhTimeTriggerEMCALBC[i]->SetYTitle("#it{t} (ns)");
998 outputContainer->Add(fhTimeTriggerEMCALBC[i]);
999
1000 fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
1001 (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
1002 Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
1003 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
1004 fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
1005 fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{t} (ns)");
1006 outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
1007
1008 fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
1009 (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
1010 Form("meson #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
1011 netabins,etamin,etamax,nphibins,phimin,phimax);
1012 fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
1013 fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
1014 outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
1015
1016 fhTimeTriggerEMCALBCUM[i] = new TH2F
1017 (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
1018 Form("meson #it{t} vs #it{E}, unmatched trigger EMCAL-BC=%d",i-5),
1019 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
1020 fhTimeTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
1021 fhTimeTriggerEMCALBCUM[i]->SetYTitle("#it{t} (ns)");
1022 outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
1023
1024 }
1025
1026 fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
1027 "cluster #it{t} vs #it{E} of clusters, no match, rematch open time",
1028 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
1029 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
1030 fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{t} (ns)");
1031 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
1032
1033
1034 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
1035 "cluster #it{t} vs #it{E} of clusters, no match, rematch with neigbour parches",
1036 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
1037 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
1038 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{t} (ns)");
1039 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
1040
1041 fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
1042 "cluster #it{t} vs #it{E} of clusters, no match, rematch open time and neigbour",
1043 nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
1044 fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
1045 fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{t} (ns)");
1046 outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
1047
1048 }
1049
1050 if(fFillHighMultHistograms)
1051 {
1052 fhPtCentrality = new TH2F("hPtCentrality","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
1053 fhPtCentrality->SetYTitle("centrality");
1054 fhPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1055 outputContainer->Add(fhPtCentrality) ;
1056
1057 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1058 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
1059 fhPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1060 outputContainer->Add(fhPtEventPlane) ;
1061 }
1062
1063 if(fAnaType == kSSCalo)
1064 {
1065 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
1066 fhPtReject->SetYTitle("#it{N}");
1067 fhPtReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1068 outputContainer->Add(fhPtReject) ;
1069
1070 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
1071 fhEReject->SetYTitle("#it{N}");
1072 fhEReject->SetXTitle("#it{E} (GeV)");
1073 outputContainer->Add(fhEReject) ;
1074
1075 fhPtPhiReject = new TH2F
1076 ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1077 fhPtPhiReject->SetYTitle("#phi (rad)");
1078 fhPtPhiReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1079 outputContainer->Add(fhPtPhiReject) ;
1080
1081 fhPtEtaReject = new TH2F
1082 ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1083 fhPtEtaReject->SetYTitle("#eta");
1084 fhPtEtaReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1085 outputContainer->Add(fhPtEtaReject) ;
1086
1087 fhEtaPhiReject = new TH2F
1088 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1089 fhEtaPhiReject->SetYTitle("#phi (rad)");
1090 fhEtaPhiReject->SetXTitle("#eta");
1091 outputContainer->Add(fhEtaPhiReject) ;
1092
1093 fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
1094 nptbins,ptmin,ptmax,20,0,20);
1095 fhNLocMaxPtReject ->SetYTitle("N maxima");
1096 fhNLocMaxPtReject ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1097 outputContainer->Add(fhNLocMaxPtReject) ;
1098
1099 }
1100
1101 fhMass = new TH2F
1102 ("hMass","all pairs #it{M}: #it{E} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1103 fhMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1104 fhMass->SetXTitle("#it{E} (GeV)");
1105 outputContainer->Add(fhMass) ;
1106
1107 fhSelectedMass = new TH2F
1108 ("hSelectedMass","Selected #pi^{0} (#eta) pairs #it{M}: E vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1109 fhSelectedMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1110 fhSelectedMass->SetXTitle("#it{E} (GeV)");
1111 outputContainer->Add(fhSelectedMass) ;
1112
1113 fhMassPt = new TH2F
1114 ("hMassPt","all pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1115 fhMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1116 fhMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1117 outputContainer->Add(fhMassPt) ;
1118
1119 fhSelectedMassPt = new TH2F
1120 ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1121 fhSelectedMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1122 fhSelectedMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1123 outputContainer->Add(fhSelectedMassPt) ;
1124
1125 if(fAnaType == kSSCalo)
1126 {
1127 fhPtLambda0NoSplitCut = new TH2F
1128 ("hPtLambda0NoSplitCut","all clusters: #it{p}_{T} vs #lambda_{0}^{2}",nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1129 fhPtLambda0NoSplitCut->SetYTitle("#lambda_{0}^{2}");
1130 fhPtLambda0NoSplitCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1131 outputContainer->Add(fhPtLambda0NoSplitCut) ;
1132
1133 for(Int_t inlm = 0; inlm < 3; inlm++)
1134 {
1135 fhMassPtLocMax[inlm] = new TH2F
1136 (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs #it{M}: #it{p}_{T} vs #it{M} and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1137 fhMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1138 fhMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1139 outputContainer->Add(fhMassPtLocMax[inlm]) ;
1140
1141 fhSelectedMassPtLocMax[inlm] = new TH2F
1142 (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1143 fhSelectedMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1144 fhSelectedMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1145 outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
1146
1147 if(fFillAllNLMHistograms)
1148 {
1149 for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1150 {
1151 fhSelectedMassPtLocMaxSM[inlm][iSM] = new TH2F
1152 (Form("hSelectedMassPtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1153 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1154 fhSelectedMassPtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1155 outputContainer->Add(fhSelectedMassPtLocMaxSM[inlm][iSM]) ;
1156
1157 fhSelectedLambda0PtLocMaxSM[inlm][iSM] = new TH2F
1158 (Form("hSelectedLambda0PtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #lambda_{0}^{2}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1159 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetYTitle("#lambda_{0}^{2}");
1160 fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1161 outputContainer->Add(fhSelectedLambda0PtLocMaxSM[inlm][iSM]) ;
1162 }
1163 }
1164
1165 if(IsDataMC())
1166 {
1167 for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
1168 {
1169 fhMCSelectedMassPtLocMax[ipart][inlm] = new TH2F
1170 (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
1171 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s, from MC %s",nlm[inlm].Data(),ptype[ipart].Data()),
1172 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1173 fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1174 fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1175 outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
1176 }
1177 }
1178 }
1179
1180 if(IsDataMC())
1181 {
1182 fhMassNoOverlap = new TH2F
1183 ("hMassNoOverlap","all pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1184 fhMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1185 fhMassNoOverlap->SetXTitle("#it{E} (GeV)");
1186 outputContainer->Add(fhMassNoOverlap) ;
1187
1188 fhSelectedMassNoOverlap = new TH2F
1189 ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1190 fhSelectedMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1191 fhSelectedMassNoOverlap->SetXTitle("#it{E} (GeV)");
1192 outputContainer->Add(fhSelectedMassNoOverlap) ;
1193
1194 fhMassPtNoOverlap = new TH2F
1195 ("hMassPtNoOverlap","all pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1196 fhMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1197 fhMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1198 outputContainer->Add(fhMassPtNoOverlap) ;
1199
1200 fhSelectedMassPtNoOverlap = new TH2F
1201 ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1202 fhSelectedMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1203 fhSelectedMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1204 outputContainer->Add(fhSelectedMassPtNoOverlap) ;
1205 }
1206 }
1207
1208 if(fAnaType != kSSCalo)
1209 {
1210 fhPtDecay = new TH1F("hPtDecay","Selected #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
1211 fhPtDecay->SetYTitle("#it{N}");
1212 fhPtDecay->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1213 outputContainer->Add(fhPtDecay) ;
1214
1215 if(IsDataMC())
1216 {
1217 for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
1218 {
1219 fhMCPtDecay[ipart] = new TH1F(Form("hPtDecay_MC%s",pname[ipart].Data()),
1220 Form("Selected #pi^{0} (#eta) decay photons, from MC %s",ptype[ipart].Data()),
1221 nptbins,ptmin,ptmax);
1222 fhMCPtDecay[ipart]->SetYTitle("#it{N}");
1223 fhMCPtDecay[ipart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1224 outputContainer->Add(fhMCPtDecay[ipart]) ;
1225 }
1226 }
1227 }
1228
1229 ////////
1230
1231 if( fFillSelectClHisto )
1232 {
1233 fhPtLambda0 = new TH2F
1234 ("hPtLambda0","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1235 fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1236 fhPtLambda0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1237 outputContainer->Add(fhPtLambda0) ;
1238
1239 fhPtLambda1 = new TH2F
1240 ("hPtLambda1","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1241 fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1242 fhPtLambda1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1243 outputContainer->Add(fhPtLambda1) ;
1244
1245 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >=0 )
1246 {
1247 fhPtLambda0NoTRD = new TH2F
1248 ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1249 fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1250 fhPtLambda0NoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1251 outputContainer->Add(fhPtLambda0NoTRD) ;
1252
1253 if(!fFillOnlySimpleSSHisto)
1254 {
1255 fhPtFracMaxCellNoTRD = new TH2F
1256 ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1257 fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1258 fhPtFracMaxCellNoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1259 outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1260 }
1261 }
1262
1263 if(!fFillOnlySimpleSSHisto)
1264 {
1265 fhPtDispersion = new TH2F
1266 ("hPtDispersion","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1267 fhPtDispersion->SetYTitle("D^{2}");
1268 fhPtDispersion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1269 outputContainer->Add(fhPtDispersion) ;
1270
1271 fhPtLambda0FracMaxCellCut = new TH2F
1272 ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1273 fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1274 fhPtLambda0FracMaxCellCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1275 outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1276
1277 fhPtFracMaxCell = new TH2F
1278 ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1279 fhPtFracMaxCell->SetYTitle("Fraction");
1280 fhPtFracMaxCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1281 outputContainer->Add(fhPtFracMaxCell) ;
1282
1283 fhPtDispEta = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1284 fhPtDispEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1285 fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1286 outputContainer->Add(fhPtDispEta);
1287
1288 fhPtDispPhi = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1289 fhPtDispPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1290 fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1291 outputContainer->Add(fhPtDispPhi);
1292
1293 fhPtSumEta = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs #it{p}_{T}", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1294 fhPtSumEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1295 fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1296 outputContainer->Add(fhPtSumEta);
1297
1298 fhPtSumPhi = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs #it{p}_{T}",
1299 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1300 fhPtSumPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1301 fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1302 outputContainer->Add(fhPtSumPhi);
1303
1304 fhPtSumEtaPhi = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",
1305 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1306 fhPtSumEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1307 fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1308 outputContainer->Add(fhPtSumEtaPhi);
1309
1310 fhPtDispEtaPhiDiff = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",
1311 nptbins,ptmin,ptmax,200, -10,10);
1312 fhPtDispEtaPhiDiff->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1313 fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1314 outputContainer->Add(fhPtDispEtaPhiDiff);
1315
1316 fhPtSphericity = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs #it{p}_{T} (GeV/#it{c})",
1317 nptbins,ptmin,ptmax, 200, -1,1);
1318 fhPtSphericity->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1319 fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1320 outputContainer->Add(fhPtSphericity);
1321
1322 for(Int_t i = 0; i < 7; i++)
1323 {
1324 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]),
1325 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1326 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1327 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1328 outputContainer->Add(fhDispEtaDispPhi[i]);
1329
1330 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]),
1331 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1332 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1333 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1334 outputContainer->Add(fhLambda0DispEta[i]);
1335
1336 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]),
1337 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1338 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1339 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1340 outputContainer->Add(fhLambda0DispPhi[i]);
1341
1342 }
1343 }
1344
1345 fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1346 nptbins,ptmin,ptmax,20,0,20);
1347 fhNLocMaxPt ->SetYTitle("N maxima");
1348 fhNLocMaxPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1349 outputContainer->Add(fhNLocMaxPt) ;
1350
1351 if(fFillAllNLMHistograms)
1352 {
1353 for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1354 {
1355 fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
1356 nptbins,ptmin,ptmax,20,0,20);
1357 fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
1358 fhNLocMaxPtSM[iSM] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1359 outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
1360 }
1361
1362 for (Int_t i = 0; i < 3; i++)
1363 {
1364 fhPtLambda0LocMax[i] = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1365 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1366 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1367 fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1368 fhPtLambda0LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1369 outputContainer->Add(fhPtLambda0LocMax[i]) ;
1370
1371 if(IsDataMC())
1372 {
1373 for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
1374 {
1375 fhMCPtLambda0LocMax[ipart][i] = new TH2F
1376 (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1377 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),ptype[ipart].Data()),
1378 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1379 fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1380 fhMCPtLambda0LocMax[ipart][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1381 outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1382 }
1383 }
1384
1385 fhPtLambda1LocMax[i] = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1386 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1387 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1388 fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1389 fhPtLambda1LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1390 outputContainer->Add(fhPtLambda1LocMax[i]) ;
1391
1392 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1393 {
1394 fhPtDispersionLocMax[i] = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1395 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1396 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1397 fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1398 fhPtDispersionLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1399 outputContainer->Add(fhPtDispersionLocMax[i]) ;
1400
1401 fhPtDispEtaLocMax[i] = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1402 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1403 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1404 fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1405 fhPtDispEtaLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1406 outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1407
1408 fhPtDispPhiLocMax[i] = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1409 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1410 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1411 fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1412 fhPtDispPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1413 outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1414
1415 fhPtSumEtaPhiLocMax[i] = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1416 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1417 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1418 fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1419 fhPtSumEtaPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1420 outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1421
1422 fhPtDispEtaPhiDiffLocMax[i] = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1423 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1424 nptbins,ptmin,ptmax,200, -10,10);
1425 fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1426 fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1427 outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1428
1429 fhPtSphericityLocMax[i] = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1430 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1431 nptbins,ptmin,ptmax,200, -1,1);
1432 fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1433 fhPtSphericityLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1434 outputContainer->Add(fhPtSphericityLocMax[i]) ;
1435 }
1436 }
1437 } // all NLM histos
1438
1439 fhPtNCells = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1440 fhPtNCells->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1441 fhPtNCells->SetYTitle("# of cells in cluster");
1442 outputContainer->Add(fhPtNCells);
1443
1444 fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1445 fhPtTime->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1446 fhPtTime->SetYTitle("t (ns)");
1447 outputContainer->Add(fhPtTime);
1448
1449 }
1450
1451 if(fAnaType != kIMCaloTracks)
1452 {
1453 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1454 fhEPairDiffTime->SetXTitle("#it{E}_{pair} (GeV)");
1455 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1456 outputContainer->Add(fhEPairDiffTime);
1457 }
1458
1459 if(fAnaType == kIMCalo)
1460 {
1461 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1462 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1463 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1464 "2 Local Maxima paired with more than 2 Local Maxima",
1465 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1466
1467 if(fFillAllNLMHistograms)
1468 {
1469 for (Int_t i = 0; i < 8 ; i++)
1470 {
1471 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1472
1473 fhMassPairLocMax[i] = new TH2F
1474 (Form("MassPairLocMax%s",combiName[i].Data()),
1475 Form("#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1476 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1477 fhMassPairLocMax[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1478 fhMassPairLocMax[i]->SetXTitle("#it{E}_{pair} (GeV)");
1479 outputContainer->Add(fhMassPairLocMax[i]) ;
1480 }
1481 }
1482 }
1483
1484 if(fFillTMHisto)
1485 {
1486 fhTrackMatchedDEta = new TH2F
1487 ("hTrackMatchedDEta",
1488 "d#eta of cluster-track vs cluster #it{p}_{T}",
1489 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1490 fhTrackMatchedDEta->SetYTitle("d#eta");
1491 fhTrackMatchedDEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1492
1493 fhTrackMatchedDPhi = new TH2F
1494 ("hTrackMatchedDPhi",
1495 "d#phi of cluster-track vs cluster #it{p}_{T}",
1496 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1497 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1498 fhTrackMatchedDPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1499
1500 fhTrackMatchedDEtaDPhi = new TH2F
1501 ("hTrackMatchedDEtaDPhi",
1502 "d#eta vs d#phi of cluster-track",
1503 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1504 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1505 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1506
1507 outputContainer->Add(fhTrackMatchedDEta) ;
1508 outputContainer->Add(fhTrackMatchedDPhi) ;
1509 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1510
1511 fhTrackMatchedDEtaPos = new TH2F
1512 ("hTrackMatchedDEtaPos",
1513 "d#eta of cluster-track vs cluster #it{p}_{T}",
1514 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1515 fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1516 fhTrackMatchedDEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1517
1518 fhTrackMatchedDPhiPos = new TH2F
1519 ("hTrackMatchedDPhiPos",
1520 "d#phi of cluster-track vs cluster #it{p}_{T}",
1521 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1522 fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1523 fhTrackMatchedDPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1524
1525 fhTrackMatchedDEtaDPhiPos = new TH2F
1526 ("hTrackMatchedDEtaDPhiPos",
1527 "d#eta vs d#phi of cluster-track",
1528 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1529 fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1530 fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1531
1532 outputContainer->Add(fhTrackMatchedDEtaPos) ;
1533 outputContainer->Add(fhTrackMatchedDPhiPos) ;
1534 outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1535
1536 fhTrackMatchedDEtaNeg = new TH2F
1537 ("hTrackMatchedDEtaNeg",
1538 "d#eta of cluster-track vs cluster #it{p}_{T}",
1539 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1540 fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1541 fhTrackMatchedDEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1542
1543 fhTrackMatchedDPhiNeg = new TH2F
1544 ("hTrackMatchedDPhiNeg",
1545 "d#phi of cluster-track vs cluster #it{p}_{T}",
1546 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1547 fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1548 fhTrackMatchedDPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1549
1550 fhTrackMatchedDEtaDPhiNeg = new TH2F
1551 ("hTrackMatchedDEtaDPhiNeg",
1552 "d#eta vs d#phi of cluster-track",
1553 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1554 fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1555 fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1556
1557 outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1558 outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1559 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1560
1561 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1562 fhdEdx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1563 fhdEdx->SetYTitle("<#it{dE}/#it{dx}>");
1564 outputContainer->Add(fhdEdx);
1565
1566 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1567 fhEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1568 fhEOverP->SetYTitle("#it{E}/#it{p}");
1569 outputContainer->Add(fhEOverP);
1570
1571 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >=0)
1572 {
1573 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1574 fhEOverPNoTRD->SetXTitle("#it{E} (GeV)");
1575 fhEOverPNoTRD->SetYTitle("#it{E}/#it{p}");
1576 outputContainer->Add(fhEOverPNoTRD);
1577 }
1578
1579 if(IsDataMC() && fFillTMHisto)
1580 {
1581 fhTrackMatchedMCParticlePt = new TH2F
1582 ("hTrackMatchedMCParticlePt",
1583 "Origin of particle vs energy",
1584 nptbins,ptmin,ptmax,8,0,8);
1585 fhTrackMatchedMCParticlePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1586 //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1587
1588 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1589 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1590 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1591 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1592 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1593 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1594 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1595 fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1596
1597 outputContainer->Add(fhTrackMatchedMCParticlePt);
1598
1599 fhTrackMatchedMCParticleDEta = new TH2F
1600 ("hTrackMatchedMCParticleDEta",
1601 "Origin of particle vs #eta residual",
1602 nresetabins,resetamin,resetamax,8,0,8);
1603 fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1604 //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1605
1606 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1607 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1608 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1609 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1610 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1611 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1612 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1613 fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1614
1615 outputContainer->Add(fhTrackMatchedMCParticleDEta);
1616
1617 fhTrackMatchedMCParticleDPhi = new TH2F
1618 ("hTrackMatchedMCParticleDPhi",
1619 "Origin of particle vs #phi residual",
1620 nresphibins,resphimin,resphimax,8,0,8);
1621 fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1622 //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1623
1624 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1625 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1626 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1627 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1628 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1629 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1630 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1631 fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1632
1633 outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1634 }
1635 }
1636
1637 if(fFillWeightHistograms)
1638 {
1639 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1640 nptbins,ptmin,ptmax, 100,0,1.);
1641 fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1642 fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
1643 outputContainer->Add(fhECellClusterRatio);
1644
1645 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1646 nptbins,ptmin,ptmax, 100,-10,0);
1647 fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1648 fhECellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1649 outputContainer->Add(fhECellClusterLogRatio);
1650
1651 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1652 nptbins,ptmin,ptmax, 100,0,1.);
1653 fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1654 fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
1655 outputContainer->Add(fhEMaxCellClusterRatio);
1656
1657 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1658 nptbins,ptmin,ptmax, 100,-10,0);
1659 fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1660 fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1661 outputContainer->Add(fhEMaxCellClusterLogRatio);
1662
1663 for(Int_t iw = 0; iw < 14; iw++)
1664 {
1665 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),
1666 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1667 fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1668 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1669 outputContainer->Add(fhLambda0ForW0[iw]);
1670
1671 // 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),
1672 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1673 // fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1674 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1675 // outputContainer->Add(fhLambda1ForW0[iw]);
1676
1677 }
1678 }
1679
1680 if(IsDataMC())
1681 {
1682 // Origin
1683
1684 fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1685 fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1686 fhMCPi0PtOrigin->SetYTitle("Origin");
1687 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1688 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1689 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1690 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1691 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1692 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1693 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1694 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1695 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1696 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1697 outputContainer->Add(fhMCPi0PtOrigin) ;
1698
1699 fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1700 fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1701 fhMCEtaPtOrigin->SetYTitle("Origin");
1702 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1703 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1704 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1705 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1706 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1707 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1708 outputContainer->Add(fhMCEtaPtOrigin) ;
1709
1710 fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1711 fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1712 fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1713 outputContainer->Add(fhMCPi0ProdVertex) ;
1714
1715 fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1716 fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1717 fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1718 outputContainer->Add(fhMCEtaProdVertex) ;
1719
1720 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1721 {
1722 fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #pi^{0} / E reco",
1723 nptbins,ptmin,ptmax,200,0,2);
1724 fhMCPi0PtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1725 fhMCPi0PtGenRecoFraction->SetYTitle("#it{E}^{#pi^{0} mother} / #it{E}^{rec}");
1726 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1727
1728 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),#it{p}_{T} versus E primary #eta / E reco",
1729 nptbins,ptmin,ptmax,200,0,2);
1730 fhMCEtaPtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1731 fhMCEtaPtGenRecoFraction->SetYTitle("#it{E}^{ #eta mother} / #it{E}^{rec}");
1732 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1733
1734 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1735 fhMCPi0DecayPt->SetYTitle("#it{N}");
1736 fhMCPi0DecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1737 outputContainer->Add(fhMCPi0DecayPt) ;
1738
1739 fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #gamma / #it{E} primary #pi^{0}",
1740 nptbins,ptmin,ptmax,100,0,1);
1741 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/#it{c})");
1742 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1743 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1744
1745 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1746 fhMCEtaDecayPt->SetYTitle("#it{N}");
1747 fhMCEtaDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1748 outputContainer->Add(fhMCEtaDecayPt) ;
1749
1750 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #gamma / E primary #eta",
1751 nptbins,ptmin,ptmax,100,0,1);
1752 fhMCEtaDecayPtFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1753 fhMCEtaDecayPtFraction->SetYTitle("#it{E}^{gen} / #it{E}^{gen-mother}");
1754 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1755
1756 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1757 fhMCOtherDecayPt->SetYTitle("#it{N}");
1758 fhMCOtherDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1759 outputContainer->Add(fhMCOtherDecayPt) ;
1760 }
1761
1762 if(fAnaType!=kSSCalo)
1763 {
1764 fhAnglePairMCPi0 = new TH2F
1765 ("AnglePairMCPi0",
1766 "Angle between decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1767 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1768 fhAnglePairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1769 outputContainer->Add(fhAnglePairMCPi0) ;
1770
1771 fhAnglePairMCEta = new TH2F
1772 ("AnglePairMCEta",
1773 "Angle between decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1774 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1775 fhAnglePairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1776 outputContainer->Add(fhAnglePairMCEta) ;
1777
1778 fhMassPairMCPi0 = new TH2F
1779 ("MassPairMCPi0",
1780 "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1781 fhMassPairMCPi0->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1782 fhMassPairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1783 outputContainer->Add(fhMassPairMCPi0) ;
1784
1785 fhMassPairMCEta = new TH2F
1786 ("MassPairMCEta",
1787 "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1788 fhMassPairMCEta->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1789 fhMassPairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1790 outputContainer->Add(fhMassPairMCEta) ;
1791 }
1792
1793 for(Int_t i = 0; i < fgkNmcTypes; i++)
1794 {
1795 fhMCE[i] = new TH1F
1796 (Form("hE_MC%s",pname[i].Data()),
1797 Form("Identified as #pi^{0} (#eta), cluster from %s",
1798 ptype[i].Data()),
1799 nptbins,ptmin,ptmax);
1800 fhMCE[i]->SetYTitle("#it{N}");
1801 fhMCE[i]->SetXTitle("#it{E} (GeV)");
1802 outputContainer->Add(fhMCE[i]) ;
1803
1804 fhMCPt[i] = new TH1F
1805 (Form("hPt_MC%s",pname[i].Data()),
1806 Form("Identified as #pi^{0} (#eta), cluster from %s",
1807 ptype[i].Data()),
1808 nptbins,ptmin,ptmax);
1809 fhMCPt[i]->SetYTitle("#it{N}");
1810 fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1811 outputContainer->Add(fhMCPt[i]) ;
1812
1813 if(fFillHighMultHistograms)
1814 {
1815 fhMCPtCentrality[i] = new TH2F
1816 (Form("hPtCentrality_MC%s",pname[i].Data()),
1817 Form("Identified as #pi^{0} (#eta), cluster from %s",
1818 ptype[i].Data()),
1819 nptbins,ptmin,ptmax, 100,0,100);
1820 fhMCPtCentrality[i]->SetYTitle("centrality");
1821 fhMCPtCentrality[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1822 outputContainer->Add(fhMCPtCentrality[i]) ;
1823 }
1824
1825 if(fAnaType == kSSCalo)
1826 {
1827 fhMCNLocMaxPt[i] = new TH2F
1828 (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1829 Form("cluster from %s, #it{p}_{T} of cluster vs NLM, accepted",ptype[i].Data()),
1830 nptbins,ptmin,ptmax,20,0,20);
1831 fhMCNLocMaxPt[i] ->SetYTitle("#it{NLM}");
1832 fhMCNLocMaxPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1833 outputContainer->Add(fhMCNLocMaxPt[i]) ;
1834
1835 fhMCNLocMaxPtReject[i] = new TH2F
1836 (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1837 Form("cluster from %s, #it{p}_{T} of cluster vs NLM, rejected",ptype[i].Data()),
1838 nptbins,ptmin,ptmax,20,0,20);
1839 fhMCNLocMaxPtReject[i] ->SetYTitle("#it{NLM}");
1840 fhMCNLocMaxPtReject[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1841 outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1842
1843 fhMCEReject[i] = new TH1F
1844 (Form("hEReject_MC%s",pname[i].Data()),
1845 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1846 ptype[i].Data()),
1847 nptbins,ptmin,ptmax);
1848 fhMCEReject[i]->SetYTitle("#it{N}");
1849 fhMCEReject[i]->SetXTitle("#it{E} (GeV)");
1850 outputContainer->Add(fhMCEReject[i]) ;
1851
1852 fhMCPtReject[i] = new TH1F
1853 (Form("hPtReject_MC%s",pname[i].Data()),
1854 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1855 ptype[i].Data()),
1856 nptbins,ptmin,ptmax);
1857 fhMCPtReject[i]->SetYTitle("#it{N}");
1858 fhMCPtReject[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1859 outputContainer->Add(fhMCPtReject[i]) ;
1860 }
1861
1862 fhMCPtPhi[i] = new TH2F
1863 (Form("hPtPhi_MC%s",pname[i].Data()),
1864 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1865 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1866 fhMCPtPhi[i]->SetYTitle("#phi");
1867 fhMCPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1868 outputContainer->Add(fhMCPtPhi[i]) ;
1869
1870 fhMCPtEta[i] = new TH2F
1871 (Form("hPtEta_MC%s",pname[i].Data()),
1872 Form("Identified as #pi^{0} (#eta), cluster from %s",
1873 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1874 fhMCPtEta[i]->SetYTitle("#eta");
1875 fhMCPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1876 outputContainer->Add(fhMCPtEta[i]) ;
1877
1878 fhMCMassPt[i] = new TH2F
1879 (Form("hMassPt_MC%s",pname[i].Data()),
1880 Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1881 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1882 fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1883 fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1884 outputContainer->Add(fhMCMassPt[i]) ;
1885
1886 fhMCSelectedMassPt[i] = new TH2F
1887 (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1888 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1889 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1890 fhMCSelectedMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1891 fhMCSelectedMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1892 outputContainer->Add(fhMCSelectedMassPt[i]) ;
1893
1894 if(fAnaType == kSSCalo)
1895 {
1896 fhMCMassPtNoOverlap[i] = new TH2F
1897 (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1898 Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1899 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1900 fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1901 fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1902 outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1903
1904 fhMCSelectedMassPtNoOverlap[i] = new TH2F
1905 (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1906 Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1907 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1908 fhMCSelectedMassPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1909 fhMCSelectedMassPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1910 outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1911 }
1912
1913 if( fFillSelectClHisto )
1914 {
1915 fhMCPtLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1916 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1917 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1918 fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1919 fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1920 outputContainer->Add(fhMCPtLambda0[i]) ;
1921
1922 fhMCPtLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1923 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1924 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1925 fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1926 fhMCPtLambda1[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1927 outputContainer->Add(fhMCPtLambda1[i]) ;
1928
1929 if(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0)
1930 {
1931 fhMCPtLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1932 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1933 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1934 fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1935 fhMCPtLambda0NoTRD[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1936 outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1937 }
1938
1939 if(!fFillOnlySimpleSSHisto)
1940 {
1941 fhMCPtDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1942 Form("Selected pair, cluster from %s : #it{p}_{T} vs dispersion^{2}",ptype[i].Data()),
1943 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1944 fhMCPtDispersion[i]->SetYTitle("#it{D}^{2}");
1945 fhMCPtDispersion[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1946 outputContainer->Add(fhMCPtDispersion[i]) ;
1947
1948 fhMCPtDispEta[i] = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1949 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1950 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1951 fhMCPtDispEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1952 fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1953 outputContainer->Add(fhMCPtDispEta[i]);
1954
1955 fhMCPtDispPhi[i] = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1956 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
1957 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1958 fhMCPtDispPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1959 fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1960 outputContainer->Add(fhMCPtDispPhi[i]);
1961
1962 fhMCPtSumEtaPhi[i] = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1963 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",ptype[i].Data()),
1964 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1965 fhMCPtSumEtaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1966 fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1967 outputContainer->Add(fhMCPtSumEtaPhi[i]);
1968
1969 fhMCPtDispEtaPhiDiff[i] = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1970 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",ptype[i].Data()),
1971 nptbins,ptmin,ptmax,200,-10,10);
1972 fhMCPtDispEtaPhiDiff[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1973 fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1974 outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1975
1976 fhMCPtSphericity[i] = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1977 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()),
1978 nptbins,ptmin,ptmax, 200,-1,1);
1979 fhMCPtSphericity[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1980 fhMCPtSphericity[i]->SetYTitle("#it{s} = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1981 outputContainer->Add(fhMCPtSphericity[i]);
1982
1983 for(Int_t ie = 0; ie < 7; ie++)
1984 {
1985 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1986 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
1987 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1988 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1989 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1990 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1991
1992 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1993 Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
1994 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1995 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1996 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1997 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1998
1999 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2000 Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
2001 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2002 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
2003 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2004 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
2005
2006 }
2007
2008 fhMCPtLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
2009 Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
2010 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2011 fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
2012 fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("#it{E} (GeV)");
2013 outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
2014
2015 fhMCPtFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
2016 Form("Selected pair, cluster from %s : #it{p}_{T} vs Max cell fraction of energy",ptype[i].Data()),
2017 nptbins,ptmin,ptmax,100,0,1);
2018 fhMCPtFracMaxCell[i]->SetYTitle("#it{Fraction}");
2019 fhMCPtFracMaxCell[i]->SetXTitle("#it{E} (GeV)");
2020 outputContainer->Add(fhMCPtFracMaxCell[i]) ;
2021
2022 }
2023 }
2024 }// MC particle loop
2025 }//Histos with MC
2026
2027 if(fAnaType==kSSCalo)
2028 {
2029 fhAsymmetry = new TH2F ("hAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
2030 nptbins,ptmin,ptmax, 200, -1,1);
2031 fhAsymmetry->SetXTitle("#it{E} (GeV)");
2032 fhAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2033 outputContainer->Add(fhAsymmetry);
2034
2035 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
2036 nptbins,ptmin,ptmax, 200, -1,1);
2037 fhSelectedAsymmetry->SetXTitle("#it{E} (GeV)");
2038 fhSelectedAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2039 outputContainer->Add(fhSelectedAsymmetry);
2040
2041 fhSplitE = new TH1F
2042 ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
2043 fhSplitE->SetYTitle("counts");
2044 fhSplitE->SetXTitle("#it{E} (GeV)");
2045 outputContainer->Add(fhSplitE) ;
2046
2047 fhSplitPt = new TH1F
2048 ("hSplitPt","Selected #pi^{0} (#eta) pairs #it{p}_{T} sum of split sub-clusters",nptbins,ptmin,ptmax);
2049 fhSplitPt->SetYTitle("counts");
2050 fhSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2051 outputContainer->Add(fhSplitPt) ;
2052
2053
2054 fhSplitPtPhi = new TH2F
2055 ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2056 fhSplitPtPhi->SetYTitle("#phi (rad)");
2057 fhSplitPtPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2058 outputContainer->Add(fhSplitPtPhi) ;
2059
2060 fhSplitPtEta = new TH2F
2061 ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2062 fhSplitPtEta->SetYTitle("#eta");
2063 fhSplitPtEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2064 outputContainer->Add(fhSplitPtEta) ;
2065
2066
2067 fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
2068 nptbins,ptmin,ptmax,20,0,20);
2069 fhNLocMaxSplitPt ->SetYTitle("#it{NLM}");
2070 fhNLocMaxSplitPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2071 outputContainer->Add(fhNLocMaxSplitPt) ;
2072
2073
2074 fhMassSplitPt = new TH2F
2075 ("hMassSplitPt","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
2076 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2077 fhMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2078 fhMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2079 outputContainer->Add(fhMassSplitPt) ;
2080
2081 fhSelectedMassSplitPt = new TH2F
2082 ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
2083 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2084 fhSelectedMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2085 fhSelectedMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2086 outputContainer->Add(fhSelectedMassSplitPt) ;
2087
2088 if(IsDataMC())
2089 {
2090 fhMassSplitPtNoOverlap = new TH2F
2091 ("hMassSplitPtNoOverlap","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2092 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2093 fhMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2094 fhMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2095 outputContainer->Add(fhMassSplitPtNoOverlap) ;
2096
2097 fhSelectedMassSplitPtNoOverlap = new TH2F
2098 ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2099 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2100 fhSelectedMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2101 fhSelectedMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2102 outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
2103
2104
2105 fhMCPi0PtRecoPtPrim = new TH2F
2106 ("hMCPi0PtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2107 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2108 fhMCPi0PtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2109 fhMCPi0PtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2110 outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
2111
2112 fhMCPi0PtRecoPtPrimNoOverlap = new TH2F
2113 ("hMCPi0PtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2114 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2115 fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2116 fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2117 outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
2118
2119 fhMCPi0SelectedPtRecoPtPrim = new TH2F
2120 ("hMCPi0SelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2121 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2122 fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2123 fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2124 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
2125
2126 fhMCPi0SelectedPtRecoPtPrimNoOverlap = new TH2F
2127 ("hMCPi0SelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2128 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2129 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2130 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2131 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
2132
2133
2134 fhMCPi0SplitPtRecoPtPrim = new TH2F
2135 ("hMCPi0SplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2136 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2137 fhMCPi0SplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2138 fhMCPi0SplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2139 outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
2140
2141 fhMCPi0SplitPtRecoPtPrimNoOverlap = new TH2F
2142 ("hMCPi0SplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2143 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2144 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2145 fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2146 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
2147
2148 fhMCPi0SelectedSplitPtRecoPtPrim = new TH2F
2149 ("hMCPi0SelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2150 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2151 fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2152 fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2153 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
2154
2155 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2156 ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2157 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2158 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2159 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2160 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
2161
2162 fhMCEtaPtRecoPtPrim = new TH2F
2163 ("hMCEtaPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2164 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2165 fhMCEtaPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2166 fhMCEtaPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2167 outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
2168
2169 fhMCEtaPtRecoPtPrimNoOverlap = new TH2F
2170 ("hMCEtaPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2171 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2172 fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2173 fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2174 outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2175
2176 fhMCEtaSelectedPtRecoPtPrim = new TH2F
2177 ("hMCEtaSelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2178 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2179 fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2180 fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2181 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2182
2183 fhMCEtaSelectedPtRecoPtPrimNoOverlap = new TH2F
2184 ("hMCEtaSelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2185 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2186 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2187 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2188 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2189
2190
2191 fhMCEtaSplitPtRecoPtPrim = new TH2F
2192 ("hMCEtaSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2193 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2194 fhMCEtaSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2195 fhMCEtaSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2196 outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2197
2198 fhMCEtaSplitPtRecoPtPrimNoOverlap = new TH2F
2199 ("hMCEtaSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2200 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2201 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2202 fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2203 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2204
2205 fhMCEtaSelectedSplitPtRecoPtPrim = new TH2F
2206 ("hMCEtaSelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2207 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2208 fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2209 fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2210 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2211
2212 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap = new TH2F
2213 ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2214 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2215 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2216 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2217 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2218
2219
2220 for(Int_t inlm = 0; inlm < 3; inlm++)
2221 {
2222 fhMCPi0PtRecoPtPrimLocMax[inlm] = new TH2F
2223 (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2224 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2225 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2226 fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2227 outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
2228
2229 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2230 (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2231 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2232 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2233 fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2234 outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
2235
2236 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] = new TH2F
2237 (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2238 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2239 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2240 fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2241 outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
2242
2243 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2244 (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2245 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2246 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2247 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2248 outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2249
2250 fhMCEtaPtRecoPtPrimLocMax[inlm] = new TH2F
2251 (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2252 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2253 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2254 fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2255 outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
2256
2257 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] = new TH2F
2258 (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2259 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2260 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2261 fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2262 outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
2263
2264 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2265 (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2266 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2267 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2268 fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2269 outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
2270
2271 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] = new TH2F
2272 (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2273 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2274 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2275 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2276 outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2277
2278 }
2279
2280 for(Int_t i = 0; i < fgkNmcTypes; i++)
2281 {
2282 fhMCPtAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2283 Form("cluster from %s : #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",ptype[i].Data()),
2284 nptbins,ptmin,ptmax, 200,-1,1);
2285 fhMCPtAsymmetry[i]->SetXTitle("#it{E} (GeV)");
2286 fhMCPtAsymmetry[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2287 outputContainer->Add(fhMCPtAsymmetry[i]);
2288
2289 fhMCSplitE[i] = new TH1F
2290 (Form("hSplitE_MC%s",pname[i].Data()),
2291 Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2292 nptbins,ptmin,ptmax);
2293 fhMCSplitE[i]->SetYTitle("counts");
2294 fhMCSplitE[i]->SetXTitle("#it{E} (GeV)");
2295 outputContainer->Add(fhMCSplitE[i]) ;
2296
2297 fhMCSplitPt[i] = new TH1F
2298 (Form("hSplitPt_MC%s",pname[i].Data()),
2299 Form("cluster from %s, #it{p}_{T} sum of split sub-clusters",ptype[i].Data()),
2300 nptbins,ptmin,ptmax);
2301 fhMCSplitPt[i]->SetYTitle("counts");
2302 fhMCSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2303 outputContainer->Add(fhMCSplitPt[i]) ;
2304
2305
2306 fhMCSplitPtPhi[i] = new TH2F
2307 (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2308 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2309 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2310 fhMCSplitPtPhi[i]->SetYTitle("#phi");
2311 fhMCSplitPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2312 outputContainer->Add(fhMCSplitPtPhi[i]) ;
2313
2314 fhMCSplitPtEta[i] = new TH2F
2315 (Form("hSplitPtEta_MC%s",pname[i].Data()),
2316 Form("Identified as #pi^{0} (#eta), cluster from %s",
2317 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2318 fhMCSplitPtEta[i]->SetYTitle("#eta");
2319 fhMCSplitPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2320 outputContainer->Add(fhMCSplitPtEta[i]) ;
2321
2322
2323 fhMCNLocMaxSplitPt[i] = new TH2F
2324 (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2325 Form("cluster from %s, #it{p}_{T} sum of split sub-clusters, for NLM",ptype[i].Data()),
2326 nptbins,ptmin,ptmax,20,0,20);
2327 fhMCNLocMaxSplitPt[i] ->SetYTitle("#it{NLM}");
2328 fhMCNLocMaxSplitPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2329 outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2330
2331 fhMCMassSplitPt[i] = new TH2F
2332 (Form("hMassSplitPt_MC%s",pname[i].Data()),
2333 Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2334 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2335 fhMCMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2336 fhMCMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2337 outputContainer->Add(fhMCMassSplitPt[i]) ;
2338
2339 fhMCSelectedMassSplitPt[i] = new TH2F
2340 (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2341 Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2342 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2343 fhMCSelectedMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2344 fhMCSelectedMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2345 outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2346
2347 fhMCMassSplitPtNoOverlap[i] = new TH2F
2348 (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2349 Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2350 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2351 fhMCMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2352 fhMCMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2353 outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2354
2355 fhMCSelectedMassSplitPtNoOverlap[i] = new TH2F
2356 (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2357 Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2358 nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2359 fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2360 fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2361 outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2362 }
2363 }
2364 }
2365
2366 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2367 {
2368 for(Int_t i = 0; i< 3; i++)
2369 {
2370 fhPtAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2371 Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ), %s",nlm[i].Data()),
2372 nptbins,ptmin,ptmax,200, -1,1);
2373 fhPtAsymmetryLocMax[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2374 fhPtAsymmetryLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2375 outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2376 }
2377
2378 for(Int_t ie = 0; ie < 7; ie++)
2379 {
2380
2381 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2382 Form("#lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2383 ssbins,ssmin,ssmax , 200,-1,1);
2384 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2385 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2386 outputContainer->Add(fhAsymmetryLambda0[ie]);
2387
2388 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2389 Form("#sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2390 ssbins,ssmin,ssmax , 200,-1,1);
2391 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2392 fhAsymmetryDispEta[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2393 outputContainer->Add(fhAsymmetryDispEta[ie]);
2394
2395 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2396 Form("#sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2397 ssbins,ssmin,ssmax , 200,-1,1);
2398 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2399 fhAsymmetryDispPhi[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2400 outputContainer->Add(fhAsymmetryDispPhi[ie]);
2401 }
2402
2403
2404 if(IsDataMC())
2405 {
2406 for(Int_t i = 0; i < fgkNmcTypes; i++)
2407 {
2408 for(Int_t ie = 0; ie < 7; ie++)
2409 {
2410 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2411 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
2412 ssbins,ssmin,ssmax , 200,-1,1);
2413 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2414 fhMCAsymmetryLambda0[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2415 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2416
2417 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2418 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
2419 ssbins,ssmin,ssmax , 200,-1,1);
2420 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2421 fhMCAsymmetryDispEta[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2422 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2423
2424 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2425 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
2426 ssbins,ssmin,ssmax , 200,-1,1);
2427 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2428 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2429 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2430 }
2431 }
2432 }
2433 }
2434
2435 if(fFillPileUpHistograms)
2436 {
2437
2438 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2439
2440 for(Int_t i = 0 ; i < 7 ; i++)
2441 {
2442 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2443 Form("Selected #pi^{0} (#eta) #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2444 fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2445 outputContainer->Add(fhPtPileUp[i]);
2446
2447 fhPtCellTimePileUp[i] = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2448 Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2449 nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2450 fhPtCellTimePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2451 fhPtCellTimePileUp[i]->SetYTitle("#it{t}_{cell} (ns)");
2452 outputContainer->Add(fhPtCellTimePileUp[i]);
2453
2454 fhPtTimeDiffPileUp[i] = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2455 Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2456 nptbins,ptmin,ptmax,400,-200,200);
2457 fhPtTimeDiffPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2458 fhPtTimeDiffPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
2459 outputContainer->Add(fhPtTimeDiffPileUp[i]);
2460
2461 }
2462
2463 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","#it{t} of cluster vs #it{E} of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2464 fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2465 fhTimePtNoCut->SetYTitle("#it{t} (ns)");
2466 outputContainer->Add(fhTimePtNoCut);
2467
2468 fhTimePtSPD = new TH2F ("hTimePt_SPD","#it{t} of cluster vs #it{E} of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2469 fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2470 fhTimePtSPD->SetYTitle("#it{t} (ns)");
2471 outputContainer->Add(fhTimePtSPD);
2472
2473 fhTimePtSPDMulti = new TH2F ("hTimePt_SPDMulti","time of cluster vs #it{E} of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2474 fhTimePtSPDMulti->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2475 fhTimePtSPDMulti->SetYTitle("#it{t} (ns)");
2476 outputContainer->Add(fhTimePtSPDMulti);
2477
2478 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","#it{t} of cluster vs #it{N} pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2479 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2480 fhTimeNPileUpVertSPD->SetXTitle("#it{t} (ns)");
2481 outputContainer->Add(fhTimeNPileUpVertSPD);
2482
2483 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","#it{t} of cluster vs #it{N} pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2484 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2485 fhTimeNPileUpVertTrack->SetXTitle("#it{t} (ns)");
2486 outputContainer->Add(fhTimeNPileUpVertTrack);
2487
2488 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","#it{t} of cluster vs #it{N} constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2489 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2490 fhTimeNPileUpVertContributors->SetXTitle("#it{t} (ns)");
2491 outputContainer->Add(fhTimeNPileUpVertContributors);
2492
2493 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
2494 fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{Z} (cm) ");
2495 fhTimePileUpMainVertexZDistance->SetXTitle("#it{t} (ns)");
2496 outputContainer->Add(fhTimePileUpMainVertexZDistance);
2497
2498 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
2499 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{Z} (cm) ");
2500 fhTimePileUpMainVertexZDiamond->SetXTitle("#it{t} (ns)");
2501 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2502
2503 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","#it{p}_{T} of cluster vs #it{N} pile-up SPD vertex",
2504 nptbins,ptmin,ptmax,20,0,20);
2505 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2506 fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2507 outputContainer->Add(fhPtNPileUpSPDVtx);
2508
2509 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","#it{p}_{T} of cluster vs #it{N} pile-up Tracks vertex",
2510 nptbins,ptmin,ptmax, 20,0,20 );
2511 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2512 fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2513 outputContainer->Add(fhPtNPileUpTrkVtx);
2514
2515 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","#it{p}_{T} of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2516 nptbins,ptmin,ptmax,20,0,20);
2517 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2518 fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2519 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2520
2521 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","#it{p}_{T} of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2522 nptbins,ptmin,ptmax, 20,0,20 );
2523 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2524 fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2525 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2526
2527 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","#it{p}_{T} of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2528 nptbins,ptmin,ptmax,20,0,20);
2529 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2530 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2531 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2532
2533 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","#it{p}_{T} of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2534 nptbins,ptmin,ptmax, 20,0,20 );
2535 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2536 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2537 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2538
2539 }
2540
2541 //Keep neutral meson selection histograms if requiered
2542 //Setting done in AliNeutralMesonSelection
2543
2544 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2545 {
2546 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2547
2548 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2549 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2550
2551 delete nmsHistos;
2552 }
2553
2554 return outputContainer ;
2555
2556}
2557
2558//_____________________________________________
2559Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2560{
2561 // Assign mc index depending on MC bit set
2562
2563 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
2564 {
2565 return kmcPi0 ;
2566 }//pi0
2567 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
2568 {
2569 return kmcEta ;
2570 }//eta
2571 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) ||
2572 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) ||
2573 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
2574 {
2575 return kmcPhoton ;
2576 }//direct photon
2577 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2578 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) )
2579 {
2580 return kmcPi0Decay ;
2581 }//decay photon from pi0
2582 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2583 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) )
2584 {
2585 return kmcEtaDecay ;
2586 }//decay photon from eta
2587 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2588 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) )
2589 {
2590 return kmcOtherDecay ;
2591 }//decay photon from other than eta or pi0
2592 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2593 {
2594 return kmcElectron ;
2595 }//electron
2596 else
2597 {
2598 return kmcHadron ;
2599 }//other particles
2600
2601}
2602
2603//__________________________________________________________________
2604void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2605 AliAODPWG4Particle * photon2,
2606 Int_t & label, Int_t & tag)
2607{
2608 // Check the labels of pare in case mother was same pi0 or eta
2609 // Set the new AOD accordingly
2610
2611 Int_t label1 = photon1->GetLabel();
2612 Int_t label2 = photon2->GetLabel();
2613
2614 if(label1 < 0 || label2 < 0 ) return ;
2615
2616 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2617 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2618 Int_t tag1 = photon1->GetTag();
2619 Int_t tag2 = photon2->GetTag();
2620
2621 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2622 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2623 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
2624 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2625 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
2626 )
2627 {
2628
2629 //Check if pi0/eta mother is the same
2630 if(GetReader()->ReadStack())
2631 {
2632 if(label1>=0)
2633 {
2634 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2635 label1 = mother1->GetFirstMother();
2636 //mother1 = GetMCStack()->Particle(label1);//pi0
2637 }
2638 if(label2>=0)
2639 {
2640 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2641 label2 = mother2->GetFirstMother();
2642 //mother2 = GetMCStack()->Particle(label2);//pi0
2643 }
2644 } // STACK
2645 else if(GetReader()->ReadAODMCParticles())
2646 {//&& (input > -1)){
2647 if(label1>=0)
2648 {
2649 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2650 label1 = mother1->GetMother();
2651 //mother1 = GetMCStack()->Particle(label1);//pi0
2652 }
2653 if(label2>=0)
2654 {
2655 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2656 label2 = mother2->GetMother();
2657 //mother2 = GetMCStack()->Particle(label2);//pi0
2658 }
2659 }// AOD
2660
2661 //printf("mother1 %d, mother2 %d\n",label1,label2);
2662 if( label1 == label2 && label1>=0 )
2663 {
2664 label = label1;
2665
2666 TLorentzVector mom1 = *(photon1->Momentum());
2667 TLorentzVector mom2 = *(photon2->Momentum());
2668
2669 Double_t angle = mom2.Angle(mom1.Vect());
2670 Double_t mass = (mom1+mom2).M();
2671 Double_t epair = (mom1+mom2).E();
2672
2673 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2674 {
2675 fhMassPairMCPi0 ->Fill(epair,mass);
2676 fhAnglePairMCPi0->Fill(epair,angle);
2677 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2678 }
2679 else
2680 {
2681 fhMassPairMCEta ->Fill(epair,mass);
2682 fhAnglePairMCEta->Fill(epair,angle);
2683 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2684 }
2685
2686 } // same label
2687 } // both from eta or pi0 decay
2688
2689}
2690
2691//____________________________________________________________________________
2692void AliAnaPi0EbE::Init()
2693{
2694 //Init
2695 //Do some checks
2696 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2697 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2698 abort();
2699 }
2700 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2701 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2702 abort();
2703 }
2704
2705}
2706
2707//____________________________________________________________________________
2708void AliAnaPi0EbE::InitParameters()
2709{
2710 //Initialize the parameters of the analysis.
2711 AddToHistogramsName("AnaPi0EbE_");
2712
2713 fInputAODGammaConvName = "PhotonsCTS" ;
2714 fAnaType = kIMCalo ;
2715 fCalorimeter = "EMCAL" ;
2716 fMinDist = 2.;
2717 fMinDist2 = 4.;
2718 fMinDist3 = 5.;
2719
2720 fNLMECutMin[0] = 10.;
2721 fNLMECutMin[1] = 6. ;
2722 fNLMECutMin[2] = 6. ;
2723}
2724
2725//__________________________________________________________________
2726void AliAnaPi0EbE::MakeAnalysisFillAOD()
2727{
2728 //Do analysis and fill aods
2729
2730 switch(fAnaType)
2731 {
2732 case kIMCalo:
2733 MakeInvMassInCalorimeter();
2734 break;
2735
2736 case kSSCalo:
2737 MakeShowerShapeIdentification();
2738 break;
2739
2740 case kIMCaloTracks:
2741 MakeInvMassInCalorimeterAndCTS();
2742 break;
2743
2744 }
2745}
2746
2747//____________________________________________
2748void AliAnaPi0EbE::MakeInvMassInCalorimeter()
2749{
2750 //Do analysis and fill aods
2751 //Search for the photon decay in calorimeters
2752 //Read photon list from AOD, produced in class AliAnaPhoton
2753 //Check if 2 photons have the mass of the pi0.
2754
2755 TLorentzVector mom1;
2756 TLorentzVector mom2;
2757 TLorentzVector mom ;
2758
2759 Int_t tag = 0;
2760 Int_t label = 0;
2761
2762 if(!GetInputAODBranch())
2763 {
2764 AliFatal(Form("No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data()));
2765 return; // coverity
2766 }
2767
2768 //Get shower shape information of clusters
2769 TObjArray *clusters = 0;
2770 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2771 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
2772
2773 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++)
2774 {
2775 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2776
2777 // Vertex cut in case of mixed events
2778 Int_t evtIndex1 = 0 ;
2779 if(GetMixedEvent())
2780 {
2781 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2782 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
2783 }
2784
2785 mom1 = *(photon1->Momentum());
2786
2787 //Get original cluster, to recover some information
2788 Int_t iclus = -1;
2789 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2790
2791 if(!cluster1)
2792 {
2793 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2794 return;
2795 }
2796
2797 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2798 {
2799 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2800
2801 // Do analysis only when one of the decays is isolated
2802 // Run AliAnaParticleIsolation before
2803 if(fSelectIsolatedDecay)
2804 {
2805 Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
2806 Bool_t isolated2 = ((AliAODPWG4ParticleCorrelation*) photon2)->IsIsolated();
2807 if(!isolated1 && !isolated2) continue;
2808 }
2809
2810 // Vertex cut in case of mixed events
2811 Int_t evtIndex2 = 0 ;
2812 if(GetMixedEvent())
2813 {
2814 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2815
2816 if(evtIndex1 == evtIndex2)
2817 continue ;
2818
2819 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
2820 }
2821
2822 mom2 = *(photon2->Momentum());
2823
2824 //Get original cluster, to recover some information
2825 Int_t iclus2 = -1;
2826 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2827 // start new loop from iclus1+1 to gain some time
2828
2829 if(!cluster2)
2830 {
2831 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2832 return;
2833 }
2834
2835 Float_t e1 = photon1->E();
2836 Float_t e2 = photon2->E();
2837
2838 //Select clusters with good time window difference
2839 Float_t tof1 = cluster1->GetTOF()*1e9;;
2840 Float_t tof2 = cluster2->GetTOF()*1e9;;
2841 Double_t t12diff = tof1-tof2;
2842 fhEPairDiffTime->Fill(e1+e2, t12diff);
2843 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2844
2845 //Play with the MC stack if available
2846 Int_t mcIndex = kmcHadron;
2847 if(IsDataMC())
2848 {
2849 HasPairSameMCMother(photon1, photon2, label, tag) ;
2850 mcIndex = GetMCIndex(tag);
2851 }
2852
2853 // Check the invariant mass for different selection on the local maxima
2854 // Name of AOD method TO BE FIXED
2855 Int_t nMaxima1 = photon1->GetFiducialArea();
2856 Int_t nMaxima2 = photon2->GetFiducialArea();
2857
2858 mom = mom1+mom2;
2859
2860 Double_t mass = mom.M();
2861 Double_t epair = mom.E();
2862 Float_t ptpair = mom.Pt();
2863
2864 if(fFillAllNLMHistograms)
2865 {
2866 if(nMaxima1==nMaxima2)
2867 {
2868 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2869 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2870 else fhMassPairLocMax[2]->Fill(epair,mass);
2871 }
2872 else if(nMaxima1==1 || nMaxima2==1)
2873 {
2874 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2875 else fhMassPairLocMax[4]->Fill(epair,mass);
2876 }
2877 else
2878 fhMassPairLocMax[5]->Fill(epair,mass);
2879
2880 // combinations with SS axis cut and NLM cut
2881 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2882 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2883 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2884 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2885 }
2886
2887 //
2888 // Skip events with too few or too many NLM
2889 //
2890 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2891
2892 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2893
2894 //Mass of all pairs
2895 fhMass ->Fill( epair,mass);
2896 fhMassPt->Fill(ptpair,mass);
2897 if(IsDataMC()) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
2898
2899 //
2900 // Select good pair (good phi, pt cuts, aperture and invariant mass)
2901 //
2902 if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)) continue;
2903
2904 if(GetDebug()>1)
2905 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",
2906 mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2907
2908 //
2909 // Tag both photons as decay if not done before
2910 // set the corresponding bit for pi0 or eta or "side" case
2911 //
2912 Int_t bit1 = photon1->GetBtag(); // temporary
2913 if( bit1 < 0 ) bit1 = 0 ; // temporary
2914 if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
2915 {
2916 if( GetDebug() > 1 )
2917 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT1 %2.2f; bit requested %d; decay bit1: In %d, ",
2918 mom1.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit1);
2919
2920 photon1->SetTagged(kTRUE); // temporary
2921 GetNeutralMesonSelection()->SetDecayBit(bit1);
2922 photon1->SetBtag(bit1); // temporary
2923
2924 if( GetDebug() > 1 )
2925 printf("Out %d \n", bit1);
2926
2927 fhPtDecay->Fill(photon1->Pt());
2928
2929 //Fill some histograms about shower shape
2930 if(fFillSelectClHisto && cluster1 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2931 FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2932
2933 if(IsDataMC())
2934 {
2935 Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
2936 fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
2937 }
2938 }
2939
2940 Int_t bit2 = photon2->GetBtag(); // temporary
2941 if( bit2 < 0 ) bit2 = 0 ; // temporary
2942 if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
2943 {
2944 if( GetDebug() > 1 )
2945 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT2 %2.2f; bit requested %d; decay bit2: In %d, ",
2946 mom2.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit2);
2947
2948 photon2->SetTagged(kTRUE); // temporary
2949 GetNeutralMesonSelection()->SetDecayBit(bit2);
2950 photon2->SetBtag(bit2); // temporary
2951
2952 if( GetDebug() > 1 )
2953 printf("Out %d \n", bit2);
2954
2955 fhPtDecay->Fill(photon2->Pt());
2956
2957 //Fill some histograms about shower shape
2958 if(fFillSelectClHisto && cluster2 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2959 FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2960
2961 if(IsDataMC())
2962 {
2963 Int_t mcIndex2 = GetMCIndex(photon2->GetTag());
2964 fhMCPtDecay[mcIndex2]->Fill(photon2->Pt());
2965 }
2966 }
2967
2968 //Mass of selected pairs
2969 fhSelectedMass ->Fill( epair,mass);
2970 fhSelectedMassPt->Fill(ptpair,mass);
2971 if(IsDataMC())fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
2972
2973 // Fill histograms to undertand pile-up before other cuts applied
2974 // Remember to relax time cuts in the reader
2975 FillPileUpHistograms(ptpair,((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2976
2977 //Create AOD for analysis
2978 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2979
2980 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2981 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2982 else
2983 {
2984 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2985 return ;
2986 }
2987 pi0.SetDetector(photon1->GetDetector());
2988
2989 // MC
2990 pi0.SetLabel(label);
2991 pi0.SetTag(tag);
2992
2993 //Set the indeces of the original caloclusters
2994 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2995 //pi0.SetInputFileIndex(input);
2996
2997 AddAODParticle(pi0);
2998
2999 }//2n photon loop
3000
3001 }//1st photon loop
3002
3003 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
3004
3005}
3006
3007//__________________________________________________
3008void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
3009{
3010 //Do analysis and fill aods
3011 //Search for the photon decay in calorimeters
3012 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
3013 //Check if 2 photons have the mass of the pi0.
3014
3015 TLorentzVector mom1;
3016 TLorentzVector mom2;
3017 TLorentzVector mom ;
3018 Int_t tag = 0;
3019 Int_t label = 0;
3020 Int_t evtIndex = 0;
3021
3022 // Check calorimeter input
3023 if(!GetInputAODBranch())
3024 {
3025 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
3026 abort();
3027 }
3028
3029 // Get the array with conversion photons
3030 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
3031 if(!inputAODGammaConv)
3032 {
3033 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
3034
3035 if(!inputAODGammaConv)
3036 {
3037 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
3038
3039 return;
3040 }
3041 }
3042
3043 //Get shower shape information of clusters
3044 TObjArray *clusters = 0;
3045 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
3046 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
3047
3048 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
3049 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
3050 if(nCTS<=0 || nCalo <=0)
3051 {
3052 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
3053 return;
3054 }
3055
3056 if(GetDebug() > 1)
3057 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
3058
3059 // Do the loop, first calo, second CTS
3060 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
3061 {
3062 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
3063 mom1 = *(photon1->Momentum());
3064
3065 // Do analysis only when one of the decays is isolated
3066 // Run AliAnaParticleIsolation before
3067 if(fSelectIsolatedDecay)
3068 {
3069 Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
3070 if(!isolated1) continue;
3071 }
3072
3073 //Get original cluster, to recover some information
3074 Int_t iclus = -1;
3075 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
3076
3077 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++)
3078 {
3079 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
3080
3081 if(GetMixedEvent())
3082 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
3083 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
3084
3085 mom2 = *(photon2->Momentum());
3086
3087 mom = mom1+mom2;
3088
3089 Double_t mass = mom.M();
3090 Double_t epair = mom.E();
3091 Float_t ptpair = mom.Pt();
3092
3093 Int_t nMaxima = photon1->GetFiducialArea();
3094 if(fFillAllNLMHistograms)
3095 {
3096 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
3097 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
3098 else fhMassPairLocMax[2]->Fill(epair,mass);
3099 }
3100
3101 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
3102 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
3103
3104 //Play with the MC stack if available
3105 Int_t mcIndex = kmcHadron;
3106 if(IsDataMC())
3107 {
3108 Int_t label2 = photon2->GetLabel();
3109 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
3110
3111 HasPairSameMCMother(photon1, photon2, label, tag) ;
3112 mcIndex = GetMCIndex(tag);
3113 }
3114
3115 //Mass of selected pairs
3116 fhMass ->Fill( epair,mass);
3117 fhMassPt->Fill(ptpair,mass);
3118 if(IsDataMC()) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
3119
3120 //
3121 // Select good pair (good phi, pt cuts, aperture and invariant mass)
3122 //
3123 if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)) continue ;
3124
3125 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",
3126 mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
3127
3128 //
3129 // Tag both photons as decay if not done before
3130 // set the corresponding bit for pi0 or eta or "side" case
3131 //
3132 Int_t bit1 = photon1->GetBtag(); // temporary
3133 if( bit1 < 0 ) bit1 = 0 ; // temporary
3134 if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
3135 {
3136 photon1->SetTagged(kTRUE); // temporary
3137 GetNeutralMesonSelection()->SetDecayBit(bit1);
3138 photon1->SetBtag(bit1); // temporary
3139 fhPtDecay->Fill(photon1->Pt());
3140
3141 //Fill some histograms about shower shape
3142 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3143 FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
3144
3145 if(IsDataMC())
3146 {
3147 Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
3148 fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
3149 }
3150 }
3151
3152 Int_t bit2 = photon2->GetBtag(); // temporary
3153 if( bit2 < 0 ) bit2 = 0 ; // temporary
3154 if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
3155 {
3156 photon2->SetTagged(kTRUE); // temporary
3157 GetNeutralMesonSelection()->SetDecayBit(bit2);
3158 photon2->SetBtag(bit2); // temporary
3159 }
3160
3161 //Mass of selected pairs
3162 fhSelectedMass ->Fill( epair,mass);
3163 fhSelectedMassPt->Fill(ptpair,mass);
3164 if(IsDataMC()) fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
3165
3166 // Fill histograms to undertand pile-up before other cuts applied
3167 // Remember to relax time cuts in the reader
3168 if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
3169
3170 //Create AOD for analysis
3171
3172 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
3173
3174 if ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
3175 else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
3176 else
3177 {
3178 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Particle type declared in AliNeutralMeson not correct, do not add \n");
3179 return ;
3180 }
3181 pi0.SetDetector(photon1->GetDetector());
3182
3183 // MC
3184 pi0.SetLabel(label);
3185 pi0.SetTag(tag);
3186
3187 //Set the indeces of the original tracks or caloclusters
3188 pi0.SetCaloLabel (photon1->GetCaloLabel(0) , -1);
3189 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
3190 //pi0.SetInputFileIndex(input);
3191
3192 AddAODParticle(pi0);
3193
3194 }//2n photon loop
3195
3196 }//1st photon loop
3197
3198 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
3199
3200}
3201
3202//_________________________________________________
3203void AliAnaPi0EbE::MakeShowerShapeIdentification()
3204{
3205 //Search for pi0 in fCalorimeter with shower shape analysis
3206
3207 TObjArray * pl = 0x0;
3208 AliVCaloCells * cells = 0x0;
3209 //Select the Calorimeter of the photon
3210 if (fCalorimeter == "EMCAL" )
3211 {
3212 pl = GetEMCALClusters();
3213 cells = GetEMCALCells();
3214 }
3215 else if (fCalorimeter == "PHOS")
3216 {
3217 AliFatal("kSSCalo case not implememted for PHOS");
3218 return; // for coverity
3219
3220 //pl = GetPHOSClusters();
3221 //cells = GetPHOSCells();
3222 }
3223
3224 if(!pl)
3225 {
3226 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3227 return;
3228 }
3229
3230 TLorentzVector mom ;
3231 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3232 {
3233 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3234
3235 Int_t evtIndex = 0 ;
3236 if (GetMixedEvent())
3237 {
3238 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3239
3240 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
3241 }
3242
3243 //Get Momentum vector,
3244 Double_t vertex[]={0,0,0};
3245 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3246 {
3247 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
3248 }//Assume that come from vertex in straight line
3249 else
3250 {
3251 calo->GetMomentum(mom,vertex) ;
3252 }
3253
3254 //If too small or big pt, skip it
3255 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
3256
3257 //Check acceptance selection
3258 if(IsFiducialCutOn())
3259 {
3260 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
3261 if(! in ) continue ;
3262 }
3263
3264 if(GetDebug() > 1)
3265 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());
3266
3267 //Play with the MC stack if available
3268 //Check origin of the candidates
3269 Int_t tag = 0 ;
3270 if(IsDataMC())
3271 {
3272 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
3273 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
3274 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
3275 }
3276
3277 //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3278
3279 //Check Distance to Bad channel, set bit.
3280 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3281 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
3282 if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3283 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3284 continue ;
3285 }
3286
3287 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
3288
3289 //If too low number of cells, skip it
3290 if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3291 {
3292 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3293 continue ;
3294 }
3295
3296 if(GetDebug() > 1)
3297 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
3298 calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
3299
3300 //.......................................
3301 // TOF cut, BE CAREFUL WITH THIS CUT
3302 Double_t tof = calo->GetTOF()*1e9;
3303 if(tof < fTimeCutMin || tof > fTimeCutMax)
3304 {
3305 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3306 continue ;
3307 }
3308
3309 //Check PID
3310 //PID selection or bit setting
3311 Int_t nMaxima = 0;
3312 Double_t mass = 0, angle = 0;
3313 Int_t absId1 =-1, absId2 =-1;
3314 Float_t distbad1 =-1, distbad2 =-1;
3315 Bool_t fidcut1 = 0, fidcut2 = 0;
3316 TLorentzVector l1, l2;
3317
3318 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3319 GetVertex(evtIndex),nMaxima,
3320 mass,angle,l1,l2,absId1,absId2,
3321 distbad1,distbad2,fidcut1,fidcut2) ;
3322
3323
3324 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
3325
3326 // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3327 if( (fCheckSplitDistToBad) &&
3328 (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3329 {
3330 if(GetDebug() > 1)
3331 Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
3332 calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
3333
3334 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3335 continue ;
3336 }
3337
3338 //Skip events with too few or too many NLM
3339 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3340 {
3341 //FillRejectedClusterHistograms(mom,tag,nMaxima);
3342 continue ;
3343 }
3344
3345 if(GetDebug() > 1)
3346 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
3347
3348 //Skip matched clusters with tracks
3349 if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3350 {
3351 FillRejectedClusterHistograms(mom,tag,nMaxima);
3352 continue ;
3353 }
3354
3355 Float_t l0 = calo->GetM02();
3356 Float_t e1 = l1.Energy();
3357 Float_t e2 = l2.Energy();
3358 TLorentzVector l12 = l1+l2;
3359 Float_t ptSplit = l12.Pt();
3360 Float_t eSplit = e1+e2;
3361
3362 //mass of all clusters
3363 fhMass ->Fill(mom.E() ,mass);
3364 fhMassPt ->Fill(mom.Pt(),mass);
3365 fhMassSplitPt->Fill(ptSplit ,mass);
3366 fhPtLambda0NoSplitCut->Fill(mom.Pt(),l0);
3367
3368 // Asymmetry of all clusters
3369 Float_t asy =-10;
3370
3371 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3372 fhAsymmetry->Fill(mom.E(),asy);
3373
3374 // Divide NLM in 3 cases, 1 local maxima, 2 local maxima, more than 2 local maxima
3375 Int_t indexMax = -1;
3376 if (nMaxima==1) indexMax = 0 ;
3377 else if(nMaxima==2) indexMax = 1 ;
3378 else indexMax = 2 ;
3379 fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3380
3381 Int_t mcIndex =-1;
3382 Int_t noverlaps = 0;
3383 Float_t ptprim = 0;
3384 if(IsDataMC())
3385 {
3386 mcIndex = GetMCIndex(tag);
3387
3388 Bool_t ok = kFALSE;
3389 Int_t mcLabel = calo->GetLabel();
3390
3391 TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3392
3393 Int_t mesonLabel = -1;
3394
3395 if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3396 {
3397 if(mcIndex == kmcPi0)
3398 {
3399 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3400 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3401 }
3402 else
3403 {
3404 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3405 if(grandmom.E() > 0 && ok) ptprim = grandmom.Pt();
3406 }
3407 }
3408
3409 const UInt_t nlabels = calo->GetNLabels();
3410 Int_t overpdg[nlabels];
3411 noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3412
3413 fhMCMassPt [mcIndex]->Fill(mom.Pt(),mass);
3414 fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3415 if(mcIndex==kmcPi0)
3416 {
3417 fhMCPi0PtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3418 fhMCPi0SplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3419 fhMCPi0PtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3420 fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3421
3422 }
3423 else if(mcIndex==kmcEta)
3424 {
3425 fhMCEtaPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3426 fhMCEtaSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3427 fhMCEtaPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3428 fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3429 }
3430
3431 if(noverlaps==0)
3432 {
3433 if(mcIndex==kmcPi0)
3434 {
3435 fhMCPi0PtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3436 fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3437 }
3438 else if(mcIndex==kmcEta)
3439 {
3440 fhMCEtaPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3441 fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3442 }
3443
3444 fhMassNoOverlap ->Fill(mom.E() ,mass);
3445 fhMassPtNoOverlap ->Fill(mom.Pt(),mass);
3446 fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3447
3448 fhMCMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3449 fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3450 }
3451
3452 fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3453 }
3454
3455 // If cluster does not pass pid, not pi0/eta, skip it.
3456 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3457 {
3458 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3459 FillRejectedClusterHistograms(mom,tag,nMaxima);
3460 continue ;
3461 }
3462
3463 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3464 {
3465 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3466 FillRejectedClusterHistograms(mom,tag,nMaxima);
3467 continue ;
3468 }
3469
3470 if(GetDebug() > 1)
3471 Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3472 mom.Pt(), idPartType);
3473
3474 //Mass and asymmetry of selected pairs
3475 fhSelectedAsymmetry ->Fill(mom.E() ,asy );
3476 fhSelectedMass ->Fill(mom.E() ,mass);
3477 fhSelectedMassPt ->Fill(mom.Pt(),mass);
3478 fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3479 if(fFillAllNLMHistograms) fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3480
3481 Int_t nSM = GetModuleNumber(calo);
3482 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0 && fFillAllNLMHistograms)
3483 {
3484 fhSelectedMassPtLocMaxSM [indexMax][nSM]->Fill(mom.Pt(),mass);
3485 fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),l0 );
3486 }
3487
3488 if(IsDataMC())
3489 {
3490 if(mcIndex==kmcPi0)
3491 {
3492 fhMCPi0SelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3493 fhMCPi0SelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3494 fhMCPi0SelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3495 fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3496 }
3497 else if(mcIndex==kmcEta)
3498 {
3499 fhMCEtaSelectedPtRecoPtPrim ->Fill(mom.Pt(),ptprim);
3500 fhMCEtaSelectedSplitPtRecoPtPrim ->Fill(ptSplit ,ptprim);
3501 fhMCEtaSelectedPtRecoPtPrimLocMax [indexMax]->Fill(mom.Pt(),ptprim);
3502 fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3503 }
3504
3505 if(noverlaps==0)
3506 {
3507 fhSelectedMassNoOverlap ->Fill(mom.E() ,mass);
3508 fhSelectedMassPtNoOverlap ->Fill(mom.Pt(),mass);
3509 fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3510
3511 if(mcIndex==kmcPi0)
3512 {
3513 fhMCPi0SelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3514 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3515 }
3516 else if(mcIndex==kmcEta)
3517 {
3518 fhMCEtaSelectedPtRecoPtPrimNoOverlap ->Fill(mom.Pt(),ptprim);
3519 fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3520 }
3521 }
3522 }
3523
3524 fhSplitE ->Fill( eSplit);
3525 fhSplitPt ->Fill(ptSplit);
3526 Float_t phi = mom.Phi();
3527 if(phi<0) phi+=TMath::TwoPi();
3528 fhSplitPtPhi ->Fill(ptSplit,phi);
3529 fhSplitPtEta ->Fill(ptSplit,mom.Eta());
3530 fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3531
3532 //Check split-clusters with good time window difference
3533 Double_t tof1 = cells->GetCellTime(absId1);
3534 GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3535 tof1*=1.e9;
3536
3537 Double_t tof2 = cells->GetCellTime(absId2);
3538 GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3539 tof2*=1.e9;
3540
3541 Double_t t12diff = tof1-tof2;
3542 fhEPairDiffTime->Fill(e1+e2, t12diff);
3543
3544 if(IsDataMC())
3545 {
3546 fhMCSplitE [mcIndex]->Fill( eSplit);
3547 fhMCSplitPt [mcIndex]->Fill(ptSplit);
3548 fhMCSplitPtPhi [mcIndex]->Fill(ptSplit,phi);
3549 fhMCSplitPtEta [mcIndex]->Fill(ptSplit,mom.Eta());
3550 fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3551 fhMCNLocMaxPt [mcIndex]->Fill(mom.Pt(),nMaxima);
3552
3553 fhMCSelectedMassPt [mcIndex]->Fill(mom.Pt(),mass);
3554 fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3555 fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3556
3557 if(noverlaps==0)
3558 {
3559 fhMCSelectedMassPtNoOverlap [mcIndex]->Fill(mom.Pt(),mass);
3560 fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3561 }
3562 }
3563
3564 // Remove clusters with NLM=x depeding on a minimim energy cut
3565 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3566 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3567 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
3568
3569 //Fill some histograms about shower shape
3570 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3571 {
3572 FillSelectedClusterHistograms(calo, mom.Pt(), nMaxima, tag, asy);
3573 }
3574
3575 // Fill histograms to undertand pile-up before other cuts applied
3576 // Remember to relax time cuts in the reader
3577 Double_t tofcluster = calo->GetTOF()*1e9;
3578
3579 FillPileUpHistograms(mom.Pt(),tofcluster,calo);
3580
3581 if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL")
3582 FillEMCALBCHistograms(mom.E(), mom.Eta(), mom.Phi(), tofcluster);
3583
3584 //-----------------------
3585 //Create AOD for analysis
3586
3587 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3588 aodpi0.SetLabel(calo->GetLabel());
3589
3590 //Set the indeces of the original caloclusters
3591 aodpi0.SetCaloLabel(calo->GetID(),-1);
3592 aodpi0.SetDetector(fCalorimeter);
3593
3594 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3595 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3596 else aodpi0.SetDistToBad(0) ;
3597
3598 // Check if cluster is pi0 via cluster splitting
3599 aodpi0.SetIdentifiedParticleType(idPartType);
3600
3601 // Add number of local maxima to AOD, method name in AOD to be FIXED
3602 aodpi0.SetFiducialArea(nMaxima);
3603
3604 aodpi0.SetTag(tag);
3605
3606 //Add AOD with pi0 object to aod branch
3607 AddAODParticle(aodpi0);
3608
3609 }//loop
3610
3611 if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3612
3613}
3614
3615//______________________________________________
3616void AliAnaPi0EbE::MakeAnalysisFillHistograms()
3617{
3618 //Do analysis and fill histograms
3619
3620 if(!GetOutputAODBranch())
3621 {
3622 AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3623 }
3624
3625 //Loop on stored AOD pi0
3626 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3627 if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3628
3629 Float_t cen = GetEventCentrality();
3630 Float_t ep = GetEventPlaneAngle();
3631
3632 for(Int_t iaod = 0; iaod < naod ; iaod++)
3633 {
3634 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3635 Int_t pdg = pi0->GetIdentifiedParticleType();
3636
3637 if( ( pdg != AliCaloPID::kPi0 && pdg != AliCaloPID::kEta ) ) continue;
3638
3639 //Fill pi0 histograms
3640 Float_t ener = pi0->E();
3641 Float_t pt = pi0->Pt();
3642 Float_t phi = pi0->Phi();
3643 if(phi < 0) phi+=TMath::TwoPi();
3644 Float_t eta = pi0->Eta();
3645
3646 fhPt ->Fill(pt );
3647 fhE ->Fill(ener);
3648
3649 fhPtEta ->Fill(pt ,eta);
3650 fhPtPhi ->Fill(pt ,phi);
3651 fhEtaPhi ->Fill(eta ,phi);
3652
3653 if(fFillHighMultHistograms)
3654 {
3655 fhPtCentrality ->Fill(pt,cen) ;
3656 fhPtEventPlane ->Fill(pt,ep ) ;
3657 }
3658
3659 if(IsDataMC())
3660 {
3661 Int_t tag = pi0->GetTag();
3662 Int_t label = pi0->GetLabel();
3663 Int_t mcIndex = GetMCIndex(tag);
3664
3665 fhMCE [mcIndex] ->Fill(ener);
3666 fhMCPt [mcIndex] ->Fill(pt);
3667 fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3668 fhMCPtEta[mcIndex] ->Fill(pt,eta);
3669
3670 if(fFillHighMultHistograms) fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3671
3672 if((mcIndex==kmcPi0Decay || mcIndex==kmcEtaDecay ||
3673 mcIndex==kmcPi0 || mcIndex==kmcEta ) &&
3674 fAnaType==kSSCalo)
3675 {
3676 Float_t efracMC = 0;
3677 Int_t momlabel = -1;
3678 Bool_t ok = kFALSE;
3679
3680 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3681 if(!ok) continue;
3682
3683 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3684 {
3685 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3686 if(grandmom.E() > 0 && ok)
3687 {
3688 efracMC = grandmom.E()/ener;
3689 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3690 }
3691 }
3692 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3693 {
3694 fhMCPi0DecayPt->Fill(pt);
3695 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3696 if(grandmom.E() > 0 && ok)
3697 {
3698 efracMC = mom.E()/grandmom.E();
3699 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3700 }
3701 }
3702 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3703 {
3704 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3705 if(grandmom.E() > 0 && ok)
3706 {
3707 efracMC = grandmom.E()/ener;
3708 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3709 }
3710 }
3711 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3712 {
3713 fhMCEtaDecayPt->Fill(pt);
3714 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3715 if(grandmom.E() > 0 && ok)
3716 {
3717 efracMC = mom.E()/grandmom.E();
3718 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3719 }
3720 }
3721 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3722 {
3723 fhMCOtherDecayPt->Fill(pt);
3724 }
3725
3726 }
3727
3728 if( mcIndex==kmcPi0 || mcIndex==kmcEta )
3729 {
3730 Float_t prodR = -1;
3731 Int_t momindex = -1;
3732 Int_t mompdg = -1;
3733 Int_t momstatus = -1;
3734
3735 if(GetReader()->ReadStack())
3736 {
3737 TParticle* ancestor = GetMCStack()->Particle(label);
3738 momindex = ancestor->GetFirstMother();
3739 if(momindex < 0) return;
3740 TParticle* mother = GetMCStack()->Particle(momindex);
3741 mompdg = TMath::Abs(mother->GetPdgCode());
3742 momstatus = mother->GetStatusCode();
3743 prodR = mother->R();
3744 }
3745 else
3746 {
3747 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3748 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
3749 momindex = ancestor->GetMother();
3750 if(momindex < 0) return;
3751 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3752 mompdg = TMath::Abs(mother->GetPdgCode());
3753 momstatus = mother->GetStatus();
3754 prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3755 }
3756
3757 if( mcIndex==kmcPi0 )
3758 {
3759 fhMCPi0ProdVertex->Fill(pt,prodR);
3760
3761 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
3762 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
3763 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
3764 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
3765 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
3766 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
3767 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
3768 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3769 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
3770 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
3771 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
3772 }
3773 else if (mcIndex==kmcEta )
3774 {
3775 fhMCEtaProdVertex->Fill(pt,prodR);
3776
3777 if (momstatus == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
3778 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
3779 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
3780 else if(mompdg == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
3781 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
3782 else if(mompdg == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
3783 else if(mompdg == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
3784 else if(mompdg >= 310 && mompdg <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3785 else if(mompdg == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
3786 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
3787 else fhMCEtaPtOrigin->Fill(pt,7.5);//other?
3788 }
3789 }
3790
3791 }//Histograms with MC
3792
3793 }// aod loop
3794
3795}
3796
3797//__________________________________________________________________
3798void AliAnaPi0EbE::Print(const Option_t * opt) const
3799{
3800 //Print some relevant parameters set for the analysis
3801 if(! opt)
3802 return;
3803
3804 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3805 AliAnaCaloTrackCorrBaseClass::Print("");
3806 printf("Analysis Type = %d \n", fAnaType) ;
3807 if(fAnaType == kSSCalo)
3808 {
3809 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
3810 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
3811 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3812 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
3813 }
3814 printf(" \n") ;
3815
3816}
3817
3818