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