Added Pi0FlowMC
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0FlowMC.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // Extension to Pi0FLOw, mimicing AliPHOSHijingEfficiency
17 // by Dmitri Peressounko, 05.02.2013
18 // Authors: Henrik Qvigstad, Dmitri Peressounko
19 // Date   : 05.04.2013
20 /* $Id$ */
21
22 #include "TChain.h"
23 #include "TTree.h"
24 #include "TObjArray.h"
25 #include "TF1.h"
26 #include "TFile.h"
27 #include "TH1F.h"
28 #include "TH2F.h"
29 #include "TH2I.h"
30 #include "TH3F.h"
31 #include "TParticle.h"
32 #include "TCanvas.h"
33 #include "TStyle.h"
34 #include "TRandom.h"
35 #include "TROOT.h"
36 #include "THashList.h"
37
38
39 #include "AliAnalysisManager.h"
40 #include "AliMCEventHandler.h"
41 #include "AliMCEvent.h"
42 #include "AliStack.h"
43 #include "AliAnalysisTaskSE.h"
44 #include "AliPHOSHijingEfficiency.h"
45 #include "AliCaloPhoton.h"
46 #include "AliPHOSGeometry.h"
47 #include "AliPHOSEsdCluster.h"
48 #include "AliPHOSCalibData.h"
49 #include "AliESDEvent.h"
50 #include "AliESDCaloCells.h"
51 #include "AliESDVertex.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliLog.h"
54 #include "AliPID.h"
55 #include "AliCDBManager.h"
56 #include "AliCentrality.h" 
57 #include "AliESDtrackCuts.h"
58 #include "AliEventplane.h"
59 #include "TProfile.h"
60 #include <TPDGCode.h>
61 #include "AliOADBContainer.h"
62
63
64 #include "AliAnalysisTaskPi0FlowMC.h"
65
66 ClassImp(AliAnalysisTaskPi0FlowMC);
67
68 //TODO: rnlin?
69
70 AliAnalysisTaskPi0FlowMC::AliAnalysisTaskPi0FlowMC(const char* name, AliAnalysisTaskPi0Flow::Period period)
71 : AliAnalysisTaskPi0Flow(name, period),
72   fStack(0x0)
73 {
74 }
75
76 AliAnalysisTaskPi0FlowMC::~AliAnalysisTaskPi0FlowMC()
77 {
78 }
79
80 void AliAnalysisTaskPi0FlowMC::MakeMCHistograms()
81 {
82   //AliAnalysisTaskPi0Flow::MakeMCHistograms();
83   
84   // MC Generated histograms
85   char key[55];
86   for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
87     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
88     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
89     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
90     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
91     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
92     fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
93     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
94     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
95     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
96     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
97     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
98     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
99     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
100     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
101     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
102     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
103     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
104     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
105     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
106     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
107     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
108     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
109     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
110     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
111   }
112   fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
113   fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
114   fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
115  
116  
117   Int_t nPt      = 200;
118   Double_t ptMin = 0;
119   Double_t ptMax = 20; 
120   fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
121   fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
122   fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
123   fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
124   fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
125   fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
126 //  fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
127 //  fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
128
129   const Int_t nM       = 500;
130   const Double_t mMin  = 0.0;
131   const Double_t mMax  = 1.0;
132
133   for(Int_t cen=0; cen < fCentEdges.GetSize()-1; cen++){
134     fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
135     fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
136     fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
137     fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
138     fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
139     fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
140     fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
141     fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
142     fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
143     fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
144     fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
145     fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
146     fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
147
148     //pairs from common parents
149     fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
150     fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
151     fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
152     fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
153     fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
154     fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
155     fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
156    
157     //common parent - pi0
158     fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
159     fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
160     fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
161     fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
162     fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
163     fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
164     fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
165     fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
166     fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
167     
168   }
169   
170   
171   //Photon contaminations
172   fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
173   fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
174   fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
175   fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
176   fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.)); 
177    
178    const Int_t nTypes=24 ;
179    char partTypes[nTypes][55] ;
180    snprintf(partTypes[0],55,"hGammaNoPrim") ; //
181    snprintf(partTypes[1],55,"hGammaPhot") ; //
182    snprintf(partTypes[2],55,"hGammaEl") ; //
183    snprintf(partTypes[3],55,"hGammaPi0") ; //
184    snprintf(partTypes[4],55,"hGammaEta") ; //
185    snprintf(partTypes[5],55,"hhGammaOmega") ; //
186    snprintf(partTypes[6],55,"hGammaPipm") ; //
187    snprintf(partTypes[7],55,"hGammaP") ; //
188    snprintf(partTypes[8],55,"hGammaPbar") ; //
189    snprintf(partTypes[9],55,"hGammaN") ; //
190    snprintf(partTypes[10],55,"hGammaNbar") ; //
191    snprintf(partTypes[11],55,"hGammaK0S") ; //
192    snprintf(partTypes[12],55,"hGammaK0L") ; //
193    snprintf(partTypes[13],55,"hGammaKpm") ; //
194    snprintf(partTypes[14],55,"hGammaKstar") ; //
195    snprintf(partTypes[15],55,"hGammaDelta") ; //
196    snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
197    snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
198    snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
199    snprintf(partTypes[19],55,"hGammaPipmEl") ; //
200    snprintf(partTypes[20],55,"hGammaPipmOther") ; //
201    snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
202    snprintf(partTypes[22],55,"hGammaPipmp") ; //
203    snprintf(partTypes[23],55,"hGammaPipmn") ; //
204  
205    const Int_t nPID=12 ;
206    char cPID[25][12] ;
207    snprintf(cPID[0],25,"All") ;
208    snprintf(cPID[1],25,"Allcore") ;
209    snprintf(cPID[2],25,"CPV") ;
210    snprintf(cPID[3],25,"CPVcore") ;
211    snprintf(cPID[4],25,"CPV2") ;
212    snprintf(cPID[5],25,"CPV2core") ;
213    snprintf(cPID[6],25,"Disp") ;
214    snprintf(cPID[7],25,"Dispcore") ;
215    snprintf(cPID[8],25,"Disp2") ;
216    snprintf(cPID[9],25,"Disp2core") ;
217    snprintf(cPID[10],25,"Both") ;
218    snprintf(cPID[11],25,"Bothcore") ;
219  
220    for(Int_t itype=0; itype<nTypes; itype++){
221      for(Int_t iPID=0; iPID<nPID; iPID++){
222        for(Int_t cen=0; cen<5; cen++){
223          fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
224        }
225      }
226    }  
227 }
228
229 void AliAnalysisTaskPi0FlowMC::DoMC()
230 {
231   //AliAnalysisTaskPi0Flow::DoMC();
232   fStack = GetMCStack();
233   
234   //TODO: Geometry IHEP?
235   //TODO: decalib.?
236   //TODO: PHOS matrix?
237   //TODO: Centrality.
238   
239   FillMCHist();
240
241   // Proccess all selected clusters
242   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
243     //Bool_t sure=  kTRUE;
244     //Int_t primary=FindPrimary(clu,sure) ;  //номер праймари частицы в стеке
245     //ph->SetPrimary(primary) ;
246     //ph->SetWeight(PrimaryWeight(primary)) ;
247   }
248   
249   FillSecondaries() ;
250 }
251
252
253
254 //___________________________________________________________________________
255 void AliAnalysisTaskPi0FlowMC::FillMCHist(){
256   //fill histograms for efficiensy etc. calculation
257
258   //---------First pi0/eta-----------------------------
259   char partName[10] ;
260   char hkey[55] ;
261
262   if(!fStack) return ;
263   for(Int_t i=0;i<fStack->GetNtrack();i++){
264      TParticle* particle =  fStack->Particle(i);
265     if(particle->GetPdgCode() == kPi0)
266       snprintf(partName,10,"pi0") ;
267     else
268       if(particle->GetPdgCode() == kEta)
269         snprintf(partName,10,"eta") ;
270       else
271         if(particle->GetPdgCode() == kGamma)
272            snprintf(partName,10,"gamma") ;
273         else
274            continue ;
275
276     //Primary particle
277     Double_t r=particle->R() ;
278     Double_t pt = particle->Pt() ;
279     //Distribution over vertex
280     FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
281     
282     if(r >kRCut)
283       continue ;
284
285     //Total number of pi0 with creation radius <1 cm
286     Double_t weight = PrimaryParticleWeight(particle) ;  
287     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCentBin) ;
288     FillHistogram(hkey,pt,weight) ;
289     if(TMath::Abs(particle->Y())<0.12){
290       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCentBin) ;
291       FillHistogram(hkey,pt,weight) ;
292     }
293
294     snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCentBin) ;
295     FillHistogram(hkey,particle->Y(),weight) ;
296     
297     Double_t phi=particle->Phi() ;
298     while(phi<0.)phi+=TMath::TwoPi() ;
299     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
300     snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCentBin) ;
301     FillHistogram(hkey,phi,weight) ;
302   }
303 }
304
305 //________________________________________________________________________
306 void AliAnalysisTaskPi0FlowMC::FillSecondaries(){
307   //Sort secondaires
308   
309   //Fill spectra of primary particles 
310   //with proper weight
311   std::cout <<  fStack << std::endl;
312   AliInfo("start");
313   for(Int_t i=0; i<fStack->GetNtrack(); i++){
314     TParticle * p = fStack->Particle(i) ;
315     if(p->R()>kRCut)
316       continue ;
317     if(TMath::Abs(p->Y())>0.5)
318       continue ;
319     Double_t w = PrimaryParticleWeight(p) ;  
320     Int_t primPdgCode=p->GetPdgCode() ;
321       switch(primPdgCode){
322         case  kGamma: FillHistogram(Form("hPrimPhot_cen%d",fCentBin),p->Pt(),w); 
323                   break ;
324         case  kElectron: 
325         case -kElectron: 
326                   FillHistogram(Form("hPrimEl_cen%d",fCentBin),p->Pt(),w); 
327                   break ;
328         case  kPi0: 
329                   FillHistogram(Form("hPrimPi0_cen%d",fCentBin),p->Pt(),w); 
330                   break ;
331         case  kEta: 
332                   FillHistogram(Form("hPrimEta_cen%d",fCentBin),p->Pt(),w); 
333                   break ;
334         case  kPiPlus: 
335         case  kPiMinus: 
336                   FillHistogram(Form("hPrimPipm_cen%d",fCentBin),p->Pt(),w); 
337                   break ;                 
338         case  kProton:  //p 
339                   FillHistogram(Form("hPrimP_cen%d",fCentBin),p->Pt(),w); 
340                   break ;                 
341         case kProtonBar:  //pbar
342                   FillHistogram(Form("hPrimPbar_cen%d",fCentBin),p->Pt(),w); 
343                   break ;                 
344         case  kNeutron:  //n 
345                   FillHistogram(Form("hPrimN_cen%d",fCentBin),p->Pt(),w); 
346                   break ;                 
347         case  kNeutronBar:  //nbar
348                   FillHistogram(Form("hPrimNbar_cen%d",fCentBin),p->Pt(),w); 
349                   break ;
350         case  310:  //nbar
351                   FillHistogram(Form("hPrimK0S_cen%d",fCentBin),p->Pt(),w); 
352                   break ;
353         case  130:  //nbar
354                   FillHistogram(Form("hPrimK0L_cen%d",fCentBin),p->Pt(),w); 
355                   break ;
356         case  321:  //K+
357         case -321:  //K-
358                   FillHistogram(Form("hPrimKpm_cen%d",fCentBin),p->Pt(),w); 
359                   break ;
360         default:           //other
361                   FillHistogram(Form("hPrimOther_cen%d",fCentBin),p->Pt(),w);    
362       }
363   }
364   AliInfo("Origins of secondary pi0s");
365   //Origins of secondary pi0s
366   for(Int_t i=0; i<fStack->GetNtrack(); i++){
367     TParticle * p = fStack->Particle(i) ;
368     if(p->GetPdgCode()!=111)
369       continue ;
370     FillHistogram("Vertex",p->Pt(),p->R());
371     if(p->R()<kRCut)
372       continue ;
373     Double_t phi=p->Phi() ;
374     while(phi<0.)phi+=TMath::TwoPi() ;
375     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
376     FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ;   
377     Double_t w = PrimaryParticleWeight(p) ;  
378     FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ;   
379   }
380
381   TLorentzVector p1;
382
383   Int_t inPHOS=fCaloPhotonsPHOS->GetEntries() ;
384   for (Int_t i1=0; i1<inPHOS-1; i1++) {
385     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
386     Double_t w1=ph1->GetWeight() ;
387     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
388       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
389       TLorentzVector p12  = *ph1  + *ph2;
390       Double_t w2=ph2->GetWeight() ;
391       Double_t w = TMath::Sqrt(w1*w2) ;
392       FillHistogram(Form("hParentAll_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
393       Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
394       if(prim>-1){
395         TParticle * particle = (TParticle *)fStack->Particle(prim);
396         FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ;
397                 
398         
399         Int_t pdgCode=particle->GetPdgCode() ;
400         if(pdgCode!=111){ //common parent - not pi0
401           if(pdgCode==22)  
402             FillHistogram(Form("hParentGamma_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
403           else{             
404             if(pdgCode==11 || pdgCode==-11){   
405               FillHistogram(Form("hParentEl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
406             }
407             else{
408               if(InPi0mass(p12.M() ,p12.Pt())){
409                 printf("Common parent: %d \n",pdgCode) ;
410               }
411               FillHistogram(Form("hParentOther_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
412             }
413           }//Not photons
414         }//Common parent not pi0
415         else{ //common parent - pi0
416           FillHistogram(Form("hParentPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; 
417           FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ;   
418           FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ;  
419           if(particle->R()<kRCut && TMath::Abs(particle->Vz())<fMaxAbsVertexZ){
420             FillHistogram(Form("hParentDirPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
421             continue ;
422           }
423           //Common particle pi0, created off-vertex
424           Int_t primPi0=particle->GetFirstMother();
425           if(primPi0==-1){
426             FillHistogram(Form("hParentPi0NoPrim_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
427           }
428           else{
429             Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ;
430             switch(primPdgCode){
431             case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //eta
432                       break ;
433             case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //omega
434                       break ;
435             case  211:  //pi+-
436             case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
437                       break ;
438             case  321:  //K+-
439             case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
440                       break ;
441             case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0S
442                       break ;
443             case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0L
444                       break ;
445             case  2212:  //p 
446             case  2112:  //n 
447                       FillHistogram(Form("hParentPi0pn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
448                       break ;
449             case -2212:  //pbar
450             case -2112:  //nbar
451                       FillHistogram(Form("hParentPi0antipn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
452                       break ;
453             default:       //other
454                       FillHistogram(Form("hParentPi0Other_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
455             }//switch     
456           }//pi0 with primary
457         }//common parent - pi0
458       }//there is common primary 
459     }//seond photon loop
460   }//first photon loop
461   
462   
463   //Now look at photon contaiminations
464   for (Int_t i1=0; i1<inPHOS-1; i1++) {
465     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
466     Int_t iprim = ph1->GetPrimary() ;
467     if(iprim<0)
468       FillAllHistograms(Form("hGammaNoPrim_cen%d",fCentBin),ph1) ; //
469     else{
470       //Find primary at vertex
471       TParticle * primPHOS = fStack->Particle(iprim) ;
472       Int_t iprimV=primPHOS->GetFirstMother();
473       TParticle * primVtx = primPHOS ;
474       while((iprimV>-1) && primVtx->R()>kRCut){
475         primVtx = fStack->Particle(iprimV) ;
476         iprimV=primVtx->GetFirstMother();
477       }
478     
479       //photon
480       Int_t primPdgCode=primVtx->GetPdgCode() ;
481       switch(primPdgCode){
482         case  22: FillAllHistograms("hGammaPhot",ph1); 
483                   break ;
484         case  11: 
485         case -11: 
486                   FillAllHistograms("hGammaEl",ph1); 
487                   break ;
488         case  111: 
489                   FillAllHistograms("hGammaPi0",ph1); 
490                   break ;
491         case  221: 
492                   FillAllHistograms("hGammaEta",ph1); 
493                   break ;
494         case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
495                   break ;
496         case  211: 
497         case -211: 
498                   FillAllHistograms("hGammaPipm",ph1); 
499                   //Find particle entered PHOS
500                   if(primVtx == primPHOS)
501                     FillAllHistograms("hGammaPipmDirect",ph1); 
502                   else{
503                     Int_t primPdgPHOS=primPHOS->GetPdgCode() ;
504                     if(primPdgPHOS==22){
505                        FillAllHistograms("hGammaPipmGamma",ph1); 
506                        FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R());
507                        FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R());
508                        break ;            
509                     }
510                     if(TMath::Abs(primPdgPHOS)==11){
511                        FillAllHistograms("hGammaPipmEl",ph1); 
512                        FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
513                        break ;            
514                     }
515                     if(TMath::Abs(primPdgPHOS)==2212){
516                        FillAllHistograms("hGammaPipmp",ph1); 
517                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
518                        break ;            
519                     }
520                     if(TMath::Abs(primPdgPHOS)==2112){
521                        FillAllHistograms("hGammaPipmn",ph1); 
522                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
523                        break ;            
524                     }
525                     FillAllHistograms("hGammaPipmOther",ph1); 
526                     FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());               
527                   }
528                   break ;                 
529         case  2212:  //p 
530                   FillAllHistograms("hGammaP",ph1); 
531                   break ;                 
532         case -2212:  //pbar
533                   FillAllHistograms("hGammaPbar",ph1); 
534                   break ;                 
535         case  2112:  //n 
536                   FillAllHistograms("hGammaN",ph1); 
537                   break ;                 
538         case -2112:  //nbar
539                   FillAllHistograms("hGammaNbar",ph1) ; // pn
540                   break ;
541         case  310:  //nbar
542                   FillAllHistograms("hGammaK0S",ph1) ; // pn
543                   break ;
544         case  130:  //nbar
545                   FillAllHistograms("hGammaK0L",ph1) ; // pn
546                   break ;
547         case  321:  //K+
548         case -321:  //K-
549                   FillAllHistograms("hGammaKpm",ph1) ; // pn
550                   break ;
551         case -323: 
552         case  323: 
553         case -313: 
554         case  313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
555                   break ;
556                   
557         case -2224 : //Deltas
558         case  2224 : //Deltas
559         case -2214 : //Deltas
560         case  2214 : //Deltas
561         case -2114 : //Deltas
562         case  2114 : //Deltas
563         case -1114 : //Deltas
564         case  1114 : //Deltas
565                   FillAllHistograms("hGammaDelta",ph1) ; // pn
566                   break ;                 
567         default:           //other
568             if(primVtx->GetPDG()->Charge())
569               FillAllHistograms("hGammaOtherCharged",ph1) ; //
570             else
571               FillAllHistograms("hGammaOtherNeutral",ph1) ; //
572       }
573     }
574   
575   }//single photons
576 }
577
578 //_____________________________________________________________________________
579 void AliAnalysisTaskPi0FlowMC::FillAllHistograms(const char * particleName,AliCaloPhoton * ph)
580 {
581   //Fill All PID histograms
582         
583   Double_t w=ph->GetWeight() ;
584   Double_t pt = ph->Pt() ;
585   Double_t ptC=ph->GetMomV2()->Pt() ;
586   FillHistogram(Form("%s_All_cen%d",particleName,fCentBin),pt,w) ;
587   FillHistogram(Form("%s_Allcore_cen%d",particleName,fCentBin),ptC,w) ;
588   if(ph->IsCPVOK()){
589     FillHistogram(Form("%s_CPV_cen%d",particleName,fCentBin),pt,w) ;
590     FillHistogram(Form("%s_CPVcore_cen%d",particleName,fCentBin),ptC,w) ;
591   }
592   if(ph->IsCPV2OK()){
593     FillHistogram(Form("%s_CPV2_cen%d",particleName,fCentBin),pt,w) ;
594     FillHistogram(Form("%s_CPV2core_cen%d",particleName,fCentBin),ptC,w) ;
595   }
596   if(ph->IsDispOK()){     
597     FillHistogram(Form("%s_Disp_cen%d",particleName,fCentBin),pt,w) ;
598     FillHistogram(Form("%s_Dispcore_cen%d",particleName,fCentBin),ptC,w) ;
599     if(ph->IsDisp2OK()){
600       FillHistogram(Form("%s_Disp2_cen%d",particleName,fCentBin),pt,w) ;
601       FillHistogram(Form("%s_Disp2core_cen%d",particleName,fCentBin),ptC,w) ;
602     }
603     if(ph->IsCPVOK()){
604       FillHistogram(Form("%s_Both_cen%d",particleName,fCentBin),pt,w) ;
605       FillHistogram(Form("%s_Bothcore_cen%d",particleName,fCentBin),ptC,w) ;
606     }
607   }  
608 }
609
610
611 //___________________________________________________________________________
612 Double_t AliAnalysisTaskPi0FlowMC::PrimaryWeight(Int_t primary){
613    //Check who is the primary and introduce weight to correct primary spectrum
614   
615   if(primary<0 || primary>=fStack->GetNtrack())
616     return 1 ;
617   //trace primaries up to IP
618   TParticle* particle =  fStack->Particle(primary);
619   Double_t r=particle->R() ;
620   Int_t mother = particle->GetFirstMother() ;
621   while(mother>-1){
622     if(r<1. && particle->GetPdgCode()==111)
623       break ;
624     particle =  fStack->Particle(mother);
625     mother = particle->GetFirstMother() ;
626     r=particle->R() ;
627   }
628
629   return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
630 }
631 //________________________________________________________________________
632 Double_t AliAnalysisTaskPi0FlowMC::PrimaryParticleWeight(TParticle * particle){
633   return 1.; //TODO: use weight.
634   
635   Int_t pdg = particle->GetPdgCode() ;
636   Int_t type=0 ;
637   if(pdg == 111 || TMath::Abs(pdg)==211){
638     type =1 ;
639   }
640   else{
641     if(TMath::Abs(pdg)<1000){ //Kaon-like
642       type =2 ;    
643     }
644     else
645       type = 3;  //baryons
646   }
647     
648   Double_t pt = particle->Pt() ;
649   if(type==1){
650    if(fCentBin==0) //0-5
651      return (1.662990+1.140890*pt-0.192088*pt*pt)/(1.-0.806630*pt+0.304771*pt*pt)+0.141690*pt ;
652    if(fCentBin==1) //5-10
653      return (1.474351+0.791492*pt-0.066369*pt*pt)/(1.-0.839338*pt+0.317312*pt*pt)+0.093289*pt ;
654    if(fCentBin==2) //10-20
655      return (1.174728+0.959681*pt-0.137695*pt*pt)/(1.-0.788873*pt+0.299538*pt*pt)+0.128759*pt ; 
656    if(fCentBin==3) //20-40
657      return (0.927335+0.475349*pt+0.004364*pt*pt)/(1.-0.817966*pt+0.309787*pt*pt)+0.086899*pt ; 
658    if(fCentBin==4) //40-60
659      return (0.676878+0.190680*pt+0.077031*pt*pt)/(1.-0.790623*pt+0.305183*pt*pt)+0.064510*pt ; 
660    if(fCentBin==5) //60-80
661      return (0.684726-0.606262*pt+0.409963*pt*pt)/(1.-1.080061*pt+0.456933*pt*pt)+0.005151*pt ; 
662   }
663   if(type==2){
664    if(fCentBin==0) //0-5
665      return (-0.417131+2.253936*pt-0.337731*pt*pt)/(1.-0.909892*pt+0.316820*pt*pt)+0.157312*pt ;
666    if(fCentBin==1) //5-10
667      return (-0.352275+1.844466*pt-0.248598*pt*pt)/(1.-0.897048*pt+0.316462*pt*pt)+0.132461*pt ; 
668    if(fCentBin==2) //10-20
669      return (-0.475481+1.975108*pt-0.336013*pt*pt)/(1.-0.801028*pt+0.276705*pt*pt)+0.188164*pt ; 
670    if(fCentBin==3) //20-40
671      return (-0.198954+1.068789*pt-0.103540*pt*pt)/(1.-0.848354*pt+0.299209*pt*pt)+0.112939*pt ; 
672    if(fCentBin==4) //40-60
673      return (-0.111052+0.664041*pt-0.019717*pt*pt)/(1.-0.804916*pt+0.300779*pt*pt)+0.095784*pt ;
674    if(fCentBin==5) //0-5
675      return (0.202788-0.439832*pt+0.564585*pt*pt)/(1.-1.254029*pt+0.679444*pt*pt)+0.016235*pt ;
676   }
677   if(type==3){
678    if(fCentBin==0) //0-5
679      return (-1.312732+2.743568*pt-0.375775*pt*pt)/(1.-0.717533*pt+0.164694*pt*pt)+0.164445*pt ;
680    if(fCentBin==1) //5-10
681      return (-1.229425+2.585889*pt-0.330164*pt*pt)/(1.-0.715892*pt+0.167386*pt*pt)+0.133085*pt ; 
682    if(fCentBin==2) //10-20
683      return (-1.135677+2.397489*pt-0.320355*pt*pt)/(1.-0.709312*pt+0.164350*pt*pt)+0.146095*pt ; 
684    if(fCentBin==3) //20-40
685      return (-0.889993+1.928263*pt-0.220785*pt*pt)/(1.-0.715991*pt+0.174729*pt*pt)+0.095098*pt ; 
686    if(fCentBin==4) //40-60
687      return (-0.539237+1.329118*pt-0.115439*pt*pt)/(1.-0.722906*pt+0.186832*pt*pt)+0.059267*pt ; 
688    if(fCentBin==5) //60-80
689      return (-0.518126+1.327628*pt-0.130881*pt*pt)/(1.-0.665649*pt+0.184300*pt*pt)+0.081701*pt ;   
690   }
691   return 1. ;  
692 }
693
694 //___________________________________________________________________________
695 AliStack* AliAnalysisTaskPi0FlowMC::GetMCStack()
696 {
697   fStack = 0;
698   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
699   if(eventHandler){
700     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
701     if( mcEventHandler)
702       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
703   }
704   return fStack;
705 }