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