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