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