PartCorr split in 2 Base and Dep; coding violations corrected; PHOS geometry can...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0EbE.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes iGetEntriesFast(s hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 /* $Id: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */\r
16 \r
17 //_________________________________________________________________________\r
18 // Class for the analysis of high pT pi0 event by event\r
19 // Pi0 identified by one of the following:\r
20 //  -Invariant mass of 2 cluster in calorimeter\r
21 //  -Shower shape analysis in calorimeter\r
22 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in TPC (in near future)\r
23 //\r
24 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)\r
25 //////////////////////////////////////////////////////////////////////////////\r
26   \r
27   \r
28 // --- ROOT system --- \r
29 #include <TList.h>\r
30 //#include "Riostream.h"\r
31 \r
32 // --- Analysis system --- \r
33 #include "AliAnaPi0EbE.h" \r
34 #include "AliLog.h"\r
35 #include "AliCaloTrackReader.h"\r
36 #include "AliIsolationCut.h"\r
37 #include "AliNeutralMesonSelection.h"\r
38 #include "AliCaloPID.h"\r
39 #include "AliAODPWG4ParticleCorrelation.h"\r
40 \r
41 ClassImp(AliAnaPi0EbE)\r
42   \r
43 //____________________________________________________________________________\r
44 AliAnaPi0EbE::AliAnaPi0EbE() : \r
45 AliAnaPartCorrBaseClass(),  fAnaType(kIMCalo),fCalorimeter(""),\r
46 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),       \r
47 fInputAODGammaConv(0x0),fInputAODGammaConvName(""),\r
48 fhPtPi0(0),fhPhiPi0(0),fhEtaPi0(0), \r
49 fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0), \r
50 fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)\r
51 {\r
52         //default ctor\r
53         \r
54         //Initialize parameters\r
55         InitParameters();\r
56         \r
57 }\r
58 \r
59 //____________________________________________________________________________\r
60 AliAnaPi0EbE::AliAnaPi0EbE(const AliAnaPi0EbE & p) : \r
61 AliAnaPartCorrBaseClass(p),  fAnaType(p.fAnaType), fCalorimeter(p.fCalorimeter),\r
62 fMinDist(p.fMinDist),fMinDist2(p.fMinDist2), fMinDist3(p.fMinDist3),\r
63 fInputAODGammaConv(new TClonesArray(*p.fInputAODGammaConv)),\r
64 fInputAODGammaConvName(p.fInputAODGammaConvName),\r
65 fhPtPi0(p.fhPtPi0),fhPhiPi0(p.fhPhiPi0),fhEtaPi0(p.fhEtaPi0), \r
66 fhPtMCNoPi0(p.fhPtMCNoPi0),fhPhiMCNoPi0(p.fhPhiMCNoPi0),fhEtaMCNoPi0(p.fhEtaMCNoPi0), \r
67 fhPtMCPi0(p.fhPtMCPi0),fhPhiMCPi0(p.fhPhiMCPi0),fhEtaMCPi0(p.fhEtaMCPi0) \r
68 {\r
69         // cpy ctor\r
70 }\r
71 \r
72 //_________________________________________________________________________\r
73 AliAnaPi0EbE & AliAnaPi0EbE::operator = (const AliAnaPi0EbE & p)\r
74 {\r
75         // assignment operator\r
76         \r
77         if(&p == this) return *this;\r
78         \r
79         fAnaType     = p.fAnaType ;\r
80         fCalorimeter = p.fCalorimeter ;\r
81         fMinDist     = p.fMinDist;\r
82         fMinDist2    = p.fMinDist2;\r
83         fMinDist3    = p.fMinDist3;\r
84         \r
85         fInputAODGammaConv     = new TClonesArray(*p.fInputAODGammaConv) ;\r
86         fInputAODGammaConvName = p.fInputAODGammaConvName ;\r
87         \r
88         fhPtPi0 = p.fhPtPi0; fhPhiPi0 = p.fhPhiPi0; fhEtaPi0 = p.fhEtaPi0;  \r
89         fhPtMCNoPi0 = p.fhPtMCNoPi0; fhPhiMCNoPi0 = p.fhPhiMCNoPi0; fhEtaMCNoPi0 = p.fhEtaMCPi0;  \r
90         fhPtMCPi0 = p.fhPtMCPi0; fhPhiMCPi0 = p.fhPhiMCPi0; fhEtaMCPi0 = p.fhEtaMCPi0;  \r
91         \r
92         return *this;\r
93         \r
94 }\r
95 \r
96 //____________________________________________________________________________\r
97 AliAnaPi0EbE::~AliAnaPi0EbE() \r
98 {\r
99         //dtor\r
100         if(fInputAODGammaConv){\r
101                 fInputAODGammaConv->Clear() ; \r
102                 delete fInputAODGammaConv ;\r
103         }\r
104 }\r
105 \r
106 //________________________________________________________________________\r
107 TList *  AliAnaPi0EbE::GetCreateOutputObjects()\r
108 {  \r
109         // Create histograms to be saved in output file and \r
110         // store them in outputContainer\r
111         TList * outputContainer = new TList() ; \r
112         outputContainer->SetName("Pi0EbEHistos") ; \r
113         \r
114         Int_t nptbins  = GetHistoNPtBins();\r
115         Int_t nphibins = GetHistoNPhiBins();\r
116         Int_t netabins = GetHistoNEtaBins();\r
117         Float_t ptmax  = GetHistoPtMax();\r
118         Float_t phimax = GetHistoPhiMax();\r
119         Float_t etamax = GetHistoEtaMax();\r
120         Float_t ptmin  = GetHistoPtMin();\r
121         Float_t phimin = GetHistoPhiMin();\r
122         Float_t etamin = GetHistoEtaMin();      \r
123         \r
124         fhPtPi0  = new TH1F("hPtPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax); \r
125         fhPtPi0->SetYTitle("N");\r
126         fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
127         outputContainer->Add(fhPtPi0) ; \r
128         \r
129         fhPhiPi0  = new TH2F\r
130     ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
131         fhPhiPi0->SetYTitle("#phi");\r
132         fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
133         outputContainer->Add(fhPhiPi0) ; \r
134         \r
135         fhEtaPi0  = new TH2F\r
136     ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
137         fhEtaPi0->SetYTitle("#eta");\r
138         fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
139         outputContainer->Add(fhEtaPi0) ;\r
140         \r
141         if(IsDataMC()){\r
142                 \r
143                 fhPtMCPi0  = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax); \r
144                 fhPtMCPi0->SetYTitle("N");\r
145                 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
146                 outputContainer->Add(fhPtMCPi0) ; \r
147                 \r
148                 fhPhiMCPi0  = new TH2F\r
149                 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
150                 fhPhiMCPi0->SetYTitle("#phi");\r
151                 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
152                 outputContainer->Add(fhPhiMCPi0) ; \r
153                 \r
154                 fhEtaMCPi0  = new TH2F\r
155                 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
156                 fhEtaMCPi0->SetYTitle("#eta");\r
157                 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
158                 outputContainer->Add(fhEtaMCPi0) ;\r
159                 \r
160                 fhPtMCNoPi0  = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax); \r
161                 fhPtMCNoPi0->SetYTitle("N");\r
162                 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
163                 outputContainer->Add(fhPtMCNoPi0) ; \r
164                 \r
165                 fhPhiMCNoPi0  = new TH2F\r
166                 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
167                 fhPhiMCNoPi0->SetYTitle("#phi");\r
168                 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
169                 outputContainer->Add(fhPhiMCNoPi0) ; \r
170                 \r
171                 fhEtaMCNoPi0  = new TH2F\r
172                 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
173                 fhEtaMCNoPi0->SetYTitle("#eta");\r
174                 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
175                 outputContainer->Add(fhEtaMCNoPi0) ;\r
176                 \r
177         }//Histos with MC\r
178     \r
179         \r
180         //Keep neutral meson selection histograms if requiered\r
181         //Setting done in AliNeutralMesonSelection\r
182         \r
183         if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){\r
184                 \r
185         TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
186         if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
187                         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
188         }\r
189         \r
190         //Save parameters used for analysis\r
191         TString parList ; //this will be list of parameters used for this analysis.\r
192         char onePar[255] ;\r
193         \r
194         sprintf(onePar,"--- AliAnaPi0EbE ---\n") ;\r
195         parList+=onePar ;       \r
196         sprintf(onePar,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;\r
197         parList+=onePar ;\r
198         \r
199         if(fAnaType == kSSCalo){\r
200                 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
201                 parList+=onePar ;\r
202                 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;\r
203                 parList+=onePar ;\r
204                 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;\r
205                 parList+=onePar ;\r
206                 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;\r
207                 parList+=onePar ;\r
208         }\r
209         \r
210         //Get parameters set in base class.\r
211         parList += GetBaseParametersList() ;\r
212         \r
213         //Get parameters set in PID class.\r
214         if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;\r
215         \r
216         TObjString *oString= new TObjString(parList) ;\r
217         outputContainer->Add(oString);\r
218         \r
219         return outputContainer ;\r
220         \r
221 }\r
222 \r
223 //__________________________________________________________________\r
224 void  AliAnaPi0EbE::MakeAnalysisFillAOD() \r
225 {\r
226         //Do analysis and fill aods\r
227         \r
228         switch(fAnaType) \r
229         {\r
230                 case kIMCalo:\r
231                         MakeInvMassInCalorimeter();\r
232                         break;\r
233                         \r
234                 case kSSCalo:\r
235                         MakeShowerShapeIdentification();\r
236                         break;\r
237                         \r
238                 case kIMCaloTracks:\r
239                         MakeInvMassInCalorimeterAndCTS();\r
240                         break;\r
241                         \r
242         }\r
243 }\r
244 \r
245 //__________________________________________________________________\r
246 void  AliAnaPi0EbE::MakeInvMassInCalorimeter() \r
247 {\r
248         //Do analysis and fill aods\r
249         //Search for the photon decay in calorimeters\r
250         //Read photon list from AOD, produced in class AliAnaPhoton\r
251         //Check if 2 photons have the mass of the pi0.\r
252         \r
253         TLorentzVector mom1;\r
254         TLorentzVector mom2;\r
255         TLorentzVector mom ;\r
256         Int_t tag1 =-1;\r
257         Int_t tag2 =-1;\r
258         Int_t tag = -1;\r
259         \r
260         if(!GetInputAODBranch())\r
261                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD with name branch < %s > \n",GetInputAODName().Data()));\r
262         \r
263         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
264                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
265                 mom1 = *(photon1->Momentum());\r
266                 \r
267                 //Play with the MC stack if available\r
268                 \r
269                 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){\r
270                         \r
271                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));\r
272                         mom2 = *(photon2->Momentum());\r
273                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
274                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))\r
275                         {\r
276                                 if(GetDebug()>1) printf("Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
277                                 \r
278                                 \r
279                                 if(IsDataMC()){\r
280                                         //Check origin of the candidates\r
281                                         tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
282                                         tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
283                                         if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
284                                         if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
285                                 }//Work with stack also   \r
286                                 \r
287                                 //Create AOD for analysis\r
288                                 mom = mom1+mom2;\r
289                                 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);\r
290                                 //pi0.SetLabel(calo->GetLabel(0));\r
291                                 pi0.SetPdg(AliCaloPID::kPi0);\r
292                                 pi0.SetDetector(photon1->GetDetector());\r
293                                 pi0.SetTag(tag);  \r
294                                 //Set the indeces of the original caloclusters  \r
295                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));\r
296                                 AddAODParticle(pi0);\r
297                         }//pi0\r
298                 }//2n photon loop\r
299                 \r
300         }//1st photon loop\r
301         \r
302         if(GetDebug() > 1) printf("End fill AODs \n");  \r
303         \r
304 }\r
305 \r
306 //__________________________________________________________________\r
307 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() \r
308 {\r
309         //Do analysis and fill aods\r
310         //Search for the photon decay in calorimeters\r
311         //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion\r
312         //Check if 2 photons have the mass of the pi0.\r
313         \r
314         TLorentzVector mom1;\r
315         TLorentzVector mom2;\r
316         TLorentzVector mom ;\r
317         Int_t tag1 =-1;\r
318         Int_t tag2 =-1;\r
319         Int_t tag = -1;\r
320         \r
321         if(!GetInputAODBranch())\r
322                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD branch with name < %s > \n",GetInputAODName().Data()));\r
323         \r
324         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
325                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
326                 mom1 = *(photon1->Momentum());\r
327                 \r
328                 //Play with the MC stack if available\r
329                 fInputAODGammaConv = (TClonesArray *) GetReader()->GetAOD()->FindListObject(fInputAODGammaConvName);\r
330                 if(!fInputAODGammaConv) \r
331                         AliFatal(Form("MakeInvMassInCaloAndCTS: No input gamma conversions in AOD branch with name < %s >",fInputAODGammaConvName.Data()));     \r
332                 \r
333                 for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){\r
334                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (fInputAODGammaConv->At(jphoton));\r
335                         mom2 = *(photon2->Momentum());\r
336                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
337                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){\r
338                                 if(GetDebug() > 1) printf("Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
339                                 \r
340                                 if(IsDataMC()){\r
341                                         //Check origin of the candidates\r
342                                         tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
343                                         tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
344                                         if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
345                                         if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
346                                 }//Work with stack also   \r
347                                 \r
348                                 //Create AOD for analysis\r
349                                 mom = mom1+mom2;\r
350                                 AliAODPWG4ParticleCorrelation pi0 = AliAODPWG4ParticleCorrelation(mom);\r
351                                 //pi0.SetLabel(calo->GetLabel(0));\r
352                                 pi0.SetPdg(AliCaloPID::kPi0);\r
353                                 pi0.SetDetector(photon1->GetDetector());\r
354                                 pi0.SetTag(tag);\r
355                                 //Set the indeces of the original tracks or caloclusters  \r
356                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);\r
357                                 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));\r
358                                 AddAODParticle(pi0);\r
359                         }//pi0\r
360                 }//2n photon loop\r
361                 \r
362         }//1st photon loop\r
363         \r
364         if(GetDebug() > 1) printf("End fill AODs \n");  \r
365         \r
366 }\r
367 \r
368 \r
369 //__________________________________________________________________\r
370 void  AliAnaPi0EbE::MakeShowerShapeIdentification() \r
371 {\r
372         //Search for pi0 in fCalorimeter with shower shape analysis \r
373         \r
374         TClonesArray * pl = new TClonesArray; \r
375         \r
376         //Get vertex for photon momentum calculation\r
377         Double_t vertex[]={0,0,0} ; //vertex ;\r
378         if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);\r
379         \r
380         //Select the Calorimeter of the photon\r
381         if(fCalorimeter == "PHOS")\r
382                 pl = GetAODPHOS();\r
383         else if (fCalorimeter == "EMCAL")\r
384                 pl = GetAODEMCAL();\r
385         //Fill AODCaloClusters and AODParticle with PHOS aods\r
386         TLorentzVector mom ;\r
387         for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){\r
388                 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));        \r
389                 \r
390                 //Cluster selection, not charged, with photon id and in fidutial cut\r
391                 //Get Momentum vector, \r
392                 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line\r
393                 //If too small or big pt, skip it\r
394                 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; \r
395                 //Check acceptance selection\r
396                 if(IsFidutialCutOn()){\r
397                         Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;\r
398                         if(! in ) continue ;\r
399                 }\r
400                 \r
401                 //Create AOD for analysis\r
402                 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);\r
403                 aodph.SetLabel(calo->GetLabel(0));\r
404                 aodph.SetDetector(fCalorimeter);\r
405                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());    \r
406                 \r
407                 //Check Distance to Bad channel, set bit.\r
408                 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel\r
409                 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;\r
410                 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)\r
411                         continue ;\r
412                 \r
413                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Bad channel cut passed %4.2f\n",distBad);\r
414                 \r
415                 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;\r
416                 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; \r
417                 else aodph.SetDistToBad(0) ;\r
418                 \r
419                 //Check PID\r
420                 //PID selection or bit setting\r
421                 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){\r
422                         //Get most probable PID, check PID weights (in MC this option is mandatory)\r
423                         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
424                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());\r
425                         //If primary is not photon, skip it.\r
426                         if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;\r
427                 }                                       \r
428                 else if(IsCaloPIDOn()){\r
429                         //Skip matched clusters with tracks\r
430                         if(calo->GetNTracksMatched() > 0) continue ;\r
431                         \r
432                         //Get most probable PID, 2 options check PID weights \r
433                         //or redo PID, recommended option for EMCal.            \r
434                         if(!IsCaloPIDRecalculationOn())\r
435                                 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
436                         else\r
437                                 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated\r
438                         \r
439                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());\r
440                         \r
441                         //If cluster does not pass pid, not photon, skip it.\r
442                         if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;                    \r
443                         \r
444                 }\r
445                 else{\r
446                         //Set PID bits for later selection (AliAnaPi0 for example)\r
447                         //GetPDG already called in SetPIDBits.\r
448                         GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);\r
449                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PID Bits set \n");            \r
450                 }\r
451                 \r
452                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());\r
453                 \r
454                 //Play with the MC stack if available\r
455                 //Check origin of the candidates\r
456                 if(IsDataMC()){\r
457                         aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));\r
458                         if(GetDebug() > 0) printf("AliAnaPi0EbE::FillAOD: Origin of candidate %d\n",aodph.GetTag());\r
459                 }//Work with stack also   \r
460                 \r
461                 //Add AOD with photon object to aod branch\r
462                 AddAODParticle(aodph);\r
463                 \r
464         }//loop\r
465         \r
466         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: End fill AODs \n");  \r
467         \r
468 }\r
469 //__________________________________________________________________\r
470 void  AliAnaPi0EbE::MakeAnalysisFillHistograms() \r
471 {\r
472         //Do analysis and fill histograms\r
473         \r
474         if(!GetOutputAODBranch())\r
475                 AliFatal(Form("MakeInvMassInCalo: No output pi0 in AOD branch with name < %s > \n",GetOutputAODName().Data()));\r
476         \r
477         //Loop on stored AOD pi0\r
478         Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
479         if(GetDebug() > 0) printf("pi0 aod branch entries %d\n", naod);\r
480         \r
481         for(Int_t iaod = 0; iaod < naod ; iaod++){\r
482                 \r
483                 AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
484                 Int_t pdg = pi0->GetPdg();\r
485                 \r
486                 if(pdg != AliCaloPID::kPi0) continue;              \r
487                 \r
488                 //Fill pi0 histograms \r
489                 Float_t pt = pi0->Pt();\r
490                 Float_t phi = pi0->Phi();\r
491                 Float_t eta = pi0->Eta();\r
492                 \r
493                 \r
494                 fhPtPi0  ->Fill(pt);\r
495                 fhPhiPi0 ->Fill(pt,phi);\r
496                 fhEtaPi0 ->Fill(pt,eta);\r
497                 \r
498                 if(IsDataMC()){\r
499                         if(pi0->GetTag()== AliCaloPID::kMCPi0Decay){\r
500                                 fhPtMCPi0  ->Fill(pt);\r
501                                 fhPhiMCPi0 ->Fill(pt,phi);\r
502                                 fhEtaMCPi0 ->Fill(pt,eta);\r
503                         }\r
504                         else{\r
505                                 fhPtMCNoPi0  ->Fill(pt);\r
506                                 fhPhiMCNoPi0 ->Fill(pt,phi);\r
507                                 fhEtaMCNoPi0 ->Fill(pt,eta);\r
508                         }\r
509                 }//Histograms with MC\r
510                 \r
511         }// aod loop\r
512         \r
513 }\r
514 \r
515 \r
516 //____________________________________________________________________________\r
517 void AliAnaPi0EbE::InitParameters()\r
518 {\r
519         \r
520         //Initialize the parameters of the analysis.\r
521         SetOutputAODClassName("AliAODPWG4Particle");\r
522         SetOutputAODName("pi0s");\r
523         fInputAODGammaConvName = "gammaconv" ;\r
524         fAnaType = kIMCalo ;\r
525         fCalorimeter = "PHOS" ;\r
526         fMinDist  = 2.;\r
527         fMinDist2 = 4.;\r
528         fMinDist3 = 5.;\r
529         \r
530 }\r
531 \r
532 //__________________________________________________________________\r
533 void AliAnaPi0EbE::Print(const Option_t * opt) const\r
534 {\r
535         \r
536         //Print some relevant parameters set for the analysis\r
537         if(! opt)\r
538                 return;\r
539         \r
540         printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
541         AliAnaPartCorrBaseClass::Print("");\r
542         printf("Analysis Type = %d \n",  fAnaType) ;\r
543         if(fAnaType == kSSCalo){     \r
544                 printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;\r
545                 printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);\r
546                 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);\r
547                 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); \r
548         } \r
549         printf("    \n") ;\r
550         \r
551\r