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