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