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