]> 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
commit for the new tag
[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)");
06373cc6 722 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
34c16486 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)");
06373cc6 728 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
34c16486 729 outputContainer->Add(fhSumPhiE);
730
06373cc6 731 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
34c16486 732 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
733 fhSumEtaPhiE->SetXTitle("E (GeV)");
06373cc6 734 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
34c16486 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}");
06373cc6 760 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
bfdcf7fb 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
af722ce4 1039 if (fAnaType!= kSSCalo)
1040 {
1041 fhAnglePairMCEta = new TH2F
1042 ("AnglePairMCEta",
1043 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1044 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1045 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1046 outputContainer->Add(fhAnglePairMCEta) ;
1047
1048 fhMassPairMCPi0 = new TH2F
1049 ("MassPairMCPi0",
1050 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1051 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1052 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1053 outputContainer->Add(fhMassPairMCPi0) ;
1054
1055 fhMassPairMCEta = new TH2F
1056 ("MassPairMCEta",
1057 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1058 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1059 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1060 outputContainer->Add(fhMassPairMCEta) ;
1061 }
e4ef72be 1062
1063 if( fFillSelectClHisto )
1064 {
1065 for(Int_t i = 0; i < 6; i++)
1066 {
1067 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1068 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1069 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1070 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1071 fhEMCLambda0[i]->SetXTitle("E (GeV)");
1072 outputContainer->Add(fhEMCLambda0[i]) ;
34c16486 1073
e4ef72be 1074 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1075 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1076 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1077 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1078 fhEMCLambda1[i]->SetXTitle("E (GeV)");
1079 outputContainer->Add(fhEMCLambda1[i]) ;
34c16486 1080
e4ef72be 1081 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1082 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1083 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1084 fhEMCDispersion[i]->SetYTitle("D^{2}");
1085 fhEMCDispersion[i]->SetXTitle("E (GeV)");
1086 outputContainer->Add(fhEMCDispersion[i]) ;
34c16486 1087
e4ef72be 1088 if(fCalorimeter=="EMCAL")
34c16486 1089 {
e4ef72be 1090 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1091 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1092 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1093 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1094 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1095 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
bfdcf7fb 1096
bfdcf7fb 1097
e4ef72be 1098 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1099 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1100 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1101 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1102 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1103 outputContainer->Add(fhMCEDispEta[i]);
bfdcf7fb 1104
e4ef72be 1105 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1106 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1107 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1108 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1109 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1110 outputContainer->Add(fhMCEDispPhi[i]);
1111
1112 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
06373cc6 1113 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()),
e4ef72be 1114 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1115 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
06373cc6 1116 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
e4ef72be 1117 outputContainer->Add(fhMCESumEtaPhi[i]);
1118
1119 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1120 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1121 nptbins,ptmin,ptmax,200,-10,10);
1122 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1123 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1124 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1125
1126 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1127 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()),
1128 nptbins,ptmin,ptmax, 200,-1,1);
1129 fhMCESphericity[i]->SetXTitle("E (GeV)");
1130 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1131 outputContainer->Add(fhMCESphericity[i]);
1132
1133 for(Int_t ie = 0; ie < 7; ie++)
1134 {
1135 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1136 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]),
1137 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1138 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1139 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1140 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1141
1142 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1143 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]),
1144 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1145 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1146 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1147 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1148
1149 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1150 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]),
1151 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1152 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1153 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1154 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1155
1156 }
1157 }
1158
1159 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1160 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1161 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1162 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1163 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1164 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1165
1166 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1167 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1168 nptbins,ptmin,ptmax,100,0,1);
1169 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1170 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1171 outputContainer->Add(fhEMCFracMaxCell[i]) ;
1172
1173 }//
1174 } // shower shape histo
34c16486 1175
521636d2 1176 } //Not MC reader
477d6cee 1177 }//Histos with MC
1178
1179
e4ef72be 1180 if(fAnaType==kSSCalo && fFillSelectClHisto )
bfdcf7fb 1181 {
1182
1183 fhAsymmetryE = new TH2F ("hAsymmetryE","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1184 nptbins,ptmin,ptmax, 200, -1,1);
1185 fhAsymmetryE->SetXTitle("E (GeV)");
1186 fhAsymmetryE->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1187 outputContainer->Add(fhAsymmetryE);
1188
1189 for(Int_t i = 0; i< 3; i++)
1190 {
1191 fhEAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1192 Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1193 nptbins,ptmin,ptmax,200, -1,1);
1194 fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1195 fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1196 outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1197 }
1198
d2655d46 1199 for(Int_t ie = 0; ie< 7; ie++)
bfdcf7fb 1200 {
1201
1202 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1203 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1204 ssbins,ssmin,ssmax , 200,-1,1);
1205 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1206 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1207 outputContainer->Add(fhAsymmetryLambda0[ie]);
1208
1209 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1210 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1211 ssbins,ssmin,ssmax , 200,-1,1);
1212 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1213 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1214 outputContainer->Add(fhAsymmetryDispEta[ie]);
1215
1216 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1217 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1218 ssbins,ssmin,ssmax , 200,-1,1);
1219 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1220 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1221 outputContainer->Add(fhAsymmetryDispPhi[ie]);
1222 }
1223
1224
1225 if(IsDataMC())
1226 {
1227 for(Int_t i = 0; i< 6; i++)
1228 {
1229 fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1230 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1231 nptbins,ptmin,ptmax, 200,-1,1);
1232 fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1233 fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1234 outputContainer->Add(fhMCEAsymmetry[i]);
1235
d2655d46 1236 for(Int_t ie = 0; ie < 7; ie++)
bfdcf7fb 1237 {
1238 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1239 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1240 ssbins,ssmin,ssmax , 200,-1,1);
1241 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1242 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1243 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1244
1245 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1246 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1247 ssbins,ssmin,ssmax , 200,-1,1);
1248 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1249 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1250 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1251
1252 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1253 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1254 ssbins,ssmin,ssmax , 200,-1,1);
1255 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1256 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1257 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1258 }
bfdcf7fb 1259 }
1260 }
bfdcf7fb 1261 }
1262
477d6cee 1263 //Keep neutral meson selection histograms if requiered
1264 //Setting done in AliNeutralMesonSelection
1265
e4ef72be 1266 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1267 {
477d6cee 1268 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
e4ef72be 1269
477d6cee 1270 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1271 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
e4ef72be 1272
5ae09196 1273 delete nmsHistos;
477d6cee 1274 }
1275
477d6cee 1276 return outputContainer ;
1277
1278}
1279
521636d2 1280//____________________________________________________________________________
1281void AliAnaPi0EbE::Init()
1282{
1283 //Init
1284 //Do some checks
1285 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1286 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1287 abort();
1288 }
1289 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1290 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1291 abort();
1292 }
1293
1294}
1295
1296//____________________________________________________________________________
1297void AliAnaPi0EbE::InitParameters()
1298{
1299 //Initialize the parameters of the analysis.
1300 AddToHistogramsName("AnaPi0EbE_");
1301
1db06135 1302 fInputAODGammaConvName = "PhotonsCTS" ;
521636d2 1303 fAnaType = kIMCalo ;
1304 fCalorimeter = "EMCAL" ;
1305 fMinDist = 2.;
1306 fMinDist2 = 4.;
1307 fMinDist3 = 5.;
1308
1309}
1310
477d6cee 1311//__________________________________________________________________
1312void AliAnaPi0EbE::MakeAnalysisFillAOD()
1313{
1314 //Do analysis and fill aods
1315
1316 switch(fAnaType)
521636d2 1317 {
477d6cee 1318 case kIMCalo:
1319 MakeInvMassInCalorimeter();
1320 break;
1321
1322 case kSSCalo:
1323 MakeShowerShapeIdentification();
1324 break;
1325
1326 case kIMCaloTracks:
1327 MakeInvMassInCalorimeterAndCTS();
1328 break;
1329
521636d2 1330 }
477d6cee 1331}
1332
42d47cb7 1333//____________________________________________
477d6cee 1334void AliAnaPi0EbE::MakeInvMassInCalorimeter()
1335{
57b97dc6 1336 //Do analysis and fill aods
1337 //Search for the photon decay in calorimeters
1338 //Read photon list from AOD, produced in class AliAnaPhoton
1339 //Check if 2 photons have the mass of the pi0.
477d6cee 1340
1341 TLorentzVector mom1;
1342 TLorentzVector mom2;
1343 TLorentzVector mom ;
b5dbb99b 1344
1345 Int_t tag = 0;
1346 Int_t label = 0;
477d6cee 1347
1348 if(!GetInputAODBranch()){
a3aebfff 1349 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
477d6cee 1350 abort();
1351 }
f8006433 1352
42d47cb7 1353 //Get shower shape information of clusters
1354 TObjArray *clusters = 0;
1355 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1356 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1357
c4a7d28a 1358 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
477d6cee 1359 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
c8fe2783 1360
c4a7d28a 1361 //Vertex cut in case of mixed events
c8fe2783 1362 Int_t evtIndex1 = 0 ;
1363 if(GetMixedEvent())
1364 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
5025c139 1365 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
477d6cee 1366 mom1 = *(photon1->Momentum());
1367
42d47cb7 1368 //Get original cluster, to recover some information
1db06135 1369 Int_t iclus = -1;
1370 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
42d47cb7 1371
1db06135 1372 if(!cluster1){
42d47cb7 1373 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1374 return;
9ab9e937 1375 }
c4a7d28a 1376
b5dbb99b 1377 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1378 {
a3aebfff 1379 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
b5dbb99b 1380
c8fe2783 1381 Int_t evtIndex2 = 0 ;
1382 if(GetMixedEvent())
1383 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
b5dbb99b 1384
c8fe2783 1385 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1386 continue ;
b5dbb99b 1387
5025c139 1388 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
b5dbb99b 1389
477d6cee 1390 mom2 = *(photon2->Momentum());
c4a7d28a 1391
1db06135 1392 //Get original cluster, to recover some information
1393 Int_t iclus2;
1394 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
42d47cb7 1395
b5dbb99b 1396 if(!cluster2)
1397 {
42d47cb7 1398 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1db06135 1399 return;
9ab9e937 1400 }
c4a7d28a 1401
42d47cb7 1402 Float_t e1 = photon1->E();
1403 Float_t e2 = photon2->E();
1404
1405 //Select clusters with good time window difference
1406 Float_t tof1 = cluster1->GetTOF()*1e9;;
1407 Float_t tof2 = cluster2->GetTOF()*1e9;;
1408 Double_t t12diff = tof1-tof2;
1409 fhEPairDiffTime->Fill(e1+e2, t12diff);
1410 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1411
b5dbb99b 1412 //Play with the MC stack if available
1413 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1414
5c46c992 1415 // Check the invariant mass for different selection on the local maxima
1416 // Name of AOD method TO BE FIXED
1417 Int_t nMaxima1 = photon1->GetFiducialArea();
1418 Int_t nMaxima2 = photon2->GetFiducialArea();
1419
1420 Double_t mass = (mom1+mom2).M();
1421 Double_t epair = (mom1+mom2).E();
1422
1423 if(nMaxima1==nMaxima2)
1424 {
1425 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1426 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1427 else fhMassPairLocMax[2]->Fill(epair,mass);
1428 }
1429 else if(nMaxima1==1 || nMaxima2==1)
1430 {
1431 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
3c1d9afb 1432 else fhMassPairLocMax[4]->Fill(epair,mass);
5c46c992 1433 }
1434 else
1435 fhMassPairLocMax[5]->Fill(epair,mass);
1436
3c1d9afb 1437 // combinations with SS axis cut and NLM cut
1438 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1439 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1440 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1441 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1442
57b97dc6 1443 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 1444 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
c8fe2783 1445 {
1446 if(GetDebug()>1)
1447 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 1448
57b97dc6 1449 //Fill some histograms about shower shape
06e81356 1450 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
5c46c992 1451 {
1452 FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1453 FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
42d47cb7 1454 }
521636d2 1455
803d06a8 1456 // Tag both photons as decay
1457 photon1->SetTagged(kTRUE);
1458 photon2->SetTagged(kTRUE);
09273901 1459
1460 fhPtDecay->Fill(photon1->Pt());
1461 fhEDecay ->Fill(photon1->E() );
1462
1463 fhPtDecay->Fill(photon2->Pt());
1464 fhEDecay ->Fill(photon2->E() );
803d06a8 1465
57b97dc6 1466 //Create AOD for analysis
c8fe2783 1467 mom = mom1+mom2;
b5dbb99b 1468
c8fe2783 1469 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1470
21a4b1c0 1471 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
c8fe2783 1472 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1473
1474 // MC
1475 pi0.SetLabel(label);
c8fe2783 1476 pi0.SetTag(tag);
b5dbb99b 1477
57b97dc6 1478 //Set the indeces of the original caloclusters
c8fe2783 1479 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
f8006433 1480 //pi0.SetInputFileIndex(input);
b5dbb99b 1481
c8fe2783 1482 AddAODParticle(pi0);
b5dbb99b 1483
c8fe2783 1484 }//pi0
57b97dc6 1485
477d6cee 1486 }//2n photon loop
1487
1488 }//1st photon loop
1489
a3aebfff 1490 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
477d6cee 1491
1492}
1493
e7fd282f 1494//__________________________________________________
477d6cee 1495void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1496{
1497 //Do analysis and fill aods
1498 //Search for the photon decay in calorimeters
1499 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1500 //Check if 2 photons have the mass of the pi0.
1501
1502 TLorentzVector mom1;
1503 TLorentzVector mom2;
1504 TLorentzVector mom ;
b5dbb99b 1505 Int_t tag = 0;
1506 Int_t label = 0;
5025c139 1507 Int_t evtIndex = 0;
1db06135 1508
1509 // Check calorimeter input
477d6cee 1510 if(!GetInputAODBranch()){
a3aebfff 1511 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
477d6cee 1512 abort();
1513 }
57b97dc6 1514
1db06135 1515 // Get the array with conversion photons
1516 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1517 if(!inputAODGammaConv) {
1518
1519 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1520
1521 if(!inputAODGammaConv) {
1522 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1523
1524 return;
1525 }
1526 }
1527
1528 //Get shower shape information of clusters
1529 TObjArray *clusters = 0;
1530 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1531 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1532
1533 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
1534 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
1535 if(nCTS<=0 || nCalo <=0) {
1536 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1537 return;
1538 }
1539
1540 if(GetDebug() > 1)
1541 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1542
1543 // Do the loop, first calo, second CTS
477d6cee 1544 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1545 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1546 mom1 = *(photon1->Momentum());
1547
1db06135 1548 //Get original cluster, to recover some information
1549 Int_t iclus = -1;
1550 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1551
1552 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1553 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
5025c139 1554 if(GetMixedEvent())
1555 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1556 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
1557
477d6cee 1558 mom2 = *(photon2->Momentum());
57b97dc6 1559
5c46c992 1560 Double_t mass = (mom1+mom2).M();
1561 Double_t epair = (mom1+mom2).E();
1562
1563 Int_t nMaxima = photon1->GetFiducialArea();
1564 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1565 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1566 else fhMassPairLocMax[2]->Fill(epair,mass);
1567
b5dbb99b 1568 //Play with the MC stack if available
1569 if(IsDataMC())
1570 {
1571 Int_t label2 = photon2->GetLabel();
1572 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1573
1574 HasPairSameMCMother(photon1, photon2, label, tag) ;
1575 }
1576
477d6cee 1577 //Select good pair (good phi, pt cuts, aperture and invariant mass)
b5dbb99b 1578 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1579 {
57b97dc6 1580 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());
1581
1db06135 1582 //Fill some histograms about shower shape
06e81356 1583 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
b5dbb99b 1584 {
5c46c992 1585 FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1db06135 1586 }
803d06a8 1587
1588 // Tag both photons as decay
1589 photon1->SetTagged(kTRUE);
1590 photon2->SetTagged(kTRUE);
1db06135 1591
09273901 1592 fhPtDecay->Fill(photon1->Pt());
1593 fhEDecay ->Fill(photon1->E() );
1594
57b97dc6 1595 //Create AOD for analysis
b5dbb99b 1596
57b97dc6 1597 mom = mom1+mom2;
b5dbb99b 1598
57b97dc6 1599 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
b5dbb99b 1600
21a4b1c0 1601 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
57b97dc6 1602 pi0.SetDetector(photon1->GetDetector());
b5dbb99b 1603
1604 // MC
1605 pi0.SetLabel(label);
57b97dc6 1606 pi0.SetTag(tag);
b5dbb99b 1607
57b97dc6 1608 //Set the indeces of the original tracks or caloclusters
1609 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1610 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
f8006433 1611 //pi0.SetInputFileIndex(input);
b5dbb99b 1612
57b97dc6 1613 AddAODParticle(pi0);
b5dbb99b 1614
477d6cee 1615 }//pi0
1616 }//2n photon loop
1617
1618 }//1st photon loop
1619
a3aebfff 1620 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
477d6cee 1621
1622}
1623
1624
e7fd282f 1625//_________________________________________________
477d6cee 1626void AliAnaPi0EbE::MakeShowerShapeIdentification()
1627{
1628 //Search for pi0 in fCalorimeter with shower shape analysis
1629
34c16486 1630 TObjArray * pl = 0x0;
1631 AliVCaloCells * cells = 0x0;
5ae09196 1632 //Select the Calorimeter of the photon
b5dbb99b 1633 if (fCalorimeter == "PHOS" )
34c16486 1634 {
1635 pl = GetPHOSClusters();
1636 cells = GetPHOSCells();
1637 }
5ae09196 1638 else if (fCalorimeter == "EMCAL")
34c16486 1639 {
1640 pl = GetEMCALClusters();
1641 cells = GetEMCALCells();
1642 }
57b97dc6 1643
34c16486 1644 if(!pl)
1645 {
5ae09196 1646 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1647 return;
1648 }
233e0df8 1649
477d6cee 1650 TLorentzVector mom ;
b5dbb99b 1651 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
1652 {
0ae57829 1653 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
477d6cee 1654
f8006433 1655 Int_t evtIndex = 0 ;
b5dbb99b 1656 if (GetMixedEvent())
1657 {
f8006433 1658 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1659 }
34c16486 1660
5025c139 1661 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
521636d2 1662
57b97dc6 1663 //Get Momentum vector,
34c16486 1664 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1665 {
1666 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
1667 }//Assume that come from vertex in straight line
1668 else
1669 {
f8006433 1670 Double_t vertex[]={0,0,0};
1671 calo->GetMomentum(mom,vertex) ;
1672 }
233e0df8 1673
57b97dc6 1674 //If too small or big pt, skip it
a426f453 1675 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
34c16486 1676
477d6cee 1677 //Check acceptance selection
b5dbb99b 1678 if(IsFiducialCutOn())
1679 {
ff45398a 1680 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 1681 if(! in ) continue ;
1682 }
1683
1684 //Create AOD for analysis
1685 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
c8fe2783 1686 aodpi0.SetLabel(calo->GetLabel());
3c1d9afb 1687
477d6cee 1688 //Set the indeces of the original caloclusters
1689 aodpi0.SetCaloLabel(calo->GetID(),-1);
1690 aodpi0.SetDetector(fCalorimeter);
1691 if(GetDebug() > 1)
ff45398a 1692 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 1693
1694 //Check Distance to Bad channel, set bit.
c8fe2783 1695 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 1696 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
1697 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
1698 continue ;
1699
34c16486 1700 //.......................................
1701 // TOF cut, BE CAREFUL WITH THIS CUT
1702 Double_t tof = calo->GetTOF()*1e9;
1703 if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
1704
1705
a3aebfff 1706 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
477d6cee 1707
3c1d9afb 1708 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
477d6cee 1709 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3c1d9afb 1710 else aodpi0.SetDistToBad(0) ;
477d6cee 1711
1712 //Check PID
1713 //PID selection or bit setting
34c16486 1714 Int_t nMaxima = 0 ;
1715 Double_t mass = 0 , angle = 0;
bfdcf7fb 1716 Double_t e1 = 0 , e2 = 0;
34c16486 1717 //Skip matched clusters with tracks
1718 if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
477d6cee 1719
34c16486 1720 // Check if cluster is pi0 via cluster splitting
1721 aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
bfdcf7fb 1722 GetVertex(evtIndex),nMaxima,
1723 mass,angle,e1,e2));
34c16486 1724
1725 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
477d6cee 1726
34c16486 1727 // If cluster does not pass pid, not pi0, skip it.
1728 // TO DO, add option for Eta ... or conversions
1729 if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
1730
1731 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",
1732 aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
477d6cee 1733
1734 //Play with the MC stack if available
1735 //Check origin of the candidates
34c16486 1736 Int_t tag = 0 ;
3c1d9afb 1737 if(IsDataMC())
1738 {
477d6cee 1739 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
26e228ff 1740 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1741 {
f8006433 1742 //aodpi0.SetInputFileIndex(input);
26e228ff 1743 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodpi0.GetInputFileIndex());
57b97dc6 1744 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
1745 aodpi0.SetTag(tag);
1746 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
477d6cee 1747 }
1748 }//Work with stack also
1749
34c16486 1750 //Fill some histograms about shower shape
1751 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1752 {
bfdcf7fb 1753 Float_t asy =-10;
1754 if(e1+e2 > 0 ) asy = (e1-e2) / (e1+e2);
1755 FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
34c16486 1756 }
1757
477d6cee 1758 //Add AOD with pi0 object to aod branch
1759 AddAODParticle(aodpi0);
1760
1761 }//loop
1762
a3aebfff 1763 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
477d6cee 1764
1765}
e7fd282f 1766//______________________________________________
477d6cee 1767void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 1768{
477d6cee 1769 //Do analysis and fill histograms
691bdd02 1770
b5dbb99b 1771 if(!GetOutputAODBranch())
1772 {
a3aebfff 1773 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
477d6cee 1774 abort();
1775 }
1776 //Loop on stored AOD pi0
1777 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff 1778 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee 1779
b5dbb99b 1780 for(Int_t iaod = 0; iaod < naod ; iaod++)
1781 {
477d6cee 1782
1783 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
21a4b1c0 1784 Int_t pdg = pi0->GetIdentifiedParticleType();
9415d854 1785
1786 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
477d6cee 1787
1788 //Fill pi0 histograms
c4a7d28a 1789 Float_t ener = pi0->E();
1790 Float_t pt = pi0->Pt();
1791 Float_t phi = pi0->Phi();
57b97dc6 1792 if(phi < 0) phi+=TMath::TwoPi();
477d6cee 1793 Float_t eta = pi0->Eta();
1794
09273901 1795 fhPt ->Fill(pt);
1796 fhE ->Fill(ener);
477d6cee 1797
09273901 1798 fhEEta ->Fill(ener,eta);
1799 fhEPhi ->Fill(ener,phi);
1800 fhEtaPhi ->Fill(eta,phi);
7a972c0c 1801
b5dbb99b 1802 if(IsDataMC())
1803 {
af722ce4 1804 if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0))
1805 {
1806 fhPtMC ->Fill(pt);
1807 fhPhiMC ->Fill(pt,phi);
1808 fhEtaMC ->Fill(pt,eta);
1809 }
1810 else
1811 {
1812 fhPtMCNo ->Fill(pt);
1813 fhPhiMCNo ->Fill(pt,phi);
1814 fhEtaMCNo ->Fill(pt,eta);
1815
477d6cee 1816 }
1817 }//Histograms with MC
1818
1819 }// aod loop
1820
1821}
1822
477d6cee 1823//__________________________________________________________________
1824void AliAnaPi0EbE::Print(const Option_t * opt) const
1825{
1826 //Print some relevant parameters set for the analysis
1827 if(! opt)
1828 return;
1829
1830 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 1831 AliAnaCaloTrackCorrBaseClass::Print("");
477d6cee 1832 printf("Analysis Type = %d \n", fAnaType) ;
1833 if(fAnaType == kSSCalo){
1834 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1835 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1836 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1837 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1838 }
1839 printf(" \n") ;
1840
1841}
78a28af3 1842
78a28af3 1843