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