Simplify AliAnaPhoton, move conversion utilities to other analysis class
[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 // Produces input for other analysis classes like AliAnaPi0, 
23 // AliAnaParticleHadronCorrelation ... 
24 //
25 // -- Author: Gustavo Conesa (LNF-INFN) 
26 //////////////////////////////////////////////////////////////////////////////
27   
28   
29 // --- ROOT system --- 
30 #include <TH2F.h>
31 #include <TH3D.h>
32 #include <TClonesArray.h>
33 #include <TObjString.h>
34 #include "TParticle.h"
35 #include "TDatabasePDG.h"
36
37 // --- Analysis system --- 
38 #include "AliAnaPhoton.h" 
39 #include "AliCaloTrackReader.h"
40 #include "AliStack.h"
41 #include "AliCaloPID.h"
42 #include "AliMCAnalysisUtils.h"
43 #include "AliFiducialCut.h"
44 #include "AliVCluster.h"
45 #include "AliAODMCParticle.h"
46 #include "AliMixedEvent.h"
47
48
49 ClassImp(AliAnaPhoton)
50   
51 //____________________________
52 AliAnaPhoton::AliAnaPhoton() : 
53     AliAnaPartCorrBaseClass(),    fCalorimeter(""), 
54     fMinDist(0.),                 fMinDist2(0.),                fMinDist3(0.), 
55     fRejectTrackMatch(0),         fTimeCutMin(-1),              fTimeCutMax(999999),         
56     fNCellsCut(0),                fFillSSHistograms(kFALSE),    
57     fNOriginHistograms(8),        fNPrimaryHistograms(4),
58
59     // Histograms
60     fhNCellsE(0),                 fhMaxCellDiffClusterE(0),
61     fhEPhoton(0),                 fhPtPhoton(0),  
62     fhPhiPhoton(0),               fhEtaPhoton(0), 
63     fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
64
65     // Shower shape histograms
66     fhDispE(0),                   fhLam0E(0),                   fhLam1E(0), 
67     fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
68
69     fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),  
70     fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
71
72     fhEtaLam0LowE(0),             fhPhiLam0LowE(0), 
73     fhEtaLam0HighE(0),            fhPhiLam0HighE(0), 
74     fhLam0DispLowE(0),            fhLam0DispHighE(0), 
75     fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0), 
76     fhDispLam1LowE(0),            fhDispLam1HighE(0),
77
78     // MC histograms
79     fhMCPhotonELambda0NoOverlap(0),     fhMCPhotonELambda0TwoOverlap(0),      fhMCPhotonELambda0NOverlap(0),
80     //Embedding
81     fhEmbeddedSignalFractionEnergy(0),     
82     fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonELambda0MostlySignal(0),  
83     fhEmbedPhotonELambda0MostlyBkg(0),  fhEmbedPhotonELambda0FullBkg(0),       
84     fhEmbedPi0ELambda0FullSignal(0),    fhEmbedPi0ELambda0MostlySignal(0),    
85     fhEmbedPi0ELambda0MostlyBkg(0),     fhEmbedPi0ELambda0FullBkg(0)         
86 {
87   //default ctor
88   
89   for(Int_t i = 0; i < 14; i++){
90     fhMCPt     [i] = 0;
91     fhMCE      [i] = 0;
92     fhMCPhi    [i] = 0;
93     fhMCEta    [i] = 0;
94     fhMCDeltaE [i] = 0;                
95     fhMCDeltaPt[i] = 0;
96     fhMC2E     [i] = 0;              
97     fhMC2Pt    [i] = 0;
98     
99   }
100   
101   for(Int_t i = 0; i < 7; i++){
102     fhPtPrimMC [i] = 0;
103     fhEPrimMC  [i] = 0;
104     fhPhiPrimMC[i] = 0;
105     fhYPrimMC  [i] = 0;
106     
107     fhPtPrimMCAcc [i] = 0;
108     fhEPrimMCAcc  [i] = 0;
109     fhPhiPrimMCAcc[i] = 0;
110     fhYPrimMCAcc  [i] = 0;
111   }  
112   
113   for(Int_t i = 0; i < 6; i++){
114     fhMCELambda0    [i]                  = 0;
115     fhMCELambda1    [i]                  = 0;
116     fhMCEDispersion [i]                  = 0;
117     fhMCNCellsE     [i]                  = 0; 
118     fhMCMaxCellDiffClusterE[i]           = 0; 
119     fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
120     fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
121     fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
122     fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
123     fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
124     fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
125   }
126   
127   //Initialize parameters
128   InitParameters();
129
130 }
131
132 //__________________________________________________________________________
133 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom) 
134 {
135   //Select clusters if they pass different cuts
136   if(GetDebug() > 2) 
137     printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
138            GetReader()->GetEventNumber(),
139            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
140   
141   //.......................................
142   //If too small or big energy, skip it
143   if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ; 
144   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
145   
146   //.......................................
147   // TOF cut, BE CAREFUL WITH THIS CUT
148   Double_t tof = calo->GetTOF()*1e9;
149   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
150   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
151   
152   //.......................................
153   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
154   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
155   
156   //.......................................
157   //Check acceptance selection
158   if(IsFiducialCutOn()){
159     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
160     if(! in ) return kFALSE ;
161   }
162   if(GetDebug() > 2) printf("Fiducial cut passed \n");
163   
164   //.......................................
165   //Skip matched clusters with tracks
166   if(fRejectTrackMatch){
167     if(IsTrackMatched(calo)) {
168       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
169       return kFALSE ;
170     }
171     else  
172       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
173   }// reject matched clusters
174   
175   //.......................................
176   //Check Distance to Bad channel, set bit.
177   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
178   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
179   if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
180     return kFALSE ;
181   }
182   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
183   //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
184
185   if(GetDebug() > 0) 
186     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
187            GetReader()->GetEventNumber(), 
188            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
189   
190   //All checks passed, cluster selected
191   return kTRUE;
192     
193 }
194
195 //_____________________________________________________________
196 void AliAnaPhoton::FillAcceptanceHistograms(){
197   //Fill acceptance histograms if MC data is available
198   
199   if(GetReader()->ReadStack()){ 
200     AliStack * stack = GetMCStack();
201     if(stack){
202       for(Int_t i=0 ; i<stack->GetNtrack(); i++){
203         TParticle * prim = stack->Particle(i) ;
204         Int_t pdg = prim->GetPdgCode();
205         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
206         //                             prim->GetName(), prim->GetPdgCode());
207         
208         if(pdg == 22){
209           
210           // Get tag of this particle photon from fragmentation, decay, prompt ...
211           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
212           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
213             //A conversion photon from a hadron, skip this kind of photon
214             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
215             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
216             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
217             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
218             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
219             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
220             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon), 
221             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
222             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton), 
223             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
224             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon), 
225             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton), 
226             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
227             
228             return;
229           }
230           
231           //Get photon kinematics
232           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception          
233           
234           Double_t photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
235           Double_t photonE   = prim->Energy() ;
236           Double_t photonPt  = prim->Pt() ;
237           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
238           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
239           Double_t photonEta = prim->Eta() ;
240           
241           
242           //Check if photons hit the Calorimeter
243           TLorentzVector lv;
244           prim->Momentum(lv);
245           Bool_t inacceptance = kFALSE;
246           if     (fCalorimeter == "PHOS"){
247             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
248               Int_t mod ;
249               Double_t x,z ;
250               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x)) 
251                 inacceptance = kTRUE;
252               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
253             }
254             else{
255               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
256                 inacceptance = kTRUE ;
257               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
258             }
259           }        
260           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
261             if(GetEMCALGeometry()){
262               
263               Int_t absID=0;
264               
265               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
266               
267               if( absID >= 0) 
268                 inacceptance = kTRUE;
269               
270               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2)) 
271               //                    inacceptance = kTRUE;
272               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
273             }
274             else{
275               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
276                 inacceptance = kTRUE ;
277               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
278             }
279           }       //In EMCAL
280           
281           //Fill histograms
282           
283           fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
284           if(TMath::Abs(photonY) < 1.0)
285           {
286             fhEPrimMC  [mcPPhoton]->Fill(photonE ) ;
287             fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
288             fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
289             fhYPrimMC[mcPPhoton]  ->Fill(photonE , photonEta) ;
290           }
291           if(inacceptance){
292             fhEPrimMCAcc[mcPPhoton]  ->Fill(photonE ) ;
293             fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
294             fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
295             fhYPrimMCAcc[mcPPhoton]  ->Fill(photonE , photonY) ;
296           }//Accepted
297           
298           //Origin of photon
299           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
300           {
301             fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
302             if(TMath::Abs(photonY) < 1.0){
303               fhEPrimMC  [mcPPrompt]->Fill(photonE ) ;
304               fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
305               fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
306               fhYPrimMC[mcPPrompt]  ->Fill(photonE , photonEta) ;
307             }   
308             if(inacceptance){
309               fhEPrimMCAcc[mcPPrompt]  ->Fill(photonE ) ;
310               fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
311               fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
312               fhYPrimMCAcc[mcPPrompt]  ->Fill(photonE , photonY) ;
313             }//Accepted
314           }
315           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation])
316           {
317             fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
318             if(TMath::Abs(photonY) < 1.0){
319               fhEPrimMC  [mcPFragmentation]->Fill(photonE ) ;
320               fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
321               fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
322               fhYPrimMC[mcPFragmentation]  ->Fill(photonE , photonEta) ;
323             }  
324             if(inacceptance){
325               fhEPrimMCAcc[mcPFragmentation]  ->Fill(photonE ) ;
326               fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
327               fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
328               fhYPrimMCAcc[mcPFragmentation]  ->Fill(photonE , photonY) ;
329             }//Accepted
330           }
331           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
332           {
333             fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
334             if(TMath::Abs(photonY) < 1.0){
335               fhEPrimMC  [mcPISR]->Fill(photonE ) ;
336               fhPtPrimMC [mcPISR]->Fill(photonPt) ;
337               fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
338               fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
339             }            
340             if(inacceptance){
341               fhEPrimMCAcc[mcPISR]  ->Fill(photonE ) ;
342               fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
343               fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
344               fhYPrimMCAcc[mcPISR]  ->Fill(photonE , photonY) ;
345             }//Accepted
346           }
347           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
348           {
349             fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
350             if(TMath::Abs(photonY) < 1.0){
351               fhEPrimMC  [mcPPi0Decay]->Fill(photonE ) ;
352               fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
353               fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
354               fhYPrimMC[mcPPi0Decay]  ->Fill(photonE , photonEta) ;
355             }     
356             if(inacceptance){
357               fhEPrimMCAcc[mcPPi0Decay]  ->Fill(photonE ) ;
358               fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
359               fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
360               fhYPrimMCAcc[mcPPi0Decay]  ->Fill(photonE , photonY) ;
361             }//Accepted
362           }
363           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || 
364                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[mcPOtherDecay])
365           {
366             fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
367             if(TMath::Abs(photonY) < 1.0){
368               fhEPrimMC  [mcPOtherDecay]->Fill(photonE ) ;
369               fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
370               fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
371               fhYPrimMC[mcPOtherDecay]  ->Fill(photonE , photonEta) ;
372             } 
373             if(inacceptance){
374               fhEPrimMCAcc[mcPOtherDecay]  ->Fill(photonE ) ;
375               fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
376               fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
377               fhYPrimMCAcc[mcPOtherDecay]  ->Fill(photonE , photonY) ;
378             }//Accepted
379           }
380           else if(fhEPrimMC[mcPOther])
381           {
382             fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
383             if(TMath::Abs(photonY) < 1.0){
384               fhEPrimMC  [mcPOther]->Fill(photonE ) ;
385               fhPtPrimMC [mcPOther]->Fill(photonPt) ;
386               fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
387               fhYPrimMC[mcPOther]  ->Fill(photonE , photonEta) ;
388             }
389             if(inacceptance){
390               fhEPrimMCAcc[mcPOther]  ->Fill(photonE ) ;
391               fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
392               fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
393               fhYPrimMCAcc[mcPOther]  ->Fill(photonE , photonY) ;
394             }//Accepted
395           }//Other origin
396         }// Primary photon 
397       }//loop on primaries      
398     }//stack exists and data is MC
399   }//read stack
400   else if(GetReader()->ReadAODMCParticles()){
401     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
402     if(mcparticles){
403       Int_t nprim = mcparticles->GetEntriesFast();
404       
405       for(Int_t i=0; i < nprim; i++)
406       {
407         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);   
408         
409         Int_t pdg = prim->GetPdgCode();
410         
411         if(pdg == 22){
412           
413           // Get tag of this particle photon from fragmentation, decay, prompt ...
414           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
415           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
416             //A conversion photon from a hadron, skip this kind of photon
417 //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
418 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
419 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
420 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
421 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
422 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
423 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon), 
424 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
425 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton), 
426 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
427 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon), 
428 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton), 
429 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
430             
431             return;
432           }
433           
434           //Get photon kinematics
435           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception      
436           
437           Double_t photonY  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
438           Double_t photonE   = prim->E() ;
439           Double_t photonPt  = prim->Pt() ;
440           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
441           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
442           Double_t photonEta = prim->Eta() ;
443           
444           //Check if photons hit the Calorimeter
445           TLorentzVector lv;
446           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
447           Bool_t inacceptance = kFALSE;
448           if     (fCalorimeter == "PHOS"){
449             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
450               Int_t mod ;
451               Double_t x,z ;
452               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
453               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
454                 inacceptance = kTRUE;
455               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
456             }
457             else{
458               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
459                 inacceptance = kTRUE ;
460               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
461             }
462           }        
463           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
464             if(GetEMCALGeometry()){
465               
466               Int_t absID=0;
467               
468               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
469               
470               if( absID >= 0) 
471                 inacceptance = kTRUE;
472               
473               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
474             }
475             else{
476               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter)) 
477                 inacceptance = kTRUE ;
478               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
479             }
480           }       //In EMCAL
481           
482           //Fill histograms
483           
484           fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
485           if(TMath::Abs(photonY) < 1.0)
486           {
487             fhEPrimMC  [mcPPhoton]->Fill(photonE ) ;
488             fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
489             fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
490             fhYPrimMC[mcPPhoton]->Fill(photonE , photonEta) ;
491           }
492           if(inacceptance){
493             fhEPrimMCAcc[mcPPhoton]  ->Fill(photonE ) ;
494             fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
495             fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
496             fhYPrimMCAcc[mcPPhoton]  ->Fill(photonE , photonY) ;
497           }//Accepted
498           
499           
500           //Origin of photon
501           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
502           {
503             fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
504             if(TMath::Abs(photonY) < 1.0){
505               fhEPrimMC  [mcPPrompt]->Fill(photonE ) ;
506               fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
507               fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
508               fhYPrimMC[mcPPrompt]->Fill(photonE , photonEta) ;
509             }   
510             if(inacceptance){
511               fhEPrimMCAcc[mcPPrompt]  ->Fill(photonE ) ;
512               fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
513               fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
514               fhYPrimMCAcc[mcPPrompt]  ->Fill(photonE , photonY) ;
515             }//Accepted
516           }
517           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation] )
518           {
519             fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
520             if(TMath::Abs(photonY) < 1.0){
521               fhEPrimMC  [mcPFragmentation]->Fill(photonE ) ;
522               fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
523               fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
524               fhYPrimMC[mcPFragmentation]->Fill(photonE , photonEta) ;
525             }  
526             if(inacceptance){
527               fhEPrimMCAcc[mcPFragmentation]  ->Fill(photonE ) ;
528               fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
529               fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
530               fhYPrimMCAcc[mcPFragmentation]  ->Fill(photonE , photonY) ;
531             }//Accepted
532           }
533           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
534           {
535             fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
536             if(TMath::Abs(photonY) < 1.0){
537               fhEPrimMC  [mcPISR]->Fill(photonE ) ;
538               fhPtPrimMC [mcPISR]->Fill(photonPt) ;
539               fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
540               fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
541             }            
542             if(inacceptance){
543               fhEPrimMCAcc[mcPISR]  ->Fill(photonE ) ;
544               fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
545               fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
546               fhYPrimMCAcc[mcPISR]  ->Fill(photonE , photonY) ;
547             }//Accepted
548           }
549           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
550           {
551             fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
552             if(TMath::Abs(photonY) < 1.0){
553               fhEPrimMC  [mcPPi0Decay]->Fill(photonE ) ;
554               fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
555               fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
556               fhYPrimMC[mcPPi0Decay]->Fill(photonE , photonEta) ;
557             }     
558             if(inacceptance){
559               fhEPrimMCAcc[mcPPi0Decay]  ->Fill(photonE ) ;
560               fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
561               fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
562               fhYPrimMCAcc[mcPPi0Decay]  ->Fill(photonE , photonY) ;
563             }//Accepted
564           }
565           else if((GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || 
566                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhEPrimMC[mcPOtherDecay])
567           {
568             fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
569             if(TMath::Abs(photonY) < 1.0){
570               fhEPrimMC  [mcPOtherDecay]->Fill(photonE ) ;
571               fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
572               fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
573               fhYPrimMC[mcPOtherDecay]->Fill(photonE , photonEta) ;
574             } 
575             if(inacceptance){
576               fhEPrimMCAcc[mcPOtherDecay]  ->Fill(photonE ) ;
577               fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
578               fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
579               fhYPrimMCAcc[mcPOtherDecay]  ->Fill(photonE , photonY) ;
580             }//Accepted
581           }
582           else if(fhEPrimMC[mcPOther])
583           {
584             fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
585             if(TMath::Abs(photonY) < 1.0){
586               fhEPrimMC  [mcPOther]->Fill(photonE ) ;
587               fhPtPrimMC [mcPOther]->Fill(photonPt) ;
588               fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
589               fhYPrimMC[mcPOther]->Fill(photonE , photonEta) ;
590             }
591             if(inacceptance){
592               fhEPrimMCAcc[mcPOther]  ->Fill(photonE ) ;
593               fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
594               fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
595               fhYPrimMCAcc[mcPOther]  ->Fill(photonE , photonY) ;
596             }//Accepted
597           }//Other origin
598         }// Primary photon 
599       }//loop on primaries      
600       
601     }//mc array exists and data is MC
602   }     // read AOD MC
603 }
604
605 //__________________________________________________________________
606 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag){
607   
608   //Fill cluster Shower Shape histograms
609   
610   if(!fFillSSHistograms || GetMixedEvent()) return;
611   
612   Float_t energy  = cluster->E();
613   Int_t   ncells  = cluster->GetNCells();
614   Float_t lambda0 = cluster->GetM02();
615   Float_t lambda1 = cluster->GetM20();
616   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
617   
618   TLorentzVector mom;
619   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
620     cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
621   else{
622     Double_t vertex[]={0,0,0};
623     cluster->GetMomentum(mom,vertex) ;
624   }
625   
626   Float_t eta = mom.Eta();
627   Float_t phi = mom.Phi();
628   if(phi < 0) phi+=TMath::TwoPi();
629   
630   fhLam0E ->Fill(energy,lambda0);
631   fhLam1E ->Fill(energy,lambda1);
632   fhDispE ->Fill(energy,disp);
633   
634   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
635     fhLam0ETRD->Fill(energy,lambda0);
636     fhLam1ETRD->Fill(energy,lambda1);
637     fhDispETRD->Fill(energy,disp);
638   }
639   
640   if(energy < 2){
641     fhNCellsLam0LowE ->Fill(ncells,lambda0);
642     fhNCellsLam1LowE ->Fill(ncells,lambda1);
643     fhNCellsDispLowE ->Fill(ncells,disp);
644     
645     fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
646     fhLam0DispLowE  ->Fill(lambda0,disp);
647     fhDispLam1LowE  ->Fill(disp,lambda1);
648     fhEtaLam0LowE   ->Fill(eta,lambda0);
649     fhPhiLam0LowE   ->Fill(phi,lambda0);  
650     
651   }
652   else {
653     fhNCellsLam0HighE ->Fill(ncells,lambda0);
654     fhNCellsLam1HighE ->Fill(ncells,lambda1);
655     fhNCellsDispHighE ->Fill(ncells,disp);
656
657     fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
658     fhLam0DispHighE  ->Fill(lambda0,disp);
659     fhDispLam1HighE  ->Fill(disp,lambda1);
660     fhEtaLam0HighE   ->Fill(eta, lambda0);
661     fhPhiLam0HighE   ->Fill(phi, lambda0);
662   }
663   
664   if(IsDataMC()){
665     
666     AliVCaloCells* cells = 0;
667     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
668     else                        cells = GetPHOSCells();
669     
670     //Fill histograms to check shape of embedded clusters
671     Float_t fraction = 0;
672     if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
673
674       Float_t clusterE = 0; // recalculate in case corrections applied.
675       Float_t cellE    = 0;
676       for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
677         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
678         clusterE+=cellE;  
679         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
680       }
681       
682       //Fraction of total energy due to the embedded signal
683       fraction/=clusterE;
684       
685       if(GetDebug() > 1 ) printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
686       
687       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
688       
689     }  // embedded fraction    
690     
691     // Get the fraction of the cluster energy that carries the cell with highest energy
692     Int_t absID             =-1 ;
693     Float_t maxCellFraction = 0.;
694     
695     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
696     
697     // Check the origin and fill histograms
698     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) && 
699        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
700        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
701        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
702       fhMCELambda0[mcssPhoton]    ->Fill(energy, lambda0);
703       fhMCELambda1[mcssPhoton]    ->Fill(energy, lambda1);
704       fhMCEDispersion[mcssPhoton] ->Fill(energy, disp);
705       fhMCNCellsE[mcssPhoton]     ->Fill(energy, ncells);
706       fhMCMaxCellDiffClusterE[mcssPhoton]->Fill(energy,maxCellFraction);
707       
708       if     (energy < 2.){
709         fhMCLambda0vsClusterMaxCellDiffE0[mcssPhoton]->Fill(lambda0, maxCellFraction);
710         fhMCNCellsvsClusterMaxCellDiffE0[mcssPhoton] ->Fill(ncells,  maxCellFraction);
711       }
712       else if(energy < 6.){
713         fhMCLambda0vsClusterMaxCellDiffE2[mcssPhoton]->Fill(lambda0, maxCellFraction);
714         fhMCNCellsvsClusterMaxCellDiffE2[mcssPhoton] ->Fill(ncells,  maxCellFraction);
715       }
716       else{
717         fhMCLambda0vsClusterMaxCellDiffE6[mcssPhoton]->Fill(lambda0, maxCellFraction);
718         fhMCNCellsvsClusterMaxCellDiffE6[mcssPhoton] ->Fill(ncells,  maxCellFraction);
719       }
720       
721       if(!GetReader()->IsEmbeddedClusterSelectionOn()){
722         //Check particle overlaps in cluster
723         
724         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
725         Int_t ancPDG = 0, ancStatus = -1;
726         TLorentzVector momentum; TVector3 prodVertex;
727         Int_t ancLabel = 0;
728         Int_t noverlaps = 1;      
729         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
730           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
731           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
732         }
733         
734         if(noverlaps == 1){
735           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
736         }
737         else if(noverlaps == 2){        
738           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
739         }
740         else if(noverlaps > 2){          
741           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
742         }
743         else {
744           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
745         }
746       }//No embedding
747       
748       //Fill histograms to check shape of embedded clusters
749       if(GetReader()->IsEmbeddedClusterSelectionOn()){
750         
751         if     (fraction > 0.9) 
752         {
753           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
754         }
755         else if(fraction > 0.5)
756         {
757           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
758         }
759         else if(fraction > 0.1)
760         { 
761           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
762         }
763         else
764         {
765           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
766         }
767       } // embedded
768       
769     }//photon   no conversion
770     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
771       fhMCELambda0[mcssElectron]    ->Fill(energy, lambda0);
772       fhMCELambda1[mcssElectron]    ->Fill(energy, lambda1);
773       fhMCEDispersion[mcssElectron] ->Fill(energy, disp);
774       fhMCNCellsE[mcssElectron]     ->Fill(energy, ncells);
775       fhMCMaxCellDiffClusterE[mcssElectron]->Fill(energy,maxCellFraction);
776
777       if     (energy < 2.){
778         fhMCLambda0vsClusterMaxCellDiffE0[mcssElectron]->Fill(lambda0, maxCellFraction);
779         fhMCNCellsvsClusterMaxCellDiffE0[mcssElectron] ->Fill(ncells,  maxCellFraction);
780       }
781       else if(energy < 6.){
782         fhMCLambda0vsClusterMaxCellDiffE2[mcssElectron]->Fill(lambda0, maxCellFraction);
783         fhMCNCellsvsClusterMaxCellDiffE2[mcssElectron] ->Fill(ncells,  maxCellFraction);
784       }
785       else{
786         fhMCLambda0vsClusterMaxCellDiffE6[mcssElectron]->Fill(lambda0, maxCellFraction);
787         fhMCNCellsvsClusterMaxCellDiffE6[mcssElectron] ->Fill(ncells,  maxCellFraction);
788       }
789     }//electron
790     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) && 
791               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
792       fhMCELambda0[mcssConversion]    ->Fill(energy, lambda0);
793       fhMCELambda1[mcssConversion]    ->Fill(energy, lambda1);
794       fhMCEDispersion[mcssConversion] ->Fill(energy, disp);
795       fhMCNCellsE[mcssConversion]     ->Fill(energy, ncells);
796       fhMCMaxCellDiffClusterE[mcssConversion]->Fill(energy,maxCellFraction);
797
798       if     (energy < 2.){
799         fhMCLambda0vsClusterMaxCellDiffE0[mcssConversion]->Fill(lambda0, maxCellFraction);
800         fhMCNCellsvsClusterMaxCellDiffE0[mcssConversion] ->Fill(ncells,  maxCellFraction);
801       }
802       else if(energy < 6.){
803         fhMCLambda0vsClusterMaxCellDiffE2[mcssConversion]->Fill(lambda0, maxCellFraction);
804         fhMCNCellsvsClusterMaxCellDiffE2[mcssConversion] ->Fill(ncells,  maxCellFraction);
805       }
806       else{
807         fhMCLambda0vsClusterMaxCellDiffE6[mcssConversion]->Fill(lambda0, maxCellFraction);
808         fhMCNCellsvsClusterMaxCellDiffE6[mcssConversion] ->Fill(ncells,  maxCellFraction);
809       }      
810       
811     }//conversion photon
812     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  ){
813       fhMCELambda0[mcssPi0]    ->Fill(energy, lambda0);
814       fhMCELambda1[mcssPi0]    ->Fill(energy, lambda1);
815       fhMCEDispersion[mcssPi0] ->Fill(energy, disp);
816       fhMCNCellsE[mcssPi0]     ->Fill(energy, ncells);
817       fhMCMaxCellDiffClusterE[mcssPi0]->Fill(energy,maxCellFraction);
818
819       if     (energy < 2.){
820         fhMCLambda0vsClusterMaxCellDiffE0[mcssPi0]->Fill(lambda0, maxCellFraction);
821         fhMCNCellsvsClusterMaxCellDiffE0[mcssPi0] ->Fill(ncells,  maxCellFraction);
822       }
823       else if(energy < 6.){
824         fhMCLambda0vsClusterMaxCellDiffE2[mcssPi0]->Fill(lambda0, maxCellFraction);
825         fhMCNCellsvsClusterMaxCellDiffE2[mcssPi0] ->Fill(ncells,  maxCellFraction);
826       }
827       else{
828         fhMCLambda0vsClusterMaxCellDiffE6[mcssPi0]->Fill(lambda0, maxCellFraction);
829         fhMCNCellsvsClusterMaxCellDiffE6[mcssPi0] ->Fill(ncells,  maxCellFraction);
830       }      
831       
832       //Fill histograms to check shape of embedded clusters
833       if(GetReader()->IsEmbeddedClusterSelectionOn()){
834         
835         if     (fraction > 0.9) 
836         {
837           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
838         }
839         else if(fraction > 0.5)
840         {
841           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
842         }
843         else if(fraction > 0.1)
844         { 
845           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
846         }
847         else
848         {
849           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
850         }
851       } // embedded      
852       
853     }//pi0
854     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ){
855       fhMCELambda0[mcssEta]    ->Fill(energy, lambda0);
856       fhMCELambda1[mcssEta]    ->Fill(energy, lambda1);
857       fhMCEDispersion[mcssEta] ->Fill(energy, disp);
858       fhMCNCellsE[mcssEta]     ->Fill(energy, ncells);
859       fhMCMaxCellDiffClusterE[mcssEta]->Fill(energy,maxCellFraction);
860
861       if     (energy < 2.){
862         fhMCLambda0vsClusterMaxCellDiffE0[mcssEta]->Fill(lambda0, maxCellFraction);
863         fhMCNCellsvsClusterMaxCellDiffE0[mcssEta] ->Fill(ncells,  maxCellFraction);
864       }
865       else if(energy < 6.){
866         fhMCLambda0vsClusterMaxCellDiffE2[mcssEta]->Fill(lambda0, maxCellFraction);
867         fhMCNCellsvsClusterMaxCellDiffE2[mcssEta] ->Fill(ncells,  maxCellFraction);
868       }
869       else{
870         fhMCLambda0vsClusterMaxCellDiffE6[mcssEta]->Fill(lambda0, maxCellFraction);
871         fhMCNCellsvsClusterMaxCellDiffE6[mcssEta] ->Fill(ncells,  maxCellFraction);
872       }
873       
874     }//eta    
875     else {
876       fhMCELambda0[mcssOther]    ->Fill(energy, lambda0);
877       fhMCELambda1[mcssOther]    ->Fill(energy, lambda1);
878       fhMCEDispersion[mcssOther] ->Fill(energy, disp);
879       fhMCNCellsE[mcssOther]     ->Fill(energy, ncells);
880       fhMCMaxCellDiffClusterE[mcssOther]->Fill(energy,maxCellFraction);
881
882       if     (energy < 2.){
883         fhMCLambda0vsClusterMaxCellDiffE0[mcssOther]->Fill(lambda0, maxCellFraction);
884         fhMCNCellsvsClusterMaxCellDiffE0[mcssOther] ->Fill(ncells,  maxCellFraction);
885       }
886       else if(energy < 6.){
887         fhMCLambda0vsClusterMaxCellDiffE2[mcssOther]->Fill(lambda0, maxCellFraction);
888         fhMCNCellsvsClusterMaxCellDiffE2[mcssOther] ->Fill(ncells,  maxCellFraction);
889       }
890       else{
891         fhMCLambda0vsClusterMaxCellDiffE6[mcssOther]->Fill(lambda0, maxCellFraction);
892         fhMCNCellsvsClusterMaxCellDiffE6[mcssOther] ->Fill(ncells,  maxCellFraction);
893       }            
894       
895     }//other particles 
896     
897   }//MC data
898   
899 }
900
901 //________________________________________________________________________
902 TObjString *  AliAnaPhoton::GetAnalysisCuts()
903 {       
904   //Save parameters used for analysis
905   TString parList ; //this will be list of parameters used for this analysis.
906   const Int_t buffersize = 255;
907   char onePar[buffersize] ;
908   
909   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
910   parList+=onePar ;     
911   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
912   parList+=onePar ;
913   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
914   parList+=onePar ;
915   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
916   parList+=onePar ;
917   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
918   parList+=onePar ;
919   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
920   parList+=onePar ;  
921   
922   //Get parameters set in base class.
923   parList += GetBaseParametersList() ;
924   
925   //Get parameters set in PID class.
926   parList += GetCaloPID()->GetPIDParametersList() ;
927   
928   //Get parameters set in FiducialCut class (not available yet)
929   //parlist += GetFidCut()->GetFidCutParametersList() 
930   
931   return new TObjString(parList) ;
932 }
933
934 //________________________________________________________________________
935 TList *  AliAnaPhoton::GetCreateOutputObjects()
936 {  
937   // Create histograms to be saved in output file and 
938   // store them in outputContainer
939   TList * outputContainer = new TList() ; 
940   outputContainer->SetName("PhotonHistos") ; 
941         
942   Int_t nptbins  = GetHistoPtBins();  Float_t ptmax  = GetHistoPtMax();  Float_t ptmin  = GetHistoPtMin(); 
943   Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin(); 
944   Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();     
945   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax = GetHistoShowerShapeMax();  Float_t ssmin = GetHistoShowerShapeMin();
946   Int_t nbins    = GetHistoNClusterCellBins(); Int_t nmax    = GetHistoNClusterCellMax(); Int_t nmin    = GetHistoNClusterCellMin(); 
947
948   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax); 
949   fhNCellsE->SetXTitle("E (GeV)");
950   fhNCellsE->SetYTitle("# of cells in cluster");
951   outputContainer->Add(fhNCellsE);  
952   
953   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
954                                     nptbins,ptmin,ptmax, 500,0,1.); 
955   fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
956   fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
957   outputContainer->Add(fhMaxCellDiffClusterE);  
958   
959   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax); 
960   fhEPhoton->SetYTitle("N");
961   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
962   outputContainer->Add(fhEPhoton) ;   
963   
964   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax); 
965   fhPtPhoton->SetYTitle("N");
966   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
967   outputContainer->Add(fhPtPhoton) ; 
968   
969   fhPhiPhoton  = new TH2F
970     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
971   fhPhiPhoton->SetYTitle("#phi (rad)");
972   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
973   outputContainer->Add(fhPhiPhoton) ; 
974   
975   fhEtaPhoton  = new TH2F
976     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
977   fhEtaPhoton->SetYTitle("#eta");
978   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
979   outputContainer->Add(fhEtaPhoton) ;
980   
981   fhEtaPhiPhoton  = new TH2F
982   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
983   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
984   fhEtaPhiPhoton->SetXTitle("#eta");
985   outputContainer->Add(fhEtaPhiPhoton) ;
986   if(GetMinPt() < 0.5){
987     fhEtaPhi05Photon  = new TH2F
988     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax); 
989     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
990     fhEtaPhi05Photon->SetXTitle("#eta");
991     outputContainer->Add(fhEtaPhi05Photon) ;
992   }
993   
994   //Shower shape
995   if(fFillSSHistograms){
996     
997     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
998     fhLam0E->SetYTitle("#lambda_{0}^{2}");
999     fhLam0E->SetXTitle("E (GeV)");
1000     outputContainer->Add(fhLam0E);  
1001     
1002     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1003     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1004     fhLam1E->SetXTitle("E (GeV)");
1005     outputContainer->Add(fhLam1E);  
1006     
1007     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1008     fhDispE->SetYTitle("D^{2}");
1009     fhDispE->SetXTitle("E (GeV) ");
1010     outputContainer->Add(fhDispE);
1011     
1012     if(fCalorimeter == "EMCAL"){
1013       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1014       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1015       fhLam0ETRD->SetXTitle("E (GeV)");
1016       outputContainer->Add(fhLam0ETRD);  
1017       
1018       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1019       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1020       fhLam1ETRD->SetXTitle("E (GeV)");
1021       outputContainer->Add(fhLam1ETRD);  
1022       
1023       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1024       fhDispETRD->SetYTitle("Dispersion^{2}");
1025       fhDispETRD->SetXTitle("E (GeV) ");
1026       outputContainer->Add(fhDispETRD);   
1027     } 
1028     
1029     fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1030     fhNCellsLam0LowE->SetXTitle("N_{cells}");
1031     fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1032     outputContainer->Add(fhNCellsLam0LowE);  
1033     
1034     fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1035     fhNCellsLam0HighE->SetXTitle("N_{cells}");
1036     fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1037     outputContainer->Add(fhNCellsLam0HighE);  
1038     
1039     fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1040     fhNCellsLam1LowE->SetXTitle("N_{cells}");
1041     fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1042     outputContainer->Add(fhNCellsLam1LowE);  
1043     
1044     fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1045     fhNCellsLam1HighE->SetXTitle("N_{cells}");
1046     fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1047     outputContainer->Add(fhNCellsLam1HighE);  
1048     
1049     fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1050     fhNCellsDispLowE->SetXTitle("N_{cells}");
1051     fhNCellsDispLowE->SetYTitle("D^{2}");
1052     outputContainer->Add(fhNCellsDispLowE);  
1053     
1054     fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax); 
1055     fhNCellsDispHighE->SetXTitle("N_{cells}");
1056     fhNCellsDispHighE->SetYTitle("D^{2}");
1057     outputContainer->Add(fhNCellsDispHighE);  
1058     
1059     fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax); 
1060     fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1061     fhEtaLam0LowE->SetXTitle("#eta");
1062     outputContainer->Add(fhEtaLam0LowE);  
1063     
1064     fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
1065     fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1066     fhPhiLam0LowE->SetXTitle("#phi");
1067     outputContainer->Add(fhPhiLam0LowE);  
1068     
1069     fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax); 
1070     fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1071     fhEtaLam0HighE->SetXTitle("#eta");
1072     outputContainer->Add(fhEtaLam0HighE);  
1073     
1074     fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
1075     fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1076     fhPhiLam0HighE->SetXTitle("#phi");
1077     outputContainer->Add(fhPhiLam0HighE);  
1078     
1079     fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1080     fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1081     fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1082     outputContainer->Add(fhLam1Lam0LowE);  
1083     
1084     fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1085     fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1086     fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1087     outputContainer->Add(fhLam1Lam0HighE);  
1088     
1089     fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1090     fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1091     fhLam0DispLowE->SetYTitle("D^{2}");
1092     outputContainer->Add(fhLam0DispLowE);  
1093     
1094     fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1095     fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1096     fhLam0DispHighE->SetYTitle("D^{2}");
1097     outputContainer->Add(fhLam0DispHighE);  
1098     
1099     fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1100     fhDispLam1LowE->SetXTitle("D^{2}");
1101     fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1102     outputContainer->Add(fhDispLam1LowE);  
1103     
1104     fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax); 
1105     fhDispLam1HighE->SetXTitle("D^{2}");
1106     fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1107     outputContainer->Add(fhDispLam1HighE);  
1108     
1109   } // Shower shape
1110   
1111   
1112   if(IsDataMC()){
1113    
1114     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
1115                         "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
1116                         "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ; 
1117     
1118     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
1119                         "Conversion", "Hadron", "AntiNeutron","AntiProton",
1120                         "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
1121     
1122     for(Int_t i = 0; i < fNOriginHistograms; i++){ 
1123       
1124       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
1125                                 Form("cluster from %s : E ",ptype[i].Data()),
1126                                 nptbins,ptmin,ptmax); 
1127       fhMCE[i]->SetXTitle("E (GeV)");
1128       outputContainer->Add(fhMCE[i]) ; 
1129       
1130       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
1131                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
1132                            nptbins,ptmin,ptmax); 
1133       fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1134       outputContainer->Add(fhMCPt[i]) ;
1135       
1136       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
1137                            Form("cluster from %s : #eta ",ptype[i].Data()),
1138                            nptbins,ptmin,ptmax,netabins,etamin,etamax); 
1139       fhMCEta[i]->SetYTitle("#eta");
1140       fhMCEta[i]->SetXTitle("E (GeV)");
1141       outputContainer->Add(fhMCEta[i]) ;
1142       
1143       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
1144                            Form("cluster from %s : #phi ",ptype[i].Data()),
1145                            nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
1146       fhMCPhi[i]->SetYTitle("#phi (rad)");
1147       fhMCPhi[i]->SetXTitle("E (GeV)");
1148       outputContainer->Add(fhMCPhi[i]) ;
1149       
1150       
1151       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
1152                                  Form("MC - Reco E from %s",pname[i].Data()), 
1153                                  nptbins,ptmin,ptmax, 200,-50,50); 
1154       fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
1155       outputContainer->Add(fhMCDeltaE[i]);
1156       
1157       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
1158                                   Form("MC - Reco p_{T} from %s",pname[i].Data()), 
1159                                   nptbins,ptmin,ptmax, 200,-50,50); 
1160       fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
1161       outputContainer->Add(fhMCDeltaPt[i]);
1162             
1163       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
1164                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()), 
1165                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
1166       fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
1167       fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
1168       outputContainer->Add(fhMC2E[i]);          
1169       
1170       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
1171                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()), 
1172                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
1173       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
1174       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
1175       outputContainer->Add(fhMC2Pt[i]);
1176       
1177       
1178     }
1179     
1180     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
1181                          "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ; 
1182     
1183     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
1184                          "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
1185     
1186     for(Int_t i = 0; i < fNPrimaryHistograms; i++){ 
1187       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
1188                            Form("primary photon %s : E ",pptype[i].Data()),
1189                            nptbins,ptmin,ptmax); 
1190       fhEPrimMC[i]->SetXTitle("E (GeV)");
1191       outputContainer->Add(fhEPrimMC[i]) ; 
1192       
1193       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
1194                             Form("primary photon %s : p_{T} ",pptype[i].Data()),
1195                             nptbins,ptmin,ptmax); 
1196       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1197       outputContainer->Add(fhPtPrimMC[i]) ;
1198       
1199       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
1200                              Form("primary photon %s : Rapidity ",pptype[i].Data()),
1201                              nptbins,ptmin,ptmax,800,-8,8); 
1202       fhYPrimMC[i]->SetYTitle("Rapidity");
1203       fhYPrimMC[i]->SetXTitle("E (GeV)");
1204       outputContainer->Add(fhYPrimMC[i]) ;
1205       
1206       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
1207                              Form("primary photon %s : #phi ",pptype[i].Data()),
1208                              nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
1209       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1210       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1211       outputContainer->Add(fhPhiPrimMC[i]) ;
1212      
1213       
1214       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
1215                                Form("primary photon %s in acceptance: E ",pptype[i].Data()),
1216                                nptbins,ptmin,ptmax); 
1217       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1218       outputContainer->Add(fhEPrimMCAcc[i]) ; 
1219       
1220       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
1221                                 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
1222                                 nptbins,ptmin,ptmax); 
1223       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1224       outputContainer->Add(fhPtPrimMCAcc[i]) ;
1225       
1226       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
1227                                  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
1228                                  nptbins,ptmin,ptmax,100,-1,1); 
1229       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1230       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1231       outputContainer->Add(fhYPrimMCAcc[i]) ;
1232       
1233       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
1234                                  Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
1235                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
1236       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1237       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1238       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1239       
1240     }
1241                 
1242     if(fFillSSHistograms){
1243       
1244       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ; 
1245       
1246       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1247       
1248       for(Int_t i = 0; i < 6; i++){ 
1249         
1250         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1251                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
1252                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1253         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1254         fhMCELambda0[i]->SetXTitle("E (GeV)");
1255         outputContainer->Add(fhMCELambda0[i]) ; 
1256         
1257         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1258                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
1259                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1260         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1261         fhMCELambda1[i]->SetXTitle("E (GeV)");
1262         outputContainer->Add(fhMCELambda1[i]) ; 
1263               
1264         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1265                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
1266                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1267         fhMCEDispersion[i]->SetYTitle("D^{2}");
1268         fhMCEDispersion[i]->SetXTitle("E (GeV)");
1269         outputContainer->Add(fhMCEDispersion[i]) ; 
1270                 
1271         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
1272                                Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()), 
1273                                     nptbins,ptmin,ptmax, nbins,nmin,nmax); 
1274         fhMCNCellsE[i]->SetXTitle("E (GeV)");
1275         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
1276         outputContainer->Add(fhMCNCellsE[i]);  
1277         
1278         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
1279                                              Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
1280                                            nptbins,ptmin,ptmax, 500,0,1.); 
1281         fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
1282         fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1283         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);  
1284         
1285         fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1286                                     Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1287                                     ssbins,ssmin,ssmax,500,0,1.); 
1288         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
1289         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1290         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ; 
1291         
1292         fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1293                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1294                                                          ssbins,ssmin,ssmax,500,0,1.); 
1295         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
1296         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1297         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ; 
1298         
1299         fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1300                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1301                                                          ssbins,ssmin,ssmax,500,0,1.); 
1302         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
1303         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1304         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ; 
1305         
1306         fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1307                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1308                                                          nbins/5,nmin,nmax/5,500,0,1.); 
1309         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
1310         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1311         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ; 
1312         
1313         fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1314                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1315                                                          nbins/5,nmin,nmax/5,500,0,1.); 
1316         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
1317         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1318         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ; 
1319         
1320         fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1321                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1322                                                          nbins/5,nmin,nmax/5,500,0,1.); 
1323         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
1324         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
1325         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ; 
1326         
1327        }// loop    
1328       
1329       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1330       {
1331         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
1332                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
1333                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1334         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1335         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1336         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ; 
1337         
1338         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1339                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
1340                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1341         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1342         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1343         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ; 
1344         
1345         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
1346                                                "cluster from Photon : E vs #lambda_{0}^{2}",
1347                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1348         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1349         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1350         outputContainer->Add(fhMCPhotonELambda0NOverlap) ; 
1351         
1352       } //No embedding     
1353       
1354       //Fill histograms to check shape of embedded clusters
1355       if(GetReader()->IsEmbeddedClusterSelectionOn())
1356       {
1357         
1358         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
1359                                                     "Energy Fraction of embedded signal versus cluster energy",
1360                                                     nptbins,ptmin,ptmax,100,0.,1.); 
1361         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1362         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1363         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ; 
1364         
1365         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
1366                                                 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1367                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1368         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1369         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1370         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ; 
1371                 
1372         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
1373                                           "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1374                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1375         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1376         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1377         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ; 
1378         
1379         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
1380                                           "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1381                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1382         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1383         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1384         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ; 
1385                 
1386         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
1387                                                    "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1388                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1389         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1390         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
1391         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ; 
1392         
1393         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
1394                                                     "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1395                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1396         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1397         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
1398         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ; 
1399                 
1400         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
1401                                                       "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1402                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1403         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1404         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
1405         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ; 
1406         
1407         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
1408                                                    "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1409                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1410         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1411         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
1412         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ; 
1413         
1414         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
1415                                                  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1416                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1417         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1418         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
1419         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ; 
1420   
1421       }// embedded histograms
1422       
1423       
1424     }// Fill SS MC histograms
1425     
1426   }//Histos with MC
1427     
1428   //Store calo PID histograms
1429   if(fRejectTrackMatch && IsCaloPIDOn()){
1430     TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
1431     for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
1432       outputContainer->Add(caloPIDHistos->At(i)) ;
1433     }
1434     delete caloPIDHistos;
1435   }
1436   
1437   return outputContainer ;
1438   
1439 }
1440
1441 //____________________________________________________________________________
1442 void AliAnaPhoton::Init()
1443 {
1444   
1445   //Init
1446   //Do some checks
1447   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1448     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1449     abort();
1450   }
1451   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1452     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1453     abort();
1454   }
1455   
1456 }
1457
1458 //____________________________________________________________________________
1459 void AliAnaPhoton::InitParameters()
1460 {
1461   
1462   //Initialize the parameters of the analysis.
1463   AddToHistogramsName("AnaPhoton_");
1464   
1465   fCalorimeter = "EMCAL" ;
1466   fMinDist     = 2.;
1467   fMinDist2    = 4.;
1468   fMinDist3    = 5.;
1469         
1470   fTimeCutMin  = -1;
1471   fTimeCutMax  = 9999999;
1472   fNCellsCut   = 0;
1473         
1474   fRejectTrackMatch       = kTRUE ;
1475         
1476 }
1477
1478 //__________________________________________________________________
1479 void  AliAnaPhoton::MakeAnalysisFillAOD() 
1480 {
1481   //Do photon analysis and fill aods
1482   
1483   //Get the vertex 
1484   Double_t v[3] = {0,0,0}; //vertex ;
1485   GetReader()->GetVertex(v);
1486   
1487   //Select the Calorimeter of the photon
1488   TObjArray * pl = 0x0; 
1489   if(fCalorimeter == "PHOS")
1490     pl = GetPHOSClusters();
1491   else if (fCalorimeter == "EMCAL")
1492     pl = GetEMCALClusters();
1493   
1494   if(!pl) {
1495     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1496     return;
1497   }
1498   
1499   //Init arrays, variables, get number of clusters
1500   TLorentzVector mom, mom2 ;
1501   Int_t nCaloClusters = pl->GetEntriesFast();
1502   
1503   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
1504   
1505   //----------------------------------------------------
1506   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1507   //----------------------------------------------------
1508   // Loop on clusters
1509   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){    
1510           
1511           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo)); 
1512     //printf("calo %d, %f\n",icalo,calo->E());
1513     
1514     //Get the index where the cluster comes, to retrieve the corresponding vertex
1515     Int_t evtIndex = 0 ; 
1516     if (GetMixedEvent()) {
1517       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
1518       //Get the vertex and check it is not too large in z
1519       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1520     }
1521     
1522     //Cluster selection, not charged, with photon id and in fiducial cut          
1523     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1524       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1525     else{
1526       Double_t vertex[]={0,0,0};
1527       calo->GetMomentum(mom,vertex) ;
1528     }
1529     
1530     //--------------------------------------
1531     // Cluster selection
1532     //--------------------------------------
1533     if(!ClusterSelected(calo,mom)) continue;
1534     
1535     //----------------------------
1536     //Create AOD for analysis
1537     //----------------------------
1538     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
1539     
1540     //...............................................
1541     //Set the indeces of the original caloclusters (MC, ID), and calorimeter  
1542     Int_t label = calo->GetLabel();
1543     aodph.SetLabel(label);
1544     aodph.SetCaloLabel(calo->GetID(),-1);
1545     aodph.SetDetector(fCalorimeter);
1546     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1547     
1548     //...............................................
1549     //Set bad channel distance bit
1550     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1551     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
1552     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
1553     else                         aodph.SetDistToBad(0) ;
1554     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
1555     
1556     //--------------------------------------------------------------------------------------
1557     //Play with the MC stack if available
1558     //--------------------------------------------------------------------------------------
1559     
1560     //Check origin of the candidates
1561     if(IsDataMC()){
1562       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
1563
1564       if(GetDebug() > 0)
1565         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
1566     }//Work with stack also   
1567     
1568     //--------------------------------------------------------------------------------------
1569     //Fill some shower shape histograms before PID is applied
1570     //--------------------------------------------------------------------------------------
1571     
1572     FillShowerShapeHistograms(calo,aodph.GetTag());
1573     
1574     //-------------------------------------
1575     //PID selection or bit setting
1576     //-------------------------------------
1577     // MC
1578     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
1579       //Get most probable PID, check PID weights (in MC this option is mandatory)
1580       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
1581       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());      
1582       //If primary is not photon, skip it.
1583       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
1584     }   
1585     //...............................................
1586     // Data, PID check on
1587     else if(IsCaloPIDOn()){
1588       //Get most probable PID, 2 options check PID weights 
1589       //or redo PID, recommended option for MCEal.              
1590       if(!IsCaloPIDRecalculationOn())
1591         aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
1592       else
1593         aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
1594       
1595       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
1596       
1597       //If cluster does not pass pid, not photon, skip it.
1598       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;                   
1599       
1600     }
1601     //...............................................
1602     // Data, PID check off
1603     else{
1604       //Set PID bits for later selection (AliAnaPi0 for example)
1605       //GetPDG already called in SetPIDBits.
1606       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
1607       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");               
1608     }
1609     
1610     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
1611         
1612     //Add AOD with photon object to aod branch
1613     AddAODParticle(aodph);
1614     
1615   }//loop
1616         
1617   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());  
1618   
1619 }
1620
1621 //__________________________________________________________________
1622 void  AliAnaPhoton::MakeAnalysisFillHistograms() 
1623 {
1624   //Fill histograms
1625   
1626   //-------------------------------------------------------------------
1627   // Access MC information in stack if requested, check that it exists. 
1628   AliStack         * stack       = 0x0;
1629   TParticle        * primary     = 0x0;   
1630   TClonesArray     * mcparticles = 0x0;
1631   AliAODMCParticle * aodprimary  = 0x0; 
1632   
1633   if(IsDataMC()){
1634     
1635     if(GetReader()->ReadStack()){
1636       stack =  GetMCStack() ;
1637       if(!stack) {
1638         printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1639         abort();
1640       }
1641       
1642     }
1643     else if(GetReader()->ReadAODMCParticles()){
1644       
1645       //Get the list of MC particles
1646       mcparticles = GetReader()->GetAODMCParticles(0);
1647       if(!mcparticles && GetDebug() > 0)        {
1648         printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
1649       } 
1650     }
1651   }// is data and MC
1652   
1653   
1654   // Get vertex
1655   Double_t v[3] = {0,0,0}; //vertex ;
1656   GetReader()->GetVertex(v);
1657   //fhVertex->Fill(v[0],v[1],v[2]);  
1658   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1659   
1660   //----------------------------------
1661   //Loop on stored AOD photons
1662   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1663   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1664   
1665   for(Int_t iaod = 0; iaod < naod ; iaod++){
1666     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1667     Int_t pdg = ph->GetIdentifiedParticleType();
1668     
1669     if(GetDebug() > 3) 
1670       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1671     
1672     //If PID used, fill histos with photons in Calorimeter fCalorimeter
1673     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue; 
1674     if(ph->GetDetector() != fCalorimeter) continue;
1675     
1676     if(GetDebug() > 2) 
1677       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1678     
1679     //................................
1680     //Fill photon histograms 
1681     Float_t ptcluster  = ph->Pt();
1682     Float_t phicluster = ph->Phi();
1683     Float_t etacluster = ph->Eta();
1684     Float_t ecluster   = ph->E();
1685     
1686     fhEPhoton   ->Fill(ecluster);
1687     fhPtPhoton  ->Fill(ptcluster);
1688     fhPhiPhoton ->Fill(ptcluster,phicluster);
1689     fhEtaPhoton ->Fill(ptcluster,etacluster);    
1690     if     (ecluster > 0.5)   fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
1691     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1692     
1693
1694     //Get original cluster, to recover some information
1695     Int_t absID             = 0; 
1696     Float_t maxCellFraction = 0;
1697     AliVCaloCells* cells    = 0;
1698     TObjArray * clusters    = 0; 
1699     if(fCalorimeter == "EMCAL"){
1700       cells    = GetEMCALCells();
1701       clusters = GetEMCALClusters();
1702     }
1703     else{
1704       cells    = GetPHOSCells();
1705       clusters = GetPHOSClusters();
1706     }
1707     
1708     Int_t iclus = -1;
1709     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus); 
1710     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
1711     
1712     fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
1713     fhNCellsE            ->Fill(ph->E(),cluster->GetNCells());
1714     
1715     //.......................................
1716     //Play with the MC data if available
1717     if(IsDataMC()){
1718       
1719       FillAcceptanceHistograms();
1720       
1721       //....................................................................
1722       // Access MC information in stack if requested, check that it exists.
1723       Int_t label =ph->GetLabel();
1724       if(label < 0) {
1725         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1726         continue;
1727       }
1728       
1729       Float_t eprim   = 0;
1730       Float_t ptprim  = 0;
1731       if(GetReader()->ReadStack()){
1732         
1733         if(label >=  stack->GetNtrack()) {
1734           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1735           continue ;
1736         }
1737         
1738         primary = stack->Particle(label);
1739         if(!primary){
1740           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1741           continue;
1742         }
1743         eprim   = primary->Energy();
1744         ptprim  = primary->Pt();                
1745         
1746       }
1747       else if(GetReader()->ReadAODMCParticles()){
1748         //Check which is the input
1749         if(ph->GetInputFileIndex() == 0){
1750           if(!mcparticles) continue;
1751           if(label >=  mcparticles->GetEntriesFast()) {
1752             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", 
1753                                        label, mcparticles->GetEntriesFast());
1754             continue ;
1755           }
1756           //Get the particle
1757           aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1758           
1759         }
1760         
1761         if(!aodprimary){
1762           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1763           continue;
1764         }
1765         
1766         eprim   = aodprimary->E();
1767         ptprim  = aodprimary->Pt();
1768         
1769       }
1770       
1771       Int_t tag =ph->GetTag();
1772       
1773       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[mcPhoton])
1774       {
1775         fhMCE  [mcPhoton] ->Fill(ecluster);
1776         fhMCPt [mcPhoton] ->Fill(ptcluster);
1777         fhMCPhi[mcPhoton] ->Fill(ecluster,phicluster);
1778         fhMCEta[mcPhoton] ->Fill(ecluster,etacluster);
1779         
1780         fhMC2E[mcPhoton]     ->Fill(ecluster, eprim);
1781         fhMC2Pt[mcPhoton]    ->Fill(ptcluster, ptprim);     
1782         fhMCDeltaE[mcPhoton] ->Fill(ecluster,eprim-ecluster);
1783         fhMCDeltaPt[mcPhoton]->Fill(ptcluster,ptprim-ptcluster);     
1784         
1785         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[mcConversion])
1786         {
1787           fhMCE  [mcConversion] ->Fill(ecluster);
1788           fhMCPt [mcConversion] ->Fill(ptcluster);
1789           fhMCPhi[mcConversion] ->Fill(ecluster,phicluster);
1790           fhMCEta[mcConversion] ->Fill(ecluster,etacluster);
1791           
1792           fhMC2E[mcConversion]     ->Fill(ecluster, eprim);
1793           fhMC2Pt[mcConversion]    ->Fill(ptcluster, ptprim);     
1794           fhMCDeltaE[mcConversion] ->Fill(ecluster,eprim-ecluster);
1795           fhMCDeltaPt[mcConversion]->Fill(ptcluster,ptprim-ptcluster);     
1796         }                       
1797         
1798         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[mcPrompt]){
1799           fhMCE  [mcPrompt] ->Fill(ecluster);
1800           fhMCPt [mcPrompt] ->Fill(ptcluster);
1801           fhMCPhi[mcPrompt] ->Fill(ecluster,phicluster);
1802           fhMCEta[mcPrompt] ->Fill(ecluster,etacluster);      
1803           
1804           fhMC2E[mcPrompt]     ->Fill(ecluster, eprim);
1805           fhMC2Pt[mcPrompt]    ->Fill(ptcluster, ptprim);     
1806           fhMCDeltaE[mcPrompt] ->Fill(ecluster,eprim-ecluster);
1807           fhMCDeltaPt[mcPrompt]->Fill(ptcluster,ptprim-ptcluster);     
1808           
1809         }
1810         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[mcFragmentation])
1811         {
1812           fhMCE  [mcFragmentation] ->Fill(ecluster);
1813           fhMCPt [mcFragmentation] ->Fill(ptcluster);
1814           fhMCPhi[mcFragmentation] ->Fill(ecluster,phicluster);
1815           fhMCEta[mcFragmentation] ->Fill(ecluster,etacluster);
1816           
1817           fhMC2E[mcFragmentation]     ->Fill(ecluster, eprim);
1818           fhMC2Pt[mcFragmentation]    ->Fill(ptcluster, ptprim);     
1819           fhMCDeltaE[mcFragmentation] ->Fill(ecluster,eprim-ecluster);
1820           fhMCDeltaPt[mcFragmentation]->Fill(ptcluster,ptprim-ptcluster);     
1821           
1822         }
1823         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[mcISR])
1824         {
1825           fhMCE  [mcISR] ->Fill(ecluster);
1826           fhMCPt [mcISR] ->Fill(ptcluster);
1827           fhMCPhi[mcISR] ->Fill(ecluster,phicluster);
1828           fhMCEta[mcISR] ->Fill(ecluster,etacluster);    
1829           
1830           fhMC2E[mcISR]     ->Fill(ecluster, eprim);
1831           fhMC2Pt[mcISR]    ->Fill(ptcluster, ptprim);     
1832           fhMCDeltaE[mcISR] ->Fill(ecluster,eprim-ecluster);
1833           fhMCDeltaPt[mcISR]->Fill(ptcluster,ptprim-ptcluster);     
1834           
1835         }
1836         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) && 
1837                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[mcPi0Decay])
1838         {
1839           fhMCE  [mcPi0Decay] ->Fill(ecluster);
1840           fhMCPt [mcPi0Decay] ->Fill(ptcluster);
1841           fhMCPhi[mcPi0Decay] ->Fill(ecluster,phicluster);
1842           fhMCEta[mcPi0Decay] ->Fill(ecluster,etacluster);
1843           
1844           fhMC2E[mcPi0Decay]     ->Fill(ecluster, eprim);
1845           fhMC2Pt[mcPi0Decay]    ->Fill(ptcluster, ptprim);     
1846           fhMCDeltaE[mcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1847           fhMCDeltaPt[mcPi0Decay]->Fill(ptcluster,ptprim-ptcluster);     
1848           
1849         }
1850         else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || 
1851                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[mcOtherDecay])
1852         {
1853           fhMCE  [mcOtherDecay] ->Fill(ecluster);
1854           fhMCPt [mcOtherDecay] ->Fill(ptcluster);
1855           fhMCPhi[mcOtherDecay] ->Fill(ecluster,phicluster);
1856           fhMCEta[mcOtherDecay] ->Fill(ecluster,etacluster);
1857           
1858           fhMC2E[mcOtherDecay]     ->Fill(ecluster, eprim);
1859           fhMC2Pt[mcOtherDecay]    ->Fill(ptcluster, ptprim);     
1860           fhMCDeltaE[mcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1861           fhMCDeltaPt[mcOtherDecay]->Fill(ptcluster,ptprim-ptcluster);     
1862           
1863         }
1864         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [mcPi0])
1865         {
1866           fhMCE  [mcPi0] ->Fill(ecluster);
1867           fhMCPt [mcPi0] ->Fill(ptcluster);
1868           fhMCPhi[mcPi0] ->Fill(ecluster,phicluster);
1869           fhMCEta[mcPi0] ->Fill(ecluster,etacluster);
1870           
1871           fhMC2E[mcPi0]     ->Fill(ecluster, eprim);
1872           fhMC2Pt[mcPi0]    ->Fill(ptcluster, ptprim);     
1873           fhMCDeltaE[mcPi0] ->Fill(ecluster,eprim-ecluster);
1874           fhMCDeltaPt[mcPi0]->Fill(ptcluster,ptprim-ptcluster);     
1875           
1876         } 
1877         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[mcEta])
1878         {
1879           fhMCE  [mcEta] ->Fill(ecluster);
1880           fhMCPt [mcEta] ->Fill(ptcluster);
1881           fhMCPhi[mcEta] ->Fill(ecluster,phicluster);
1882           fhMCEta[mcEta] ->Fill(ecluster,etacluster);
1883           
1884           fhMC2E[mcEta]     ->Fill(ecluster, eprim);
1885           fhMC2Pt[mcEta]    ->Fill(ptcluster, ptprim);     
1886           fhMCDeltaE[mcEta] ->Fill(ecluster,eprim-ecluster);
1887           fhMCDeltaPt[mcEta]->Fill(ptcluster,ptprim-ptcluster);     
1888           
1889         }      
1890       }
1891       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[mcAntiNeutron])
1892       {
1893         fhMCE  [mcAntiNeutron] ->Fill(ecluster);
1894         fhMCPt [mcAntiNeutron] ->Fill(ptcluster);
1895         fhMCPhi[mcAntiNeutron] ->Fill(ecluster,phicluster);
1896         fhMCEta[mcAntiNeutron] ->Fill(ecluster,etacluster);
1897         
1898         fhMC2E[mcAntiNeutron]     ->Fill(ecluster, eprim);
1899         fhMC2Pt[mcAntiNeutron]    ->Fill(ptcluster, ptprim);     
1900         fhMCDeltaE[mcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1901         fhMCDeltaPt[mcAntiNeutron]->Fill(ptcluster,ptprim-ptcluster);     
1902         
1903       }
1904       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[mcAntiProton])
1905       {
1906         fhMCE  [mcAntiProton] ->Fill(ecluster);
1907         fhMCPt [mcAntiProton] ->Fill(ptcluster);
1908         fhMCPhi[mcAntiProton] ->Fill(ecluster,phicluster);
1909         fhMCEta[mcAntiProton] ->Fill(ecluster,etacluster);
1910         
1911         fhMC2E[mcAntiProton]     ->Fill(ecluster, eprim);
1912         fhMC2Pt[mcAntiProton]    ->Fill(ptcluster, ptprim);     
1913         fhMCDeltaE[mcAntiProton] ->Fill(ecluster,eprim-ecluster);
1914         fhMCDeltaPt[mcAntiProton]->Fill(ecluster,ptprim-ptcluster);     
1915         
1916       } 
1917       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[mcElectron])
1918       {
1919         fhMCE  [mcElectron] ->Fill(ecluster);
1920         fhMCPt [mcElectron] ->Fill(ptcluster);
1921         fhMCPhi[mcElectron] ->Fill(ecluster,phicluster);
1922         fhMCEta[mcElectron] ->Fill(ecluster,etacluster);
1923         
1924         fhMC2E[mcElectron]     ->Fill(ecluster, eprim);
1925         fhMC2Pt[mcElectron]    ->Fill(ptcluster, ptprim);     
1926         fhMCDeltaE[mcElectron] ->Fill(ecluster,eprim-ecluster);
1927         fhMCDeltaPt[mcElectron]->Fill(ecluster,ptprim-ptcluster);             
1928       }     
1929       else if( fhMCE[mcOther]){
1930         fhMCE  [mcOther] ->Fill(ecluster);
1931         fhMCPt [mcOther] ->Fill(ptcluster);
1932         fhMCPhi[mcOther] ->Fill(ecluster,phicluster);
1933         fhMCEta[mcOther] ->Fill(ecluster,etacluster);
1934         
1935         fhMC2E[mcOther]     ->Fill(ecluster, eprim);
1936         fhMC2Pt[mcOther]    ->Fill(ptcluster, ptprim);     
1937         fhMCDeltaE[mcOther] ->Fill(ecluster,eprim-ecluster);
1938         fhMCDeltaPt[mcOther]->Fill(ecluster,ptprim-ptcluster);     
1939         
1940         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1941         //                                      ph->GetLabel(),ph->Pt());
1942         //                for(Int_t i = 0; i < 20; i++) {
1943         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1944         //                }
1945         //                printf("\n");
1946         
1947       }
1948       
1949     }//Histograms with MC
1950     
1951   }// aod loop
1952   
1953 }
1954
1955
1956 //__________________________________________________________________
1957 void AliAnaPhoton::Print(const Option_t * opt) const
1958 {
1959   //Print some relevant parameters set for the analysis
1960   
1961   if(! opt)
1962     return;
1963   
1964   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1965   AliAnaPartCorrBaseClass::Print(" ");
1966
1967   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1968   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1969   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1970   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1971   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1972   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1973   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
1974   printf("    \n") ;
1975         
1976