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