rearrange the x-y input of histograms for CORRFW requierements
[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 #include <TClonesArray.h>
30 #include <TObjString.h>
31 //#include <Riostream.h>
32 #include "TParticle.h"
33
34 // --- Analysis system --- 
35 #include "AliAnaPhoton.h" 
36 #include "AliCaloTrackReader.h"
37 #include "AliStack.h"
38 #include "AliCaloPID.h"
39 #include "AliMCAnalysisUtils.h"
40 #include "AliFidutialCut.h"
41 #include "AliAODCaloCluster.h"
42 #include "AliAODMCParticle.h"
43
44 ClassImp(AliAnaPhoton)
45   
46 //____________________________________________________________________________
47   AliAnaPhoton::AliAnaPhoton() : 
48     AliAnaPartCorrBaseClass(), fCalorimeter(""), 
49     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
50         fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
51     //MC
52     fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
53         fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0), 
54     fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0), 
55     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0), 
56     fhPtISR(0),fhPhiISR(0),fhEtaISR(0), 
57     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0), 
58     fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0), 
59     fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0), 
60     fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
61 {
62   //default ctor
63   
64   //Initialize parameters
65   InitParameters();
66
67 }
68
69 //____________________________________________________________________________
70 AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) : 
71   AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
72    fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
73    fRejectTrackMatch(g.fRejectTrackMatch),
74    fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton), 
75   //MC
76   fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt), 
77   fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
78   fh2E(g.fh2E), fh2Pt(g.fh2Pt), 
79   fhPtMCPhoton(g.fhPtMCPhoton),fhPhiMCPhoton(g.fhPhiMCPhoton),fhEtaMCPhoton(g.fhEtaMCPhoton),
80   fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt), 
81   fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation), 
82   fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR), 
83   fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay), 
84   fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay), 
85   fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion), 
86   fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
87   
88 {
89   // cpy ctor
90   
91 }
92
93 //_________________________________________________________________________
94 AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
95 {
96   // assignment operator
97   
98   if(&g == this) return *this;
99
100   fCalorimeter = g.fCalorimeter ;
101   fMinDist  = g.fMinDist;
102   fMinDist2 = g.fMinDist2;
103   fMinDist3 = g.fMinDist3;
104   fRejectTrackMatch = g.fRejectTrackMatch;
105       
106   fhPtPhoton  = g.fhPtPhoton ; 
107   fhPhiPhoton = g.fhPhiPhoton ;
108   fhEtaPhoton = g.fhEtaPhoton ;
109
110   fhDeltaE   = g.fhDeltaE;  
111   fhDeltaPt  = g.fhDeltaPt;
112   fhRatioE   = g.fhRatioE;  
113   fhRatioPt  = g.fhRatioPt;
114   fh2E   = g.fh2E;  
115   fh2Pt  = g.fh2Pt;
116  
117   fhPtMCPhoton  = g.fhPtMCPhoton ; 
118   fhPhiMCPhoton = g.fhPhiMCPhoton ;
119   fhEtaMCPhoton = g.fhEtaMCPhoton ;     
120   fhPtPrompt  = g.fhPtPrompt;
121   fhPhiPrompt = g.fhPhiPrompt;
122   fhEtaPrompt = g.fhEtaPrompt; 
123   fhPtFragmentation  = g.fhPtFragmentation;
124   fhPhiFragmentation = g.fhPhiFragmentation;
125   fhEtaFragmentation = g.fhEtaFragmentation; 
126   fhPtISR  = g.fhPtISR;
127   fhPhiISR = g.fhPhiISR; 
128   fhEtaISR = g.fhEtaISR; 
129   fhPtPi0Decay  = g.fhPtPi0Decay;
130   fhPhiPi0Decay = g.fhPhiPi0Decay;
131   fhEtaPi0Decay = g.fhEtaPi0Decay; 
132   fhPtOtherDecay  = g.fhPtOtherDecay;
133   fhPhiOtherDecay = g.fhPhiOtherDecay;
134   fhEtaOtherDecay = g.fhEtaOtherDecay; 
135   fhPtConversion  = g. fhPtConversion;
136   fhPhiConversion = g.fhPhiConversion;
137   fhEtaConversion = g.fhEtaConversion; 
138   fhPtUnknown  = g.fhPtUnknown;
139   fhPhiUnknown = g.fhPhiUnknown;
140   fhEtaUnknown = g.fhEtaUnknown; 
141
142   return *this;
143   
144 }
145
146 //____________________________________________________________________________
147 AliAnaPhoton::~AliAnaPhoton() 
148 {
149   //dtor
150
151 }
152
153
154 //________________________________________________________________________
155 TList *  AliAnaPhoton::GetCreateOutputObjects()
156 {  
157   // Create histograms to be saved in output file and 
158   // store them in outputContainer
159   TList * outputContainer = new TList() ; 
160   outputContainer->SetName("PhotonHistos") ; 
161   
162   Int_t nptbins  = GetHistoNPtBins();
163   Int_t nphibins = GetHistoNPhiBins();
164   Int_t netabins = GetHistoNEtaBins();
165   Float_t ptmax  = GetHistoPtMax();
166   Float_t phimax = GetHistoPhiMax();
167   Float_t etamax = GetHistoEtaMax();
168   Float_t ptmin  = GetHistoPtMin();
169   Float_t phimin = GetHistoPhiMin();
170   Float_t etamin = GetHistoEtaMin();    
171   
172   //Histograms of highest Photon identified in Event
173   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
174   fhPtPhoton->SetYTitle("N");
175   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
176   outputContainer->Add(fhPtPhoton) ; 
177   
178   fhPhiPhoton  = new TH2F
179     ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
180   fhPhiPhoton->SetYTitle("#phi");
181   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
182   outputContainer->Add(fhPhiPhoton) ; 
183   
184   fhEtaPhoton  = new TH2F
185     ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
186   fhEtaPhoton->SetYTitle("#eta");
187   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
188   outputContainer->Add(fhEtaPhoton) ;
189   
190   if(IsDataMC()){
191     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50); 
192     fhDeltaE->SetXTitle("#Delta E (GeV)");
193     outputContainer->Add(fhDeltaE);
194                 
195     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50); 
196     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
197     outputContainer->Add(fhDeltaPt);
198
199     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 200,0,2); 
200     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
201     outputContainer->Add(fhRatioE);
202     
203     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2); 
204     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
205     outputContainer->Add(fhRatioPt);    
206
207     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
208     fh2E->SetXTitle("E_{rec} (GeV)");
209     fh2E->SetYTitle("E_{gen} (GeV)");
210     outputContainer->Add(fh2E);          
211     
212     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
213     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
214     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
215     outputContainer->Add(fh2Pt);
216    
217         fhPtMCPhoton  = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
218         fhPtMCPhoton->SetYTitle("N");
219         fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
220         outputContainer->Add(fhPtMCPhoton) ; 
221           
222         fhPhiMCPhoton  = new TH2F
223         ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
224         fhPhiMCPhoton->SetYTitle("#phi");
225         fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
226         outputContainer->Add(fhPhiMCPhoton) ; 
227           
228         fhEtaMCPhoton  = new TH2F
229         ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
230         fhEtaMCPhoton->SetYTitle("#eta");
231         fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
232         outputContainer->Add(fhEtaMCPhoton) ;
233           
234     fhPtPrompt  = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax); 
235     fhPtPrompt->SetYTitle("N");
236     fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
237     outputContainer->Add(fhPtPrompt) ; 
238     
239     fhPhiPrompt  = new TH2F
240       ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
241     fhPhiPrompt->SetYTitle("#phi");
242     fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
243     outputContainer->Add(fhPhiPrompt) ; 
244     
245     fhEtaPrompt  = new TH2F
246       ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
247     fhEtaPrompt->SetYTitle("#eta");
248     fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
249     outputContainer->Add(fhEtaPrompt) ;
250     
251     fhPtFragmentation  = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax); 
252     fhPtFragmentation->SetYTitle("N");
253     fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
254     outputContainer->Add(fhPtFragmentation) ; 
255     
256     fhPhiFragmentation  = new TH2F
257       ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
258     fhPhiFragmentation->SetYTitle("#phi");
259     fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
260     outputContainer->Add(fhPhiFragmentation) ; 
261     
262     fhEtaFragmentation  = new TH2F
263       ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
264     fhEtaFragmentation->SetYTitle("#eta");
265     fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
266     outputContainer->Add(fhEtaFragmentation) ;
267     
268     fhPtISR  = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax); 
269     fhPtISR->SetYTitle("N");
270     fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
271     outputContainer->Add(fhPtISR) ; 
272     
273     fhPhiISR  = new TH2F
274       ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
275     fhPhiISR->SetYTitle("#phi");
276     fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
277     outputContainer->Add(fhPhiISR) ; 
278     
279     fhEtaISR  = new TH2F
280       ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
281     fhEtaISR->SetYTitle("#eta");
282     fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
283     outputContainer->Add(fhEtaISR) ;
284     
285     fhPtPi0Decay  = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
286     fhPtPi0Decay->SetYTitle("N");
287     fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
288     outputContainer->Add(fhPtPi0Decay) ; 
289     
290     fhPhiPi0Decay  = new TH2F
291       ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
292     fhPhiPi0Decay->SetYTitle("#phi");
293     fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
294     outputContainer->Add(fhPhiPi0Decay) ; 
295     
296     fhEtaPi0Decay  = new TH2F
297       ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
298     fhEtaPi0Decay->SetYTitle("#eta");
299     fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
300     outputContainer->Add(fhEtaPi0Decay) ;
301     
302     fhPtOtherDecay  = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
303     fhPtOtherDecay->SetYTitle("N");
304     fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
305     outputContainer->Add(fhPtOtherDecay) ; 
306     
307     fhPhiOtherDecay  = new TH2F
308       ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
309     fhPhiOtherDecay->SetYTitle("#phi");
310     fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
311     outputContainer->Add(fhPhiOtherDecay) ; 
312     
313     fhEtaOtherDecay  = new TH2F
314       ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
315     fhEtaOtherDecay->SetYTitle("#eta");
316     fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
317     outputContainer->Add(fhEtaOtherDecay) ;
318     
319     fhPtConversion  = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
320     fhPtConversion->SetYTitle("N");
321     fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
322     outputContainer->Add(fhPtConversion) ; 
323     
324     fhPhiConversion  = new TH2F
325       ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
326     fhPhiConversion->SetYTitle("#phi");
327     fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
328     outputContainer->Add(fhPhiConversion) ; 
329     
330     fhEtaConversion  = new TH2F
331       ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
332     fhEtaConversion->SetYTitle("#eta");
333     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
334     outputContainer->Add(fhEtaConversion) ;
335     
336     fhPtUnknown  = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
337     fhPtUnknown->SetYTitle("N");
338     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
339     outputContainer->Add(fhPtUnknown) ; 
340     
341     fhPhiUnknown  = new TH2F
342       ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
343     fhPhiUnknown->SetYTitle("#phi");
344     fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
345     outputContainer->Add(fhPhiUnknown) ; 
346     
347     fhEtaUnknown  = new TH2F
348       ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
349     fhEtaUnknown->SetYTitle("#eta");
350     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
351     outputContainer->Add(fhEtaUnknown) ;
352         
353   }//Histos with MC
354   
355   //Save parameters used for analysis
356   TString parList ; //this will be list of parameters used for this analysis.
357   char onePar[255] ;
358   
359   sprintf(onePar,"--- AliAnaPhoton ---\n") ;
360   parList+=onePar ;     
361   sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
362   parList+=onePar ;
363   sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
364   parList+=onePar ;
365   sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
366   parList+=onePar ;
367   sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
368   parList+=onePar ;
369   sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
370   parList+=onePar ;  
371   
372   //Get parameters set in base class.
373   parList += GetBaseParametersList() ;
374   
375   //Get parameters set in PID class.
376   parList += GetCaloPID()->GetPIDParametersList() ;
377   
378   //Get parameters set in FidutialCut class (not available yet)
379   //parlist += GetFidCut()->GetFidCutParametersList() 
380   
381   TObjString *oString= new TObjString(parList) ;
382   outputContainer->Add(oString);
383   
384   return outputContainer ;
385   
386 }
387
388 //____________________________________________________________________________
389 void AliAnaPhoton::Init()
390 {
391   
392   //Init
393   //Do some checks
394   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
395     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
396     abort();
397   }
398   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
399     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
400     abort();
401   }
402   
403 }
404
405
406 //____________________________________________________________________________
407 void AliAnaPhoton::InitParameters()
408 {
409   
410   //Initialize the parameters of the analysis.
411   SetOutputAODClassName("AliAODPWG4Particle");
412   SetOutputAODName("PWG4Particle");
413
414   AddToHistogramsName("AnaPhoton_");
415
416   fCalorimeter = "PHOS" ;
417   fMinDist  = 2.;
418   fMinDist2 = 4.;
419   fMinDist3 = 5.;
420   fRejectTrackMatch = kTRUE ;
421 }
422
423 //__________________________________________________________________
424 void  AliAnaPhoton::MakeAnalysisFillAOD() 
425 {
426   //Do analysis and fill aods
427   //Search for photons in fCalorimeter 
428   TObjArray * pl = new TObjArray(); 
429   
430   //Get vertex for photon momentum calculation
431   Double_t vertex[]={0,0,0} ; //vertex ;
432   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
433   //Select the Calorimeter of the photon
434   if(fCalorimeter == "PHOS")
435     pl = GetAODPHOS();
436   else if (fCalorimeter == "EMCAL")
437     pl = GetAODEMCAL();
438   
439   //Fill AODCaloClusters and AODParticle with PHOS aods
440   TLorentzVector mom ;
441   
442   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
443     AliAODCaloCluster * calo =  (AliAODCaloCluster*) (pl->At(icalo));   
444     //Cluster selection, not charged, with photon id and in fidutial cut
445     //Get Momentum vector, 
446     calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
447     //If too small or big pt, skip it
448     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
449           
450         //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %f, pT %f, Ecl %f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), mom.E(), mom.Pt(),calo->E());
451
452     //Check acceptance selection
453     if(IsFidutialCutOn()){
454       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
455       if(! in ) continue ;
456     }
457         
458         //Input from second AOD?
459         Int_t input = 0;
460         if     (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1;
461         else if(fCalorimeter == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= icalo) input = 1;
462
463         //Create AOD for analysis
464     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
465         Int_t label = calo->GetLabel(0);
466         if(input == 1) label+=GetLabelShift();
467     aodph.SetLabel(label);
468     //printf("Index %d, Id %d\n",icalo, calo->GetID());
469     //Set the indeces of the original caloclusters  
470     aodph.SetCaloLabel(calo->GetID(),-1);
471     aodph.SetDetector(fCalorimeter);
472     if(GetDebug() > 1) 
473       printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());  
474     
475     //Check Distance to Bad channel, set bit.
476     Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
477     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
478     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
479       continue ;
480     
481     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
482     
483     if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
484     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
485     else aodph.SetDistToBad(0) ;
486     
487     //Check PID
488     //PID selection or bit setting
489     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
490       //Get most probable PID, check PID weights (in MC this option is mandatory)
491       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
492       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());         
493       //If primary is not photon, skip it.
494       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
495     }                                   
496     else if(IsCaloPIDOn()){
497       //Skip matched clusters with tracks
498       if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
499       
500       //Get most probable PID, 2 options check PID weights 
501       //or redo PID, recommended option for EMCal.              
502       if(!IsCaloPIDRecalculationOn())
503         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
504       else
505         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
506       
507       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
508       
509       //If cluster does not pass pid, not photon, skip it.
510       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;                      
511       
512     }
513     else{
514       //Set PID bits for later selection (AliAnaPi0 for example)
515       //GetPDG already called in SetPIDBits.
516       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
517       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");               
518     }
519     
520     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
521     
522     //Play with the MC stack if available
523     //Check origin of the candidates
524     if(IsDataMC()){
525                 
526       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetReader(), input));
527       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
528     }//Work with stack also   
529     
530     //Add AOD with photon object to aod branch
531     AddAODParticle(aodph);
532     
533   }//loop
534   
535   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs \n");  
536   
537 }
538
539 //__________________________________________________________________
540 void  AliAnaPhoton::MakeAnalysisFillHistograms() 
541 {
542     //Do analysis and fill histograms
543             
544         // Access MC information in stack if requested, check that it exists.   
545         AliStack * stack = 0x0;
546         TParticle * primary = 0x0;   
547         TClonesArray * mcparticles0 = 0x0;
548         TClonesArray * mcparticles1 = 0x0;
549         AliAODMCParticle * aodprimary = 0x0; 
550         if(IsDataMC()){
551                 
552                 if(GetReader()->ReadStack()){
553                         stack =  GetMCStack() ;
554                         if(!stack) {
555                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
556                                 abort();
557                         }
558                 
559                 }
560                 else if(GetReader()->ReadAODMCParticles()){
561         
562                         //Get the list of MC particles
563                         mcparticles0 = GetReader()->GetAODMCParticles(0);
564                         if(!mcparticles0 && GetDebug() > 0)     {
565                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
566                         }       
567                         if(GetReader()->GetSecondInputAODTree()){
568                                 mcparticles1 = GetReader()->GetAODMCParticles(1);
569                                 if(!mcparticles1 && GetDebug() > 0)     {
570                                         printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
571                                 }
572                         }               
573                         
574                 }
575         }// is data and MC
576         
577         //Loop on stored AOD photons
578         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
579         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
580         
581         for(Int_t iaod = 0; iaod < naod ; iaod++){
582     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
583     Int_t pdg = ph->GetPdg();
584     
585     if(GetDebug() > 3) 
586       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
587     
588     //If PID used, fill histos with photons in Calorimeter fCalorimeter
589     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue; 
590     if(ph->GetDetector() != fCalorimeter) continue;
591     
592     if(GetDebug() > 2) 
593       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
594     
595     //Fill photon histograms 
596     Float_t ptcluster  = ph->Pt();
597     Float_t phicluster = ph->Phi();
598     Float_t etacluster = ph->Eta();
599     Float_t ecluster   = ph->E();
600    
601     fhPtPhoton  ->Fill(ptcluster);
602     fhPhiPhoton ->Fill(ptcluster,phicluster);
603     fhEtaPhoton ->Fill(ptcluster,etacluster);
604
605     //Play with the MC data if available
606     if(IsDataMC()){
607                 
608       Int_t tag =ph->GetTag();
609                 
610         if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
611                 {
612                         fhPtMCPhoton  ->Fill(ptcluster);
613                         fhPhiMCPhoton ->Fill(ptcluster,phicluster);
614                         fhEtaMCPhoton ->Fill(ptcluster,etacluster);
615                                 
616                         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
617                         {
618                                 fhPtConversion  ->Fill(ptcluster);
619                                 fhPhiConversion ->Fill(ptcluster,phicluster);
620                                 fhEtaConversion ->Fill(ptcluster,etacluster);
621                         }                       
622                         
623                         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
624                                 fhPtPrompt  ->Fill(ptcluster);
625                                 fhPhiPrompt ->Fill(ptcluster,phicluster);
626                                 fhEtaPrompt ->Fill(ptcluster,etacluster);
627                         }
628                         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
629                         {
630                                 fhPtFragmentation  ->Fill(ptcluster);
631                                 fhPhiFragmentation ->Fill(ptcluster,phicluster);
632                                 fhEtaFragmentation ->Fill(ptcluster,etacluster);
633                         }
634                         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
635                         {
636                                 fhPtISR  ->Fill(ptcluster);
637                                 fhPhiISR ->Fill(ptcluster,phicluster);
638                                 fhEtaISR ->Fill(ptcluster,etacluster);
639                         }
640                         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
641                         {
642                                 fhPtPi0Decay  ->Fill(ptcluster);
643                                 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
644                                 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
645                         }
646                         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
647                         {
648                                 fhPtOtherDecay  ->Fill(ptcluster);
649                                 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
650                                 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
651                         }
652                 }
653       else{
654                   fhPtUnknown  ->Fill(ptcluster);
655                   fhPhiUnknown ->Fill(ptcluster,phicluster);
656                   fhEtaUnknown ->Fill(ptcluster,etacluster);
657
658 //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
659 //                                      ph->GetLabel(),ph->Pt());
660 //                for(Int_t i = 0; i < 20; i++) {
661 //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
662 //                }
663 //                printf("\n");
664         
665       }
666         
667                 
668         // Access MC information in stack if requested, check that it exists.
669         Int_t label =ph->GetLabel();
670         if(label < 0) {
671                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
672                 continue;
673         }
674                 
675         Float_t eprim   = 0;
676         Float_t ptprim  = 0;
677         if(GetReader()->ReadStack()){
678                                 
679                 if(label >=  stack->GetNtrack()) {
680                         if(GetDebug() > 2)  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
681                         continue ;
682                 }
683                 
684                 primary = stack->Particle(label);
685                 if(!primary){
686                         printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
687                         continue;
688                 }
689                 eprim   = primary->Energy();
690                 ptprim  = primary->Pt();                
691                 
692         }
693         else if(GetReader()->ReadAODMCParticles()){
694                 //Check which is the input
695                 if(label < GetLabelShift()){
696                         if(!mcparticles0) continue;
697                         if(label >=  mcparticles0->GetEntriesFast()) {
698                                 if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", 
699                                                                                    label, mcparticles0->GetEntriesFast());
700                                 continue ;
701                         }
702                         //Get the particle
703                         aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
704
705                 }
706                 else {//Second input
707                         if(!mcparticles1) continue;
708                         label-=GetLabelShift();
709                         if(label >=  mcparticles1->GetEntriesFast()) {
710                                 if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", 
711                                                                                    label, mcparticles1->GetEntriesFast());
712                                 continue ;
713                         }
714                         //Get the particle
715                         aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
716
717                 }//second input
718
719                 if(!aodprimary){
720                         printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
721                         continue;
722                 }
723                 
724                 eprim   = aodprimary->E();
725                 ptprim  = aodprimary->Pt();
726                 
727         }
728         
729         fh2E     ->Fill(ecluster, eprim);
730         fh2Pt    ->Fill(ptcluster, ptprim);     
731         fhDeltaE ->Fill(eprim-ecluster);
732         fhDeltaPt->Fill(ptprim-ptcluster);     
733         if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
734         if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);              
735                 
736     }//Histograms with MC
737     
738   }// aod loop
739   
740 }
741
742
743 //__________________________________________________________________
744 void AliAnaPhoton::Print(const Option_t * opt) const
745 {
746   //Print some relevant parameters set for the analysis
747   
748   if(! opt)
749     return;
750   
751   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
752   AliAnaPartCorrBaseClass::Print(" ");
753
754   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
755   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
756   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
757   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
758   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
759
760   printf("    \n") ;
761         
762