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