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