]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPi0EbE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add cluster disitribution histograms as a function of eta/phi and depending on V0...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / 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 **************************************************************************/
15/* $Id: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18// Class for the analysis of high pT pi0 event by event
19// Pi0 identified by one of the following:
20// -Invariant mass of 2 cluster in calorimeter
21// -Shower shape analysis in calorimeter
21a4b1c0 22// -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
477d6cee 23//
24// -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TList.h>
30#include <TClonesArray.h>
0c1383b5 31#include <TObjString.h>
477d6cee 32
33// --- Analysis system ---
34#include "AliAnaPi0EbE.h"
35#include "AliCaloTrackReader.h"
36#include "AliIsolationCut.h"
37#include "AliNeutralMesonSelection.h"
38#include "AliCaloPID.h"
39#include "AliMCAnalysisUtils.h"
477d6cee 40#include "AliStack.h"
ff45398a 41#include "AliFiducialCut.h"
477d6cee 42#include "TParticle.h"
0ae57829 43#include "AliVCluster.h"
477d6cee 44#include "AliAODEvent.h"
591cc579 45#include "AliAODMCParticle.h"
477d6cee 46
47ClassImp(AliAnaPi0EbE)
48
78a28af3 49//____________________________
477d6cee 50AliAnaPi0EbE::AliAnaPi0EbE() :
521636d2 51 AliAnaPartCorrBaseClass(), fAnaType(kIMCalo), fCalorimeter(""),
78a28af3 52 fMinDist(0.),fMinDist2(0.), fMinDist3(0.), fFillWeightHistograms(kFALSE),
1db06135 53 fInputAODGammaConvName(""),
521636d2 54 //Histograms
7a972c0c 55 fhPtPi0(0), fhEPi0(0),
56 fhEEtaPi0(0), fhEPhiPi0(0), fhEtaPhiPi0(0),
521636d2 57 //Shower shape histos
3bfcb597 58 fhEDispersion(0), fhELambda0(0), fhELambda1(0),
59 fhELambda0NoTRD(0), fhELambda0FracMaxCellCut(0),
60 fhEFracMaxCell(0), fhEFracMaxCellNoTRD(0),
42d47cb7 61 fhENCells(0), fhETime(0), fhEPairDiffTime(0),
521636d2 62 //MC histos
63 fhPtMCNoPi0(0), fhPhiMCNoPi0(0), fhEtaMCNoPi0(0),
78a28af3 64 fhPtMCPi0(0), fhPhiMCPi0(0), fhEtaMCPi0(0),
65 // Weight studies
66 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
67 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0)
477d6cee 68{
69 //default ctor
70
6db946bd 71 for(Int_t i = 0; i < 6; i++){
521636d2 72 fhEMCLambda0[i] = 0;
3bfcb597 73 fhEMCLambda0NoTRD[i]= 0;
74 fhEMCLambda0FracMaxCellCut[i]= 0;
75 fhEMCFracMaxCell[i] = 0;
521636d2 76 fhEMCLambda1[i] = 0;
77 fhEMCDispersion[i] = 0;
521636d2 78 }
79
78a28af3 80 //Weight studies
1a72f6c5 81 for(Int_t i =0; i < 14; i++){
78a28af3 82 fhLambda0ForW0[i] = 0;
1a72f6c5 83 //fhLambda1ForW0[i] = 0;
78a28af3 84 }
85
477d6cee 86 //Initialize parameters
87 InitParameters();
88
89}
477d6cee 90
42d47cb7 91//_____________________________________________________________________________________
92void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, const Int_t tag){
93
94 // Fill shower shape, timing and other histograms for selected clusters from decay
95
96 Float_t e = cluster->E();
97 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
98 Float_t l0 = cluster->GetM02();
99 Float_t l1 = cluster->GetM20();
100 Int_t nSM = GetModuleNumber(cluster);
101
102 AliVCaloCells * cell = 0x0;
103 if(fCalorimeter == "PHOS")
104 cell = GetPHOSCells();
105 else
106 cell = GetEMCALCells();
107
108 Float_t maxCellFraction = 0;
109 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
110 fhEFracMaxCell->Fill(e,maxCellFraction);
111
112 FillWeightHistograms(cluster);
113
114 fhEDispersion->Fill(e, disp);
115 fhELambda0 ->Fill(e, l0 );
116 fhELambda1 ->Fill(e, l1 );
117
118 if(fCalorimeter=="EMCAL" && nSM < 6) {
119 fhELambda0NoTRD->Fill(e, l0 );
120 fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
121 }
122
123 if(maxCellFraction < 0.5)
124 fhELambda0FracMaxCellCut->Fill(e, l0 );
125
126 fhETime ->Fill(e, cluster->GetTOF()*1.e9);
127 fhENCells->Fill(e, cluster->GetNCells());
128
129 if(IsDataMC()) {
130 //Photon1
131 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ){
c5693f62 132 fhEMCLambda0[kmcPi0] ->Fill(e, l0);
133 fhEMCLambda1[kmcPi0] ->Fill(e, l1);
134 fhEMCDispersion[kmcPi0] ->Fill(e, disp);
42d47cb7 135
c5693f62 136 fhEMCFracMaxCell[kmcPi0]->Fill(e,maxCellFraction);
42d47cb7 137 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 138 fhEMCLambda0NoTRD[kmcPi0]->Fill(e, l0 );
42d47cb7 139 if(maxCellFraction < 0.5)
c5693f62 140 fhEMCLambda0FracMaxCellCut[kmcPi0]->Fill(e, l0 );
42d47cb7 141
142 }//pi0
143 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ){
c5693f62 144 fhEMCLambda0[kmcEta] ->Fill(e, l0);
145 fhEMCLambda1[kmcEta] ->Fill(e, l1);
146 fhEMCDispersion[kmcEta] ->Fill(e, disp);
147 fhEMCFracMaxCell[kmcEta]->Fill(e,maxCellFraction);
42d47cb7 148 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 149 fhEMCLambda0NoTRD[kmcEta]->Fill(e, l0 );
42d47cb7 150 if(maxCellFraction < 0.5)
c5693f62 151 fhEMCLambda0FracMaxCellCut[kmcEta]->Fill(e, l0 );
42d47cb7 152 }//eta
153 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
154 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) ){
c5693f62 155 fhEMCLambda0[kmcConversion] ->Fill(e, l0);
156 fhEMCLambda1[kmcConversion] ->Fill(e, l1);
157 fhEMCDispersion[kmcConversion] ->Fill(e, disp);
158 fhEMCFracMaxCell[kmcConversion]->Fill(e,maxCellFraction);
42d47cb7 159 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 160 fhEMCLambda0NoTRD[kmcConversion]->Fill(e, l0 );
42d47cb7 161 if(maxCellFraction < 0.5)
c5693f62 162 fhEMCLambda0FracMaxCellCut[kmcConversion]->Fill(e, l0 );
42d47cb7 163 }//conversion photon
164 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ){
c5693f62 165 fhEMCLambda0[kmcPhoton] ->Fill(e, l0);
166 fhEMCLambda1[kmcPhoton] ->Fill(e, l1);
167 fhEMCDispersion[kmcPhoton] ->Fill(e, disp);
168 fhEMCFracMaxCell[kmcPhoton]->Fill(e,maxCellFraction);
42d47cb7 169 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 170 fhEMCLambda0NoTRD[kmcPhoton]->Fill(e, l0 );
42d47cb7 171 if(maxCellFraction < 0.5)
c5693f62 172 fhEMCLambda0FracMaxCellCut[kmcPhoton]->Fill(e, l0 );
42d47cb7 173 }//photon no conversion
174 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)){
c5693f62 175 fhEMCLambda0[kmcElectron] ->Fill(e, l0);
176 fhEMCLambda1[kmcElectron] ->Fill(e, l1);
177 fhEMCDispersion[kmcElectron] ->Fill(e, disp);
178 fhEMCFracMaxCell[kmcElectron]->Fill(e,maxCellFraction);
42d47cb7 179 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 180 fhEMCLambda0NoTRD[kmcElectron]->Fill(e, l0 );
42d47cb7 181 if(maxCellFraction < 0.5)
c5693f62 182 fhEMCLambda0FracMaxCellCut[kmcElectron]->Fill(e, l0 );
42d47cb7 183 }//electron
184 else {
c5693f62 185 fhEMCLambda0[kmcHadron] ->Fill(e, l0);
186 fhEMCLambda1[kmcHadron] ->Fill(e, l1);
187 fhEMCDispersion[kmcHadron] ->Fill(e, disp);
188 fhEMCFracMaxCell[kmcHadron]->Fill(e,maxCellFraction);
42d47cb7 189 if(fCalorimeter=="EMCAL" && nSM < 6)
c5693f62 190 fhEMCLambda0NoTRD[kmcHadron]->Fill(e, l0 );
42d47cb7 191 if(maxCellFraction < 0.5)
c5693f62 192 fhEMCLambda0FracMaxCellCut[kmcHadron]->Fill(e, l0 );
42d47cb7 193 }//other particles
194 }//MC
195}
196
197//________________________________________________________
198void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
199{
200 // Calculate weights and fill histograms
201
202 if(!fFillWeightHistograms || GetMixedEvent()) return;
203
204 AliVCaloCells* cells = 0;
205 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
206 else cells = GetPHOSCells();
207
208 // First recalculate energy in case non linearity was applied
209 Float_t energy = 0;
210 Float_t ampMax = 0;
211 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
212
213 Int_t id = clus->GetCellsAbsId()[ipos];
214
215 //Recalibrate cell energy if needed
216 Float_t amp = cells->GetCellAmplitude(id);
217 RecalibrateCellAmplitude(amp,id);
218
219 energy += amp;
220
221 if(amp> ampMax)
222 ampMax = amp;
223
224 } // energy loop
225
226 if(energy <=0 ) {
227 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
228 return;
229 }
230
231 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
232 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
233
234 //Get the ratio and log ratio to all cells in cluster
235 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
236 Int_t id = clus->GetCellsAbsId()[ipos];
237
238 //Recalibrate cell energy if needed
239 Float_t amp = cells->GetCellAmplitude(id);
240 RecalibrateCellAmplitude(amp,id);
241
242 fhECellClusterRatio ->Fill(energy,amp/energy);
243 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
244 }
245
246 //Recalculate shower shape for different W0
247 if(fCalorimeter=="EMCAL"){
248
249 Float_t l0org = clus->GetM02();
250 Float_t l1org = clus->GetM20();
251 Float_t dorg = clus->GetDispersion();
252
1a72f6c5 253 for(Int_t iw = 0; iw < 14; iw++){
254 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
42d47cb7 255 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
256
257 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 258 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
42d47cb7 259
260 } // w0 loop
261
262 // Set the original values back
263 clus->SetM02(l0org);
264 clus->SetM20(l1org);
265 clus->SetDispersion(dorg);
266
267 }// EMCAL
268}
269
78a28af3 270//___________________________________________
0c1383b5 271TObjString * AliAnaPi0EbE::GetAnalysisCuts()
272{
273 //Save parameters used for analysis
521636d2 274 TString parList ; //this will be list of parameters used for this analysis.
275 const Int_t buffersize = 255;
276 char onePar[buffersize] ;
277
278 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
279 parList+=onePar ;
280 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
281 parList+=onePar ;
282
283 if(fAnaType == kSSCalo){
284 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
285 parList+=onePar ;
286 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
287 parList+=onePar ;
288 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
289 parList+=onePar ;
290 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
291 parList+=onePar ;
292 }
293
294 //Get parameters set in base class.
295 parList += GetBaseParametersList() ;
296
297 //Get parameters set in PID class.
298 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
299
300 return new TObjString(parList) ;
0c1383b5 301}
302
78a28af3 303//_____________________________________________
477d6cee 304TList * AliAnaPi0EbE::GetCreateOutputObjects()
305{
306 // Create histograms to be saved in output file and
307 // store them in outputContainer
308 TList * outputContainer = new TList() ;
309 outputContainer->SetName("Pi0EbEHistos") ;
310
57b97dc6 311 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
312 Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
313 Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
314 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
521636d2 315 Int_t tdbins = GetHistoDiffTimeBins() ; Float_t tdmax = GetHistoDiffTimeMax(); Float_t tdmin = GetHistoDiffTimeMin();
42d47cb7 316 Int_t tbins = GetHistoTimeBins() ; Float_t tmax = GetHistoTimeMax(); Float_t tmin = GetHistoTimeMin();
317 Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
318
477d6cee 319 fhPtPi0 = new TH1F("hPtPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
320 fhPtPi0->SetYTitle("N");
321 fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
322 outputContainer->Add(fhPtPi0) ;
521636d2 323
c4a7d28a 324 fhEPi0 = new TH1F("hEPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
325 fhEPi0->SetYTitle("N");
326 fhEPi0->SetXTitle("E #pi^{0}(GeV)");
327 outputContainer->Add(fhEPi0) ;
477d6cee 328
7a972c0c 329 fhEPhiPi0 = new TH2F
330 ("hEPhiPi0","Selected #pi^{0} pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
331 fhEPhiPi0->SetYTitle("#phi #pi^{0}(GeV)");
332 fhEPhiPi0->SetXTitle("E (GeV) #pi^{0}(GeV)");
333 outputContainer->Add(fhEPhiPi0) ;
334
335 fhEEtaPi0 = new TH2F
336 ("hEEtaPi0","Selected #pi^{0} pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
337 fhEEtaPi0->SetYTitle("#eta #pi^{0}(GeV)");
338 fhEEtaPi0->SetXTitle("E (GeV) #pi^{0}(GeV)");
339 outputContainer->Add(fhEEtaPi0) ;
340
341 fhEtaPhiPi0 = new TH2F
342 ("hEtaPhiPi0","Selected #pi^{0} pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
343 fhEtaPhiPi0->SetYTitle("#phi #pi^{0}(GeV)");
344 fhEtaPhiPi0->SetXTitle("#eta #pi^{0}(GeV)");
345 outputContainer->Add(fhEtaPhiPi0) ;
57b97dc6 346
c4a7d28a 347 ////////
57b97dc6 348
1db06135 349 if(fAnaType == kIMCalo || fAnaType == kIMCaloTracks ){
c4a7d28a 350
521636d2 351 fhEDispersion = new TH2F
352 ("hEDispersion","Selected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
353 fhEDispersion->SetYTitle("D^{2}");
354 fhEDispersion->SetXTitle("E (GeV)");
355 outputContainer->Add(fhEDispersion) ;
356
357 fhELambda0 = new TH2F
358 ("hELambda0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
359 fhELambda0->SetYTitle("#lambda_{0}^{2}");
360 fhELambda0->SetXTitle("E (GeV)");
361 outputContainer->Add(fhELambda0) ;
3bfcb597 362
42d47cb7 363 fhELambda1 = new TH2F
364 ("hELambda1","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
365 fhELambda1->SetYTitle("#lambda_{1}^{2}");
366 fhELambda1->SetXTitle("E (GeV)");
367 outputContainer->Add(fhELambda1) ;
368
3bfcb597 369 fhELambda0FracMaxCellCut = new TH2F
370 ("hELambda0FracMaxCellCut","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
371 fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
372 fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
373 outputContainer->Add(fhELambda0FracMaxCellCut) ;
374
375 fhEFracMaxCell = new TH2F
376 ("hEFracMaxCell","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
377 fhEFracMaxCell->SetYTitle("Fraction");
378 fhEFracMaxCell->SetXTitle("E (GeV)");
379 outputContainer->Add(fhEFracMaxCell) ;
380
381 if(fCalorimeter=="EMCAL"){
382 fhELambda0NoTRD = new TH2F
383 ("hELambda0NoTRD","Selected #pi^{0} pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
384 fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
385 fhELambda0NoTRD->SetXTitle("E (GeV)");
386 outputContainer->Add(fhELambda0NoTRD) ;
387
388 fhEFracMaxCellNoTRD = new TH2F
389 ("hEFracMaxCellNoTRD","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
390 fhEFracMaxCellNoTRD->SetYTitle("Fraction");
391 fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
392 outputContainer->Add(fhEFracMaxCellNoTRD) ;
393 }
521636d2 394
42d47cb7 395 fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
396 fhENCells->SetXTitle("E (GeV)");
397 fhENCells->SetYTitle("# of cells in cluster");
398 outputContainer->Add(fhENCells);
399
400 fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
401 fhETime->SetXTitle("E (GeV)");
402 fhETime->SetYTitle(" t (ns)");
403 outputContainer->Add(fhETime);
521636d2 404
e7fd282f 405 }// Invariant mass analysis in calorimeters and calorimeter + conversion photons
406
407 if(fAnaType == kIMCalo){
42d47cb7 408 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
409 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
410 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
411 outputContainer->Add(fhEPairDiffTime);
e7fd282f 412 }
477d6cee 413
78a28af3 414 if(fFillWeightHistograms){
415
416 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
417 nptbins,ptmin,ptmax, 100,0,1.);
418 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
419 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
420 outputContainer->Add(fhECellClusterRatio);
421
422 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 423 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 424 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 425 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 426 outputContainer->Add(fhECellClusterLogRatio);
427
428 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
429 nptbins,ptmin,ptmax, 100,0,1.);
430 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
431 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
432 outputContainer->Add(fhEMaxCellClusterRatio);
433
434 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1a72f6c5 435 nptbins,ptmin,ptmax, 100,-10,0);
78a28af3 436 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 437 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
78a28af3 438 outputContainer->Add(fhEMaxCellClusterLogRatio);
439
1a72f6c5 440 for(Int_t iw = 0; iw < 14; iw++){
441 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 442 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
443 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
444 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
445 outputContainer->Add(fhLambda0ForW0[iw]);
446
1a72f6c5 447// 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),
448// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
449// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
450// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
451// outputContainer->Add(fhLambda1ForW0[iw]);
78a28af3 452
453 }
454 }
455
477d6cee 456 if(IsDataMC()) {
457 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
458 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
459
460 fhPtMCPi0 = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax);
461 fhPtMCPi0->SetYTitle("N");
462 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
463 outputContainer->Add(fhPtMCPi0) ;
464
465 fhPhiMCPi0 = new TH2F
5ae09196 466 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 467 fhPhiMCPi0->SetYTitle("#phi");
468 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
469 outputContainer->Add(fhPhiMCPi0) ;
470
471 fhEtaMCPi0 = new TH2F
5ae09196 472 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 473 fhEtaMCPi0->SetYTitle("#eta");
474 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
475 outputContainer->Add(fhEtaMCPi0) ;
476
477 fhPtMCNoPi0 = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax);
478 fhPtMCNoPi0->SetYTitle("N");
479 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
480 outputContainer->Add(fhPtMCNoPi0) ;
481
482 fhPhiMCNoPi0 = new TH2F
5ae09196 483 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 484 fhPhiMCNoPi0->SetYTitle("#phi");
485 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
486 outputContainer->Add(fhPhiMCNoPi0) ;
487
488 fhEtaMCNoPi0 = new TH2F
5ae09196 489 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 490 fhEtaMCNoPi0->SetYTitle("#eta");
491 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
492 outputContainer->Add(fhEtaMCNoPi0) ;
493
c4a7d28a 494 if(fAnaType == kIMCalo){
6db946bd 495 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
496 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
497 for(Int_t i = 0; i < 6; i++){
521636d2 498
499 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
500 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
501 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
502 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
503 fhEMCLambda0[i]->SetXTitle("E (GeV)");
504 outputContainer->Add(fhEMCLambda0[i]) ;
505
3bfcb597 506 if(fCalorimeter=="EMCAL"){
507 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
508 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
509 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
510 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
511 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
512 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
513 }
514
515 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
516 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
517 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
518 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
519 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
520 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
521
522 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
523 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
524 nptbins,ptmin,ptmax,100,0,1);
525 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
526 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
527 outputContainer->Add(fhEMCFracMaxCell[i]) ;
528
521636d2 529 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
530 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
531 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
532 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
533 fhEMCLambda1[i]->SetXTitle("E (GeV)");
534 outputContainer->Add(fhEMCLambda1[i]) ;
7c65ad18 535
521636d2 536 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
537 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
538 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
539 fhEMCDispersion[i]->SetYTitle("D^{2}");
540 fhEMCDispersion[i]->SetXTitle("E (GeV)");
541 outputContainer->Add(fhEMCDispersion[i]) ;
7c65ad18 542
521636d2 543 }//
c4a7d28a 544
545 }//kIMCalo
521636d2 546 } //Not MC reader
477d6cee 547 }//Histos with MC
548
549
550 //Keep neutral meson selection histograms if requiered
551 //Setting done in AliNeutralMesonSelection
552
553 if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
554
555 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
556 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
557 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
5ae09196 558 delete nmsHistos;
a14fd526 559
477d6cee 560 }
561
477d6cee 562 return outputContainer ;
563
564}
565
521636d2 566//____________________________________________________________________________
567void AliAnaPi0EbE::Init()
568{
569 //Init
570 //Do some checks
571 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
572 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
573 abort();
574 }
575 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
576 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
577 abort();
578 }
579
580}
581
582//____________________________________________________________________________
583void AliAnaPi0EbE::InitParameters()
584{
585 //Initialize the parameters of the analysis.
586 AddToHistogramsName("AnaPi0EbE_");
587
1db06135 588 fInputAODGammaConvName = "PhotonsCTS" ;
521636d2 589 fAnaType = kIMCalo ;
590 fCalorimeter = "EMCAL" ;
591 fMinDist = 2.;
592 fMinDist2 = 4.;
593 fMinDist3 = 5.;
594
595}
596
477d6cee 597//__________________________________________________________________
598void AliAnaPi0EbE::MakeAnalysisFillAOD()
599{
600 //Do analysis and fill aods
601
602 switch(fAnaType)
521636d2 603 {
477d6cee 604 case kIMCalo:
605 MakeInvMassInCalorimeter();
606 break;
607
608 case kSSCalo:
609 MakeShowerShapeIdentification();
610 break;
611
612 case kIMCaloTracks:
613 MakeInvMassInCalorimeterAndCTS();
614 break;
615
521636d2 616 }
477d6cee 617}
618
42d47cb7 619//____________________________________________
477d6cee 620void AliAnaPi0EbE::MakeInvMassInCalorimeter()
621{
57b97dc6 622 //Do analysis and fill aods
623 //Search for the photon decay in calorimeters
624 //Read photon list from AOD, produced in class AliAnaPhoton
625 //Check if 2 photons have the mass of the pi0.
477d6cee 626
627 TLorentzVector mom1;
628 TLorentzVector mom2;
629 TLorentzVector mom ;
591cc579 630 Int_t tag1 = 0;
631 Int_t tag2 = 0;
632 Int_t tag = 0;
477d6cee 633
634 if(!GetInputAODBranch()){
a3aebfff 635 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
477d6cee 636 abort();
637 }
f8006433 638
42d47cb7 639 //Get shower shape information of clusters
640 TObjArray *clusters = 0;
641 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
642 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
643
c4a7d28a 644 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
477d6cee 645 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
c8fe2783 646
c4a7d28a 647 //Vertex cut in case of mixed events
c8fe2783 648 Int_t evtIndex1 = 0 ;
649 if(GetMixedEvent())
650 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
5025c139 651 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
477d6cee 652 mom1 = *(photon1->Momentum());
653
42d47cb7 654 //Get original cluster, to recover some information
1db06135 655 Int_t iclus = -1;
656 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
42d47cb7 657
1db06135 658 if(!cluster1){
42d47cb7 659 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
660 return;
9ab9e937 661 }
c4a7d28a 662
663 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++){
477d6cee 664
a3aebfff 665 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
c8fe2783 666 Int_t evtIndex2 = 0 ;
667 if(GetMixedEvent())
668 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
669 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
670 continue ;
5025c139 671 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
477d6cee 672 mom2 = *(photon2->Momentum());
c4a7d28a 673
1db06135 674 //Get original cluster, to recover some information
675 Int_t iclus2;
676 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
42d47cb7 677
1db06135 678 if(!cluster2){
42d47cb7 679 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1db06135 680 return;
9ab9e937 681 }
c4a7d28a 682
42d47cb7 683 Float_t e1 = photon1->E();
684 Float_t e2 = photon2->E();
685
686 //Select clusters with good time window difference
687 Float_t tof1 = cluster1->GetTOF()*1e9;;
688 Float_t tof2 = cluster2->GetTOF()*1e9;;
689 Double_t t12diff = tof1-tof2;
690 fhEPairDiffTime->Fill(e1+e2, t12diff);
691 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
692
57b97dc6 693 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 694 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
c8fe2783 695 {
696 if(GetDebug()>1)
697 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 698
57b97dc6 699 //Play with the MC stack if available
c8fe2783 700 if(IsDataMC()){
57b97dc6 701 //Check origin of the candidates
c8fe2783 702 Int_t label1 = photon1->GetLabel();
703 Int_t label2 = photon2->GetLabel();
0db79843 704 if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
705 if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
c8fe2783 706
707 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
0db79843 708 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
709 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
c8fe2783 710
57b97dc6 711 //Check if pi0 mother is the same
c8fe2783 712 if(GetReader()->ReadStack()){
0db79843 713 if(label1>=0){
714 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
715 label1 = mother1->GetFirstMother();
716 //mother1 = GetMCStack()->Particle(label1);//pi0
717 }
718 if(label2>=0){
719 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
720 label2 = mother2->GetFirstMother();
721 //mother2 = GetMCStack()->Particle(label2);//pi0
521636d2 722 }
c8fe2783 723 }
f8006433 724 else if(GetReader()->ReadAODMCParticles()){//&& (input > -1)){
0db79843 725 if(label1>=0){
726 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
727 label1 = mother1->GetMother();
728 //mother1 = GetMCStack()->Particle(label1);//pi0
521636d2 729 }
0db79843 730 if(label2>=0){
731 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
732 label2 = mother2->GetMother();
733 //mother2 = GetMCStack()->Particle(label2);//pi0
734 }
c8fe2783 735 }
736
57b97dc6 737 //printf("mother1 %d, mother2 %d\n",label1,label2);
0db79843 738 if(label1 == label2 && label1>=0)
c8fe2783 739 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
740 }
741 }//Work with stack also
742
78a28af3 743
57b97dc6 744 //Fill some histograms about shower shape
9ab9e937 745 if(clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
42d47cb7 746 FillSelectedClusterHistograms(cluster1, tag1);
747 FillSelectedClusterHistograms(cluster2, tag2);
748 }
521636d2 749
803d06a8 750 // Tag both photons as decay
751 photon1->SetTagged(kTRUE);
752 photon2->SetTagged(kTRUE);
753
57b97dc6 754 //Create AOD for analysis
c8fe2783 755 mom = mom1+mom2;
756 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
57b97dc6 757 //pi0.SetLabel(calo->GetLabel());
21a4b1c0 758 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
c8fe2783 759 pi0.SetDetector(photon1->GetDetector());
760 pi0.SetTag(tag);
57b97dc6 761 //Set the indeces of the original caloclusters
c8fe2783 762 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
f8006433 763 //pi0.SetInputFileIndex(input);
c8fe2783 764 AddAODParticle(pi0);
765 }//pi0
57b97dc6 766
477d6cee 767 }//2n photon loop
768
769 }//1st photon loop
770
a3aebfff 771 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
477d6cee 772
773}
774
e7fd282f 775//__________________________________________________
477d6cee 776void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
777{
778 //Do analysis and fill aods
779 //Search for the photon decay in calorimeters
780 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
781 //Check if 2 photons have the mass of the pi0.
782
783 TLorentzVector mom1;
784 TLorentzVector mom2;
785 TLorentzVector mom ;
591cc579 786 Int_t tag1 = 0;
787 Int_t tag2 = 0;
788 Int_t tag = 0;
5025c139 789 Int_t evtIndex = 0;
1db06135 790
791 // Check calorimeter input
477d6cee 792 if(!GetInputAODBranch()){
a3aebfff 793 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
477d6cee 794 abort();
795 }
57b97dc6 796
1db06135 797 // Get the array with conversion photons
798 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
799 if(!inputAODGammaConv) {
800
801 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
802
803 if(!inputAODGammaConv) {
804 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
805
806 return;
807 }
808 }
809
810 //Get shower shape information of clusters
811 TObjArray *clusters = 0;
812 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
813 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
814
815 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
816 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
817 if(nCTS<=0 || nCalo <=0) {
818 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
819 return;
820 }
821
822 if(GetDebug() > 1)
823 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
824
825 // Do the loop, first calo, second CTS
477d6cee 826 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
827 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
828 mom1 = *(photon1->Momentum());
829
1db06135 830 //Get original cluster, to recover some information
831 Int_t iclus = -1;
832 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
833
834 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
835 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
5025c139 836 if(GetMixedEvent())
837 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
838 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
839
477d6cee 840 mom2 = *(photon2->Momentum());
57b97dc6 841
477d6cee 842 //Select good pair (good phi, pt cuts, aperture and invariant mass)
3bfcb597 843 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)){
57b97dc6 844 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());
845
846 if(IsDataMC()){
847 Int_t label1 = photon1->GetLabel();
848 Int_t label2 = photon2->GetLabel();
0db79843 849 if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
850 if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
57b97dc6 851 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
0db79843 852 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
853 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
57b97dc6 854 //Check if pi0 mother is the same
855
856 if(GetReader()->ReadStack()){
0db79843 857 if(label1>=0){
858 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
859 label1 = mother1->GetFirstMother();
860 //mother1 = GetMCStack()->Particle(label1);//pi0
861 }
862 if(label2>=0){
863 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
864 label2 = mother2->GetFirstMother();
865 //mother2 = GetMCStack()->Particle(label2);//pi0
866 }
57b97dc6 867 }
0db79843 868 else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
869 if(label1>=0){
870 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
871 label1 = mother1->GetMother();
872 //mother1 = GetMCStack()->Particle(label1);//pi0
873 }
874 if(label2>=0){
875 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
876 label2 = mother2->GetMother();
877 //mother2 = GetMCStack()->Particle(label2);//pi0
878 }
57b97dc6 879 }
880
881 //printf("mother1 %d, mother2 %d\n",label1,label2);
0db79843 882 if(label1 == label2 && label1>=0)
57b97dc6 883 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
884 }
885 }//Work with stack also
886
1db06135 887 //Fill some histograms about shower shape
888 if(cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
889 FillSelectedClusterHistograms(cluster, tag1);
890 }
803d06a8 891
892 // Tag both photons as decay
893 photon1->SetTagged(kTRUE);
894 photon2->SetTagged(kTRUE);
1db06135 895
57b97dc6 896 //Create AOD for analysis
897 mom = mom1+mom2;
898 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
899 //pi0.SetLabel(calo->GetLabel());
21a4b1c0 900 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
57b97dc6 901 pi0.SetDetector(photon1->GetDetector());
902 pi0.SetTag(tag);
903 //Set the indeces of the original tracks or caloclusters
904 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
905 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
f8006433 906 //pi0.SetInputFileIndex(input);
57b97dc6 907 AddAODParticle(pi0);
477d6cee 908 }//pi0
909 }//2n photon loop
910
911 }//1st photon loop
912
a3aebfff 913 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
477d6cee 914
915}
916
917
e7fd282f 918//_________________________________________________
477d6cee 919void AliAnaPi0EbE::MakeShowerShapeIdentification()
920{
921 //Search for pi0 in fCalorimeter with shower shape analysis
922
4a745797 923 TObjArray * pl = 0x0;
5ae09196 924 //Select the Calorimeter of the photon
925 if(fCalorimeter == "PHOS")
be518ab0 926 pl = GetPHOSClusters();
5ae09196 927 else if (fCalorimeter == "EMCAL")
be518ab0 928 pl = GetEMCALClusters();
57b97dc6 929
5ae09196 930 if(!pl) {
931 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
932 return;
933 }
233e0df8 934
477d6cee 935 TLorentzVector mom ;
936 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
0ae57829 937 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
477d6cee 938
f8006433 939 Int_t evtIndex = 0 ;
940 if (GetMixedEvent()) {
941 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
942 }
5025c139 943 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
521636d2 944
57b97dc6 945 //Get Momentum vector,
f8006433 946 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
947 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
948 else{
949 Double_t vertex[]={0,0,0};
950 calo->GetMomentum(mom,vertex) ;
951 }
233e0df8 952
57b97dc6 953 //If too small or big pt, skip it
477d6cee 954 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
955 //Check acceptance selection
ff45398a 956 if(IsFiducialCutOn()){
957 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 958 if(! in ) continue ;
959 }
960
961 //Create AOD for analysis
962 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
c8fe2783 963 aodpi0.SetLabel(calo->GetLabel());
477d6cee 964 //Set the indeces of the original caloclusters
965 aodpi0.SetCaloLabel(calo->GetID(),-1);
966 aodpi0.SetDetector(fCalorimeter);
967 if(GetDebug() > 1)
ff45398a 968 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 969
970 //Check Distance to Bad channel, set bit.
c8fe2783 971 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 972 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
973 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
974 continue ;
975
a3aebfff 976 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
477d6cee 977
978 if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
979 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
980 else aodpi0.SetDistToBad(0) ;
981
982 //Check PID
983 //PID selection or bit setting
49b5c49b 984 if(IsCaloPIDOn()){
477d6cee 985 //Skip matched clusters with tracks
49b5c49b 986 if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
477d6cee 987
49b5c49b 988 // Get most probable PID, 2 options check bayesian PID weights or redo PID
989 // By default, redo PID
990
991 aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
477d6cee 992
21a4b1c0 993 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
477d6cee 994
995 //If cluster does not pass pid, not pi0, skip it.
21a4b1c0 996 if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
477d6cee 997
998 }
999 else{
49b5c49b 1000 //Set PID bits for later selection
477d6cee 1001 //GetPDG already called in SetPIDBits.
49b5c49b 1002 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils(), GetReader()->GetInputEvent());
a3aebfff 1003 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
477d6cee 1004 }
1005
21a4b1c0 1006 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
477d6cee 1007
1008 //Play with the MC stack if available
1009 //Check origin of the candidates
1010 if(IsDataMC()){
1011 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
57b97dc6 1012 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
f8006433 1013 //aodpi0.SetInputFileIndex(input);
57b97dc6 1014 Int_t tag =0;
1015 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
1016 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
1017 aodpi0.SetTag(tag);
1018 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
477d6cee 1019 }
1020 }//Work with stack also
1021
1022 //Add AOD with pi0 object to aod branch
1023 AddAODParticle(aodpi0);
1024
1025 }//loop
1026
a3aebfff 1027 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
477d6cee 1028
1029}
e7fd282f 1030//______________________________________________
477d6cee 1031void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 1032{
477d6cee 1033 //Do analysis and fill histograms
691bdd02 1034
477d6cee 1035 if(!GetOutputAODBranch()){
a3aebfff 1036 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
477d6cee 1037 abort();
1038 }
1039 //Loop on stored AOD pi0
1040 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff 1041 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee 1042
1043 for(Int_t iaod = 0; iaod < naod ; iaod++){
1044
1045 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
21a4b1c0 1046 Int_t pdg = pi0->GetIdentifiedParticleType();
9415d854 1047
1048 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
477d6cee 1049
1050 //Fill pi0 histograms
c4a7d28a 1051 Float_t ener = pi0->E();
1052 Float_t pt = pi0->Pt();
1053 Float_t phi = pi0->Phi();
57b97dc6 1054 if(phi < 0) phi+=TMath::TwoPi();
477d6cee 1055 Float_t eta = pi0->Eta();
1056
c4a7d28a 1057 fhPtPi0 ->Fill(pt);
1058 fhEPi0 ->Fill(ener);
477d6cee 1059
7a972c0c 1060 fhEEtaPi0 ->Fill(ener,eta);
1061 fhEPhiPi0 ->Fill(ener,phi);
1062 fhEtaPhiPi0 ->Fill(eta,phi);
1063
477d6cee 1064 if(IsDataMC()){
1065 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
57b97dc6 1066 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1067 if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
1068 fhPtMCPi0 ->Fill(pt);
1069 fhPhiMCPi0 ->Fill(pt,phi);
1070 fhEtaMCPi0 ->Fill(pt,eta);
1071 }
1072 else{
1073 fhPtMCNoPi0 ->Fill(pt);
1074 fhPhiMCNoPi0 ->Fill(pt,phi);
1075 fhEtaMCNoPi0 ->Fill(pt,eta);
1076 }
477d6cee 1077 }
1078 }//Histograms with MC
1079
1080 }// aod loop
1081
1082}
1083
477d6cee 1084//__________________________________________________________________
1085void AliAnaPi0EbE::Print(const Option_t * opt) const
1086{
1087 //Print some relevant parameters set for the analysis
1088 if(! opt)
1089 return;
1090
1091 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1092 AliAnaPartCorrBaseClass::Print("");
1093 printf("Analysis Type = %d \n", fAnaType) ;
1094 if(fAnaType == kSSCalo){
1095 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1096 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1097 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1098 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1099 }
1100 printf(" \n") ;
1101
1102}
78a28af3 1103
1104//___________________________________________________________________________________
1105void AliAnaPi0EbE::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
1106{
1107 //Recaculate cell energy if recalibration factor
1108
1109 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1110 Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
1111
1112 if (GetCaloUtils()->IsRecalibrationOn()) {
1113 if(fCalorimeter == "PHOS") {
1114 amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1115 }
1116 else {
1117 amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1118 }
1119 }
1120}
1121
78a28af3 1122