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