1)AliAODParticleCorrelation.h: Data member fIsolation changed from float to bool
[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 "AliMCAnalysisUtils.h"\r
40 #include "AliAODPWG4ParticleCorrelation.h"\r
41     #include "AliStack.h"\r
42 ClassImp(AliAnaPi0EbE)\r
43   \r
44 //____________________________________________________________________________\r
45 AliAnaPi0EbE::AliAnaPi0EbE() : \r
46 AliAnaPartCorrBaseClass(),  fAnaType(kIMCalo),fCalorimeter(""),\r
47 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),       \r
48 fInputAODGammaConv(0x0),fInputAODGammaConvName(""),\r
49 fhPtPi0(0),fhPhiPi0(0),fhEtaPi0(0), \r
50 fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0), \r
51 fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)\r
52 {\r
53         //default ctor\r
54         \r
55         //Initialize parameters\r
56         InitParameters();\r
57         \r
58 }\r
59 \r
60 //____________________________________________________________________________\r
61 AliAnaPi0EbE::AliAnaPi0EbE(const AliAnaPi0EbE & p) : \r
62 AliAnaPartCorrBaseClass(p),  fAnaType(p.fAnaType), fCalorimeter(p.fCalorimeter),\r
63 fMinDist(p.fMinDist),fMinDist2(p.fMinDist2), fMinDist3(p.fMinDist3),\r
64 fInputAODGammaConv(new TClonesArray(*p.fInputAODGammaConv)),\r
65 fInputAODGammaConvName(p.fInputAODGammaConvName),\r
66 fhPtPi0(p.fhPtPi0),fhPhiPi0(p.fhPhiPi0),fhEtaPi0(p.fhEtaPi0), \r
67 fhPtMCNoPi0(p.fhPtMCNoPi0),fhPhiMCNoPi0(p.fhPhiMCNoPi0),fhEtaMCNoPi0(p.fhEtaMCNoPi0), \r
68 fhPtMCPi0(p.fhPtMCPi0),fhPhiMCPi0(p.fhPhiMCPi0),fhEtaMCPi0(p.fhEtaMCPi0) \r
69 {\r
70         // cpy ctor\r
71 }\r
72 \r
73 //_________________________________________________________________________\r
74 AliAnaPi0EbE & AliAnaPi0EbE::operator = (const AliAnaPi0EbE & p)\r
75 {\r
76         // assignment operator\r
77         \r
78         if(&p == this) return *this;\r
79         \r
80         fAnaType     = p.fAnaType ;\r
81         fCalorimeter = p.fCalorimeter ;\r
82         fMinDist     = p.fMinDist;\r
83         fMinDist2    = p.fMinDist2;\r
84         fMinDist3    = p.fMinDist3;\r
85         \r
86         fInputAODGammaConv     = new TClonesArray(*p.fInputAODGammaConv) ;\r
87         fInputAODGammaConvName = p.fInputAODGammaConvName ;\r
88         \r
89         fhPtPi0 = p.fhPtPi0; fhPhiPi0 = p.fhPhiPi0; fhEtaPi0 = p.fhEtaPi0;  \r
90         fhPtMCNoPi0 = p.fhPtMCNoPi0; fhPhiMCNoPi0 = p.fhPhiMCNoPi0; fhEtaMCNoPi0 = p.fhEtaMCPi0;  \r
91         fhPtMCPi0 = p.fhPtMCPi0; fhPhiMCPi0 = p.fhPhiMCPi0; fhEtaMCPi0 = p.fhEtaMCPi0;  \r
92         \r
93         return *this;\r
94         \r
95 }\r
96 \r
97 //____________________________________________________________________________\r
98 AliAnaPi0EbE::~AliAnaPi0EbE() \r
99 {\r
100         //dtor\r
101         if(fInputAODGammaConv){\r
102                 fInputAODGammaConv->Clear() ; \r
103                 delete fInputAODGammaConv ;\r
104         }\r
105 }\r
106 \r
107 //________________________________________________________________________\r
108 TList *  AliAnaPi0EbE::GetCreateOutputObjects()\r
109 {  \r
110         // Create histograms to be saved in output file and \r
111         // store them in outputContainer\r
112         TList * outputContainer = new TList() ; \r
113         outputContainer->SetName("Pi0EbEHistos") ; \r
114         \r
115         Int_t nptbins  = GetHistoNPtBins();\r
116         Int_t nphibins = GetHistoNPhiBins();\r
117         Int_t netabins = GetHistoNEtaBins();\r
118         Float_t ptmax  = GetHistoPtMax();\r
119         Float_t phimax = GetHistoPhiMax();\r
120         Float_t etamax = GetHistoEtaMax();\r
121         Float_t ptmin  = GetHistoPtMin();\r
122         Float_t phimin = GetHistoPhiMin();\r
123         Float_t etamin = GetHistoEtaMin();      \r
124         \r
125         fhPtPi0  = new TH1F("hPtPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax); \r
126         fhPtPi0->SetYTitle("N");\r
127         fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
128         outputContainer->Add(fhPtPi0) ; \r
129         \r
130         fhPhiPi0  = new TH2F\r
131     ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
132         fhPhiPi0->SetYTitle("#phi");\r
133         fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
134         outputContainer->Add(fhPhiPi0) ; \r
135         \r
136         fhEtaPi0  = new TH2F\r
137     ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
138         fhEtaPi0->SetYTitle("#eta");\r
139         fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
140         outputContainer->Add(fhEtaPi0) ;\r
141         \r
142         if(IsDataMC()) {
143           if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
144              GetReader()->GetDataType() != AliCaloTrackReader::kMC){\r
145                 \r
146                 fhPtMCPi0  = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax); \r
147                 fhPtMCPi0->SetYTitle("N");\r
148                 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
149                 outputContainer->Add(fhPtMCPi0) ; \r
150                 \r
151                 fhPhiMCPi0  = new TH2F\r
152                 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
153                 fhPhiMCPi0->SetYTitle("#phi");\r
154                 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
155                 outputContainer->Add(fhPhiMCPi0) ; \r
156                 \r
157                 fhEtaMCPi0  = new TH2F\r
158                 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
159                 fhEtaMCPi0->SetYTitle("#eta");\r
160                 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
161                 outputContainer->Add(fhEtaMCPi0) ;\r
162                 \r
163                 fhPtMCNoPi0  = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax); \r
164                 fhPtMCNoPi0->SetYTitle("N");\r
165                 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
166                 outputContainer->Add(fhPtMCNoPi0) ; \r
167                 \r
168                 fhPhiMCNoPi0  = new TH2F\r
169                 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
170                 fhPhiMCNoPi0->SetYTitle("#phi");\r
171                 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
172                 outputContainer->Add(fhPhiMCNoPi0) ; \r
173                 \r
174                 fhEtaMCNoPi0  = new TH2F\r
175                 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
176                 fhEtaMCNoPi0->SetYTitle("#eta");\r
177                 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
178                 outputContainer->Add(fhEtaMCNoPi0) ;\r
179                 \r
180         }
181         }//Histos with MC\r
182     \r
183         \r
184         //Keep neutral meson selection histograms if requiered\r
185         //Setting done in AliNeutralMesonSelection\r
186         \r
187         if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){\r
188                 \r
189         TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
190         if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
191                         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
192         }\r
193         \r
194         //Save parameters used for analysis\r
195         TString parList ; //this will be list of parameters used for this analysis.\r
196         char onePar[255] ;\r
197         \r
198         sprintf(onePar,"--- AliAnaPi0EbE ---\n") ;\r
199         parList+=onePar ;       \r
200         sprintf(onePar,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;\r
201         parList+=onePar ;\r
202         \r
203         if(fAnaType == kSSCalo){\r
204                 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
205                 parList+=onePar ;\r
206                 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;\r
207                 parList+=onePar ;\r
208                 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;\r
209                 parList+=onePar ;\r
210                 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;\r
211                 parList+=onePar ;\r
212         }\r
213         \r
214         //Get parameters set in base class.\r
215         parList += GetBaseParametersList() ;\r
216         \r
217         //Get parameters set in PID class.\r
218         if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;\r
219         \r
220         TObjString *oString= new TObjString(parList) ;\r
221         outputContainer->Add(oString);\r
222         \r
223         return outputContainer ;\r
224         \r
225 }\r
226 \r
227 //__________________________________________________________________\r
228 void  AliAnaPi0EbE::MakeAnalysisFillAOD() \r
229 {\r
230         //Do analysis and fill aods\r
231         \r
232         switch(fAnaType) \r
233         {\r
234                 case kIMCalo:\r
235                         MakeInvMassInCalorimeter();\r
236                         break;\r
237                         \r
238                 case kSSCalo:\r
239                         MakeShowerShapeIdentification();\r
240                         break;\r
241                         \r
242                 case kIMCaloTracks:\r
243                         MakeInvMassInCalorimeterAndCTS();\r
244                         break;\r
245                         \r
246         }\r
247 }\r
248 \r
249 //__________________________________________________________________\r
250 void  AliAnaPi0EbE::MakeInvMassInCalorimeter() \r
251 {\r
252         //Do analysis and fill aods\r
253         //Search for the photon decay in calorimeters\r
254         //Read photon list from AOD, produced in class AliAnaPhoton\r
255         //Check if 2 photons have the mass of the pi0.\r
256         \r
257         TLorentzVector mom1;\r
258         TLorentzVector mom2;\r
259         TLorentzVector mom ;\r
260         Int_t tag1 =-1;\r
261         Int_t tag2 =-1;\r
262         Int_t tag = -1;\r
263         \r
264         if(!GetInputAODBranch())\r
265                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD with name branch < %s > \n",GetInputAODName().Data()));\r
266         \r
267         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
268                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
269                 mom1 = *(photon1->Momentum());\r
270                 \r
271         \r
272                 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){\r
273                         \r
274                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));\r
275                         mom2 = *(photon2->Momentum());\r
276                 \r
277                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
278                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))\r
279                         {\r
280                                 if(GetDebug()>1) \r
281                                         printf("AliAnaPi0EbE::kIMCalo::Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
282                                 \r
283                                                 //Play with the MC stack if available\r
284                                 if(IsDataMC()){\r
285                                         //Check origin of the candidates\r
286                                         tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
287                                         tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
288                                         \r
289                                         if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCalo::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
290                                         if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){\r
291                                         \r
292                                          //Check if pi0 mother is the same\r
293                                         Int_t label1 = photon1->GetLabel();\r
294                                         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
295                                         label1 = mother1->GetFirstMother();\r
296                                         //mother1 = GetMCStack()->Particle(label1);//pi0\r
297                                         \r
298                                         Int_t label2 = photon2->GetLabel();\r
299                                         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
300                                         label2 = mother2->GetFirstMother();\r
301                                         //mother2 = GetMCStack()->Particle(label2);//pi0\r
302                                         \r
303                                         //printf("mother1 %d, mother2 %d\n",label1,label2);\r
304                                         if(label1 == label2)\r
305                                                 tag = AliMCAnalysisUtils::kMCPi0;\r
306                                         }\r
307                                 }//Work with stack also   \r
308                                 \r
309                                 //Create AOD for analysis\r
310                                 mom = mom1+mom2;\r
311                                 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);\r
312                                 //pi0.SetLabel(calo->GetLabel(0));\r
313                                 pi0.SetPdg(AliCaloPID::kPi0);\r
314                                 pi0.SetDetector(photon1->GetDetector());\r
315                                 pi0.SetTag(tag);  \r
316                                 //Set the indeces of the original caloclusters  \r
317                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));\r
318                                 AddAODParticle(pi0);\r
319                         }//pi0\r
320                 }//2n photon loop\r
321                 \r
322         }//1st photon loop\r
323         \r
324         if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCalo::End fill AODs \n");  \r
325         \r
326 }\r
327 \r
328 //__________________________________________________________________\r
329 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() \r
330 {\r
331         //Do analysis and fill aods\r
332         //Search for the photon decay in calorimeters\r
333         //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion\r
334         //Check if 2 photons have the mass of the pi0.\r
335         \r
336         TLorentzVector mom1;\r
337         TLorentzVector mom2;\r
338         TLorentzVector mom ;\r
339         Int_t tag1 =-1;\r
340         Int_t tag2 =-1;\r
341         Int_t tag = -1;\r
342         \r
343         if(!GetInputAODBranch())\r
344                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD branch with name < %s > \n",GetInputAODName().Data()));\r
345         \r
346         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
347                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
348                 mom1 = *(photon1->Momentum());\r
349                 \r
350                 //Play with the MC stack if available\r
351                 fInputAODGammaConv = (TClonesArray *) GetReader()->GetAOD()->FindListObject(fInputAODGammaConvName);\r
352                 if(!fInputAODGammaConv) \r
353                         AliFatal(Form("MakeInvMassInCaloAndCTS: No input gamma conversions in AOD branch with name < %s >",fInputAODGammaConvName.Data()));     \r
354                 \r
355                 for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){\r
356                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (fInputAODGammaConv->At(jphoton));\r
357                         mom2 = *(photon2->Momentum());\r
358                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
359                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){\r
360                                 if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCaloCTS::Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
361                                 \r
362                                 if(IsDataMC()){\r
363                                         //Check origin of the candidates\r
364                                         tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
365                                         tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
366                                         if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCaloCTS::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
367                                         if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){\r
368                                                 //Check if pi0 mother is the same\r
369                                                 Int_t label1 = photon1->GetLabel();\r
370                                                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
371                                                 label1 = mother1->GetFirstMother();\r
372                                                 //mother1 = GetMCStack()->Particle(label1);//pi0\r
373                                         \r
374                                                 Int_t label2 = photon2->GetLabel();\r
375                                                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
376                                                 label2 = mother2->GetFirstMother();\r
377                                                 //mother2 = GetMCStack()->Particle(label2);//pi0\r
378                                         \r
379                                                 //printf("mother1 %d, mother2 %d\n",label1,label2);\r
380                                                 if(label1 == label2)\r
381                                                         tag = AliMCAnalysisUtils::kMCPi0;\r
382                                         }\r
383                                 }//Work with stack also   \r
384                                 \r
385                                 //Create AOD for analysis\r
386                                 mom = mom1+mom2;\r
387                                 AliAODPWG4ParticleCorrelation pi0 = AliAODPWG4ParticleCorrelation(mom);\r
388                                 //pi0.SetLabel(calo->GetLabel(0));\r
389                                 pi0.SetPdg(AliCaloPID::kPi0);\r
390                                 pi0.SetDetector(photon1->GetDetector());\r
391                                 pi0.SetTag(tag);\r
392                                 //Set the indeces of the original tracks or caloclusters  \r
393                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);\r
394                                 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));\r
395                                 AddAODParticle(pi0);\r
396                         }//pi0\r
397                 }//2n photon loop\r
398                 \r
399         }//1st photon loop\r
400         \r
401         if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCaloCTS::End fill AODs \n");  \r
402         \r
403 }\r
404 \r
405 \r
406 //__________________________________________________________________\r
407 void  AliAnaPi0EbE::MakeShowerShapeIdentification() \r
408 {\r
409         //Search for pi0 in fCalorimeter with shower shape analysis \r
410         \r
411         TClonesArray * pl = new TClonesArray; \r
412         \r
413         //Get vertex for photon momentum calculation\r
414         Double_t vertex[]={0,0,0} ; //vertex ;\r
415         if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);\r
416         \r
417         //Select the Calorimeter of the photon\r
418         if(fCalorimeter == "PHOS")\r
419                 pl = GetAODPHOS();\r
420         else if (fCalorimeter == "EMCAL")\r
421                 pl = GetAODEMCAL();\r
422         //Fill AODCaloClusters and AODParticle with PHOS aods\r
423         TLorentzVector mom ;\r
424         for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){\r
425                 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));        \r
426                 \r
427                 //Cluster selection, not charged, with pi0 id and in fidutial cut\r
428                 //Get Momentum vector, \r
429                 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line\r
430                 //If too small or big pt, skip it\r
431                 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; \r
432                 //Check acceptance selection\r
433                 if(IsFidutialCutOn()){\r
434                         Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;\r
435                         if(! in ) continue ;\r
436                 }\r
437                 \r
438                 //Create AOD for analysis\r
439                 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);\r
440                 aodpi0.SetLabel(calo->GetLabel(0));\r
441                 //Set the indeces of the original caloclusters  
442                 aodpi0.SetCaloLabel(calo->GetID(),-1);
443                 aodpi0.SetDetector(fCalorimeter);\r
444                 if(GetDebug() > 1) 
445                   printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());       \r
446                 \r
447                 //Check Distance to Bad channel, set bit.\r
448                 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel\r
449                 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;\r
450                 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)\r
451                         continue ;\r
452                 \r
453                 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Bad channel cut passed %4.2f\n",distBad);\r
454                 \r
455                 if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;\r
456                 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ; \r
457                 else aodpi0.SetDistToBad(0) ;\r
458                 \r
459                 //Check PID\r
460                 //PID selection or bit setting\r
461                 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){\r
462                         //Get most probable PID, check PID weights (in MC this option is mandatory)\r
463                         aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
464                         if(GetDebug() > 1) 
465                           printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());\r
466                         //If primary is not pi0, skip it.\r
467                         if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;\r
468                 }                                       \r
469                 else if(IsCaloPIDOn()){\r
470                         //Skip matched clusters with tracks\r
471                         if(calo->GetNTracksMatched() > 0) continue ;\r
472                         \r
473                         //Get most probable PID, 2 options check PID weights \r
474                         //or redo PID, recommended option for EMCal.            \r
475                         if(!IsCaloPIDRecalculationOn())\r
476                                 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
477                         else\r
478                                 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated\r
479                         \r
480                         if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());\r
481                         \r
482                         //If cluster does not pass pid, not pi0, skip it.\r
483                         if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;                      \r
484                         \r
485                 }\r
486                 else{\r
487                         //Set PID bits for later selection (AliAnaPi0 for example)\r
488                         //GetPDG already called in SetPIDBits.\r
489                         GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0);\r
490                         if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PID Bits set \n");         \r
491                 }\r
492                 \r
493                 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());\r
494                 \r
495                 //Play with the MC stack if available\r
496                 //Check origin of the candidates\r
497                 if(IsDataMC()){\r
498                   if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
499                      GetReader()->GetDataType() != AliCaloTrackReader::kMC){\r
500                         aodpi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));\r
501                         if(GetDebug() > 0) printf("AliAnaPi0EbE::kSSCalo::EbE::FillAOD: Origin of candidate %d\n",aodpi0.GetTag());\r
502                   }
503                 }//Work with stack also   \r
504                 \r
505                 //Add AOD with pi0 object to aod branch\r
506                 AddAODParticle(aodpi0);\r
507                 \r
508         }//loop\r
509         \r
510         if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: End fill AODs \n");  \r
511         \r
512 }\r
513 //__________________________________________________________________\r
514 void  AliAnaPi0EbE::MakeAnalysisFillHistograms() \r
515 {\r
516         //Do analysis and fill histograms\r
517         \r
518         if(!GetOutputAODBranch())\r
519                 AliFatal(Form("MakeInvMassInCalo: No output pi0 in AOD branch with name < %s > \n",GetOutputAODName().Data()));\r
520         \r
521         //Loop on stored AOD pi0\r
522         Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
523         if(GetDebug() > 0) printf("AliAnaPi0EbE::Histo::pi0 aod branch entries %d\n", naod);\r
524         \r
525         for(Int_t iaod = 0; iaod < naod ; iaod++){\r
526                 \r
527                 AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
528                 Int_t pdg = pi0->GetPdg();\r
529                 \r
530                 if(pdg != AliCaloPID::kPi0) continue;              \r
531                 \r
532                 //Fill pi0 histograms \r
533                 Float_t pt = pi0->Pt();\r
534                 Float_t phi = pi0->Phi();\r
535                 Float_t eta = pi0->Eta();\r
536                 \r
537                 \r
538                 fhPtPi0  ->Fill(pt);\r
539                 fhPhiPi0 ->Fill(pt,phi);\r
540                 fhEtaPi0 ->Fill(pt,eta);\r
541                 \r
542                 if(IsDataMC()){\r
543                   if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
544                      GetReader()->GetDataType() != AliCaloTrackReader::kMC){\r
545                         if(pi0->GetTag()== AliMCAnalysisUtils::kMCPi0){\r
546                                 fhPtMCPi0  ->Fill(pt);\r
547                                 fhPhiMCPi0 ->Fill(pt,phi);\r
548                                 fhEtaMCPi0 ->Fill(pt,eta);\r
549                         }\r
550                         else{\r
551                                 fhPtMCNoPi0  ->Fill(pt);\r
552                                 fhPhiMCNoPi0 ->Fill(pt,phi);\r
553                                 fhEtaMCNoPi0 ->Fill(pt,eta);\r
554                         }\r
555                   }
556                 }//Histograms with MC\r
557                 \r
558         }// aod loop\r
559         \r
560 }\r
561 \r
562 \r
563 //____________________________________________________________________________
564 void AliAnaPi0EbE::Init()
565 {
566   
567   //Init
568   //Do some checks
569   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
570     printf("!!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
571     abort();
572   }
573   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
574     printf("!!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
575     abort();
576   }
577   
578 }
579
580 //____________________________________________________________________________\r
581 void AliAnaPi0EbE::InitParameters()\r
582 {\r
583         \r
584         //Initialize the parameters of the analysis.\r
585         SetOutputAODClassName("AliAODPWG4Particle");\r
586         SetOutputAODName("pi0s");\r
587         fInputAODGammaConvName = "gammaconv" ;\r
588         fAnaType = kIMCalo ;\r
589         fCalorimeter = "EMCAL" ;\r
590         fMinDist  = 2.;\r
591         fMinDist2 = 4.;\r
592         fMinDist3 = 5.;\r
593         \r
594 }\r
595 \r
596 //__________________________________________________________________\r
597 void AliAnaPi0EbE::Print(const Option_t * opt) const\r
598 {\r
599         \r
600         //Print some relevant parameters set for the analysis\r
601         if(! opt)\r
602                 return;\r
603         \r
604         printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
605         AliAnaPartCorrBaseClass::Print("");\r
606         printf("Analysis Type = %d \n",  fAnaType) ;\r
607         if(fAnaType == kSSCalo){     \r
608                 printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;\r
609                 printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);\r
610                 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);\r
611                 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); \r
612         } \r
613         printf("    \n") ;\r
614         \r
615\r