]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
correcting cone exess
[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//////////////////////////////////////////////////////////////////////////////
25
26
27// --- ROOT system ---
28#include <TList.h>
29#include <TClonesArray.h>
0c1383b5 30#include <TObjString.h>
477d6cee 31
32// --- Analysis system ---
33#include "AliAnaPi0EbE.h"
34#include "AliCaloTrackReader.h"
35#include "AliIsolationCut.h"
36#include "AliNeutralMesonSelection.h"
37#include "AliCaloPID.h"
38#include "AliMCAnalysisUtils.h"
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)
48
78a28af3 49//____________________________
477d6cee 50AliAnaPi0EbE::AliAnaPi0EbE() :
09273901 51 AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo), fCalorimeter(""),
34c16486 52 fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
e671adc2 53 fNLMCutMin(-1), fNLMCutMax(10),
e671adc2 54 fTimeCutMin(-10000), fTimeCutMax(10000),
2ad19c3d 55 fFillPileUpHistograms(0),
764ab1f4 56 fFillWeightHistograms(kFALSE), fFillTMHisto(0),
57 fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1),
1db06135 58 fInputAODGammaConvName(""),
5c46c992 59 // Histograms
09273901 60 fhPt(0), fhE(0),
61 fhEEta(0), fhEPhi(0), fhEtaPhi(0),
c8710850 62 fhPtCentrality(), fhPtEventPlane(0),
17f5b4b6 63 fhPtReject(0), fhEReject(0),
40d3ce60 64 fhEEtaReject(0), fhEPhiReject(0), fhEtaPhiReject(0),
4650f5cf 65 fhMass(0), fhAsymmetry(0),
66 fhSelectedMass(0), fhSelectedAsymmetry(0),
09273901 67 fhPtDecay(0), fhEDecay(0),
5c46c992 68 // Shower shape histos
09273901 69 fhEDispersion(0), fhELambda0(0), fhELambda1(0),
70 fhELambda0NoTRD(0), fhELambda0FracMaxCellCut(0),
71 fhEFracMaxCell(0), fhEFracMaxCellNoTRD(0),
34c16486 72 fhENCells(0), fhETime(0), fhEPairDiffTime(0),
73 fhDispEtaE(0), fhDispPhiE(0),
74 fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
4650f5cf 75 fhDispEtaPhiDiffE(0), fhSphericityE(0),
34c16486 76
5c46c992 77 // MC histos
40d3ce60 78 fhMCE(), fhMCPt(),
79 fhMCPhi(), fhMCEta(),
17f5b4b6 80 fhMCEReject(), fhMCPtReject(),
81 fhMCPtCentrality(),
883411b2 82 fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
51a0ace5 83 fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
3455f821 84 fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
85 fhMCOtherDecayPt(0),
b5dbb99b 86 fhMassPairMCPi0(0), fhMassPairMCEta(0),
87 fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
78a28af3 88 // Weight studies
09273901 89 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
90 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
31ae6d59 91 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
5c46c992 92 fhTrackMatchedMCParticle(0), fhdEdx(0),
93 fhEOverP(0), fhEOverPNoTRD(0),
94 // Number of local maxima in cluster
2ad19c3d 95 fhNLocMax(0),
96 // PileUp
97 fhTimeENoCut(0), fhTimeESPD(0), fhTimeESPDMulti(0),
98 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
99 fhTimeNPileUpVertContributors(0),
100 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
477d6cee 101{
102 //default ctor
103
34c16486 104 for(Int_t i = 0; i < 6; i++)
105 {
40d3ce60 106 fhMCE [i] = 0;
3455f821 107 fhMCPt [i] = 0;
108 fhMCPhi [i] = 0;
109 fhMCEta [i] = 0;
17f5b4b6 110 fhMCPtCentrality [i] = 0;
111
34c16486 112 fhEMCLambda0 [i] = 0;
113 fhEMCLambda0NoTRD [i] = 0;
3bfcb597 114 fhEMCLambda0FracMaxCellCut[i]= 0;
34c16486 115 fhEMCFracMaxCell [i] = 0;
116 fhEMCLambda1 [i] = 0;
117 fhEMCDispersion [i] = 0;
118
bfdcf7fb 119 fhMCEDispEta [i] = 0;
120 fhMCEDispPhi [i] = 0;
121 fhMCESumEtaPhi [i] = 0;
122 fhMCEDispEtaPhiDiff[i] = 0;
123 fhMCESphericity [i] = 0;
124 fhMCEAsymmetry [i] = 0;
125
d2655d46 126 for(Int_t j = 0; j < 7; j++)
bfdcf7fb 127 {
128 fhMCLambda0DispEta [j][i] = 0;
129 fhMCLambda0DispPhi [j][i] = 0;
130 fhMCDispEtaDispPhi [j][i] = 0;
131 fhMCAsymmetryLambda0 [j][i] = 0;
132 fhMCAsymmetryDispEta [j][i] = 0;
133 fhMCAsymmetryDispPhi [j][i] = 0;
134 }
34c16486 135 }
136
d2655d46 137 for(Int_t j = 0; j < 7; j++)
bfdcf7fb 138 {
139 fhLambda0DispEta [j] = 0;
140 fhLambda0DispPhi [j] = 0;
141 fhDispEtaDispPhi [j] = 0;
142 fhAsymmetryLambda0 [j] = 0;
143 fhAsymmetryDispEta [j] = 0;
144 fhAsymmetryDispPhi [j] = 0;
5e5e056f 145
146 fhPtPi0PileUp [j] = 0;
147 }
bfdcf7fb 148
34c16486 149 for(Int_t i = 0; i < 3; i++)
150 {
151 fhELambda0LocMax [i] = 0;
152 fhELambda1LocMax [i] = 0;
153 fhEDispersionLocMax [i] = 0;
154 fhEDispEtaLocMax [i] = 0;
155 fhEDispPhiLocMax [i] = 0;
156 fhESumEtaPhiLocMax [i] = 0;
157 fhEDispEtaPhiDiffLocMax[i] = 0;
158 fhESphericityLocMax [i] = 0;
bfdcf7fb 159 fhEAsymmetryLocMax [i] = 0;
521636d2 160 }
161
78a28af3 162 //Weight studies
1a72f6c5 163 for(Int_t i =0; i < 14; i++){
78a28af3 164 fhLambda0ForW0[i] = 0;
1a72f6c5 165 //fhLambda1ForW0[i] = 0;
3c1d9afb 166 if(i<8)fhMassPairLocMax[i] = 0;
78a28af3 167 }
168
477d6cee 169 //Initialize parameters
170 InitParameters();
171
172}
477d6cee 173
40d3ce60 174//_______________________________________________________________________________
175void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time)
2ad19c3d 176{
177 // Fill some histograms to understand pile-up
178 if(!fFillPileUpHistograms) return;
179
180 //printf("E %f, time %f\n",energy,time);
181 AliVEvent * event = GetReader()->GetInputEvent();
182
183 fhTimeENoCut->Fill(energy,time);
184 if(GetReader()->IsPileUpFromSPD()) fhTimeESPD ->Fill(energy,time);
185 if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
186
de101942 187 if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
2ad19c3d 188
189 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
190 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
191
192 // N pile up vertices
193 Int_t nVerticesSPD = -1;
194 Int_t nVerticesTracks = -1;
195
196 if (esdEv)
197 {
198 nVerticesSPD = esdEv->GetNumberOfPileupVerticesSPD();
199 nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
200
201 }//ESD
202 else if (aodEv)
203 {
204 nVerticesSPD = aodEv->GetNumberOfPileupVerticesSPD();
205 nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
206 }//AOD
207
208 fhTimeNPileUpVertSPD ->Fill(time,nVerticesSPD);
209 fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
210
211 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
212 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
213
214 Int_t ncont = -1;
5559f30a 215 Float_t z1 = -1, z2 = -1;
2ad19c3d 216 Float_t diamZ = -1;
217 for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
218 {
219 if (esdEv)
220 {
221 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
222 ncont=pv->GetNContributors();
223 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
224 z2 = pv->GetZ();
225 diamZ = esdEv->GetDiamondZ();
226 }//ESD
227 else if (aodEv)
228 {
229 AliAODVertex *pv=aodEv->GetVertex(iVert);
230 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
231 ncont=pv->GetNContributors();
232 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
233 z2=pv->GetZ();
234 diamZ = aodEv->GetDiamondZ();
235 }// AOD
236
237 Double_t distZ = TMath::Abs(z2-z1);
238 diamZ = TMath::Abs(z2-diamZ);
239
240 fhTimeNPileUpVertContributors ->Fill(time,ncont);
241 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
242 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
243
244 }// loop
245}
246
40d3ce60 247
248//___________________________________________________________________________________________
249void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
250{
251 // Fill histograms that do not pass the identification (SS case only)
252
253 Float_t ener = mom.E();
254 Float_t pt = mom.Pt();
255 Float_t phi = mom.Phi();
256 if(phi < 0) phi+=TMath::TwoPi();
257 Float_t eta = mom.Eta();
258
259 fhPtReject ->Fill(pt);
260 fhEReject ->Fill(ener);
261
262 fhEEtaReject ->Fill(ener,eta);
263 fhEPhiReject ->Fill(ener,phi);
264 fhEtaPhiReject ->Fill(eta,phi);
265
266 if(IsDataMC())
267 {
268 Int_t mcIndex = GetMCIndex(mctag);
269 fhMCEReject [mcIndex] ->Fill(ener);
270 fhMCPtReject [mcIndex] ->Fill(pt);
271 }
272}
273
42d47cb7 274//_____________________________________________________________________________________
5c46c992 275void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
276 const Int_t nMaxima,
bfdcf7fb 277 const Int_t tag,
278 const Float_t asy)
5c46c992 279{
42d47cb7 280 // Fill shower shape, timing and other histograms for selected clusters from decay
281
282 Float_t e = cluster->E();
283 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
284 Float_t l0 = cluster->GetM02();
285 Float_t l1 = cluster->GetM20();
286 Int_t nSM = GetModuleNumber(cluster);
09273901 287
bfdcf7fb 288 Int_t ebin = -1;
289 if (e < 2 ) ebin = 0;
290 else if (e < 4 ) ebin = 1;
291 else if (e < 6 ) ebin = 2;
292 else if (e < 10) ebin = 3;
d2655d46 293 else if (e < 15) ebin = 4;
294 else if (e < 20) ebin = 5;
295 else ebin = 6;
296
bfdcf7fb 297 Int_t indexMax = -1;
298 if (nMaxima==1) indexMax = 0 ;
299 else if(nMaxima==2) indexMax = 1 ;
300 else indexMax = 2 ;
301
302
42d47cb7 303 AliVCaloCells * cell = 0x0;
304 if(fCalorimeter == "PHOS")
305 cell = GetPHOSCells();
306 else
307 cell = GetEMCALCells();
308
309 Float_t maxCellFraction = 0;
310 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
311 fhEFracMaxCell->Fill(e,maxCellFraction);
312
313 FillWeightHistograms(cluster);
314
315 fhEDispersion->Fill(e, disp);
316 fhELambda0 ->Fill(e, l0 );
317 fhELambda1 ->Fill(e, l1 );
318
34c16486 319 Float_t ll0 = 0., ll1 = 0.;
320 Float_t dispp= 0., dEta = 0., dPhi = 0.;
321 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
764ab1f4 322 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
34c16486 323 {
324 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
325 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
764ab1f4 326
34c16486 327 fhDispEtaE -> Fill(e,dEta);
328 fhDispPhiE -> Fill(e,dPhi);
329 fhSumEtaE -> Fill(e,sEta);
330 fhSumPhiE -> Fill(e,sPhi);
331 fhSumEtaPhiE -> Fill(e,sEtaPhi);
332 fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
26e228ff 333 if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
34c16486 334
bfdcf7fb 335 fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
336 fhLambda0DispEta[ebin]->Fill(l0 ,dEta);
337 fhLambda0DispPhi[ebin]->Fill(l0 ,dPhi);
34c16486 338
bfdcf7fb 339 if (fAnaType==kSSCalo)
340 {
341 // Asymmetry histograms
bfdcf7fb 342 fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
343 fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
344 fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
345 }
34c16486 346 }
347
5c46c992 348 fhNLocMax->Fill(e,nMaxima);
5c46c992 349
34c16486 350 fhELambda0LocMax [indexMax]->Fill(e,l0);
351 fhELambda1LocMax [indexMax]->Fill(e,l1);
352 fhEDispersionLocMax[indexMax]->Fill(e,disp);
764ab1f4 353
354 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
34c16486 355 {
356 fhEDispEtaLocMax [indexMax]-> Fill(e,dEta);
357 fhEDispPhiLocMax [indexMax]-> Fill(e,dPhi);
358 fhESumEtaPhiLocMax [indexMax]-> Fill(e,sEtaPhi);
359 fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
bfdcf7fb 360 if(dEta+dPhi>0) fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
361 if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e ,asy);
362
34c16486 363 }
364
b5dbb99b 365 if(fCalorimeter=="EMCAL" && nSM < 6)
366 {
42d47cb7 367 fhELambda0NoTRD->Fill(e, l0 );
368 fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
369 }
370
371 if(maxCellFraction < 0.5)
372 fhELambda0FracMaxCellCut->Fill(e, l0 );
373
374 fhETime ->Fill(e, cluster->GetTOF()*1.e9);
375 fhENCells->Fill(e, cluster->GetNCells());
376
09273901 377 // Fill Track matching control histograms
b5dbb99b 378 if(fFillTMHisto)
379 {
09273901 380 Float_t dZ = cluster->GetTrackDz();
381 Float_t dR = cluster->GetTrackDx();
382
b5dbb99b 383 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
384 {
09273901 385 dR = 2000., dZ = 2000.;
31ae6d59 386 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
09273901 387 }
388 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
389
b5dbb99b 390 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
391 {
09273901 392 fhTrackMatchedDEta->Fill(e,dZ);
393 fhTrackMatchedDPhi->Fill(e,dR);
b5dbb99b 394 if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
09273901 395 }
31ae6d59 396
397 // Check dEdx and E/p of matched clusters
398
399 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
400 {
4bfeae64 401 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
31ae6d59 402
34c16486 403 if(track)
404 {
31ae6d59 405 Float_t dEdx = track->GetTPCsignal();
406 fhdEdx->Fill(e, dEdx);
407
408 Float_t eOverp = e/track->P();
409 fhEOverP->Fill(e, eOverp);
4bfeae64 410
b5dbb99b 411 if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e, eOverp);
412
31ae6d59 413 }
4bfeae64 414 //else
415 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
416
417
31ae6d59 418
b5dbb99b 419 if(IsDataMC())
420 {
31ae6d59 421 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
422 {
423 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
424 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle->Fill(e, 2.5 );
425 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle->Fill(e, 0.5 );
426 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle->Fill(e, 1.5 );
427 else fhTrackMatchedMCParticle->Fill(e, 3.5 );
428
429 }
430 else
431 {
432 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
433 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle->Fill(e, 6.5 );
434 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle->Fill(e, 4.5 );
435 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle->Fill(e, 5.5 );
436 else fhTrackMatchedMCParticle->Fill(e, 7.5 );
437 }
438 } // MC
439 }
09273901 440 }// Track matching histograms
441
b5dbb99b 442 if(IsDataMC())
443 {
3455f821 444 Int_t mcIndex = GetMCIndex(tag);
34c16486 445
446 fhEMCLambda0[mcIndex] ->Fill(e, l0);
447 fhEMCLambda1[mcIndex] ->Fill(e, l1);
448 fhEMCDispersion[mcIndex] ->Fill(e, disp);
449 fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
450
451 if(fCalorimeter=="EMCAL" && nSM < 6)
452 fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0 );
764ab1f4 453
34c16486 454 if(maxCellFraction < 0.5)
455 fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0 );
456
764ab1f4 457 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
34c16486 458 {
459 fhMCEDispEta [mcIndex]-> Fill(e,dEta);
460 fhMCEDispPhi [mcIndex]-> Fill(e,dPhi);
461 fhMCESumEtaPhi [mcIndex]-> Fill(e,sEtaPhi);
462 fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
26e228ff 463 if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
464
bfdcf7fb 465 if (fAnaType==kSSCalo)
466 {
bfdcf7fb 467 fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
468 fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
469 fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
470 }
471
472 fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
473 fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0 ,dEta);
474 fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0 ,dPhi);
34c16486 475
476 }
477
42d47cb7 478 }//MC
bfdcf7fb 479
42d47cb7 480}
481
482//________________________________________________________
483void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
484{
485 // Calculate weights and fill histograms
486
487 if(!fFillWeightHistograms || GetMixedEvent()) return;
488
489 AliVCaloCells* cells = 0;
490 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
491 else cells = GetPHOSCells();
492
493 // First recalculate energy in case non linearity was applied
494 Float_t energy = 0;
495 Float_t ampMax = 0;
b5dbb99b 496 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
497 {
42d47cb7 498
499 Int_t id = clus->GetCellsAbsId()[ipos];
500
501 //Recalibrate cell energy if needed
502 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 503 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
42d47cb7 504
505 energy += amp;
506
507 if(amp> ampMax)
508 ampMax = amp;
509
510 } // energy loop
511
b5dbb99b 512 if(energy <=0 )
513 {
42d47cb7 514 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
515 return;
516 }
517
518 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
519 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
520
521 //Get the ratio and log ratio to all cells in cluster
b5dbb99b 522 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
523 {
42d47cb7 524 Int_t id = clus->GetCellsAbsId()[ipos];
525
526 //Recalibrate cell energy if needed
527 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 528 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
42d47cb7 529
530 fhECellClusterRatio ->Fill(energy,amp/energy);
531 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
532 }
533
534 //Recalculate shower shape for different W0
535 if(fCalorimeter=="EMCAL"){
536
537 Float_t l0org = clus->GetM02();
538 Float_t l1org = clus->GetM20();
539 Float_t dorg = clus->GetDispersion();
540
b5dbb99b 541 for(Int_t iw = 0; iw < 14; iw++)
542 {
1a72f6c5 543 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
42d47cb7 544 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
545
546 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 547 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
42d47cb7 548
549 } // w0 loop
550
551 // Set the original values back
552 clus->SetM02(l0org);
553 clus->SetM20(l1org);
554 clus->SetDispersion(dorg);
555
556 }// EMCAL
557}
558
b5dbb99b 559//__________________________________________
560TObjString * AliAnaPi0EbE::GetAnalysisCuts()
0c1383b5 561{
562 //Save parameters used for analysis
521636d2 563 TString parList ; //this will be list of parameters used for this analysis.
564 const Int_t buffersize = 255;
565 char onePar[buffersize] ;
566
567 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
568 parList+=onePar ;
569 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
570 parList+=onePar ;
571
b5dbb99b 572 if(fAnaType == kSSCalo)
573 {
521636d2 574 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
575 parList+=onePar ;
576 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
577 parList+=onePar ;
578 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
579 parList+=onePar ;
580 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
581 parList+=onePar ;
582 }
583
584 //Get parameters set in base class.
585 parList += GetBaseParametersList() ;
586
587 //Get parameters set in PID class.
588 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
589
590 return new TObjString(parList) ;
0c1383b5 591}
592
78a28af3 593//_____________________________________________
477d6cee 594TList * AliAnaPi0EbE::GetCreateOutputObjects()
595{
596 // Create histograms to be saved in output file and
597 // store them in outputContainer
598 TList * outputContainer = new TList() ;
599 outputContainer->SetName("Pi0EbEHistos") ;
600
745913ae 601 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
602 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
603 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
604 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
605 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
606 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
607 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
42d47cb7 608
b5dbb99b 609 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
610 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
611 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
612
09273901 613 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
614 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
615 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
616 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
617 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
618 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
619
31ae6d59 620 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
621 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
622 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
623 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
624 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
625 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
626
2ad19c3d 627 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();
628 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
629 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
630
bfdcf7fb 631 TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
632 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
633 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
d2655d46 634 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
bfdcf7fb 635
b9947879 636 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
09273901 637 fhPt->SetYTitle("N");
9fb80477 638 fhPt->SetXTitle("p_{T} (GeV/c)");
09273901 639 outputContainer->Add(fhPt) ;
640
b9947879 641 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
09273901 642 fhE->SetYTitle("N");
b9947879 643 fhE->SetXTitle("E (GeV)");
09273901 644 outputContainer->Add(fhE) ;
645
646 fhEPhi = new TH2F
b9947879 647 ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
648 fhEPhi->SetYTitle("#phi (rad)");
649 fhEPhi->SetXTitle("E (GeV)");
09273901 650 outputContainer->Add(fhEPhi) ;
651
652 fhEEta = new TH2F
b9947879 653 ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
654 fhEEta->SetYTitle("#eta");
9fb80477 655 fhEEta->SetXTitle("E (GeV)");
09273901 656 outputContainer->Add(fhEEta) ;
657
658 fhEtaPhi = new TH2F
b9947879 659 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
660 fhEtaPhi->SetYTitle("#phi (rad)");
661 fhEtaPhi->SetXTitle("#eta");
09273901 662 outputContainer->Add(fhEtaPhi) ;
663
c8710850 664 fhPtCentrality = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
665 fhPtCentrality->SetYTitle("centrality");
666 fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
667 outputContainer->Add(fhPtCentrality) ;
668
669 fhPtEventPlane = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
670 fhPtEventPlane->SetYTitle("Event plane angle (rad)");
671 fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
672 outputContainer->Add(fhPtEventPlane) ;
673
40d3ce60 674 if(fAnaType == kSSCalo)
675 {
676 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
677 fhPtReject->SetYTitle("N");
678 fhPtReject->SetXTitle("p_{T} (GeV/c)");
679 outputContainer->Add(fhPtReject) ;
680
681 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
682 fhEReject->SetYTitle("N");
683 fhEReject->SetXTitle("E (GeV)");
684 outputContainer->Add(fhEReject) ;
685
686 fhEPhiReject = new TH2F
687 ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
688 fhEPhiReject->SetYTitle("#phi (rad)");
689 fhEPhiReject->SetXTitle("E (GeV)");
690 outputContainer->Add(fhEPhiReject) ;
691
692 fhEEtaReject = new TH2F
693 ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
694 fhEEtaReject->SetYTitle("#eta");
695 fhEEtaReject->SetXTitle("E (GeV)");
696 outputContainer->Add(fhEEtaReject) ;
697
698 fhEtaPhiReject = new TH2F
699 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
700 fhEtaPhiReject->SetYTitle("#phi (rad)");
701 fhEtaPhiReject->SetXTitle("#eta");
702 outputContainer->Add(fhEtaPhiReject) ;
703 }
704
f02db2c0 705 fhMass = new TH2F
706 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
707 fhMass->SetYTitle("mass (GeV/c^{2})");
708 fhMass->SetXTitle("E (GeV)");
709 outputContainer->Add(fhMass) ;
710
711 fhSelectedMass = new TH2F
712 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
713 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
714 fhSelectedMass->SetXTitle("E (GeV)");
715 outputContainer->Add(fhSelectedMass) ;
716
34c16486 717 if(fAnaType != kSSCalo)
718 {
719 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
720 fhPtDecay->SetYTitle("N");
721 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
722 outputContainer->Add(fhPtDecay) ;
723
724 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
725 fhEDecay->SetYTitle("N");
726 fhEDecay->SetXTitle("E (GeV)");
727 outputContainer->Add(fhEDecay) ;
728 }
57b97dc6 729
c4a7d28a 730 ////////
57b97dc6 731
34c16486 732 if( fFillSelectClHisto )
b5dbb99b 733 {
c4a7d28a 734
521636d2 735 fhEDispersion = new TH2F
b9947879 736 ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 737 fhEDispersion->SetYTitle("D^{2}");
738 fhEDispersion->SetXTitle("E (GeV)");
739 outputContainer->Add(fhEDispersion) ;
740
741 fhELambda0 = new TH2F
b9947879 742 ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 743 fhELambda0->SetYTitle("#lambda_{0}^{2}");
744 fhELambda0->SetXTitle("E (GeV)");
745 outputContainer->Add(fhELambda0) ;
3bfcb597 746
42d47cb7 747 fhELambda1 = new TH2F
b9947879 748 ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
42d47cb7 749 fhELambda1->SetYTitle("#lambda_{1}^{2}");
750 fhELambda1->SetXTitle("E (GeV)");
751 outputContainer->Add(fhELambda1) ;
752
3bfcb597 753 fhELambda0FracMaxCellCut = new TH2F
b9947879 754 ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3bfcb597 755 fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
756 fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
757 outputContainer->Add(fhELambda0FracMaxCellCut) ;
758
759 fhEFracMaxCell = new TH2F
b9947879 760 ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
3bfcb597 761 fhEFracMaxCell->SetYTitle("Fraction");
762 fhEFracMaxCell->SetXTitle("E (GeV)");
763 outputContainer->Add(fhEFracMaxCell) ;
5c46c992 764
06e81356 765 if(fCalorimeter=="EMCAL")
766 {
3bfcb597 767 fhELambda0NoTRD = new TH2F
b9947879 768 ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3bfcb597 769 fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
770 fhELambda0NoTRD->SetXTitle("E (GeV)");
771 outputContainer->Add(fhELambda0NoTRD) ;
772
773 fhEFracMaxCellNoTRD = new TH2F
b9947879 774 ("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 775 fhEFracMaxCellNoTRD->SetYTitle("Fraction");
776 fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
777 outputContainer->Add(fhEFracMaxCellNoTRD) ;
34c16486 778
764ab1f4 779 if(!fFillOnlySimpleSSHisto)
34c16486 780 {
764ab1f4 781 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);
782 fhDispEtaE->SetXTitle("E (GeV)");
783 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
784 outputContainer->Add(fhDispEtaE);
785
786 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);
787 fhDispPhiE->SetXTitle("E (GeV)");
788 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
789 outputContainer->Add(fhDispPhiE);
790
791 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);
792 fhSumEtaE->SetXTitle("E (GeV)");
793 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
794 outputContainer->Add(fhSumEtaE);
795
796 fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
797 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
798 fhSumPhiE->SetXTitle("E (GeV)");
799 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
800 outputContainer->Add(fhSumPhiE);
801
802 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
803 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
804 fhSumEtaPhiE->SetXTitle("E (GeV)");
805 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
806 outputContainer->Add(fhSumEtaPhiE);
bfdcf7fb 807
764ab1f4 808 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
809 nptbins,ptmin,ptmax,200, -10,10);
810 fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
811 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
812 outputContainer->Add(fhDispEtaPhiDiffE);
bfdcf7fb 813
764ab1f4 814 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
815 nptbins,ptmin,ptmax, 200, -1,1);
816 fhSphericityE->SetXTitle("E (GeV)");
817 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
818 outputContainer->Add(fhSphericityE);
bfdcf7fb 819
764ab1f4 820 for(Int_t i = 0; i < 7; i++)
821 {
822 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]),
823 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
824 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
825 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
826 outputContainer->Add(fhDispEtaDispPhi[i]);
827
828 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]),
829 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
830 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
831 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
832 outputContainer->Add(fhLambda0DispEta[i]);
833
834 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]),
835 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
836 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
837 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
838 outputContainer->Add(fhLambda0DispPhi[i]);
839
840 }
34c16486 841 }
842 }
843
844 fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
845 nptbins,ptmin,ptmax,10,0,10);
846 fhNLocMax ->SetYTitle("N maxima");
847 fhNLocMax ->SetXTitle("E (GeV)");
848 outputContainer->Add(fhNLocMax) ;
521636d2 849
34c16486 850 for (Int_t i = 0; i < 3; i++)
851 {
852 fhELambda0LocMax[i] = new TH2F(Form("hELambda0LocMax%d",i+1),
853 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
854 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
855 fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
856 fhELambda0LocMax[i]->SetXTitle("E (GeV)");
857 outputContainer->Add(fhELambda0LocMax[i]) ;
858
859 fhELambda1LocMax[i] = new TH2F(Form("hELambda1LocMax%d",i+1),
860 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
861 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
862 fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
863 fhELambda1LocMax[i]->SetXTitle("E (GeV)");
864 outputContainer->Add(fhELambda1LocMax[i]) ;
865
866 fhEDispersionLocMax[i] = new TH2F(Form("hEDispersionLocMax%d",i+1),
867 Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
868 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
869 fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
870 fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
871 outputContainer->Add(fhEDispersionLocMax[i]) ;
872
764ab1f4 873 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
34c16486 874 {
875 fhEDispEtaLocMax[i] = new TH2F(Form("hEDispEtaLocMax%d",i+1),
876 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
877 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
878 fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
879 fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
880 outputContainer->Add(fhEDispEtaLocMax[i]) ;
881
882 fhEDispPhiLocMax[i] = new TH2F(Form("hEDispPhiLocMax%d",i+1),
883 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
884 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
885 fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
886 fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
887 outputContainer->Add(fhEDispPhiLocMax[i]) ;
888
889 fhESumEtaPhiLocMax[i] = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
890 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
891 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
892 fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
893 fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
894 outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
895
896 fhEDispEtaPhiDiffLocMax[i] = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
897 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
898 nptbins,ptmin,ptmax,200, -10,10);
899 fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
900 fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
901 outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
902
903 fhESphericityLocMax[i] = new TH2F(Form("hESphericityLocMax%d",i+1),
904 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
905 nptbins,ptmin,ptmax,200, -1,1);
906 fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
907 fhESphericityLocMax[i]->SetXTitle("E (GeV)");
908 outputContainer->Add(fhESphericityLocMax[i]) ;
909 }
910
911 }
912
42d47cb7 913 fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
914 fhENCells->SetXTitle("E (GeV)");
915 fhENCells->SetYTitle("# of cells in cluster");
916 outputContainer->Add(fhENCells);
917
918 fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
919 fhETime->SetXTitle("E (GeV)");
9fb80477 920 fhETime->SetYTitle("t (ns)");
42d47cb7 921 outputContainer->Add(fhETime);
521636d2 922
764ab1f4 923 }
e7fd282f 924
06e81356 925 if(fAnaType == kIMCalo)
926 {
42d47cb7 927 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
928 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
929 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
930 outputContainer->Add(fhEPairDiffTime);
5c46c992 931
3c1d9afb 932 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
5c46c992 933 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
934 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
3c1d9afb 935 "2 Local Maxima paired with more than 2 Local Maxima",
936 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
5c46c992 937
3c1d9afb 938 for (Int_t i = 0; i < 8 ; i++)
5c46c992 939 {
940
941 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
942
943 fhMassPairLocMax[i] = new TH2F
944 (Form("MassPairLocMax%s",combiName[i].Data()),
945 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
946 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
947 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
948 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
949 outputContainer->Add(fhMassPairLocMax[i]) ;
950 }
e7fd282f 951 }
477d6cee 952
b5dbb99b 953 if(fFillTMHisto)
954 {
09273901 955 fhTrackMatchedDEta = new TH2F
31ae6d59 956 ("hTrackMatchedDEta",
09273901 957 "d#eta of cluster-track vs cluster energy",
958 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
959 fhTrackMatchedDEta->SetYTitle("d#eta");
960 fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
961
962 fhTrackMatchedDPhi = new TH2F
31ae6d59 963 ("hTrackMatchedDPhi",
09273901 964 "d#phi of cluster-track vs cluster energy",
965 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
966 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
967 fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
968
969 fhTrackMatchedDEtaDPhi = new TH2F
31ae6d59 970 ("hTrackMatchedDEtaDPhi",
09273901 971 "d#eta vs d#phi of cluster-track vs cluster energy",
972 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
973 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
974 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
975
976 outputContainer->Add(fhTrackMatchedDEta) ;
977 outputContainer->Add(fhTrackMatchedDPhi) ;
978 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
31ae6d59 979
980 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
981 fhdEdx->SetXTitle("E (GeV)");
982 fhdEdx->SetYTitle("<dE/dx>");
983 outputContainer->Add(fhdEdx);
984
985 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
986 fhEOverP->SetXTitle("E (GeV)");
987 fhEOverP->SetYTitle("E/p");
b5dbb99b 988 outputContainer->Add(fhEOverP);
989
990 if(fCalorimeter=="EMCAL")
991 {
992 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
993 fhEOverPNoTRD->SetXTitle("E (GeV)");
994 fhEOverPNoTRD->SetYTitle("E/p");
995 outputContainer->Add(fhEOverPNoTRD);
996 }
31ae6d59 997
764ab1f4 998 if(IsDataMC() && fFillTMHisto)
31ae6d59 999 {
1000 fhTrackMatchedMCParticle = new TH2F
1001 ("hTrackMatchedMCParticle",
1002 "Origin of particle vs energy",
1003 nptbins,ptmin,ptmax,8,0,8);
1004 fhTrackMatchedMCParticle->SetXTitle("E (GeV)");
1005 //fhTrackMatchedMCParticle->SetYTitle("Particle type");
1006
1007 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
1008 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
1009 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1010 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
1011 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1012 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1013 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1014 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1015
1016 outputContainer->Add(fhTrackMatchedMCParticle);
1017 }
09273901 1018 }
1019
b5dbb99b 1020 if(fFillWeightHistograms)
1021 {
78a28af3 1022 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1023 nptbins,ptmin,ptmax, 100,0,1.);
1024 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1025 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1026 outputContainer->Add(fhECellClusterRatio);
1027
1028 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 1029 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 1030 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1031 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 1032 outputContainer->Add(fhECellClusterLogRatio);
1033
1034 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1035 nptbins,ptmin,ptmax, 100,0,1.);
1036 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1037 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1038 outputContainer->Add(fhEMaxCellClusterRatio);
1039
1040 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 1041 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 1042 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1043 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 1044 outputContainer->Add(fhEMaxCellClusterLogRatio);
1045
b5dbb99b 1046 for(Int_t iw = 0; iw < 14; iw++)
1047 {
1a72f6c5 1048 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),
78a28af3 1049 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1050 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1051 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1052 outputContainer->Add(fhLambda0ForW0[iw]);
1053
1a72f6c5 1054// 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),
1055// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1056// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1057// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1058// outputContainer->Add(fhLambda1ForW0[iw]);
78a28af3 1059
1060 }
1061 }
1062
b5dbb99b 1063 if(IsDataMC())
1064 {
3455f821 1065 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1066 {
883411b2 1067 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",
1068 nptbins,ptmin,ptmax,200,0,2);
1069 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1070 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1071 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
51a0ace5 1072
883411b2 1073 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1074 nptbins,ptmin,ptmax,200,0,2);
1075 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1076 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1077 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
51a0ace5 1078
3455f821 1079 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1080 fhMCPi0DecayPt->SetYTitle("N");
1081 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1082 outputContainer->Add(fhMCPi0DecayPt) ;
1083
883411b2 1084 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}",
3455f821 1085 nptbins,ptmin,ptmax,100,0,1);
1086 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
883411b2 1087 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
3455f821 1088 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1089
51a0ace5 1090 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
3455f821 1091 fhMCEtaDecayPt->SetYTitle("N");
1092 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1093 outputContainer->Add(fhMCEtaDecayPt) ;
1094
883411b2 1095 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
3455f821 1096 nptbins,ptmin,ptmax,100,0,1);
1097 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
883411b2 1098 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
3455f821 1099 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1100
51a0ace5 1101 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
3455f821 1102 fhMCOtherDecayPt->SetYTitle("N");
1103 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1104 outputContainer->Add(fhMCOtherDecayPt) ;
1105
1106 }
1107
477d6cee 1108 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
b5dbb99b 1109 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1110 {
477d6cee 1111
b5dbb99b 1112 fhAnglePairMCPi0 = new TH2F
1113 ("AnglePairMCPi0",
1114 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1115 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1116 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1117 outputContainer->Add(fhAnglePairMCPi0) ;
1118
af722ce4 1119 if (fAnaType!= kSSCalo)
1120 {
1121 fhAnglePairMCEta = new TH2F
1122 ("AnglePairMCEta",
1123 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1124 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1125 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1126 outputContainer->Add(fhAnglePairMCEta) ;
1127
1128 fhMassPairMCPi0 = new TH2F
1129 ("MassPairMCPi0",
1130 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1131 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1132 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1133 outputContainer->Add(fhMassPairMCPi0) ;
1134
1135 fhMassPairMCEta = new TH2F
1136 ("MassPairMCEta",
1137 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1138 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1139 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1140 outputContainer->Add(fhMassPairMCEta) ;
1141 }
e4ef72be 1142
3455f821 1143 for(Int_t i = 0; i < 6; i++)
1144 {
1145
40d3ce60 1146 fhMCE[i] = new TH1F
1147 (Form("hE_MC%s",pname[i].Data()),
1148 Form("Identified as #pi^{0} (#eta), cluster from %s",
1149 ptype[i].Data()),
1150 nptbins,ptmin,ptmax);
1151 fhMCE[i]->SetYTitle("N");
1152 fhMCE[i]->SetXTitle("E (GeV)");
1153 outputContainer->Add(fhMCE[i]) ;
1154
3455f821 1155 fhMCPt[i] = new TH1F
1156 (Form("hPt_MC%s",pname[i].Data()),
1157 Form("Identified as #pi^{0} (#eta), cluster from %s",
1158 ptype[i].Data()),
1159 nptbins,ptmin,ptmax);
1160 fhMCPt[i]->SetYTitle("N");
1161 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1162 outputContainer->Add(fhMCPt[i]) ;
1163
17f5b4b6 1164 fhMCPtCentrality[i] = new TH2F
1165 (Form("hPtCentrality_MC%s",pname[i].Data()),
1166 Form("Identified as #pi^{0} (#eta), cluster from %s",
1167 ptype[i].Data()),
1168 nptbins,ptmin,ptmax, 100,0,100);
1169 fhMCPtCentrality[i]->SetYTitle("centrality");
1170 fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1171 outputContainer->Add(fhMCPtCentrality[i]) ;
1172
40d3ce60 1173 if(fAnaType == kSSCalo)
1174 {
1175 fhMCEReject[i] = new TH1F
1176 (Form("hEReject_MC%s",pname[i].Data()),
1177 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1178 ptype[i].Data()),
1179 nptbins,ptmin,ptmax);
1180 fhMCEReject[i]->SetYTitle("N");
1181 fhMCEReject[i]->SetXTitle("E (GeV)");
1182 outputContainer->Add(fhMCEReject[i]) ;
1183
1184 fhMCPtReject[i] = new TH1F
1185 (Form("hPtReject_MC%s",pname[i].Data()),
1186 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1187 ptype[i].Data()),
1188 nptbins,ptmin,ptmax);
1189 fhMCPtReject[i]->SetYTitle("N");
1190 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1191 outputContainer->Add(fhMCPtReject[i]) ;
1192 }
1193
3455f821 1194 fhMCPhi[i] = new TH2F
1195 (Form("hPhi_MC%s",pname[i].Data()),
1196 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1197 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1198 fhMCPhi[i]->SetYTitle("#phi");
1199 fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1200 outputContainer->Add(fhMCPhi[i]) ;
1201
1202 fhMCEta[i] = new TH2F
1203 (Form("hEta_MC%s",pname[i].Data()),
1204 Form("Identified as #pi^{0} (#eta), cluster from %s",
1205 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1206 fhMCEta[i]->SetYTitle("#eta");
1207 fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1208 outputContainer->Add(fhMCEta[i]) ;
1209
1210
1211 if( fFillSelectClHisto )
1212 {
e4ef72be 1213 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1214 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1215 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1216 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1217 fhEMCLambda0[i]->SetXTitle("E (GeV)");
1218 outputContainer->Add(fhEMCLambda0[i]) ;
34c16486 1219
e4ef72be 1220 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1221 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1222 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1223 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1224 fhEMCLambda1[i]->SetXTitle("E (GeV)");
1225 outputContainer->Add(fhEMCLambda1[i]) ;
34c16486 1226
e4ef72be 1227 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1228 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1229 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1230 fhEMCDispersion[i]->SetYTitle("D^{2}");
1231 fhEMCDispersion[i]->SetXTitle("E (GeV)");
1232 outputContainer->Add(fhEMCDispersion[i]) ;
34c16486 1233
e4ef72be 1234 if(fCalorimeter=="EMCAL")
34c16486 1235 {
e4ef72be 1236 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1237 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1238 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1239 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1240 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1241 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
bfdcf7fb 1242
764ab1f4 1243 if(!fFillOnlySimpleSSHisto)
e4ef72be 1244 {
764ab1f4 1245 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1246 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1247 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1248 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1249 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1250 outputContainer->Add(fhMCEDispEta[i]);
1251
1252 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1253 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1254 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1255 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1256 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1257 outputContainer->Add(fhMCEDispPhi[i]);
1258
1259 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1260 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()),
1261 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1262 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1263 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1264 outputContainer->Add(fhMCESumEtaPhi[i]);
e4ef72be 1265
764ab1f4 1266 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1267 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1268 nptbins,ptmin,ptmax,200,-10,10);
1269 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1270 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1271 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
e4ef72be 1272
764ab1f4 1273 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1274 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()),
1275 nptbins,ptmin,ptmax, 200,-1,1);
1276 fhMCESphericity[i]->SetXTitle("E (GeV)");
1277 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1278 outputContainer->Add(fhMCESphericity[i]);
e4ef72be 1279
764ab1f4 1280 for(Int_t ie = 0; ie < 7; ie++)
1281 {
1282 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1283 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]),
1284 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1285 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1286 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1287 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1288
1289 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1290 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]),
1291 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1292 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1293 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1294 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1295
1296 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1297 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]),
1298 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1299 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1300 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1301 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1302
1303 }
1304 }
e4ef72be 1305 }
1306
1307 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1308 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1309 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1310 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1311 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1312 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1313
1314 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1315 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1316 nptbins,ptmin,ptmax,100,0,1);
1317 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1318 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1319 outputContainer->Add(fhEMCFracMaxCell[i]) ;
1320
1321 }//
1322 } // shower shape histo
34c16486 1323
521636d2 1324 } //Not MC reader
477d6cee 1325 }//Histos with MC
1326
4650f5cf 1327 if(fAnaType==kSSCalo)
1328 {
1329 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1330 nptbins,ptmin,ptmax, 200, -1,1);
1331 fhAsymmetry->SetXTitle("E (GeV)");
1332 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1333 outputContainer->Add(fhAsymmetry);
1334
1335 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1336 nptbins,ptmin,ptmax, 200, -1,1);
1337 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1338 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1339 outputContainer->Add(fhSelectedAsymmetry);
1340
1341 if(IsDataMC())
1342 {
1343 for(Int_t i = 0; i< 6; i++)
1344 {
1345 fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1346 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1347 nptbins,ptmin,ptmax, 200,-1,1);
1348 fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1349 fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1350 outputContainer->Add(fhMCEAsymmetry[i]);
1351 }
1352 }
1353 }
477d6cee 1354
764ab1f4 1355 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
bfdcf7fb 1356 {
1357
bfdcf7fb 1358
1359 for(Int_t i = 0; i< 3; i++)
1360 {
1361 fhEAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1362 Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1363 nptbins,ptmin,ptmax,200, -1,1);
1364 fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1365 fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1366 outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1367 }
1368
d2655d46 1369 for(Int_t ie = 0; ie< 7; ie++)
bfdcf7fb 1370 {
1371
1372 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1373 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1374 ssbins,ssmin,ssmax , 200,-1,1);
1375 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1376 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1377 outputContainer->Add(fhAsymmetryLambda0[ie]);
1378
1379 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1380 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1381 ssbins,ssmin,ssmax , 200,-1,1);
1382 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1383 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1384 outputContainer->Add(fhAsymmetryDispEta[ie]);
1385
1386 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1387 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1388 ssbins,ssmin,ssmax , 200,-1,1);
1389 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1390 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1391 outputContainer->Add(fhAsymmetryDispPhi[ie]);
1392 }
1393
1394
1395 if(IsDataMC())
1396 {
1397 for(Int_t i = 0; i< 6; i++)
1398 {
d2655d46 1399 for(Int_t ie = 0; ie < 7; ie++)
bfdcf7fb 1400 {
1401 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1402 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1403 ssbins,ssmin,ssmax , 200,-1,1);
1404 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1405 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1406 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1407
1408 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1409 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1410 ssbins,ssmin,ssmax , 200,-1,1);
1411 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1412 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1413 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1414
1415 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1416 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1417 ssbins,ssmin,ssmax , 200,-1,1);
1418 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1419 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1420 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1421 }
bfdcf7fb 1422 }
1423 }
bfdcf7fb 1424 }
1425
2ad19c3d 1426 if(fFillPileUpHistograms)
1427 {
5e5e056f 1428
1429 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1430
1431 for(Int_t i = 0 ; i < 7 ; i++)
1432 {
1433 fhPtPi0PileUp[i] = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
1434 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1435 fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
1436 outputContainer->Add(fhPtPi0PileUp[i]);
1437 }
1438
2ad19c3d 1439 fhTimeENoCut = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1440 fhTimeENoCut->SetXTitle("E (GeV)");
1441 fhTimeENoCut->SetYTitle("time (ns)");
1442 outputContainer->Add(fhTimeENoCut);
1443
1444 fhTimeESPD = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1445 fhTimeESPD->SetXTitle("E (GeV)");
1446 fhTimeESPD->SetYTitle("time (ns)");
1447 outputContainer->Add(fhTimeESPD);
1448
1449 fhTimeESPDMulti = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1450 fhTimeESPDMulti->SetXTitle("E (GeV)");
1451 fhTimeESPDMulti->SetYTitle("time (ns)");
1452 outputContainer->Add(fhTimeESPDMulti);
1453
1454 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1455 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1456 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1457 outputContainer->Add(fhTimeNPileUpVertSPD);
1458
1459 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1460 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1461 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1462 outputContainer->Add(fhTimeNPileUpVertTrack);
1463
1464 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1465 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1466 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1467 outputContainer->Add(fhTimeNPileUpVertContributors);
1468
1469 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);
1470 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1471 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1472 outputContainer->Add(fhTimePileUpMainVertexZDistance);
1473
1474 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1475 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1476 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1477 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1478
1479 }
1480
477d6cee 1481 //Keep neutral meson selection histograms if requiered
1482 //Setting done in AliNeutralMesonSelection
1483
e4ef72be 1484 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1485 {
477d6cee 1486 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
e4ef72be 1487
477d6cee 1488 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1489 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
e4ef72be 1490
5ae09196 1491 delete nmsHistos;
477d6cee 1492 }
1493
477d6cee 1494 return outputContainer ;
1495
1496}
1497
3455f821 1498//_____________________________________________
1499Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
1500{
1501
1502 // Assign mc index depending on MC bit set
1503
1504 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
1505 {
1506 return kmcPi0 ;
1507 }//pi0
1508 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
1509 {
1510 return kmcEta ;
1511 }//eta
1512 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
51a0ace5 1513 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
3455f821 1514 {
1515 return kmcConversion ;
1516 }//conversion photon
1517 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
1518 {
1519 return kmcPhoton ;
1520 }//photon no conversion
1521 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
1522 {
1523 return kmcElectron ;
1524 }//electron
1525 else
1526 {
1527 return kmcHadron ;
1528 }//other particles
1529
1530}
1531
1532//__________________________________________________________________
1533void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
1534 AliAODPWG4Particle * photon2,
1535 Int_t & label, Int_t & tag)
1536{
1537 // Check the labels of pare in case mother was same pi0 or eta
1538 // Set the new AOD accordingly
1539
1540 Int_t label1 = photon1->GetLabel();
1541 Int_t label2 = photon2->GetLabel();
1542
1543 if(label1 < 0 || label2 < 0 ) return ;
1544
2644ead9 1545 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
1546 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
3455f821 1547 Int_t tag1 = photon1->GetTag();
1548 Int_t tag2 = photon2->GetTag();
1549
1550 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1551 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
1552 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
1553 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
1554 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
1555 )
1556 {
1557
1558 //Check if pi0/eta mother is the same
1559 if(GetReader()->ReadStack())
1560 {
1561 if(label1>=0)
1562 {
1563 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1564 label1 = mother1->GetFirstMother();
1565 //mother1 = GetMCStack()->Particle(label1);//pi0
1566 }
1567 if(label2>=0)
1568 {
1569 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1570 label2 = mother2->GetFirstMother();
1571 //mother2 = GetMCStack()->Particle(label2);//pi0
1572 }
1573 } // STACK
1574 else if(GetReader()->ReadAODMCParticles())
1575 {//&& (input > -1)){
1576 if(label1>=0)
1577 {
2644ead9 1578 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
3455f821 1579 label1 = mother1->GetMother();
1580 //mother1 = GetMCStack()->Particle(label1);//pi0
1581 }
1582 if(label2>=0)
1583 {
2644ead9 1584 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
3455f821 1585 label2 = mother2->GetMother();
1586 //mother2 = GetMCStack()->Particle(label2);//pi0
1587 }
1588 }// AOD
1589
1590 //printf("mother1 %d, mother2 %d\n",label1,label2);
1591 if( label1 == label2 && label1>=0 )
1592 {
1593
1594 label = label1;
1595
1596 TLorentzVector mom1 = *(photon1->Momentum());
1597 TLorentzVector mom2 = *(photon2->Momentum());
1598
1599 Double_t angle = mom2.Angle(mom1.Vect());
1600 Double_t mass = (mom1+mom2).M();
1601 Double_t epair = (mom1+mom2).E();
1602
1603 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
1604 {
1605 fhMassPairMCPi0 ->Fill(epair,mass);
1606 fhAnglePairMCPi0->Fill(epair,angle);
1607 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1608 }
1609 else
1610 {
1611 fhMassPairMCEta ->Fill(epair,mass);
1612 fhAnglePairMCEta->Fill(epair,angle);
1613 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
1614 }
1615
1616 } // same label
1617 } // both from eta or pi0 decay
1618
1619}
1620
521636d2 1621//____________________________________________________________________________
1622void AliAnaPi0EbE::Init()
1623{
1624 //Init
1625 //Do some checks
1626 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1627 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1628 abort();
1629 }
1630 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1631 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1632 abort();
1633 }
1634
1635}
1636
1637//____________________________________________________________________________
1638void AliAnaPi0EbE::InitParameters()
1639{
1640 //Initialize the parameters of the analysis.
1641 AddToHistogramsName("AnaPi0EbE_");
1642
1db06135 1643 fInputAODGammaConvName = "PhotonsCTS" ;
521636d2 1644 fAnaType = kIMCalo ;
1645 fCalorimeter = "EMCAL" ;
1646 fMinDist = 2.;
1647 fMinDist2 = 4.;
1648 fMinDist3 = 5.;
1649
4d97a954 1650 fNLMECutMin[0] = 10.;
1651 fNLMECutMin[1] = 6. ;
1652 fNLMECutMin[2] = 6. ;
1653
521636d2 1654}
1655
477d6cee 1656//__________________________________________________________________
1657void AliAnaPi0EbE::MakeAnalysisFillAOD()
1658{
1659 //Do analysis and fill aods
1660
1661 switch(fAnaType)
521636d2 1662 {
477d6cee 1663 case kIMCalo:
1664 MakeInvMassInCalorimeter();
1665 break;
1666
1667 case kSSCalo:
1668 MakeShowerShapeIdentification();
1669 break;
1670
1671 case kIMCaloTracks:
1672 MakeInvMassInCalorimeterAndCTS();
1673 break;
1674
521636d2 1675 }
477d6cee 1676}
1677
42d47cb7 1678//____________________________________________
477d6cee 1679void AliAnaPi0EbE::MakeInvMassInCalorimeter()
1680{
57b97dc6 1681 //Do analysis and fill aods
1682 //Search for the photon decay in calorimeters
1683 //Read photon list from AOD, produced in class AliAnaPhoton
1684 //Check if 2 photons have the mass of the pi0.
477d6cee 1685
1686 TLorentzVector mom1;
1687 TLorentzVector mom2;
1688 TLorentzVector mom ;
b5dbb99b 1689
1690 Int_t tag = 0;
1691 Int_t label = 0;
477d6cee 1692
1693 if(!GetInputAODBranch()){
a3aebfff 1694 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
477d6cee 1695 abort();
1696 }
f8006433 1697
42d47cb7 1698 //Get shower shape information of clusters
1699 TObjArray *clusters = 0;
1700 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1701 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1702
c4a7d28a 1703 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
477d6cee 1704 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
c8fe2783 1705
c4a7d28a 1706 //Vertex cut in case of mixed events
c8fe2783 1707 Int_t evtIndex1 = 0 ;
1708 if(GetMixedEvent())
1709 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
5025c139 1710 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
477d6cee 1711 mom1 = *(photon1->Momentum());
1712
42d47cb7 1713 //Get original cluster, to recover some information
1db06135 1714 Int_t iclus = -1;
1715 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
42d47cb7 1716
1db06135 1717 if(!cluster1){
42d47cb7 1718 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1719 return;
9ab9e937 1720 }
c4a7d28a 1721
b5dbb99b 1722 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1723 {
a3aebfff 1724 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
b5dbb99b 1725
c8fe2783 1726 Int_t evtIndex2 = 0 ;
1727 if(GetMixedEvent())
1728 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
b5dbb99b 1729
c8fe2783 1730 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1731 continue ;
b5dbb99b 1732
5025c139 1733 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
b5dbb99b 1734
477d6cee 1735 mom2 = *(photon2->Momentum());
c4a7d28a 1736
1db06135 1737 //Get original cluster, to recover some information
1738 Int_t iclus2;
1739 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
42d47cb7 1740
b5dbb99b 1741 if(!cluster2)
1742 {
42d47cb7 1743 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1db06135 1744 return;
9ab9e937 1745 }
c4a7d28a 1746
42d47cb7 1747 Float_t e1 = photon1->E();
1748 Float_t e2 = photon2->E();
1749
1750 //Select clusters with good time window difference
1751 Float_t tof1 = cluster1->GetTOF()*1e9;;
1752 Float_t tof2 = cluster2->GetTOF()*1e9;;
1753 Double_t t12diff = tof1-tof2;
1754 fhEPairDiffTime->Fill(e1+e2, t12diff);
1755 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1756
b5dbb99b 1757 //Play with the MC stack if available
1758 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1759
5c46c992 1760 // Check the invariant mass for different selection on the local maxima
1761 // Name of AOD method TO BE FIXED
1762 Int_t nMaxima1 = photon1->GetFiducialArea();
1763 Int_t nMaxima2 = photon2->GetFiducialArea();
1764
1765 Double_t mass = (mom1+mom2).M();
1766 Double_t epair = (mom1+mom2).E();
1767
1768 if(nMaxima1==nMaxima2)
1769 {
1770 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1771 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1772 else fhMassPairLocMax[2]->Fill(epair,mass);
1773 }
1774 else if(nMaxima1==1 || nMaxima2==1)
1775 {
1776 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
3c1d9afb 1777 else fhMassPairLocMax[4]->Fill(epair,mass);
5c46c992 1778 }
1779 else
1780 fhMassPairLocMax[5]->Fill(epair,mass);
1781
3c1d9afb 1782 // combinations with SS axis cut and NLM cut
1783 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1784 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1785 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1786 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1787
a6e83e39 1788 //Skip events with too few or too many NLM
fb51265c 1789 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
a6e83e39 1790
1791 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
1792
f02db2c0 1793 //Mass of all pairs
1794 fhMass->Fill(epair,(mom1+mom2).M());
a6e83e39 1795
57b97dc6 1796 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 1797 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
c8fe2783 1798 {
1799 if(GetDebug()>1)
1800 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 1801
57b97dc6 1802 //Fill some histograms about shower shape
06e81356 1803 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
5c46c992 1804 {
1805 FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1806 FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
42d47cb7 1807 }
521636d2 1808
803d06a8 1809 // Tag both photons as decay
1810 photon1->SetTagged(kTRUE);
1811 photon2->SetTagged(kTRUE);
09273901 1812
1813 fhPtDecay->Fill(photon1->Pt());
1814 fhEDecay ->Fill(photon1->E() );
1815
1816 fhPtDecay->Fill(photon2->Pt());
1817 fhEDecay ->Fill(photon2->E() );
2ad19c3d 1818
57b97dc6 1819 //Create AOD for analysis
c8fe2783 1820 mom = mom1+mom2;
a6e83e39 1821
f02db2c0 1822 //Mass of selected pairs
1823 fhSelectedMass->Fill(epair,mom.M());
1824
2ad19c3d 1825 // Fill histograms to undertand pile-up before other cuts applied
1826 // Remember to relax time cuts in the reader
1827 FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);
1828
c8fe2783 1829 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1830
21a4b1c0 1831 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
c8fe2783 1832 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1833
1834 // MC
1835 pi0.SetLabel(label);
c8fe2783 1836 pi0.SetTag(tag);
b5dbb99b 1837
57b97dc6 1838 //Set the indeces of the original caloclusters
c8fe2783 1839 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
f8006433 1840 //pi0.SetInputFileIndex(input);
b5dbb99b 1841
c8fe2783 1842 AddAODParticle(pi0);
b5dbb99b 1843
c8fe2783 1844 }//pi0
57b97dc6 1845
477d6cee 1846 }//2n photon loop
1847
1848 }//1st photon loop
1849
a3aebfff 1850 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
477d6cee 1851
1852}
1853
e7fd282f 1854//__________________________________________________
477d6cee 1855void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1856{
1857 //Do analysis and fill aods
1858 //Search for the photon decay in calorimeters
1859 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1860 //Check if 2 photons have the mass of the pi0.
1861
1862 TLorentzVector mom1;
1863 TLorentzVector mom2;
1864 TLorentzVector mom ;
b5dbb99b 1865 Int_t tag = 0;
1866 Int_t label = 0;
5025c139 1867 Int_t evtIndex = 0;
1db06135 1868
1869 // Check calorimeter input
477d6cee 1870 if(!GetInputAODBranch()){
a3aebfff 1871 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
477d6cee 1872 abort();
1873 }
57b97dc6 1874
1db06135 1875 // Get the array with conversion photons
1876 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1877 if(!inputAODGammaConv) {
1878
1879 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1880
1881 if(!inputAODGammaConv) {
1882 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1883
1884 return;
1885 }
1886 }
1887
1888 //Get shower shape information of clusters
1889 TObjArray *clusters = 0;
1890 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1891 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1892
1893 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
1894 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
a6e83e39 1895 if(nCTS<=0 || nCalo <=0)
1896 {
1db06135 1897 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1898 return;
1899 }
1900
1901 if(GetDebug() > 1)
1902 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1903
1904 // Do the loop, first calo, second CTS
477d6cee 1905 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1906 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1907 mom1 = *(photon1->Momentum());
1908
1db06135 1909 //Get original cluster, to recover some information
1910 Int_t iclus = -1;
1911 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1912
1913 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1914 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
5025c139 1915 if(GetMixedEvent())
1916 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1917 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
1918
477d6cee 1919 mom2 = *(photon2->Momentum());
57b97dc6 1920
5c46c992 1921 Double_t mass = (mom1+mom2).M();
1922 Double_t epair = (mom1+mom2).E();
1923
1924 Int_t nMaxima = photon1->GetFiducialArea();
1925 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1926 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1927 else fhMassPairLocMax[2]->Fill(epair,mass);
1928
a6e83e39 1929 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
1930 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
1931
b5dbb99b 1932 //Play with the MC stack if available
1933 if(IsDataMC())
1934 {
1935 Int_t label2 = photon2->GetLabel();
2644ead9 1936 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
b5dbb99b 1937
1938 HasPairSameMCMother(photon1, photon2, label, tag) ;
1939 }
1940
f02db2c0 1941 //Mass of selected pairs
e671adc2 1942 fhMass->Fill(epair,(mom1+mom2).M());
f02db2c0 1943
477d6cee 1944 //Select good pair (good phi, pt cuts, aperture and invariant mass)
b5dbb99b 1945 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1946 {
57b97dc6 1947 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());
1948
1db06135 1949 //Fill some histograms about shower shape
06e81356 1950 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
b5dbb99b 1951 {
5c46c992 1952 FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1db06135 1953 }
803d06a8 1954
1955 // Tag both photons as decay
1956 photon1->SetTagged(kTRUE);
1957 photon2->SetTagged(kTRUE);
1db06135 1958
09273901 1959 fhPtDecay->Fill(photon1->Pt());
1960 fhEDecay ->Fill(photon1->E() );
1961
57b97dc6 1962 //Create AOD for analysis
b5dbb99b 1963
57b97dc6 1964 mom = mom1+mom2;
b5dbb99b 1965
f02db2c0 1966 //Mass of selected pairs
1967 fhSelectedMass->Fill(epair,mom.M());
1968
2ad19c3d 1969 // Fill histograms to undertand pile-up before other cuts applied
1970 // Remember to relax time cuts in the reader
41ebe136 1971 if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
2ad19c3d 1972
57b97dc6 1973 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1974
21a4b1c0 1975 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
57b97dc6 1976 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1977
1978 // MC
1979 pi0.SetLabel(label);
57b97dc6 1980 pi0.SetTag(tag);
b5dbb99b 1981
57b97dc6 1982 //Set the indeces of the original tracks or caloclusters
1983 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1984 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
f8006433 1985 //pi0.SetInputFileIndex(input);
b5dbb99b 1986
57b97dc6 1987 AddAODParticle(pi0);
b5dbb99b 1988
477d6cee 1989 }//pi0
1990 }//2n photon loop
1991
1992 }//1st photon loop
1993
a3aebfff 1994 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
477d6cee 1995
1996}
1997
1998
e7fd282f 1999//_________________________________________________
477d6cee 2000void AliAnaPi0EbE::MakeShowerShapeIdentification()
2001{
2002 //Search for pi0 in fCalorimeter with shower shape analysis
2003
34c16486 2004 TObjArray * pl = 0x0;
2005 AliVCaloCells * cells = 0x0;
5ae09196 2006 //Select the Calorimeter of the photon
b5dbb99b 2007 if (fCalorimeter == "PHOS" )
34c16486 2008 {
2009 pl = GetPHOSClusters();
2010 cells = GetPHOSCells();
2011 }
5ae09196 2012 else if (fCalorimeter == "EMCAL")
34c16486 2013 {
2014 pl = GetEMCALClusters();
2015 cells = GetEMCALCells();
2016 }
57b97dc6 2017
34c16486 2018 if(!pl)
2019 {
5ae09196 2020 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2021 return;
2022 }
233e0df8 2023
477d6cee 2024 TLorentzVector mom ;
b5dbb99b 2025 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2026 {
0ae57829 2027 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
477d6cee 2028
f8006433 2029 Int_t evtIndex = 0 ;
b5dbb99b 2030 if (GetMixedEvent())
2031 {
f8006433 2032 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2033 }
34c16486 2034
5025c139 2035 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
521636d2 2036
57b97dc6 2037 //Get Momentum vector,
a6e83e39 2038 Double_t vertex[]={0,0,0};
34c16486 2039 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2040 {
2041 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2042 }//Assume that come from vertex in straight line
2043 else
2044 {
f8006433 2045 calo->GetMomentum(mom,vertex) ;
2046 }
233e0df8 2047
57b97dc6 2048 //If too small or big pt, skip it
a426f453 2049 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
34c16486 2050
477d6cee 2051 //Check acceptance selection
b5dbb99b 2052 if(IsFiducialCutOn())
2053 {
ff45398a 2054 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 2055 if(! in ) continue ;
2056 }
2057
477d6cee 2058 if(GetDebug() > 1)
a6e83e39 2059 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());
2060
477d6cee 2061 //Check Distance to Bad channel, set bit.
c8fe2783 2062 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 2063 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2064 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2065 continue ;
a6e83e39 2066
2067 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2068
34c16486 2069 //.......................................
2070 // TOF cut, BE CAREFUL WITH THIS CUT
2071 Double_t tof = calo->GetTOF()*1e9;
2072 if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2073
40d3ce60 2074 //Play with the MC stack if available
2075 //Check origin of the candidates
2076 Int_t tag = 0 ;
2077 if(IsDataMC())
2078 {
2644ead9 2079 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
40d3ce60 2080 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2081 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2082 }
2083
a6e83e39 2084 //Skip matched clusters with tracks
40d3ce60 2085 if(IsTrackMatched(calo, GetReader()->GetInputEvent()))
2086 {
2087 FillRejectedClusterHistograms(mom,tag);
2088 continue ;
2089 }
2090
477d6cee 2091
2092 //Check PID
2093 //PID selection or bit setting
34c16486 2094 Int_t nMaxima = 0 ;
2095 Double_t mass = 0 , angle = 0;
bfdcf7fb 2096 Double_t e1 = 0 , e2 = 0;
19391b8c 2097 Int_t absId1 = -1; Int_t absId2 = -1;
2098
a6e83e39 2099 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2100 GetVertex(evtIndex),nMaxima,
19391b8c 2101 mass,angle,e1,e2,absId1,absId2) ;
34c16486 2102
a6e83e39 2103 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
40d3ce60 2104
2105
a6e83e39 2106 //Skip events with too few or too many NLM
40d3ce60 2107 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2108 {
2109 FillRejectedClusterHistograms(mom,tag);
2110 continue ;
2111 }
2112
4d97a954 2113 if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
2114 if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
2115 if(nMaxima > 2 && fNLMECutMin[2] > mom.E()) continue;
2116
bb2d339b 2117 if(GetDebug() > 1)
2118 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
a6e83e39 2119
f02db2c0 2120 //mass of all clusters
2121 fhMass->Fill(mom.E(),mass);
4650f5cf 2122
2123 // Asymmetry of all clusters
2124 Float_t asy =-10;
2125 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2126 fhAsymmetry->Fill(mom.E(),asy);
2127
4650f5cf 2128 if(IsDataMC())
2129 {
4650f5cf 2130 Int_t mcIndex = GetMCIndex(tag);
2131 fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
40d3ce60 2132 }
f02db2c0 2133
a6e83e39 2134 // If cluster does not pass pid, not pi0/eta, skip it.
bb2d339b 2135 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
2136 {
2137 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
40d3ce60 2138 FillRejectedClusterHistograms(mom,tag);
bb2d339b 2139 continue ;
2140 }
2141
2142 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
2143 {
2144 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
40d3ce60 2145 FillRejectedClusterHistograms(mom,tag);
bb2d339b 2146 continue ;
2147 }
a6e83e39 2148
bb2d339b 2149 if(GetDebug() > 1)
2150 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
a6e83e39 2151 mom.Pt(), idPartType);
34c16486 2152
667432ef 2153 //Mass and asymmetry of selected pairs
e671adc2 2154 fhSelectedAsymmetry->Fill(mom.E(),asy);
4650f5cf 2155 fhSelectedMass ->Fill(mom.E(),mass);
f02db2c0 2156
a6e83e39 2157 //-----------------------
2158 //Create AOD for analysis
477d6cee 2159
a6e83e39 2160 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2161 aodpi0.SetLabel(calo->GetLabel());
2162
2163 //Set the indeces of the original caloclusters
2164 aodpi0.SetCaloLabel(calo->GetID(),-1);
2165 aodpi0.SetDetector(fCalorimeter);
2166
2167 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2168 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
2169 else aodpi0.SetDistToBad(0) ;
2170
2171 // Check if cluster is pi0 via cluster splitting
2172 aodpi0.SetIdentifiedParticleType(idPartType);
2173
8736d400 2174 // Add number of local maxima to AOD, method name in AOD to be FIXED
2175 aodpi0.SetFiducialArea(nMaxima);
2176
4650f5cf 2177 aodpi0.SetTag(tag);
477d6cee 2178
34c16486 2179 //Fill some histograms about shower shape
2180 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2181 {
bfdcf7fb 2182 FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2ad19c3d 2183 }
2184
2185 // Fill histograms to undertand pile-up before other cuts applied
2186 // Remember to relax time cuts in the reader
2187 FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
2188
477d6cee 2189 //Add AOD with pi0 object to aod branch
2190 AddAODParticle(aodpi0);
2191
2192 }//loop
2193
a3aebfff 2194 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
477d6cee 2195
2196}
e7fd282f 2197//______________________________________________
477d6cee 2198void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 2199{
477d6cee 2200 //Do analysis and fill histograms
691bdd02 2201
b5dbb99b 2202 if(!GetOutputAODBranch())
2203 {
a3aebfff 2204 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
477d6cee 2205 abort();
2206 }
2207 //Loop on stored AOD pi0
2208 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff 2209 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee 2210
c8710850 2211 Float_t cen = GetEventCentrality();
2212 Float_t ep = GetEventPlaneAngle();
2213
b5dbb99b 2214 for(Int_t iaod = 0; iaod < naod ; iaod++)
2215 {
477d6cee 2216
2217 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
21a4b1c0 2218 Int_t pdg = pi0->GetIdentifiedParticleType();
9415d854 2219
2220 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
477d6cee 2221
2222 //Fill pi0 histograms
c4a7d28a 2223 Float_t ener = pi0->E();
2224 Float_t pt = pi0->Pt();
2225 Float_t phi = pi0->Phi();
57b97dc6 2226 if(phi < 0) phi+=TMath::TwoPi();
477d6cee 2227 Float_t eta = pi0->Eta();
2228
c8710850 2229 fhPt ->Fill(pt );
09273901 2230 fhE ->Fill(ener);
477d6cee 2231
09273901 2232 fhEEta ->Fill(ener,eta);
2233 fhEPhi ->Fill(ener,phi);
c8710850 2234 fhEtaPhi ->Fill(eta ,phi);
7a972c0c 2235
c8710850 2236 fhPtCentrality ->Fill(pt,cen) ;
2237 fhPtEventPlane ->Fill(pt,ep ) ;
2238
5e5e056f 2239 if(fFillPileUpHistograms)
2240 {
2241 if(GetReader()->IsPileUpFromSPD()) fhPtPi0PileUp[0]->Fill(pt);
2242 if(GetReader()->IsPileUpFromEMCal()) fhPtPi0PileUp[1]->Fill(pt);
2243 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPi0PileUp[2]->Fill(pt);
2244 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPi0PileUp[3]->Fill(pt);
2245 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPi0PileUp[4]->Fill(pt);
2246 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPi0PileUp[5]->Fill(pt);
2247 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
2248 }
2249
2250
b5dbb99b 2251 if(IsDataMC())
2252 {
3455f821 2253 Int_t tag = pi0->GetTag();
2254 Int_t mcIndex = GetMCIndex(tag);
2255
40d3ce60 2256 fhMCE [mcIndex] ->Fill(ener);
3455f821 2257 fhMCPt [mcIndex] ->Fill(pt);
2258 fhMCPhi[mcIndex] ->Fill(pt,phi);
2259 fhMCEta[mcIndex] ->Fill(pt,eta);
2260
17f5b4b6 2261 fhMCPtCentrality[mcIndex]->Fill(pt,cen);
2262
883411b2 2263 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
af722ce4 2264 {
3455f821 2265 Float_t efracMC = 0;
2266 Int_t label = pi0->GetLabel();
51a0ace5 2267
2268 Bool_t ok = kFALSE;
2269 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2270 if(!ok) continue;
2271
2272 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2273 {
2274 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2275 if(grandmom.E() > 0 && ok)
2276 {
883411b2 2277 efracMC = grandmom.E()/ener;
2278 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
51a0ace5 2279 }
2280 }
2281 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3455f821 2282 {
2283 fhMCPi0DecayPt->Fill(pt);
883411b2 2284 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
51a0ace5 2285 if(grandmom.E() > 0 && ok)
2286 {
2287 efracMC = mom.E()/grandmom.E();
2288 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2289 }
3455f821 2290 }
51a0ace5 2291 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2292 {
2293 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2294 if(grandmom.E() > 0 && ok)
2295 {
883411b2 2296 efracMC = grandmom.E()/ener;
2297 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
51a0ace5 2298 }
2299 }
3455f821 2300 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2301 {
2302 fhMCEtaDecayPt->Fill(pt);
51a0ace5 2303 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2304 if(grandmom.E() > 0 && ok)
2305 {
2306 efracMC = mom.E()/grandmom.E();
2307 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2308 }
3455f821 2309 }
2310 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2311 {
2312 fhMCOtherDecayPt->Fill(pt);
2313 }
af722ce4 2314
477d6cee 2315 }
3455f821 2316
477d6cee 2317 }//Histograms with MC
2318
2319 }// aod loop
2320
2321}
2322
477d6cee 2323//__________________________________________________________________
2324void AliAnaPi0EbE::Print(const Option_t * opt) const
2325{
2326 //Print some relevant parameters set for the analysis
2327 if(! opt)
2328 return;
2329
2330 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2331 AliAnaCaloTrackCorrBaseClass::Print("");
477d6cee 2332 printf("Analysis Type = %d \n", fAnaType) ;
2333 if(fAnaType == kSSCalo){
2334 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
2335 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2336 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2337 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2338 }
2339 printf(" \n") ;
2340
2341}
78a28af3 2342
78a28af3 2343