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