Move histogram setting ranges from analysis modules to base class (to be moved to...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0EbE.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 is 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: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17 //_________________________________________________________________________
18 // Class for the analysis of high pT pi0 event by event
19 // Pi0 identified by one of the following:
20 //  -Invariant mass of 2 cluster in calorimeter
21 //  -Shower shape analysis in calorimeter
22 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
23 //
24 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
25 //////////////////////////////////////////////////////////////////////////////
26   
27   
28 // --- ROOT system --- 
29 #include <TList.h>
30 #include <TClonesArray.h>
31 #include <TObjString.h>
32 #include <TH3F.h>
33 //#include "Riostream.h"
34
35 // --- Analysis system --- 
36 #include "AliAnaPi0EbE.h" 
37 #include "AliCaloTrackReader.h"
38 #include "AliIsolationCut.h"
39 #include "AliNeutralMesonSelection.h"
40 #include "AliCaloPID.h"
41 #include "AliMCAnalysisUtils.h"
42 #include "AliStack.h"
43 #include "AliFiducialCut.h"
44 #include "TParticle.h"
45 #include "AliVCluster.h"
46 #include "AliAODEvent.h"
47 #include "AliAODMCParticle.h"
48
49 ClassImp(AliAnaPi0EbE)
50   
51 //____________________________________________________________________________
52 AliAnaPi0EbE::AliAnaPi0EbE() : 
53     AliAnaPartCorrBaseClass(),    fAnaType(kIMCalo),            fCalorimeter(""),
54     fMinDist(0.),fMinDist2(0.),   fMinDist3(0.),        
55     fInputAODGammaConv(0x0),      fInputAODGammaConvName(""),
56     //Histograms
57     fhPtPi0(0),                   fhEPi0(0),                    fhEEtaPhiPi0(0),
58     //Shower shape histos
59     fhEDispersion(0),             fhELambda0(0),                fhELambda1(0),             
60     fhEdDispersion(0),            fhEdLambda0(0),               fhEdLambda1(0),            
61     //Time histograms
62     fhClusterPairDiffTimeE(0),    fhClusterPairDiffTimeAsy(0),
63     //MC histos
64     fhPtMCNoPi0(0),               fhPhiMCNoPi0(0),              fhEtaMCNoPi0(0), 
65     fhPtMCPi0(0),                 fhPhiMCPi0(0),                fhEtaMCPi0(0)
66 {
67   //default ctor
68   
69   for(Int_t i = 0; i < 5; i++){
70     fhEMCLambda0[i]     = 0;
71     fhEMCLambda1[i]     = 0;
72     fhEMCDispersion[i]  = 0;
73     fhEMCdLambda0[i]    = 0;
74     fhEMCdLambda1[i]    = 0;
75     fhEMCdDispersion[i] = 0;
76   }
77   
78   //Initialize parameters
79   InitParameters();
80   
81 }
82
83 //____________________________________________________________________________
84 AliAnaPi0EbE::~AliAnaPi0EbE() 
85 {
86   //dtor
87   if(fInputAODGammaConv){
88     fInputAODGammaConv->Clear() ; 
89     delete fInputAODGammaConv ;
90   }
91 }
92
93 //________________________________________________________________________
94 TObjString *  AliAnaPi0EbE::GetAnalysisCuts()
95 {       
96         //Save parameters used for analysis
97   TString parList ; //this will be list of parameters used for this analysis.
98   const Int_t buffersize = 255;
99   char onePar[buffersize] ;
100   
101   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
102   parList+=onePar ;     
103   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
104   parList+=onePar ;
105   
106   if(fAnaType == kSSCalo){
107     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
108     parList+=onePar ;
109     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
110     parList+=onePar ;
111     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
112     parList+=onePar ;
113     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
114     parList+=onePar ;
115   }
116   
117   //Get parameters set in base class.
118   parList += GetBaseParametersList() ;
119   
120   //Get parameters set in PID class.
121   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
122   
123   return new TObjString(parList) ;
124 }
125
126 //________________________________________________________________________
127 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
128 {  
129   // Create histograms to be saved in output file and 
130   // store them in outputContainer
131   TList * outputContainer = new TList() ; 
132   outputContainer->SetName("Pi0EbEHistos") ; 
133   
134   Int_t nptbins  = GetHistoPtBins();           Float_t ptmax  = GetHistoPtMax();           Float_t ptmin  = GetHistoPtMin();
135   Int_t nphibins = GetHistoPhiBins();          Float_t phimax = GetHistoPhiMax();          Float_t phimin = GetHistoPhiMin();
136   Int_t netabins = GetHistoEtaBins();          Float_t etamax = GetHistoEtaMax();          Float_t etamin = GetHistoEtaMin();
137   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax  = GetHistoShowerShapeMax();  Float_t ssmin  = GetHistoShowerShapeMin();
138   Int_t tdbins   = GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistoDiffTimeMax();     Float_t tdmin  = GetHistoDiffTimeMin();
139   
140   fhPtPi0  = new TH1F("hPtPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax); 
141   fhPtPi0->SetYTitle("N");
142   fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
143   outputContainer->Add(fhPtPi0) ; 
144   
145   fhEPi0  = new TH1F("hEPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax); 
146   fhEPi0->SetYTitle("N");
147   fhEPi0->SetXTitle("E  #pi^{0}(GeV)");
148   outputContainer->Add(fhEPi0) ; 
149   
150   fhEEtaPhiPi0  = new TH3F
151   ("hEEtaPhiPi0","Selected #pi^{0} pairs: E vs #eta vs #phi",nptbins,ptmin,ptmax,netabins,etamin,etamax, nphibins,phimin,phimax); 
152   fhEEtaPhiPi0->SetZTitle("#phi");
153   fhEEtaPhiPi0->SetYTitle("#eta");
154   fhEEtaPhiPi0->SetXTitle("E (GeV)");
155   outputContainer->Add(fhEEtaPhiPi0) ; 
156   
157   ////////
158   
159   if(fAnaType == kIMCalo){
160     
161     fhEDispersion  = new TH2F
162     ("hEDispersion","Selected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
163     fhEDispersion->SetYTitle("D^{2}");
164     fhEDispersion->SetXTitle("E (GeV)");
165     outputContainer->Add(fhEDispersion) ; 
166     
167     fhELambda0  = new TH2F
168     ("hELambda0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
169     fhELambda0->SetYTitle("#lambda_{0}^{2}");
170     fhELambda0->SetXTitle("E (GeV)");
171     outputContainer->Add(fhELambda0) ; 
172     
173     fhELambda1  = new TH2F
174     ("hELambda1","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
175     fhELambda1->SetYTitle("#lambda_{1}^{2}");
176     fhELambda1->SetXTitle("E (GeV)");
177     outputContainer->Add(fhELambda1) ; 
178     
179     
180     fhEdDispersion  = new TH2F
181     ("hEdDispersion","Selected #pi^{0} pairs: E vs dispersion^{2}/N_{cells}^{2}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
182     fhEdDispersion->SetYTitle("D^{2}");
183     fhEdDispersion->SetXTitle("E (GeV)");
184     outputContainer->Add(fhEdDispersion) ; 
185     
186     fhEdLambda0  = new TH2F
187     ("hEdLambda0","Selected #pi^{0} pairs: E vs #lambda_{0}^{2}/N_{cells}^{2}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
188     fhEdLambda0->SetYTitle("#lambda_{0}^{2}");
189     fhEdLambda0->SetXTitle("E (GeV)");
190     outputContainer->Add(fhEdLambda0) ; 
191     
192     fhEdLambda1  = new TH2F
193     ("hEdLambda1","Selected #pi^{0} pairs: E vs #lambda_{1}^{2}/N_{cells}^{2}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
194     fhEdLambda1->SetYTitle("#lambda_{1}^{2}");
195     fhEdLambda1->SetXTitle("E (GeV)");
196     outputContainer->Add(fhEdLambda1) ; 
197     
198     
199     fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
200     fhClusterPairDiffTimeE->SetXTitle("E_{pair} (GeV)");
201     fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
202     outputContainer->Add(fhClusterPairDiffTimeE);
203     
204     fhClusterPairDiffTimeAsy = new TH2F("hClusterPairDiffTime","cluster pair time difference vs pair asymmetry",100,0,1, tdbins,tdmin,tdmax);
205     fhClusterPairDiffTimeAsy->SetXTitle("Asymmetry");
206     fhClusterPairDiffTimeAsy->SetYTitle("#Delta t (ns)");
207     outputContainer->Add(fhClusterPairDiffTimeAsy);    
208     
209   }// Invariant mass analysis in calorimeters only
210   
211   if(IsDataMC()) {
212     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
213        GetReader()->GetDataType() != AliCaloTrackReader::kMC){
214       
215       fhPtMCPi0  = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax); 
216       fhPtMCPi0->SetYTitle("N");
217       fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
218       outputContainer->Add(fhPtMCPi0) ; 
219       
220       fhPhiMCPi0  = new TH2F
221       ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
222       fhPhiMCPi0->SetYTitle("#phi");
223       fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
224       outputContainer->Add(fhPhiMCPi0) ; 
225       
226       fhEtaMCPi0  = new TH2F
227       ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
228       fhEtaMCPi0->SetYTitle("#eta");
229       fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
230       outputContainer->Add(fhEtaMCPi0) ;
231       
232       fhPtMCNoPi0  = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax); 
233       fhPtMCNoPi0->SetYTitle("N");
234       fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
235       outputContainer->Add(fhPtMCNoPi0) ; 
236       
237       fhPhiMCNoPi0  = new TH2F
238       ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
239       fhPhiMCNoPi0->SetYTitle("#phi");
240       fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
241       outputContainer->Add(fhPhiMCNoPi0) ; 
242       
243       fhEtaMCNoPi0  = new TH2F
244       ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
245       fhEtaMCNoPi0->SetYTitle("#eta");
246       fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
247       outputContainer->Add(fhEtaMCNoPi0) ;
248       
249       if(fAnaType == kIMCalo){
250         TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","e^{#pm}", "hadron"}; 
251         TString pname[] ={"Photon","Conversion",     "Pi0",    "Electron","Hadron"};
252         for(Int_t i = 0; i < 5; i++){ 
253           
254           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
255                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
256                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
257           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
258           fhEMCLambda0[i]->SetXTitle("E (GeV)");
259           outputContainer->Add(fhEMCLambda0[i]) ; 
260           
261           fhEMCdLambda0[i]  = new TH2F(Form("hEdLambda0_MC%s",pname[i].Data()),
262                                        Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}/N_{cells}^{2}",ptype[i].Data()),
263                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
264           fhEMCdLambda0[i]->SetYTitle("d#lambda_{0}^{2}");
265           fhEMCdLambda0[i]->SetXTitle("E (GeV)");
266           outputContainer->Add(fhEMCdLambda0[i]) ; 
267           
268           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
269                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
270                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
271           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
272           fhEMCLambda1[i]->SetXTitle("E (GeV)");
273           outputContainer->Add(fhEMCLambda1[i]) ; 
274           
275           fhEMCdLambda1[i]  = new TH2F(Form("hEdLambda1_MC%s",pname[i].Data()),
276                                        Form("Selected pair, cluster from %s : E vs d#lambda_{1}^{2}/N_{cells}^{2}",ptype[i].Data()),
277                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
278           fhEMCdLambda1[i]->SetYTitle("d#lambda_{1}^{2}");
279           fhEMCdLambda1[i]->SetXTitle("E (GeV)");
280           outputContainer->Add(fhEMCdLambda1[i]) ; 
281           
282           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
283                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
284                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
285           fhEMCDispersion[i]->SetYTitle("D^{2}");
286           fhEMCDispersion[i]->SetXTitle("E (GeV)");
287           outputContainer->Add(fhEMCDispersion[i]) ; 
288           
289           fhEMCdDispersion[i]  = new TH2F(Form("hEdDispersion_MC%s",pname[i].Data()),
290                                           Form("Selected pair, cluster from %s : E vs dispersion^{2}/N_{cells}^{2}",ptype[i].Data()),
291                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50); 
292           fhEMCdDispersion[i]->SetYTitle("dD^{2}");
293           fhEMCdDispersion[i]->SetXTitle("E (GeV)");
294           outputContainer->Add(fhEMCdDispersion[i]) ;   
295           
296         }//
297         
298       }//kIMCalo
299     } //Not MC reader
300   }//Histos with MC
301   
302   
303   //Keep neutral meson selection histograms if requiered
304   //Setting done in AliNeutralMesonSelection
305   
306   if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
307     
308     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
309     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
310       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
311     delete nmsHistos;
312           
313   }
314   
315   return outputContainer ;
316   
317 }
318
319 //____________________________________________________________________________
320 void AliAnaPi0EbE::Init()
321
322   //Init
323   //Do some checks
324   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
325     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
326     abort();
327   }
328   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
329     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
330     abort();
331   }
332   
333 }
334
335 //____________________________________________________________________________
336 void AliAnaPi0EbE::InitParameters()
337 {
338   //Initialize the parameters of the analysis.  
339   AddToHistogramsName("AnaPi0EbE_");
340   
341   fInputAODGammaConvName = "gammaconv" ;
342   fAnaType = kIMCalo ;
343   fCalorimeter = "EMCAL" ;
344   fMinDist  = 2.;
345   fMinDist2 = 4.;
346   fMinDist3 = 5.;
347   
348 }
349
350 //__________________________________________________________________
351 void  AliAnaPi0EbE::MakeAnalysisFillAOD() 
352 {
353   //Do analysis and fill aods
354   
355   switch(fAnaType) 
356   {
357     case kIMCalo:
358       MakeInvMassInCalorimeter();
359       break;
360       
361     case kSSCalo:
362       MakeShowerShapeIdentification();
363       break;
364       
365     case kIMCaloTracks:
366       MakeInvMassInCalorimeterAndCTS();
367       break;
368       
369   }
370 }
371
372 //__________________________________________________________________
373 void  AliAnaPi0EbE::MakeInvMassInCalorimeter() 
374 {
375   //Do analysis and fill aods
376   //Search for the photon decay in calorimeters
377   //Read photon list from AOD, produced in class AliAnaPhoton
378   //Check if 2 photons have the mass of the pi0.
379   
380   TLorentzVector mom1;
381   TLorentzVector mom2;
382   TLorentzVector mom ;
383   Int_t tag1 = 0;
384   Int_t tag2 = 0;
385   Int_t tag  = 0;
386   
387   if(!GetInputAODBranch()){
388     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
389     abort();
390   }
391   
392   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
393     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
394     
395     //Vertex cut in case of mixed events
396     Int_t evtIndex1 = 0 ; 
397     if(GetMixedEvent())
398       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
399     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
400     mom1 = *(photon1->Momentum());
401     
402     //Get shower shape information of clusters
403     TObjArray *clusters = 0;
404     if     (fCalorimeter="EMCAL") clusters = GetEMCALClusters();
405     else if(fCalorimeter="PHOS" ) clusters = GetPHOSClusters() ;
406     
407     Bool_t bFound1        = kFALSE;
408     Int_t  caloLabel1     = photon1->GetCaloLabel(0);
409     Bool_t iclus1         = -1;    
410     AliVCluster *cluster1 = 0; 
411     if(clusters) {
412       //Get original cluster, to recover some information
413       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
414         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
415         if(cluster){
416           if     (cluster->GetID()==caloLabel1) {
417             bFound1  = kTRUE  ;
418             cluster1 = cluster;
419             iclus1   = iclus;
420           }
421         }      
422         if(bFound1) break;
423       }// calorimeter clusters loop
424     }
425     
426     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++){
427       
428       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
429       Int_t evtIndex2 = 0 ; 
430       if(GetMixedEvent())
431         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
432       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
433         continue ; 
434       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
435       mom2 = *(photon2->Momentum());
436       
437       //Photon1
438       Float_t e1    = photon1->E();
439       Float_t tof1  = -1;
440       Float_t disp1 = -1;
441       Float_t l01   = -1;
442       Float_t l11   = -1; 
443       Int_t ncells1 = 0;
444       
445       //Photon2
446       Float_t e2    = photon2->E();
447       Float_t disp2 = -1;
448       Float_t tof2  = -1;
449       Float_t l02   = -1;
450       Float_t l12   = -1; 
451       Int_t ncells2 = 0;
452       Bool_t bFound2        = kFALSE;
453       Int_t  caloLabel2     = photon2->GetCaloLabel(0);
454       AliVCluster *cluster2 = 0; 
455       if(clusters){      
456         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
457           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
458           if(cluster){
459             if(cluster->GetID()==caloLabel2) {
460               bFound2  = kTRUE  ;
461               cluster2 = cluster;
462             }          
463           }      
464           if(bFound2) break;
465         }// calorimeter clusters loop
466         
467         //Photon/Cluster 1
468         if(cluster1 && bFound1){
469           disp1   = cluster1->GetDispersion()*cluster1->GetDispersion();
470           l11     = cluster1->GetM20();
471           l01     = cluster1->GetM02();
472           tof1    = cluster1->GetTOF()*1e9;
473           ncells1 = cluster1->GetNCells()*cluster1->GetNCells();
474         }
475         //        else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
476         //                     photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
477         
478         //Photon/Cluster 2
479         if(cluster2 && bFound2){
480           disp2   = cluster2->GetDispersion()*cluster2->GetDispersion();
481           l12     = cluster2->GetM20();
482           l02     = cluster2->GetM02();
483           tof2    = cluster2->GetTOF()*1e9;
484           ncells2 = cluster2->GetNCells()*cluster2->GetNCells();
485         }
486         //        else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
487         //                     photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());    
488         
489         //Select clusters with good time window difference
490         Double_t t12diff = tof1-tof2;
491         Float_t asymmetry = TMath::Abs(e1-e2)/(e1+e2);
492         fhClusterPairDiffTimeE  ->Fill(e1+e2,    t12diff);
493         fhClusterPairDiffTimeAsy->Fill(asymmetry,t12diff);
494         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
495       }
496       
497       //Select good pair (good phi, pt cuts, aperture and invariant mass)
498       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
499       {
500         if(GetDebug()>1) 
501           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
502         
503         //Play with the MC stack if available
504         if(IsDataMC()){
505           //Check origin of the candidates
506           Int_t  label1 = photon1->GetLabel();
507           Int_t  label2 = photon2->GetLabel();
508           if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
509           if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
510           
511           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
512           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && 
513              GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
514             
515             //Check if pi0 mother is the same
516             if(GetReader()->ReadStack()){ 
517               if(label1>=0){
518                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
519                 label1 = mother1->GetFirstMother();
520                 //mother1 = GetMCStack()->Particle(label1);//pi0
521               }
522               if(label2>=0){
523                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
524                 label2 = mother2->GetFirstMother();
525                 //mother2 = GetMCStack()->Particle(label2);//pi0
526               }
527             }
528             else if(GetReader()->ReadAODMCParticles()){//&& (input > -1)){
529               if(label1>=0){
530                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
531                 label1 = mother1->GetMother();
532                 //mother1 = GetMCStack()->Particle(label1);//pi0
533               }
534               if(label2>=0){
535                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
536                 label2 = mother2->GetMother();
537                 //mother2 = GetMCStack()->Particle(label2);//pi0
538               }
539             }
540             
541             //printf("mother1 %d, mother2 %d\n",label1,label2);
542             if(label1 == label2 && label1>=0)
543               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
544           }
545         }//Work with stack also   
546         
547         //Fill some histograms about shower shape
548         if(clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
549           //Photon1 
550           
551           //printf("Signal Cl1: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",
552           //       e,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
553           
554           fhEDispersion   ->Fill(e1, disp1);   
555           fhELambda0->Fill(e1, l01  );  
556           fhELambda1->Fill(e1, l11  );  
557           
558           fhEdDispersion   ->Fill(e1, disp1/ncells1);   
559           fhEdLambda0->Fill(e1, l01/ncells1  );  
560           fhEdLambda1->Fill(e1, l11/ncells1  );  
561           
562           //Photon2
563           //printf("Signal Cl2: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",e
564           //     ,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
565           
566           fhEDispersion   ->Fill(e2, disp2);   
567           fhELambda0->Fill(e2, l02  ); 
568           fhELambda1->Fill(e2, l12  ); 
569           
570           fhEdDispersion   ->Fill(e2, disp2/ncells2);   
571           fhEdLambda0->Fill(e2, l02/ncells2  ); 
572           fhEdLambda1->Fill(e2, l12/ncells2  ); 
573           
574           if(IsDataMC()) {
575             //Photon1
576             if( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPhoton) && 
577                !GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
578               fhEMCLambda0[0]    ->Fill(e1, l01);
579               fhEMCdLambda0[0]   ->Fill(e1, l01/ncells1);
580               fhEMCLambda1[0]    ->Fill(e1, l11);
581               fhEMCdLambda1[0]   ->Fill(e1, l11/ncells1);
582               fhEMCDispersion[0] ->Fill(e1, disp1);
583               fhEMCdDispersion[0]->Fill(e1, disp1/ncells1);                
584             }//photon   no conversion
585             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCElectron)){
586               fhEMCLambda0[3]    ->Fill(e1, l01);
587               fhEMCdLambda0[3]   ->Fill(e1, l01/ncells1);
588               fhEMCLambda1[3]    ->Fill(e1, l11);
589               fhEMCdLambda1[3]   ->Fill(e1, l11/ncells1);
590               fhEMCDispersion[3] ->Fill(e1, disp1);
591               fhEMCdDispersion[3]->Fill(e1, disp1/ncells1);
592             }//electron
593             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) && 
594                       GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
595               fhEMCLambda0[1]    ->Fill(e1, l01);
596               fhEMCdLambda0[1]   ->Fill(e1, l01/ncells1);
597               fhEMCLambda1[1]    ->Fill(e1, l11);
598               fhEMCdLambda1[1]   ->Fill(e1, l11/ncells1);
599               fhEMCDispersion[1] ->Fill(e1, disp1);
600               fhEMCdDispersion[1]->Fill(e1, disp1/ncells1);           
601             }//conversion photon
602             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0)  ){
603               fhEMCLambda0[2]    ->Fill(e1, l01);
604               fhEMCdLambda0[2]   ->Fill(e1, l01/ncells1);
605               fhEMCLambda1[2]    ->Fill(e1, l11);
606               fhEMCdLambda1[2]   ->Fill(e1, l11/ncells1);
607               fhEMCDispersion[2] ->Fill(e1, disp1);
608               fhEMCdDispersion[2]->Fill(e1, disp1/ncells1);     
609             }//pi0
610             else {
611               fhEMCLambda0[4]    ->Fill(e1, l01);
612               fhEMCdLambda0[4]   ->Fill(e1, l01/ncells1);
613               fhEMCLambda1[4]    ->Fill(e1, l11);
614               fhEMCdLambda1[4]   ->Fill(e1, l11/ncells1);
615               fhEMCDispersion[4] ->Fill(e1, disp1);
616               fhEMCdDispersion[4]->Fill(e1, disp1/ncells1);    
617             }//other particles 
618             
619             //Photon 2
620             if( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPhoton) && 
621                !GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
622               fhEMCLambda0[0]    ->Fill(e2, l02);
623               fhEMCdLambda0[0]   ->Fill(e2, l02/ncells2);
624               fhEMCLambda1[0]    ->Fill(e2, l12);
625               fhEMCdLambda1[0]   ->Fill(e2, l12/ncells2);
626               fhEMCDispersion[0] ->Fill(e2, disp2);
627               fhEMCdDispersion[0]->Fill(e2, disp2/ncells2);            
628             }//photon   no conversion
629             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCElectron)){
630               fhEMCLambda0[3]    ->Fill(e2, l02);
631               fhEMCdLambda0[3]   ->Fill(e2, l02/ncells2);
632               fhEMCLambda1[3]    ->Fill(e2, l12);
633               fhEMCdLambda1[3]   ->Fill(e2, l12/ncells2);
634               fhEMCDispersion[3] ->Fill(e2, disp2);
635               fhEMCdDispersion[3]->Fill(e2, disp2/ncells2);    
636             }//electron
637             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && 
638                       GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
639               fhEMCLambda0[1]    ->Fill(e2, l02);
640               fhEMCdLambda0[1]   ->Fill(e2, l02/ncells2);
641               fhEMCLambda1[1]    ->Fill(e2, l12);
642               fhEMCdLambda1[1]   ->Fill(e2, l12/ncells2);
643               fhEMCDispersion[1] ->Fill(e2, disp2);
644               fhEMCdDispersion[1]->Fill(e2, disp2/ncells2);                
645             }//conversion photon
646             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0)  ){
647               fhEMCLambda0[2]    ->Fill(e2, l02);
648               fhEMCdLambda0[2]   ->Fill(e2, l02/ncells2);
649               fhEMCLambda1[2]    ->Fill(e2, l12);
650               fhEMCdLambda1[2]   ->Fill(e2, l12/ncells2);
651               fhEMCDispersion[2] ->Fill(e2, disp2);
652               fhEMCdDispersion[2]->Fill(e2, disp2/ncells2);                
653             }//pi0
654             else {
655               fhEMCLambda0[4]    ->Fill(e2, l02);
656               fhEMCdLambda0[4]   ->Fill(e2, l02/ncells2);
657               fhEMCLambda1[4]    ->Fill(e2, l12);
658               fhEMCdLambda1[4]   ->Fill(e2, l12/ncells2);
659               fhEMCDispersion[4] ->Fill(e2, disp2);
660               fhEMCdDispersion[4]->Fill(e2, disp2/ncells2);                
661             }//other particles 
662           }//is datamc
663         }//MC histograms
664         
665         //Create AOD for analysis
666         mom = mom1+mom2;
667         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
668         //pi0.SetLabel(calo->GetLabel());
669         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
670         pi0.SetDetector(photon1->GetDetector());
671         pi0.SetTag(tag);  
672         //Set the indeces of the original caloclusters  
673         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
674         //pi0.SetInputFileIndex(input);
675         AddAODParticle(pi0);
676       }//pi0
677       
678     }//2n photon loop
679     
680   }//1st photon loop
681   
682   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");  
683   
684 }
685
686 //__________________________________________________________________
687 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() 
688 {
689   //Do analysis and fill aods
690   //Search for the photon decay in calorimeters
691   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
692   //Check if 2 photons have the mass of the pi0.
693   
694   TLorentzVector mom1;
695   TLorentzVector mom2;
696   TLorentzVector mom ;
697   Int_t tag1 = 0;
698   Int_t tag2 = 0;
699   Int_t tag  = 0;
700   Int_t evtIndex = 0;
701   if(!GetInputAODBranch()){
702     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
703     abort();
704   }
705   
706   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
707     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
708     mom1 = *(photon1->Momentum());
709     
710     //Play with the MC stack if available
711     fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
712     if(!fInputAODGammaConv) {
713       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >, STOP\n",fInputAODGammaConvName.Data());
714       abort();  
715     }
716     for(Int_t jphoton = 0; jphoton < fInputAODGammaConv->GetEntriesFast(); jphoton++){
717       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (fInputAODGammaConv->At(jphoton));
718       if(GetMixedEvent())
719         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
720       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
721       
722       mom2 = *(photon2->Momentum());
723       
724       //Int_t input = -1;       //if -1 photons come from different files, not a pi0
725       //if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
726       
727       //Select good pair (good phi, pt cuts, aperture and invariant mass)
728       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
729         if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
730         
731         if(IsDataMC()){
732           Int_t label1 = photon1->GetLabel();
733           Int_t label2 = photon2->GetLabel();
734           if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
735           if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
736           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
737           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && 
738              GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
739             //Check if pi0 mother is the same
740             
741             if(GetReader()->ReadStack()){ 
742               if(label1>=0){
743                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
744                 label1 = mother1->GetFirstMother();
745                 //mother1 = GetMCStack()->Particle(label1);//pi0
746               }
747               if(label2>=0){
748                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
749                 label2 = mother2->GetFirstMother();
750                 //mother2 = GetMCStack()->Particle(label2);//pi0
751               }
752             }
753             else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
754               if(label1>=0){
755                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
756                 label1 = mother1->GetMother();
757                 //mother1 = GetMCStack()->Particle(label1);//pi0
758               }
759               if(label2>=0){
760                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
761                 label2 = mother2->GetMother();
762                 //mother2 = GetMCStack()->Particle(label2);//pi0
763               }
764             }
765             
766             //printf("mother1 %d, mother2 %d\n",label1,label2);
767             if(label1 == label2 && label1>=0)
768               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
769           }
770         }//Work with stack also   
771         
772         //Create AOD for analysis
773         mom = mom1+mom2;
774         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
775         //pi0.SetLabel(calo->GetLabel());
776         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
777         pi0.SetDetector(photon1->GetDetector());
778         pi0.SetTag(tag);
779         //Set the indeces of the original tracks or caloclusters  
780         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
781         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
782         //pi0.SetInputFileIndex(input);
783         AddAODParticle(pi0);
784       }//pi0
785     }//2n photon loop
786     
787   }//1st photon loop
788   
789   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");  
790   
791 }
792
793
794 //__________________________________________________________________
795 void  AliAnaPi0EbE::MakeShowerShapeIdentification() 
796 {
797   //Search for pi0 in fCalorimeter with shower shape analysis 
798   
799   TObjArray * pl = 0x0; 
800   //Select the Calorimeter of the photon
801   if(fCalorimeter == "PHOS")
802     pl = GetPHOSClusters();
803   else if (fCalorimeter == "EMCAL")
804     pl = GetEMCALClusters();
805   
806   if(!pl) {
807     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
808     return;
809   }  
810         
811   TLorentzVector mom ;
812   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
813     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));        
814     
815     Int_t evtIndex = 0 ; 
816     if (GetMixedEvent()) {
817       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
818     }
819     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
820     
821     //Get Momentum vector, 
822     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
823       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
824     else{
825       Double_t vertex[]={0,0,0};
826       calo->GetMomentum(mom,vertex) ;
827     }
828           
829     //If too small or big pt, skip it
830     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
831     //Check acceptance selection
832     if(IsFiducialCutOn()){
833       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
834       if(! in ) continue ;
835     }
836     
837     //Create AOD for analysis
838     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
839     aodpi0.SetLabel(calo->GetLabel());
840     //Set the indeces of the original caloclusters  
841     aodpi0.SetCaloLabel(calo->GetID(),-1);
842     aodpi0.SetDetector(fCalorimeter);
843     if(GetDebug() > 1) 
844       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());    
845     
846     //Check Distance to Bad channel, set bit.
847     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
848     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
849     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
850       continue ;
851     
852     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
853     
854     if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
855     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ; 
856     else aodpi0.SetDistToBad(0) ;
857     
858     //Check PID
859     //PID selection or bit setting
860     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
861       //Get most probable PID, check PID weights (in MC this option is mandatory)
862       aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
863       if(GetDebug() > 1) 
864         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
865       //If primary is not pi0, skip it.
866       if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
867     }                                   
868     else if(IsCaloPIDOn()){
869       //Skip matched clusters with tracks
870       if(IsTrackMatched(calo)) continue ;
871       
872       //Get most probable PID, 2 options check PID weights 
873       //or redo PID, recommended option for EMCal.              
874       if(!IsCaloPIDRecalculationOn())
875         aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
876       else
877         aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
878       
879       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
880       
881       //If cluster does not pass pid, not pi0, skip it.
882       if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;                     
883       
884     }
885     else{
886       //Set PID bits for later selection (AliAnaPi0 for example)
887       //GetPDG already called in SetPIDBits.
888       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils());
889       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");             
890     }
891     
892     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
893     
894     //Play with the MC stack if available
895     //Check origin of the candidates
896     if(IsDataMC()){
897       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
898          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
899         //aodpi0.SetInputFileIndex(input);
900         Int_t tag       =0;
901         tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
902         //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
903         aodpi0.SetTag(tag);
904         if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
905       }
906     }//Work with stack also   
907     
908     //Add AOD with pi0 object to aod branch
909     AddAODParticle(aodpi0);
910     
911   }//loop
912   
913   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");  
914   
915 }
916 //__________________________________________________________________
917 void  AliAnaPi0EbE::MakeAnalysisFillHistograms() 
918 {
919   //Do analysis and fill histograms
920   
921   if(!GetOutputAODBranch()){
922     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
923     abort();
924   }
925   //Loop on stored AOD pi0
926   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
927   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
928   
929   for(Int_t iaod = 0; iaod < naod ; iaod++){
930     
931     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
932     Int_t pdg = pi0->GetIdentifiedParticleType();
933           
934     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;              
935     
936     //Fill pi0 histograms 
937     Float_t ener  = pi0->E();
938     Float_t pt    = pi0->Pt();
939     Float_t phi   = pi0->Phi();
940     if(phi < 0) phi+=TMath::TwoPi();
941     Float_t eta = pi0->Eta();
942     
943     fhPtPi0      ->Fill(pt);
944     fhEPi0       ->Fill(ener);
945     fhEEtaPhiPi0 ->Fill(ener,eta,phi);
946     
947     if(IsDataMC()){
948       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
949          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
950         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
951           fhPtMCPi0  ->Fill(pt);
952           fhPhiMCPi0 ->Fill(pt,phi);
953           fhEtaMCPi0 ->Fill(pt,eta);
954         }
955         else{
956           fhPtMCNoPi0  ->Fill(pt);
957           fhPhiMCNoPi0 ->Fill(pt,phi);
958           fhEtaMCNoPi0 ->Fill(pt,eta);
959         }
960       }
961     }//Histograms with MC
962     
963   }// aod loop
964   
965 }
966
967 //__________________________________________________________________
968 void AliAnaPi0EbE::Print(const Option_t * opt) const
969 {
970   //Print some relevant parameters set for the analysis
971   if(! opt)
972     return;
973   
974   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
975   AliAnaPartCorrBaseClass::Print("");
976   printf("Analysis Type = %d \n",  fAnaType) ;
977   if(fAnaType == kSSCalo){     
978     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
979     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
980     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
981     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); 
982   } 
983   printf("    \n") ;
984   
985