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