]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
Bug fix
[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),
995c6150 54 fUseSplitAsyCut(kFALSE),
e671adc2 55 fTimeCutMin(-10000), fTimeCutMax(10000),
2ad19c3d 56 fFillPileUpHistograms(0),
764ab1f4 57 fFillWeightHistograms(kFALSE), fFillTMHisto(0),
58 fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1),
1db06135 59 fInputAODGammaConvName(""),
5c46c992 60 // Histograms
09273901 61 fhPt(0), fhE(0),
62 fhEEta(0), fhEPhi(0), fhEtaPhi(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
40d3ce60 661 if(fAnaType == kSSCalo)
662 {
663 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
664 fhPtReject->SetYTitle("N");
665 fhPtReject->SetXTitle("p_{T} (GeV/c)");
666 outputContainer->Add(fhPtReject) ;
667
668 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
669 fhEReject->SetYTitle("N");
670 fhEReject->SetXTitle("E (GeV)");
671 outputContainer->Add(fhEReject) ;
672
673 fhEPhiReject = new TH2F
674 ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
675 fhEPhiReject->SetYTitle("#phi (rad)");
676 fhEPhiReject->SetXTitle("E (GeV)");
677 outputContainer->Add(fhEPhiReject) ;
678
679 fhEEtaReject = new TH2F
680 ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
681 fhEEtaReject->SetYTitle("#eta");
682 fhEEtaReject->SetXTitle("E (GeV)");
683 outputContainer->Add(fhEEtaReject) ;
684
685 fhEtaPhiReject = new TH2F
686 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
687 fhEtaPhiReject->SetYTitle("#phi (rad)");
688 fhEtaPhiReject->SetXTitle("#eta");
689 outputContainer->Add(fhEtaPhiReject) ;
690 }
691
f02db2c0 692 fhMass = new TH2F
693 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
694 fhMass->SetYTitle("mass (GeV/c^{2})");
695 fhMass->SetXTitle("E (GeV)");
696 outputContainer->Add(fhMass) ;
697
698 fhSelectedMass = new TH2F
699 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
700 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
701 fhSelectedMass->SetXTitle("E (GeV)");
702 outputContainer->Add(fhSelectedMass) ;
703
34c16486 704 if(fAnaType != kSSCalo)
705 {
706 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
707 fhPtDecay->SetYTitle("N");
708 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
709 outputContainer->Add(fhPtDecay) ;
710
711 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
712 fhEDecay->SetYTitle("N");
713 fhEDecay->SetXTitle("E (GeV)");
714 outputContainer->Add(fhEDecay) ;
715 }
57b97dc6 716
c4a7d28a 717 ////////
57b97dc6 718
34c16486 719 if( fFillSelectClHisto )
b5dbb99b 720 {
c4a7d28a 721
521636d2 722 fhEDispersion = new TH2F
b9947879 723 ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 724 fhEDispersion->SetYTitle("D^{2}");
725 fhEDispersion->SetXTitle("E (GeV)");
726 outputContainer->Add(fhEDispersion) ;
727
728 fhELambda0 = new TH2F
b9947879 729 ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 730 fhELambda0->SetYTitle("#lambda_{0}^{2}");
731 fhELambda0->SetXTitle("E (GeV)");
732 outputContainer->Add(fhELambda0) ;
3bfcb597 733
42d47cb7 734 fhELambda1 = new TH2F
b9947879 735 ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
42d47cb7 736 fhELambda1->SetYTitle("#lambda_{1}^{2}");
737 fhELambda1->SetXTitle("E (GeV)");
738 outputContainer->Add(fhELambda1) ;
739
3bfcb597 740 fhELambda0FracMaxCellCut = new TH2F
b9947879 741 ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3bfcb597 742 fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
743 fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
744 outputContainer->Add(fhELambda0FracMaxCellCut) ;
745
746 fhEFracMaxCell = new TH2F
b9947879 747 ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
3bfcb597 748 fhEFracMaxCell->SetYTitle("Fraction");
749 fhEFracMaxCell->SetXTitle("E (GeV)");
750 outputContainer->Add(fhEFracMaxCell) ;
5c46c992 751
06e81356 752 if(fCalorimeter=="EMCAL")
753 {
3bfcb597 754 fhELambda0NoTRD = new TH2F
b9947879 755 ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3bfcb597 756 fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
757 fhELambda0NoTRD->SetXTitle("E (GeV)");
758 outputContainer->Add(fhELambda0NoTRD) ;
759
760 fhEFracMaxCellNoTRD = new TH2F
b9947879 761 ("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 762 fhEFracMaxCellNoTRD->SetYTitle("Fraction");
763 fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
764 outputContainer->Add(fhEFracMaxCellNoTRD) ;
34c16486 765
764ab1f4 766 if(!fFillOnlySimpleSSHisto)
34c16486 767 {
764ab1f4 768 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);
769 fhDispEtaE->SetXTitle("E (GeV)");
770 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
771 outputContainer->Add(fhDispEtaE);
772
773 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);
774 fhDispPhiE->SetXTitle("E (GeV)");
775 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
776 outputContainer->Add(fhDispPhiE);
777
778 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);
779 fhSumEtaE->SetXTitle("E (GeV)");
780 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
781 outputContainer->Add(fhSumEtaE);
782
783 fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
784 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
785 fhSumPhiE->SetXTitle("E (GeV)");
786 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
787 outputContainer->Add(fhSumPhiE);
788
789 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
790 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
791 fhSumEtaPhiE->SetXTitle("E (GeV)");
792 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
793 outputContainer->Add(fhSumEtaPhiE);
bfdcf7fb 794
764ab1f4 795 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
796 nptbins,ptmin,ptmax,200, -10,10);
797 fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
798 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
799 outputContainer->Add(fhDispEtaPhiDiffE);
bfdcf7fb 800
764ab1f4 801 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
802 nptbins,ptmin,ptmax, 200, -1,1);
803 fhSphericityE->SetXTitle("E (GeV)");
804 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
805 outputContainer->Add(fhSphericityE);
bfdcf7fb 806
764ab1f4 807 for(Int_t i = 0; i < 7; i++)
808 {
809 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]),
810 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
811 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
812 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
813 outputContainer->Add(fhDispEtaDispPhi[i]);
814
815 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]),
816 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
817 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
818 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
819 outputContainer->Add(fhLambda0DispEta[i]);
820
821 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]),
822 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
823 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
824 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
825 outputContainer->Add(fhLambda0DispPhi[i]);
826
827 }
34c16486 828 }
829 }
830
831 fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
832 nptbins,ptmin,ptmax,10,0,10);
833 fhNLocMax ->SetYTitle("N maxima");
834 fhNLocMax ->SetXTitle("E (GeV)");
835 outputContainer->Add(fhNLocMax) ;
521636d2 836
34c16486 837 for (Int_t i = 0; i < 3; i++)
838 {
839 fhELambda0LocMax[i] = new TH2F(Form("hELambda0LocMax%d",i+1),
840 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
841 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
842 fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
843 fhELambda0LocMax[i]->SetXTitle("E (GeV)");
844 outputContainer->Add(fhELambda0LocMax[i]) ;
845
846 fhELambda1LocMax[i] = new TH2F(Form("hELambda1LocMax%d",i+1),
847 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
848 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
849 fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
850 fhELambda1LocMax[i]->SetXTitle("E (GeV)");
851 outputContainer->Add(fhELambda1LocMax[i]) ;
852
853 fhEDispersionLocMax[i] = new TH2F(Form("hEDispersionLocMax%d",i+1),
854 Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
855 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
856 fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
857 fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
858 outputContainer->Add(fhEDispersionLocMax[i]) ;
859
764ab1f4 860 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
34c16486 861 {
862 fhEDispEtaLocMax[i] = new TH2F(Form("hEDispEtaLocMax%d",i+1),
863 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
864 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
865 fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
866 fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
867 outputContainer->Add(fhEDispEtaLocMax[i]) ;
868
869 fhEDispPhiLocMax[i] = new TH2F(Form("hEDispPhiLocMax%d",i+1),
870 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
871 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
872 fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
873 fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
874 outputContainer->Add(fhEDispPhiLocMax[i]) ;
875
876 fhESumEtaPhiLocMax[i] = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
877 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
878 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
879 fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
880 fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
881 outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
882
883 fhEDispEtaPhiDiffLocMax[i] = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
884 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
885 nptbins,ptmin,ptmax,200, -10,10);
886 fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
887 fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
888 outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
889
890 fhESphericityLocMax[i] = new TH2F(Form("hESphericityLocMax%d",i+1),
891 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
892 nptbins,ptmin,ptmax,200, -1,1);
893 fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
894 fhESphericityLocMax[i]->SetXTitle("E (GeV)");
895 outputContainer->Add(fhESphericityLocMax[i]) ;
896 }
897
898 }
899
42d47cb7 900 fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
901 fhENCells->SetXTitle("E (GeV)");
902 fhENCells->SetYTitle("# of cells in cluster");
903 outputContainer->Add(fhENCells);
904
905 fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
906 fhETime->SetXTitle("E (GeV)");
9fb80477 907 fhETime->SetYTitle("t (ns)");
42d47cb7 908 outputContainer->Add(fhETime);
521636d2 909
764ab1f4 910 }
e7fd282f 911
06e81356 912 if(fAnaType == kIMCalo)
913 {
42d47cb7 914 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
915 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
916 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
917 outputContainer->Add(fhEPairDiffTime);
5c46c992 918
3c1d9afb 919 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
5c46c992 920 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
921 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
3c1d9afb 922 "2 Local Maxima paired with more than 2 Local Maxima",
923 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
5c46c992 924
3c1d9afb 925 for (Int_t i = 0; i < 8 ; i++)
5c46c992 926 {
927
928 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
929
930 fhMassPairLocMax[i] = new TH2F
931 (Form("MassPairLocMax%s",combiName[i].Data()),
932 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
933 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
934 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
935 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
936 outputContainer->Add(fhMassPairLocMax[i]) ;
937 }
e7fd282f 938 }
477d6cee 939
b5dbb99b 940 if(fFillTMHisto)
941 {
09273901 942 fhTrackMatchedDEta = new TH2F
31ae6d59 943 ("hTrackMatchedDEta",
09273901 944 "d#eta of cluster-track vs cluster energy",
945 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
946 fhTrackMatchedDEta->SetYTitle("d#eta");
947 fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
948
949 fhTrackMatchedDPhi = new TH2F
31ae6d59 950 ("hTrackMatchedDPhi",
09273901 951 "d#phi of cluster-track vs cluster energy",
952 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
953 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
954 fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
955
956 fhTrackMatchedDEtaDPhi = new TH2F
31ae6d59 957 ("hTrackMatchedDEtaDPhi",
09273901 958 "d#eta vs d#phi of cluster-track vs cluster energy",
959 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
960 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
961 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
962
963 outputContainer->Add(fhTrackMatchedDEta) ;
964 outputContainer->Add(fhTrackMatchedDPhi) ;
965 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
31ae6d59 966
967 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
968 fhdEdx->SetXTitle("E (GeV)");
969 fhdEdx->SetYTitle("<dE/dx>");
970 outputContainer->Add(fhdEdx);
971
972 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
973 fhEOverP->SetXTitle("E (GeV)");
974 fhEOverP->SetYTitle("E/p");
b5dbb99b 975 outputContainer->Add(fhEOverP);
976
977 if(fCalorimeter=="EMCAL")
978 {
979 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
980 fhEOverPNoTRD->SetXTitle("E (GeV)");
981 fhEOverPNoTRD->SetYTitle("E/p");
982 outputContainer->Add(fhEOverPNoTRD);
983 }
31ae6d59 984
764ab1f4 985 if(IsDataMC() && fFillTMHisto)
31ae6d59 986 {
987 fhTrackMatchedMCParticle = new TH2F
988 ("hTrackMatchedMCParticle",
989 "Origin of particle vs energy",
990 nptbins,ptmin,ptmax,8,0,8);
991 fhTrackMatchedMCParticle->SetXTitle("E (GeV)");
992 //fhTrackMatchedMCParticle->SetYTitle("Particle type");
993
994 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
995 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
996 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
997 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
998 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
999 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1000 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1001 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1002
1003 outputContainer->Add(fhTrackMatchedMCParticle);
1004 }
09273901 1005 }
1006
b5dbb99b 1007 if(fFillWeightHistograms)
1008 {
78a28af3 1009 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1010 nptbins,ptmin,ptmax, 100,0,1.);
1011 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1012 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1013 outputContainer->Add(fhECellClusterRatio);
1014
1015 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 1016 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 1017 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1018 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 1019 outputContainer->Add(fhECellClusterLogRatio);
1020
1021 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1022 nptbins,ptmin,ptmax, 100,0,1.);
1023 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1024 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1025 outputContainer->Add(fhEMaxCellClusterRatio);
1026
1027 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 1028 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 1029 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 1030 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 1031 outputContainer->Add(fhEMaxCellClusterLogRatio);
1032
b5dbb99b 1033 for(Int_t iw = 0; iw < 14; iw++)
1034 {
1a72f6c5 1035 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 1036 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1037 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1038 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1039 outputContainer->Add(fhLambda0ForW0[iw]);
1040
1a72f6c5 1041// 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),
1042// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1043// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1044// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1045// outputContainer->Add(fhLambda1ForW0[iw]);
78a28af3 1046
1047 }
1048 }
1049
b5dbb99b 1050 if(IsDataMC())
1051 {
3455f821 1052 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1053 {
883411b2 1054 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",
1055 nptbins,ptmin,ptmax,200,0,2);
1056 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1057 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1058 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
51a0ace5 1059
883411b2 1060 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1061 nptbins,ptmin,ptmax,200,0,2);
1062 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1063 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1064 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
51a0ace5 1065
3455f821 1066 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1067 fhMCPi0DecayPt->SetYTitle("N");
1068 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1069 outputContainer->Add(fhMCPi0DecayPt) ;
1070
883411b2 1071 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 1072 nptbins,ptmin,ptmax,100,0,1);
1073 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
883411b2 1074 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
3455f821 1075 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1076
51a0ace5 1077 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
3455f821 1078 fhMCEtaDecayPt->SetYTitle("N");
1079 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1080 outputContainer->Add(fhMCEtaDecayPt) ;
1081
883411b2 1082 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
3455f821 1083 nptbins,ptmin,ptmax,100,0,1);
1084 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
883411b2 1085 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
3455f821 1086 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1087
51a0ace5 1088 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
3455f821 1089 fhMCOtherDecayPt->SetYTitle("N");
1090 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1091 outputContainer->Add(fhMCOtherDecayPt) ;
1092
1093 }
1094
477d6cee 1095 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
b5dbb99b 1096 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1097 {
477d6cee 1098
b5dbb99b 1099 fhAnglePairMCPi0 = new TH2F
1100 ("AnglePairMCPi0",
1101 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1102 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1103 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1104 outputContainer->Add(fhAnglePairMCPi0) ;
1105
af722ce4 1106 if (fAnaType!= kSSCalo)
1107 {
1108 fhAnglePairMCEta = new TH2F
1109 ("AnglePairMCEta",
1110 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1111 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1112 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1113 outputContainer->Add(fhAnglePairMCEta) ;
1114
1115 fhMassPairMCPi0 = new TH2F
1116 ("MassPairMCPi0",
1117 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1118 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1119 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1120 outputContainer->Add(fhMassPairMCPi0) ;
1121
1122 fhMassPairMCEta = new TH2F
1123 ("MassPairMCEta",
1124 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1125 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1126 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1127 outputContainer->Add(fhMassPairMCEta) ;
1128 }
e4ef72be 1129
3455f821 1130 for(Int_t i = 0; i < 6; i++)
1131 {
1132
40d3ce60 1133 fhMCE[i] = new TH1F
1134 (Form("hE_MC%s",pname[i].Data()),
1135 Form("Identified as #pi^{0} (#eta), cluster from %s",
1136 ptype[i].Data()),
1137 nptbins,ptmin,ptmax);
1138 fhMCE[i]->SetYTitle("N");
1139 fhMCE[i]->SetXTitle("E (GeV)");
1140 outputContainer->Add(fhMCE[i]) ;
1141
3455f821 1142 fhMCPt[i] = new TH1F
1143 (Form("hPt_MC%s",pname[i].Data()),
1144 Form("Identified as #pi^{0} (#eta), cluster from %s",
1145 ptype[i].Data()),
1146 nptbins,ptmin,ptmax);
1147 fhMCPt[i]->SetYTitle("N");
1148 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1149 outputContainer->Add(fhMCPt[i]) ;
1150
40d3ce60 1151 if(fAnaType == kSSCalo)
1152 {
1153 fhMCEReject[i] = new TH1F
1154 (Form("hEReject_MC%s",pname[i].Data()),
1155 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1156 ptype[i].Data()),
1157 nptbins,ptmin,ptmax);
1158 fhMCEReject[i]->SetYTitle("N");
1159 fhMCEReject[i]->SetXTitle("E (GeV)");
1160 outputContainer->Add(fhMCEReject[i]) ;
1161
1162 fhMCPtReject[i] = new TH1F
1163 (Form("hPtReject_MC%s",pname[i].Data()),
1164 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1165 ptype[i].Data()),
1166 nptbins,ptmin,ptmax);
1167 fhMCPtReject[i]->SetYTitle("N");
1168 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1169 outputContainer->Add(fhMCPtReject[i]) ;
1170 }
1171
3455f821 1172 fhMCPhi[i] = new TH2F
1173 (Form("hPhi_MC%s",pname[i].Data()),
1174 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1175 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1176 fhMCPhi[i]->SetYTitle("#phi");
1177 fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1178 outputContainer->Add(fhMCPhi[i]) ;
1179
1180 fhMCEta[i] = new TH2F
1181 (Form("hEta_MC%s",pname[i].Data()),
1182 Form("Identified as #pi^{0} (#eta), cluster from %s",
1183 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1184 fhMCEta[i]->SetYTitle("#eta");
1185 fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1186 outputContainer->Add(fhMCEta[i]) ;
1187
1188
1189 if( fFillSelectClHisto )
1190 {
e4ef72be 1191 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1192 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1193 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1194 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1195 fhEMCLambda0[i]->SetXTitle("E (GeV)");
1196 outputContainer->Add(fhEMCLambda0[i]) ;
34c16486 1197
e4ef72be 1198 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1199 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1200 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1201 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1202 fhEMCLambda1[i]->SetXTitle("E (GeV)");
1203 outputContainer->Add(fhEMCLambda1[i]) ;
34c16486 1204
e4ef72be 1205 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1206 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1207 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1208 fhEMCDispersion[i]->SetYTitle("D^{2}");
1209 fhEMCDispersion[i]->SetXTitle("E (GeV)");
1210 outputContainer->Add(fhEMCDispersion[i]) ;
34c16486 1211
e4ef72be 1212 if(fCalorimeter=="EMCAL")
34c16486 1213 {
e4ef72be 1214 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1215 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1216 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1217 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1218 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1219 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
bfdcf7fb 1220
764ab1f4 1221 if(!fFillOnlySimpleSSHisto)
e4ef72be 1222 {
764ab1f4 1223 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1224 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1225 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1226 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1227 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1228 outputContainer->Add(fhMCEDispEta[i]);
1229
1230 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1231 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1232 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1233 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1234 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1235 outputContainer->Add(fhMCEDispPhi[i]);
1236
1237 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1238 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()),
1239 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1240 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1241 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1242 outputContainer->Add(fhMCESumEtaPhi[i]);
e4ef72be 1243
764ab1f4 1244 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1245 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1246 nptbins,ptmin,ptmax,200,-10,10);
1247 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1248 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1249 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
e4ef72be 1250
764ab1f4 1251 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1252 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()),
1253 nptbins,ptmin,ptmax, 200,-1,1);
1254 fhMCESphericity[i]->SetXTitle("E (GeV)");
1255 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1256 outputContainer->Add(fhMCESphericity[i]);
e4ef72be 1257
764ab1f4 1258 for(Int_t ie = 0; ie < 7; ie++)
1259 {
1260 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1261 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]),
1262 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1263 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1264 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1265 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1266
1267 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1268 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]),
1269 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1270 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1271 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1272 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1273
1274 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1275 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]),
1276 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1277 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1278 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1279 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1280
1281 }
1282 }
e4ef72be 1283 }
1284
1285 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1286 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1287 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1288 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1289 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1290 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1291
1292 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1293 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1294 nptbins,ptmin,ptmax,100,0,1);
1295 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1296 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1297 outputContainer->Add(fhEMCFracMaxCell[i]) ;
1298
1299 }//
1300 } // shower shape histo
34c16486 1301
521636d2 1302 } //Not MC reader
477d6cee 1303 }//Histos with MC
1304
4650f5cf 1305 if(fAnaType==kSSCalo)
1306 {
1307 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1308 nptbins,ptmin,ptmax, 200, -1,1);
1309 fhAsymmetry->SetXTitle("E (GeV)");
1310 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1311 outputContainer->Add(fhAsymmetry);
1312
1313 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1314 nptbins,ptmin,ptmax, 200, -1,1);
1315 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1316 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1317 outputContainer->Add(fhSelectedAsymmetry);
1318
1319 if(IsDataMC())
1320 {
1321 for(Int_t i = 0; i< 6; i++)
1322 {
1323 fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1324 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1325 nptbins,ptmin,ptmax, 200,-1,1);
1326 fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1327 fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1328 outputContainer->Add(fhMCEAsymmetry[i]);
1329 }
1330 }
1331 }
477d6cee 1332
764ab1f4 1333 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
bfdcf7fb 1334 {
1335
bfdcf7fb 1336
1337 for(Int_t i = 0; i< 3; i++)
1338 {
1339 fhEAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1340 Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1341 nptbins,ptmin,ptmax,200, -1,1);
1342 fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1343 fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1344 outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1345 }
1346
d2655d46 1347 for(Int_t ie = 0; ie< 7; ie++)
bfdcf7fb 1348 {
1349
1350 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1351 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1352 ssbins,ssmin,ssmax , 200,-1,1);
1353 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1354 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1355 outputContainer->Add(fhAsymmetryLambda0[ie]);
1356
1357 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1358 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1359 ssbins,ssmin,ssmax , 200,-1,1);
1360 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1361 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1362 outputContainer->Add(fhAsymmetryDispEta[ie]);
1363
1364 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1365 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1366 ssbins,ssmin,ssmax , 200,-1,1);
1367 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1368 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1369 outputContainer->Add(fhAsymmetryDispPhi[ie]);
1370 }
1371
1372
1373 if(IsDataMC())
1374 {
1375 for(Int_t i = 0; i< 6; i++)
1376 {
d2655d46 1377 for(Int_t ie = 0; ie < 7; ie++)
bfdcf7fb 1378 {
1379 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1380 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1381 ssbins,ssmin,ssmax , 200,-1,1);
1382 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1383 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1384 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1385
1386 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1387 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1388 ssbins,ssmin,ssmax , 200,-1,1);
1389 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1390 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1391 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1392
1393 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1394 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1395 ssbins,ssmin,ssmax , 200,-1,1);
1396 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1397 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1398 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1399 }
bfdcf7fb 1400 }
1401 }
bfdcf7fb 1402 }
1403
2ad19c3d 1404 if(fFillPileUpHistograms)
1405 {
5e5e056f 1406
1407 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1408
1409 for(Int_t i = 0 ; i < 7 ; i++)
1410 {
1411 fhPtPi0PileUp[i] = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
1412 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1413 fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
1414 outputContainer->Add(fhPtPi0PileUp[i]);
1415 }
1416
2ad19c3d 1417 fhTimeENoCut = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1418 fhTimeENoCut->SetXTitle("E (GeV)");
1419 fhTimeENoCut->SetYTitle("time (ns)");
1420 outputContainer->Add(fhTimeENoCut);
1421
1422 fhTimeESPD = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1423 fhTimeESPD->SetXTitle("E (GeV)");
1424 fhTimeESPD->SetYTitle("time (ns)");
1425 outputContainer->Add(fhTimeESPD);
1426
1427 fhTimeESPDMulti = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1428 fhTimeESPDMulti->SetXTitle("E (GeV)");
1429 fhTimeESPDMulti->SetYTitle("time (ns)");
1430 outputContainer->Add(fhTimeESPDMulti);
1431
1432 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1433 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1434 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1435 outputContainer->Add(fhTimeNPileUpVertSPD);
1436
1437 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1438 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1439 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1440 outputContainer->Add(fhTimeNPileUpVertTrack);
1441
1442 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1443 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1444 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1445 outputContainer->Add(fhTimeNPileUpVertContributors);
1446
1447 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);
1448 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1449 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1450 outputContainer->Add(fhTimePileUpMainVertexZDistance);
1451
1452 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1453 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1454 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1455 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1456
1457 }
1458
477d6cee 1459 //Keep neutral meson selection histograms if requiered
1460 //Setting done in AliNeutralMesonSelection
1461
e4ef72be 1462 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1463 {
477d6cee 1464 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
e4ef72be 1465
477d6cee 1466 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1467 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
e4ef72be 1468
5ae09196 1469 delete nmsHistos;
477d6cee 1470 }
1471
477d6cee 1472 return outputContainer ;
1473
1474}
1475
3455f821 1476//_____________________________________________
1477Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
1478{
1479
1480 // Assign mc index depending on MC bit set
1481
1482 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
1483 {
1484 return kmcPi0 ;
1485 }//pi0
1486 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
1487 {
1488 return kmcEta ;
1489 }//eta
1490 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
51a0ace5 1491 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
3455f821 1492 {
1493 return kmcConversion ;
1494 }//conversion photon
1495 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
1496 {
1497 return kmcPhoton ;
1498 }//photon no conversion
1499 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
1500 {
1501 return kmcElectron ;
1502 }//electron
1503 else
1504 {
1505 return kmcHadron ;
1506 }//other particles
1507
1508}
1509
1510//__________________________________________________________________
1511void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
1512 AliAODPWG4Particle * photon2,
1513 Int_t & label, Int_t & tag)
1514{
1515 // Check the labels of pare in case mother was same pi0 or eta
1516 // Set the new AOD accordingly
1517
1518 Int_t label1 = photon1->GetLabel();
1519 Int_t label2 = photon2->GetLabel();
1520
1521 if(label1 < 0 || label2 < 0 ) return ;
1522
1523 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
1524 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
1525 Int_t tag1 = photon1->GetTag();
1526 Int_t tag2 = photon2->GetTag();
1527
1528 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1529 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
1530 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
1531 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
1532 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
1533 )
1534 {
1535
1536 //Check if pi0/eta mother is the same
1537 if(GetReader()->ReadStack())
1538 {
1539 if(label1>=0)
1540 {
1541 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1542 label1 = mother1->GetFirstMother();
1543 //mother1 = GetMCStack()->Particle(label1);//pi0
1544 }
1545 if(label2>=0)
1546 {
1547 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1548 label2 = mother2->GetFirstMother();
1549 //mother2 = GetMCStack()->Particle(label2);//pi0
1550 }
1551 } // STACK
1552 else if(GetReader()->ReadAODMCParticles())
1553 {//&& (input > -1)){
1554 if(label1>=0)
1555 {
1556 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
1557 label1 = mother1->GetMother();
1558 //mother1 = GetMCStack()->Particle(label1);//pi0
1559 }
1560 if(label2>=0)
1561 {
1562 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
1563 label2 = mother2->GetMother();
1564 //mother2 = GetMCStack()->Particle(label2);//pi0
1565 }
1566 }// AOD
1567
1568 //printf("mother1 %d, mother2 %d\n",label1,label2);
1569 if( label1 == label2 && label1>=0 )
1570 {
1571
1572 label = label1;
1573
1574 TLorentzVector mom1 = *(photon1->Momentum());
1575 TLorentzVector mom2 = *(photon2->Momentum());
1576
1577 Double_t angle = mom2.Angle(mom1.Vect());
1578 Double_t mass = (mom1+mom2).M();
1579 Double_t epair = (mom1+mom2).E();
1580
1581 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
1582 {
1583 fhMassPairMCPi0 ->Fill(epair,mass);
1584 fhAnglePairMCPi0->Fill(epair,angle);
1585 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1586 }
1587 else
1588 {
1589 fhMassPairMCEta ->Fill(epair,mass);
1590 fhAnglePairMCEta->Fill(epair,angle);
1591 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
1592 }
1593
1594 } // same label
1595 } // both from eta or pi0 decay
1596
1597}
1598
521636d2 1599//____________________________________________________________________________
1600void AliAnaPi0EbE::Init()
1601{
1602 //Init
1603 //Do some checks
1604 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1605 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1606 abort();
1607 }
1608 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1609 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1610 abort();
1611 }
1612
1613}
1614
1615//____________________________________________________________________________
1616void AliAnaPi0EbE::InitParameters()
1617{
1618 //Initialize the parameters of the analysis.
1619 AddToHistogramsName("AnaPi0EbE_");
1620
1db06135 1621 fInputAODGammaConvName = "PhotonsCTS" ;
521636d2 1622 fAnaType = kIMCalo ;
1623 fCalorimeter = "EMCAL" ;
1624 fMinDist = 2.;
1625 fMinDist2 = 4.;
1626 fMinDist3 = 5.;
1627
1628}
1629
477d6cee 1630//__________________________________________________________________
1631void AliAnaPi0EbE::MakeAnalysisFillAOD()
1632{
1633 //Do analysis and fill aods
1634
1635 switch(fAnaType)
521636d2 1636 {
477d6cee 1637 case kIMCalo:
1638 MakeInvMassInCalorimeter();
1639 break;
1640
1641 case kSSCalo:
1642 MakeShowerShapeIdentification();
1643 break;
1644
1645 case kIMCaloTracks:
1646 MakeInvMassInCalorimeterAndCTS();
1647 break;
1648
521636d2 1649 }
477d6cee 1650}
1651
42d47cb7 1652//____________________________________________
477d6cee 1653void AliAnaPi0EbE::MakeInvMassInCalorimeter()
1654{
57b97dc6 1655 //Do analysis and fill aods
1656 //Search for the photon decay in calorimeters
1657 //Read photon list from AOD, produced in class AliAnaPhoton
1658 //Check if 2 photons have the mass of the pi0.
477d6cee 1659
1660 TLorentzVector mom1;
1661 TLorentzVector mom2;
1662 TLorentzVector mom ;
b5dbb99b 1663
1664 Int_t tag = 0;
1665 Int_t label = 0;
477d6cee 1666
1667 if(!GetInputAODBranch()){
a3aebfff 1668 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
477d6cee 1669 abort();
1670 }
f8006433 1671
42d47cb7 1672 //Get shower shape information of clusters
1673 TObjArray *clusters = 0;
1674 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1675 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1676
c4a7d28a 1677 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
477d6cee 1678 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
c8fe2783 1679
c4a7d28a 1680 //Vertex cut in case of mixed events
c8fe2783 1681 Int_t evtIndex1 = 0 ;
1682 if(GetMixedEvent())
1683 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
5025c139 1684 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
477d6cee 1685 mom1 = *(photon1->Momentum());
1686
42d47cb7 1687 //Get original cluster, to recover some information
1db06135 1688 Int_t iclus = -1;
1689 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
42d47cb7 1690
1db06135 1691 if(!cluster1){
42d47cb7 1692 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1693 return;
9ab9e937 1694 }
c4a7d28a 1695
b5dbb99b 1696 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1697 {
a3aebfff 1698 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
b5dbb99b 1699
c8fe2783 1700 Int_t evtIndex2 = 0 ;
1701 if(GetMixedEvent())
1702 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
b5dbb99b 1703
c8fe2783 1704 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1705 continue ;
b5dbb99b 1706
5025c139 1707 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
b5dbb99b 1708
477d6cee 1709 mom2 = *(photon2->Momentum());
c4a7d28a 1710
1db06135 1711 //Get original cluster, to recover some information
1712 Int_t iclus2;
1713 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
42d47cb7 1714
b5dbb99b 1715 if(!cluster2)
1716 {
42d47cb7 1717 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1db06135 1718 return;
9ab9e937 1719 }
c4a7d28a 1720
42d47cb7 1721 Float_t e1 = photon1->E();
1722 Float_t e2 = photon2->E();
1723
1724 //Select clusters with good time window difference
1725 Float_t tof1 = cluster1->GetTOF()*1e9;;
1726 Float_t tof2 = cluster2->GetTOF()*1e9;;
1727 Double_t t12diff = tof1-tof2;
1728 fhEPairDiffTime->Fill(e1+e2, t12diff);
1729 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1730
b5dbb99b 1731 //Play with the MC stack if available
1732 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1733
5c46c992 1734 // Check the invariant mass for different selection on the local maxima
1735 // Name of AOD method TO BE FIXED
1736 Int_t nMaxima1 = photon1->GetFiducialArea();
1737 Int_t nMaxima2 = photon2->GetFiducialArea();
1738
1739 Double_t mass = (mom1+mom2).M();
1740 Double_t epair = (mom1+mom2).E();
1741
1742 if(nMaxima1==nMaxima2)
1743 {
1744 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1745 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1746 else fhMassPairLocMax[2]->Fill(epair,mass);
1747 }
1748 else if(nMaxima1==1 || nMaxima2==1)
1749 {
1750 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
3c1d9afb 1751 else fhMassPairLocMax[4]->Fill(epair,mass);
5c46c992 1752 }
1753 else
1754 fhMassPairLocMax[5]->Fill(epair,mass);
1755
3c1d9afb 1756 // combinations with SS axis cut and NLM cut
1757 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1758 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1759 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1760 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1761
a6e83e39 1762 //Skip events with too few or too many NLM
fb51265c 1763 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
a6e83e39 1764
1765 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
1766
f02db2c0 1767 //Mass of all pairs
1768 fhMass->Fill(epair,(mom1+mom2).M());
a6e83e39 1769
57b97dc6 1770 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 1771 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
c8fe2783 1772 {
1773 if(GetDebug()>1)
1774 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 1775
57b97dc6 1776 //Fill some histograms about shower shape
06e81356 1777 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
5c46c992 1778 {
1779 FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1780 FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
42d47cb7 1781 }
521636d2 1782
803d06a8 1783 // Tag both photons as decay
1784 photon1->SetTagged(kTRUE);
1785 photon2->SetTagged(kTRUE);
09273901 1786
1787 fhPtDecay->Fill(photon1->Pt());
1788 fhEDecay ->Fill(photon1->E() );
1789
1790 fhPtDecay->Fill(photon2->Pt());
1791 fhEDecay ->Fill(photon2->E() );
2ad19c3d 1792
57b97dc6 1793 //Create AOD for analysis
c8fe2783 1794 mom = mom1+mom2;
a6e83e39 1795
f02db2c0 1796 //Mass of selected pairs
1797 fhSelectedMass->Fill(epair,mom.M());
1798
2ad19c3d 1799 // Fill histograms to undertand pile-up before other cuts applied
1800 // Remember to relax time cuts in the reader
1801 FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);
1802
c8fe2783 1803 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1804
21a4b1c0 1805 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
c8fe2783 1806 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1807
1808 // MC
1809 pi0.SetLabel(label);
c8fe2783 1810 pi0.SetTag(tag);
b5dbb99b 1811
57b97dc6 1812 //Set the indeces of the original caloclusters
c8fe2783 1813 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
f8006433 1814 //pi0.SetInputFileIndex(input);
b5dbb99b 1815
c8fe2783 1816 AddAODParticle(pi0);
b5dbb99b 1817
c8fe2783 1818 }//pi0
57b97dc6 1819
477d6cee 1820 }//2n photon loop
1821
1822 }//1st photon loop
1823
a3aebfff 1824 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
477d6cee 1825
1826}
1827
e7fd282f 1828//__________________________________________________
477d6cee 1829void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1830{
1831 //Do analysis and fill aods
1832 //Search for the photon decay in calorimeters
1833 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1834 //Check if 2 photons have the mass of the pi0.
1835
1836 TLorentzVector mom1;
1837 TLorentzVector mom2;
1838 TLorentzVector mom ;
b5dbb99b 1839 Int_t tag = 0;
1840 Int_t label = 0;
5025c139 1841 Int_t evtIndex = 0;
1db06135 1842
1843 // Check calorimeter input
477d6cee 1844 if(!GetInputAODBranch()){
a3aebfff 1845 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
477d6cee 1846 abort();
1847 }
57b97dc6 1848
1db06135 1849 // Get the array with conversion photons
1850 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1851 if(!inputAODGammaConv) {
1852
1853 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1854
1855 if(!inputAODGammaConv) {
1856 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1857
1858 return;
1859 }
1860 }
1861
1862 //Get shower shape information of clusters
1863 TObjArray *clusters = 0;
1864 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1865 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1866
1867 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
1868 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
a6e83e39 1869 if(nCTS<=0 || nCalo <=0)
1870 {
1db06135 1871 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1872 return;
1873 }
1874
1875 if(GetDebug() > 1)
1876 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1877
1878 // Do the loop, first calo, second CTS
477d6cee 1879 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1880 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1881 mom1 = *(photon1->Momentum());
1882
1db06135 1883 //Get original cluster, to recover some information
1884 Int_t iclus = -1;
1885 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1886
1887 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1888 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
5025c139 1889 if(GetMixedEvent())
1890 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1891 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
1892
477d6cee 1893 mom2 = *(photon2->Momentum());
57b97dc6 1894
5c46c992 1895 Double_t mass = (mom1+mom2).M();
1896 Double_t epair = (mom1+mom2).E();
1897
1898 Int_t nMaxima = photon1->GetFiducialArea();
1899 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1900 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1901 else fhMassPairLocMax[2]->Fill(epair,mass);
1902
a6e83e39 1903 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
1904 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
1905
b5dbb99b 1906 //Play with the MC stack if available
1907 if(IsDataMC())
1908 {
1909 Int_t label2 = photon2->GetLabel();
1910 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1911
1912 HasPairSameMCMother(photon1, photon2, label, tag) ;
1913 }
1914
f02db2c0 1915 //Mass of selected pairs
e671adc2 1916 fhMass->Fill(epair,(mom1+mom2).M());
f02db2c0 1917
477d6cee 1918 //Select good pair (good phi, pt cuts, aperture and invariant mass)
b5dbb99b 1919 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1920 {
57b97dc6 1921 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());
1922
1db06135 1923 //Fill some histograms about shower shape
06e81356 1924 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
b5dbb99b 1925 {
5c46c992 1926 FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1db06135 1927 }
803d06a8 1928
1929 // Tag both photons as decay
1930 photon1->SetTagged(kTRUE);
1931 photon2->SetTagged(kTRUE);
1db06135 1932
09273901 1933 fhPtDecay->Fill(photon1->Pt());
1934 fhEDecay ->Fill(photon1->E() );
1935
57b97dc6 1936 //Create AOD for analysis
b5dbb99b 1937
57b97dc6 1938 mom = mom1+mom2;
b5dbb99b 1939
f02db2c0 1940 //Mass of selected pairs
1941 fhSelectedMass->Fill(epair,mom.M());
1942
2ad19c3d 1943 // Fill histograms to undertand pile-up before other cuts applied
1944 // Remember to relax time cuts in the reader
41ebe136 1945 if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
2ad19c3d 1946
57b97dc6 1947 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1948
21a4b1c0 1949 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
57b97dc6 1950 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1951
1952 // MC
1953 pi0.SetLabel(label);
57b97dc6 1954 pi0.SetTag(tag);
b5dbb99b 1955
57b97dc6 1956 //Set the indeces of the original tracks or caloclusters
1957 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1958 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
f8006433 1959 //pi0.SetInputFileIndex(input);
b5dbb99b 1960
57b97dc6 1961 AddAODParticle(pi0);
b5dbb99b 1962
477d6cee 1963 }//pi0
1964 }//2n photon loop
1965
1966 }//1st photon loop
1967
a3aebfff 1968 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
477d6cee 1969
1970}
1971
1972
e7fd282f 1973//_________________________________________________
477d6cee 1974void AliAnaPi0EbE::MakeShowerShapeIdentification()
1975{
1976 //Search for pi0 in fCalorimeter with shower shape analysis
1977
34c16486 1978 TObjArray * pl = 0x0;
1979 AliVCaloCells * cells = 0x0;
5ae09196 1980 //Select the Calorimeter of the photon
b5dbb99b 1981 if (fCalorimeter == "PHOS" )
34c16486 1982 {
1983 pl = GetPHOSClusters();
1984 cells = GetPHOSCells();
1985 }
5ae09196 1986 else if (fCalorimeter == "EMCAL")
34c16486 1987 {
1988 pl = GetEMCALClusters();
1989 cells = GetEMCALCells();
1990 }
57b97dc6 1991
34c16486 1992 if(!pl)
1993 {
5ae09196 1994 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1995 return;
1996 }
233e0df8 1997
477d6cee 1998 TLorentzVector mom ;
b5dbb99b 1999 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2000 {
0ae57829 2001 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
477d6cee 2002
f8006433 2003 Int_t evtIndex = 0 ;
b5dbb99b 2004 if (GetMixedEvent())
2005 {
f8006433 2006 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2007 }
34c16486 2008
5025c139 2009 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
521636d2 2010
57b97dc6 2011 //Get Momentum vector,
a6e83e39 2012 Double_t vertex[]={0,0,0};
34c16486 2013 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2014 {
2015 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2016 }//Assume that come from vertex in straight line
2017 else
2018 {
f8006433 2019 calo->GetMomentum(mom,vertex) ;
2020 }
233e0df8 2021
57b97dc6 2022 //If too small or big pt, skip it
a426f453 2023 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
34c16486 2024
477d6cee 2025 //Check acceptance selection
b5dbb99b 2026 if(IsFiducialCutOn())
2027 {
ff45398a 2028 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 2029 if(! in ) continue ;
2030 }
2031
477d6cee 2032 if(GetDebug() > 1)
a6e83e39 2033 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());
2034
477d6cee 2035 //Check Distance to Bad channel, set bit.
c8fe2783 2036 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 2037 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2038 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2039 continue ;
a6e83e39 2040
2041 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2042
34c16486 2043 //.......................................
2044 // TOF cut, BE CAREFUL WITH THIS CUT
2045 Double_t tof = calo->GetTOF()*1e9;
2046 if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2047
40d3ce60 2048 //Play with the MC stack if available
2049 //Check origin of the candidates
2050 Int_t tag = 0 ;
2051 if(IsDataMC())
2052 {
2053 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),0);
2054 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2055 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2056 }
2057
a6e83e39 2058 //Skip matched clusters with tracks
40d3ce60 2059 if(IsTrackMatched(calo, GetReader()->GetInputEvent()))
2060 {
2061 FillRejectedClusterHistograms(mom,tag);
2062 continue ;
2063 }
2064
477d6cee 2065
2066 //Check PID
2067 //PID selection or bit setting
34c16486 2068 Int_t nMaxima = 0 ;
2069 Double_t mass = 0 , angle = 0;
bfdcf7fb 2070 Double_t e1 = 0 , e2 = 0;
a6e83e39 2071 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2072 GetVertex(evtIndex),nMaxima,
2073 mass,angle,e1,e2) ;
34c16486 2074
a6e83e39 2075 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
40d3ce60 2076
2077
a6e83e39 2078 //Skip events with too few or too many NLM
40d3ce60 2079 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2080 {
2081 FillRejectedClusterHistograms(mom,tag);
2082 continue ;
2083 }
2084
bb2d339b 2085 if(GetDebug() > 1)
2086 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
a6e83e39 2087
f02db2c0 2088 //mass of all clusters
2089 fhMass->Fill(mom.E(),mass);
4650f5cf 2090
2091 // Asymmetry of all clusters
2092 Float_t asy =-10;
2093 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2094 fhAsymmetry->Fill(mom.E(),asy);
2095
4650f5cf 2096 if(IsDataMC())
2097 {
4650f5cf 2098 Int_t mcIndex = GetMCIndex(tag);
2099 fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
40d3ce60 2100 }
f02db2c0 2101
a6e83e39 2102 // If cluster does not pass pid, not pi0/eta, skip it.
bb2d339b 2103 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
2104 {
2105 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
40d3ce60 2106 FillRejectedClusterHistograms(mom,tag);
bb2d339b 2107 continue ;
2108 }
2109
2110 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
2111 {
2112 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
40d3ce60 2113 FillRejectedClusterHistograms(mom,tag);
bb2d339b 2114 continue ;
2115 }
a6e83e39 2116
bb2d339b 2117 if(GetDebug() > 1)
2118 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
a6e83e39 2119 mom.Pt(), idPartType);
34c16486 2120
e671adc2 2121 fhSelectedAsymmetry->Fill(mom.E(),asy);
2122
167f2534 2123 if( fUseSplitAsyCut && !GetCaloPID()->IsInPi0SplitAsymmetryRange(mom.E(),asy,nMaxima) )
e671adc2 2124 {
2125 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Too large asymmetry\n");
2126 FillRejectedClusterHistograms(mom,tag);
2127 continue ;
2128 }
2129
f02db2c0 2130 //Mass of selected pairs
4650f5cf 2131 fhSelectedMass ->Fill(mom.E(),mass);
f02db2c0 2132
a6e83e39 2133 //-----------------------
2134 //Create AOD for analysis
477d6cee 2135
a6e83e39 2136 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2137 aodpi0.SetLabel(calo->GetLabel());
2138
2139 //Set the indeces of the original caloclusters
2140 aodpi0.SetCaloLabel(calo->GetID(),-1);
2141 aodpi0.SetDetector(fCalorimeter);
2142
2143 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2144 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
2145 else aodpi0.SetDistToBad(0) ;
2146
2147 // Check if cluster is pi0 via cluster splitting
2148 aodpi0.SetIdentifiedParticleType(idPartType);
2149
8736d400 2150 // Add number of local maxima to AOD, method name in AOD to be FIXED
2151 aodpi0.SetFiducialArea(nMaxima);
2152
4650f5cf 2153 aodpi0.SetTag(tag);
477d6cee 2154
34c16486 2155 //Fill some histograms about shower shape
2156 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2157 {
bfdcf7fb 2158 FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2ad19c3d 2159 }
2160
2161 // Fill histograms to undertand pile-up before other cuts applied
2162 // Remember to relax time cuts in the reader
2163 FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
2164
477d6cee 2165 //Add AOD with pi0 object to aod branch
2166 AddAODParticle(aodpi0);
2167
2168 }//loop
2169
a3aebfff 2170 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
477d6cee 2171
2172}
e7fd282f 2173//______________________________________________
477d6cee 2174void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 2175{
477d6cee 2176 //Do analysis and fill histograms
691bdd02 2177
b5dbb99b 2178 if(!GetOutputAODBranch())
2179 {
a3aebfff 2180 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
477d6cee 2181 abort();
2182 }
2183 //Loop on stored AOD pi0
2184 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff 2185 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee 2186
b5dbb99b 2187 for(Int_t iaod = 0; iaod < naod ; iaod++)
2188 {
477d6cee 2189
2190 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
21a4b1c0 2191 Int_t pdg = pi0->GetIdentifiedParticleType();
9415d854 2192
2193 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
477d6cee 2194
2195 //Fill pi0 histograms
c4a7d28a 2196 Float_t ener = pi0->E();
2197 Float_t pt = pi0->Pt();
2198 Float_t phi = pi0->Phi();
57b97dc6 2199 if(phi < 0) phi+=TMath::TwoPi();
477d6cee 2200 Float_t eta = pi0->Eta();
2201
09273901 2202 fhPt ->Fill(pt);
2203 fhE ->Fill(ener);
477d6cee 2204
09273901 2205 fhEEta ->Fill(ener,eta);
2206 fhEPhi ->Fill(ener,phi);
2207 fhEtaPhi ->Fill(eta,phi);
7a972c0c 2208
5e5e056f 2209 if(fFillPileUpHistograms)
2210 {
2211 if(GetReader()->IsPileUpFromSPD()) fhPtPi0PileUp[0]->Fill(pt);
2212 if(GetReader()->IsPileUpFromEMCal()) fhPtPi0PileUp[1]->Fill(pt);
2213 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPi0PileUp[2]->Fill(pt);
2214 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPi0PileUp[3]->Fill(pt);
2215 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPi0PileUp[4]->Fill(pt);
2216 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPi0PileUp[5]->Fill(pt);
2217 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
2218 }
2219
2220
b5dbb99b 2221 if(IsDataMC())
2222 {
3455f821 2223 Int_t tag = pi0->GetTag();
2224 Int_t mcIndex = GetMCIndex(tag);
2225
40d3ce60 2226 fhMCE [mcIndex] ->Fill(ener);
3455f821 2227 fhMCPt [mcIndex] ->Fill(pt);
2228 fhMCPhi[mcIndex] ->Fill(pt,phi);
2229 fhMCEta[mcIndex] ->Fill(pt,eta);
2230
883411b2 2231 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
af722ce4 2232 {
3455f821 2233 Float_t efracMC = 0;
2234 Int_t label = pi0->GetLabel();
51a0ace5 2235
2236 Bool_t ok = kFALSE;
2237 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2238 if(!ok) continue;
2239
2240 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2241 {
2242 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2243 if(grandmom.E() > 0 && ok)
2244 {
883411b2 2245 efracMC = grandmom.E()/ener;
2246 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
51a0ace5 2247 }
2248 }
2249 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3455f821 2250 {
2251 fhMCPi0DecayPt->Fill(pt);
883411b2 2252 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
51a0ace5 2253 if(grandmom.E() > 0 && ok)
2254 {
2255 efracMC = mom.E()/grandmom.E();
2256 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2257 }
3455f821 2258 }
51a0ace5 2259 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2260 {
2261 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2262 if(grandmom.E() > 0 && ok)
2263 {
883411b2 2264 efracMC = grandmom.E()/ener;
2265 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
51a0ace5 2266 }
2267 }
3455f821 2268 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2269 {
2270 fhMCEtaDecayPt->Fill(pt);
51a0ace5 2271 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2272 if(grandmom.E() > 0 && ok)
2273 {
2274 efracMC = mom.E()/grandmom.E();
2275 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2276 }
3455f821 2277 }
2278 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2279 {
2280 fhMCOtherDecayPt->Fill(pt);
2281 }
af722ce4 2282
477d6cee 2283 }
3455f821 2284
477d6cee 2285 }//Histograms with MC
2286
2287 }// aod loop
2288
2289}
2290
477d6cee 2291//__________________________________________________________________
2292void AliAnaPi0EbE::Print(const Option_t * opt) const
2293{
2294 //Print some relevant parameters set for the analysis
2295 if(! opt)
2296 return;
2297
2298 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2299 AliAnaCaloTrackCorrBaseClass::Print("");
477d6cee 2300 printf("Analysis Type = %d \n", fAnaType) ;
2301 if(fAnaType == kSSCalo){
2302 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
2303 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2304 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2305 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2306 }
2307 printf(" \n") ;
2308
2309}
78a28af3 2310
78a28af3 2311