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