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