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