New class for MC AOD analysis (B.Polishchuk)
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0FlowMCAOD.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 // Adopted for AOD analysis by Paul Batzing and Boris Polishchuk (10.03.2014).
21 /* $Id$ */
22
23 #include "TChain.h"
24 #include "TTree.h"
25 #include "TObjArray.h"
26 #include "TF1.h"
27 #include "TFile.h"
28 #include "TH1F.h"
29 #include "TH2F.h"
30 #include "TH2I.h"
31 #include "TH3F.h"
32 #include "TParticle.h"
33 #include "TCanvas.h"
34 #include "TStyle.h"
35 #include "TRandom.h"
36 #include "TROOT.h"
37 #include "THashList.h"
38 #include "TArray.h"
39 #include "TArrayD.h"
40
41
42 #include "AliAnalysisManager.h"
43 #include "AliMCEventHandler.h"
44 #include "AliMCEvent.h"
45 #include "AliStack.h"
46 #include "AliAnalysisTaskSE.h"
47 #include "AliPHOSHijingEfficiency.h"
48 #include "AliCaloPhoton.h"
49 #include "AliPHOSGeometry.h"
50 #include "AliPHOSEsdCluster.h"
51 #include "AliPHOSCalibData.h"
52 #include "AliLog.h"
53 #include "AliPID.h"
54 #include "AliCDBManager.h"
55 #include "AliCentrality.h" 
56 #include "AliEventplane.h"
57 #include "TProfile.h"
58 #include <TPDGCode.h>
59 #include "AliOADBContainer.h"
60 #include "AliAODInputHandler.h"
61 #include "AliAODMCParticle.h"
62
63 #include "AliAnalysisTaskPi0FlowMCAOD.h"
64 ClassImp(AliAnalysisTaskPi0FlowMCAOD);
65
66 const Double_t AliAnalysisTaskPi0FlowMCAOD::kRCut = 1.;
67
68 AliAnalysisTaskPi0FlowMCAOD::AliAnalysisTaskPi0FlowMCAOD(const char* name, AliAnalysisTaskPi0Flow::Period period)
69 : AliAnalysisTaskPi0Flow(name, period),fMcArray(0x0)
70 {
71 }
72
73 AliAnalysisTaskPi0FlowMCAOD::~AliAnalysisTaskPi0FlowMCAOD()
74 {
75 }
76
77 void AliAnalysisTaskPi0FlowMCAOD::UserCreateOutputObjects()
78 {
79   // Do standard Pi0Flow CreateOuputObjects
80   AliAnalysisTaskPi0Flow::UserCreateOutputObjects();
81   
82   // MC Generated histograms
83   char key[55];
84   for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
85     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
86     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
87     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
88     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
89     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
90     fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
91     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
92     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
93     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
94     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
95     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
96     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
97     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
98     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
99     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
100     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
101     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
102     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
103     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
104     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
105     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
106     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
107     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
108     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
109
110     snprintf(key,55,"hMC_all_K0S_cen%d",cent) ;
111     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
112
113     snprintf(key,55,"hMC_unitEta_K0S_cen%d",cent) ;
114     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
115
116   }
117   fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
118   fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
119   fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
120  
121  
122   Int_t nPt      = 200;
123   Double_t ptMin = 0;
124   Double_t ptMax = 20; 
125   fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
126   fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
127   fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
128   fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
129   fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
130   fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
131 //  fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
132 //  fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
133
134   const Int_t nM       = 500;
135   const Double_t mMin  = 0.0;
136   const Double_t mMax  = 1.0;
137
138   for(Int_t cen=0; cen < fCentEdges.GetSize()-1; cen++){
139     fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
140     fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
141     fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
142     fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
143     fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
144     fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
145     fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
146     fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
147     fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
148     fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
149     fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
150     fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
151     fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
152
153     //pairs from common parents
154     fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
155     fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
156     fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
157     fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
158     fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
159     fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
160     fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
161    
162     //common parent - pi0
163     fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
164     fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
165     fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
166     fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
167     fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
168     fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
169     fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
170     fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
171     fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
172     
173   }
174   
175   
176   //Photon contaminations
177   fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
178   fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
179   fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
180   fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
181   fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.)); 
182    
183    const Int_t nTypes=24 ;
184    char partTypes[nTypes][55] ;
185    snprintf(partTypes[0],55,"hGammaNoPrim") ; //
186    snprintf(partTypes[1],55,"hGammaPhot") ; //
187    snprintf(partTypes[2],55,"hGammaEl") ; //
188    snprintf(partTypes[3],55,"hGammaPi0") ; //
189    snprintf(partTypes[4],55,"hGammaEta") ; //
190    snprintf(partTypes[5],55,"hhGammaOmega") ; //
191    snprintf(partTypes[6],55,"hGammaPipm") ; //
192    snprintf(partTypes[7],55,"hGammaP") ; //
193    snprintf(partTypes[8],55,"hGammaPbar") ; //
194    snprintf(partTypes[9],55,"hGammaN") ; //
195    snprintf(partTypes[10],55,"hGammaNbar") ; //
196    snprintf(partTypes[11],55,"hGammaK0S") ; //
197    snprintf(partTypes[12],55,"hGammaK0L") ; //
198    snprintf(partTypes[13],55,"hGammaKpm") ; //
199    snprintf(partTypes[14],55,"hGammaKstar") ; //
200    snprintf(partTypes[15],55,"hGammaDelta") ; //
201    snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
202    snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
203    snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
204    snprintf(partTypes[19],55,"hGammaPipmEl") ; //
205    snprintf(partTypes[20],55,"hGammaPipmOther") ; //
206    snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
207    snprintf(partTypes[22],55,"hGammaPipmp") ; //
208    snprintf(partTypes[23],55,"hGammaPipmn") ; //
209  
210    const Int_t nPID=12 ;
211    char cPID[12][25] ;
212    snprintf(cPID[0],25,"All") ;
213    snprintf(cPID[1],25,"Allcore") ;
214    snprintf(cPID[2],25,"CPV") ;
215    snprintf(cPID[3],25,"CPVcore") ;
216    snprintf(cPID[4],25,"CPV2") ;
217    snprintf(cPID[5],25,"CPV2core") ;
218    snprintf(cPID[6],25,"Disp") ;
219    snprintf(cPID[7],25,"Dispcore") ;
220    snprintf(cPID[8],25,"Disp2") ;
221    snprintf(cPID[9],25,"Disp2core") ;
222    snprintf(cPID[10],25,"Both") ;
223    snprintf(cPID[11],25,"Bothcore") ;
224  
225    for(Int_t itype=0; itype<nTypes; itype++){
226      for(Int_t iPID=0; iPID<nPID; iPID++){
227        for(Int_t cen=0; cen < fCentEdges.GetSize()-1; cen++){
228          fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
229        }
230      }
231    }
232
233   PostData(1, fOutputContainer);
234 }
235
236 void AliAnalysisTaskPi0FlowMCAOD::UserExec(Option_t* option)
237 {
238   fMcArray = GetMCArray();
239   if(!fMcArray) { PostData(1, fOutputContainer); return; }
240  
241   AliAnalysisTaskPi0Flow::UserExec(option);
242 }
243
244
245 void AliAnalysisTaskPi0FlowMCAOD::SelectPhotonClusters()
246 {
247   AliAnalysisTaskPi0Flow::SelectPhotonClusters();
248   
249   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
250     AliCaloPhoton * photon = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1);
251     AliVCluster* cluster = photon->GetCluster();
252     Bool_t sure=  kTRUE;
253     Int_t primary=FindPrimary(cluster,sure) ;
254     photon->SetPrimary(primary);
255     photon->SetWeight(PrimaryWeight(primary)) ;
256   }
257
258   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
259     AliCaloPhoton * photon = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1);
260     Int_t primary = photon->GetPrimary();
261     AliAODMCParticle* p = GetParticle(primary);
262     if(R(p) >kRCut) {
263       if(p->PdgCode()==11 || p->PdgCode()==-11) continue;
264       else { fCaloPhotonsPHOS->Remove(photon); fCaloPhotonsPHOS->Compress(); }
265     }
266   }
267     
268 }
269
270 void AliAnalysisTaskPi0FlowMCAOD::FillSelectedClusterHistograms()
271 {
272   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
273     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
274
275     Double_t dphiA=ph1->Phi()-fRPV0A ;
276     while(dphiA<0)dphiA+=TMath::Pi() ;
277     while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
278
279     Double_t dphiC=ph1->Phi()-fRPV0C ;
280     while(dphiC<0)dphiC+=TMath::Pi() ;
281     while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
282
283     Double_t dphiT=ph1->Phi()-fRP ;
284     while(dphiT<0)dphiT+=TMath::Pi() ;
285     while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
286     
287     Double_t pt = ph1->Pt() ;
288     Double_t ptcore = ph1->GetMomV2()->Pt() ;
289     Double_t w = ph1->GetWeight();
290
291     FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),pt,dphiA, w) ;
292     FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),pt,dphiC, w) ;
293     if(fHaveTPCRP)
294       FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),pt,dphiT, w) ;
295     FillHistogram(Form("hPhotPhiV0AAllcore_cen%d",fCentBin),ptcore,dphiA, w) ;
296     FillHistogram(Form("hPhotPhiV0CAllcore_cen%d",fCentBin),ptcore,dphiC, w) ;
297     if(fHaveTPCRP)
298       FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCentBin),ptcore,dphiT, w) ;
299
300     FillHistogram(Form("hPhotAll_cen%d",fCentBin),pt, w) ;
301     FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ptcore, w) ;
302     if(ph1->IsntUnfolded()){
303       FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),pt, w) ;
304       FillHistogram(Form("hPhotPhiV0AAllwou_cen%d",fCentBin),pt,dphiA, w) ;
305       FillHistogram(Form("hPhotPhiV0CAllwou_cen%d",fCentBin),pt,dphiC, w) ;
306       if(fHaveTPCRP)
307         FillHistogram(Form("hPhotPhiTPCAllwou_cen%d",fCentBin),pt,dphiT, w) ;
308     }
309     if(ph1->IsCPVOK()){
310       FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),pt,dphiA, w) ;
311       FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),pt,dphiC, w) ;
312       if(fHaveTPCRP)
313         FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),pt,dphiT, w) ;
314
315       FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ptcore,dphiA, w) ;
316       FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ptcore,dphiC, w) ;
317       if(fHaveTPCRP)
318         FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ptcore,dphiT, w) ;
319
320       FillHistogram(Form("hPhotCPV_cen%d",fCentBin),pt, w) ;
321       FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ptcore, w) ;
322     }
323     if(ph1->IsCPV2OK()){
324       FillHistogram(Form("hPhotPhiV0ACPV2_cen%d",fCentBin),pt,dphiA, w) ;
325       FillHistogram(Form("hPhotPhiV0CCPV2_cen%d",fCentBin),pt,dphiC, w) ;
326       if(fHaveTPCRP)
327         FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCentBin),pt,dphiT, w) ;
328
329       FillHistogram(Form("hPhotPhiV0ACPV2core_cen%d",fCentBin),ptcore,dphiA, w) ;
330       FillHistogram(Form("hPhotPhiV0CCPV2core_cen%d",fCentBin),ptcore,dphiC, w) ;
331       if(fHaveTPCRP)
332         FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCentBin),ptcore,dphiT, w) ;
333       FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),pt, w) ;
334       FillHistogram(Form("hPhotCPV2core_cen%d",fCentBin),ptcore, w) ;
335     }
336     if(ph1->IsDispOK()){
337       FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),pt,dphiA, w) ;
338       FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),pt,dphiC, w) ;
339       if(fHaveTPCRP)
340         FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),pt,dphiT, w) ;
341
342       FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ptcore,dphiA, w) ;
343       FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ptcore,dphiC, w) ;
344       if(fHaveTPCRP)
345         FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ptcore,dphiT, w) ;
346
347       if(ph1->IsntUnfolded()){
348         FillHistogram(Form("hPhotPhiV0ADispwou_cen%d",fCentBin),pt,dphiA, w) ;
349         FillHistogram(Form("hPhotPhiV0CDispwou_cen%d",fCentBin),pt,dphiC, w) ;
350         if(fHaveTPCRP)
351           FillHistogram(Form("hPhotPhiTPCDispwou_cen%d",fCentBin),pt,dphiT, w) ;
352
353       }
354       FillHistogram(Form("hPhotDisp_cen%d",fCentBin),pt, w) ;
355       FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ptcore, w) ;
356       if(ph1->IsntUnfolded()){
357         FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),pt, w) ;
358       }
359       if(ph1->IsCPVOK()){
360         FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),pt,dphiA, w) ;
361         FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),pt,dphiC, w) ;
362         if(fHaveTPCRP)
363           FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),pt,dphiT, w) ;
364
365         FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ptcore,dphiA, w) ;
366         FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ptcore,dphiC, w) ;
367         if(fHaveTPCRP)
368           FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ptcore,dphiT, w) ;
369
370         FillHistogram(Form("hPhotBoth_cen%d",fCentBin),pt, w) ;
371         FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ptcore, w) ;
372       }
373     }
374     if(ph1->IsDisp2OK()){
375       FillHistogram(Form("hPhotPhiV0ADisp2_cen%d",fCentBin),pt,dphiA, w) ;
376       FillHistogram(Form("hPhotPhiV0CDisp2_cen%d",fCentBin),pt,dphiC, w) ;
377       if(fHaveTPCRP)
378         FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCentBin),pt,dphiT, w) ;
379       FillHistogram(Form("hPhotPhiV0ADisp2core_cen%d",fCentBin),ptcore,dphiA, w) ;
380       FillHistogram(Form("hPhotPhiV0CDisp2core_cen%d",fCentBin),ptcore,dphiC, w) ;
381       if(fHaveTPCRP)
382         FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCentBin),ptcore,dphiT, w) ;
383
384       FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),pt, w) ;
385       FillHistogram(Form("hPhotDisp2core_cen%d",fCentBin),ptcore, w) ;
386       if(ph1->IsCPVOK()){
387         FillHistogram(Form("hPhotPhiV0ABoth2_cen%d",fCentBin),pt,dphiA, w) ;
388         FillHistogram(Form("hPhotPhiV0CBoth2_cen%d",fCentBin),pt,dphiC, w) ;
389         if(fHaveTPCRP)
390           FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCentBin),pt,dphiT, w) ;
391
392         FillHistogram(Form("hPhotPhiV0ABoth2core_cen%d",fCentBin),ptcore,dphiA, w) ;
393         FillHistogram(Form("hPhotPhiV0CBoth2core_cen%d",fCentBin),ptcore,dphiC, w) ;
394         if(fHaveTPCRP)
395           FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCentBin),ptcore,dphiT, w) ;
396
397         FillHistogram(Form("hPhotBoth2_cen%d",fCentBin),pt, w) ;
398         FillHistogram(Form("hPhotBoth2core_cen%d",fCentBin),ptcore, w) ;
399       }
400     }
401   }
402 }
403
404 void AliAnalysisTaskPi0FlowMCAOD::ConsiderPi0s()
405 {
406   char key[55];
407   for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
408     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
409     const Double_t w1 = ph1->GetWeight();
410     for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
411       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
412       TLorentzVector p12  = *ph1  + *ph2;
413       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
414       
415       const Double_t w2 = ph2->GetWeight();
416       Double_t w = TMath::Sqrt(w1*w2);
417       
418       FillHistogram("hPHOSphi",fCentrality,p12.Pt(),p12.Phi(), w) ;
419       Double_t dphiA=p12.Phi()-fRPV0A ;
420       while(dphiA<0)dphiA+=TMath::Pi() ;
421       while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
422
423       Double_t dphiC=p12.Phi()-fRPV0C ;
424       while(dphiC<0)dphiC+=TMath::Pi() ;
425       while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
426
427       Double_t dphiT=p12.Phi()-fRP ;
428       while(dphiT<0)dphiT+=TMath::Pi() ;
429       while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
430
431       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
432       Double_t m=p12.M() ;
433       Double_t mcore=pv12.M() ;
434       Double_t pt=p12.Pt() ;
435       Double_t ptcore=pv12.Pt() ;
436       Double_t pt1=ph1->Pt() ;
437       Double_t pt2=ph2->Pt() ;
438       Double_t ptcore1=ph1->GetMomV2()->Pt() ;
439       Double_t ptcore2=ph2->GetMomV2()->Pt() ;
440
441       FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA, w) ;
442       FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC, w) ;
443       if(fHaveTPCRP)
444         FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT, w) ;
445
446       FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
447       FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
448       if(fHaveTPCRP)
449         FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
450
451
452       FillHistogram(Form("hPi0All_cen%d",fCentBin),m,pt, w) ;
453       FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),mcore,ptcore, w) ;
454       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
455         FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),m,pt, w) ;
456         FillHistogram(Form("hMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA, w) ;
457         FillHistogram(Form("hMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC, w) ;
458         if(fHaveTPCRP)
459           FillHistogram(Form("hMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT, w) ;
460       }
461
462       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt1, w) ;
463       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt2, w) ;
464       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore1, w) ;
465       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore2, w) ;
466       if(ph1->IsntUnfolded())
467         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt1, w) ;
468       if(ph2->IsntUnfolded())
469         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt2, w) ;
470       if(ph1->IsCPVOK()){
471         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt1, w) ;
472         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore1, w) ;
473       }
474       if(ph2->IsCPVOK()){
475         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt2, w) ;
476         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore2, w) ;
477       }
478       if(ph1->IsCPV2OK()){
479         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt1, w) ;
480         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
481       }
482       if(ph2->IsCPV2OK()){
483         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt2, w) ;
484         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
485       }
486       if(ph1->IsDispOK()){
487         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt1, w) ;
488         if(ph1->IsntUnfolded()){
489           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt1, w) ;
490         }
491         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore1, w) ;
492       }
493       if(ph2->IsDispOK()){
494         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt2, w) ;
495         if(ph1->IsntUnfolded()){
496           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt2, w) ;
497         }
498         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore2, w) ;
499       }
500       if(ph1->IsDisp2OK()){
501         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt1, w) ;
502         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
503       }
504       if(ph2->IsDisp2OK()){
505         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt2, w) ;
506         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
507       }
508       if(ph1->IsDispOK() && ph1->IsCPVOK()){
509         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt1, w) ;
510         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore1, w) ;
511       }
512       if(ph2->IsDispOK() && ph2->IsCPVOK()){
513         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt2, w) ;
514         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore2, w) ;
515       }
516       if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
517         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt1, w) ;
518         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore1, w) ;
519       }
520       if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
521         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt2, w) ;
522         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore2, w) ;
523       }
524
525
526       if(a<kAlphaCut){
527         FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),m,pt, w) ;
528       }
529
530       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
531         snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
532         FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA, w) ;
533         FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC, w) ;
534         if(fHaveTPCRP)
535           FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT, w) ;
536
537         FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
538         FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
539         if(fHaveTPCRP)
540           FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
541
542         FillHistogram(Form("hPi0CPV_cen%d",fCentBin),m,pt, w) ;
543         FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),mcore, ptcore, w) ;
544
545         if(a<kAlphaCut){
546           FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),m,pt, w) ;
547         }
548       }
549       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
550         FillHistogram(Form("hMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA, w) ;
551         FillHistogram(Form("hMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC, w) ;
552         if(fHaveTPCRP)
553           FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT, w) ;
554         FillHistogram(Form("hMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
555         FillHistogram(Form("hMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
556         if(fHaveTPCRP)
557           FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
558         
559         FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),m,pt, w) ;
560         FillHistogram(Form("hPi0CPV2core_cen%d",fCentBin),mcore, ptcore, w) ;
561         if(a<kAlphaCut){
562           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),m,pt, w) ;
563         }
564       }
565       if(ph1->IsDispOK() && ph2->IsDispOK()){
566         snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
567         FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA, w) ;
568         FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC, w) ;
569         if(fHaveTPCRP)
570           FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT, w) ;
571         
572         FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
573         FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
574         if(fHaveTPCRP)
575           FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
576
577         FillHistogram(Form("hPi0Disp_cen%d",fCentBin),m,pt, w) ;
578         FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),mcore, ptcore, w) ;
579         
580         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
581           FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),m,pt, w) ;
582
583           FillHistogram(Form("hMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA, w) ;
584           FillHistogram(Form("hMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC, w) ;
585           if(fHaveTPCRP)
586             FillHistogram(Form("hMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT, w) ;
587         }
588
589         if(a<kAlphaCut){
590           FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),m,pt, w) ;
591         }
592         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
593           FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA, w) ;
594           FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC, w) ;
595           if(fHaveTPCRP)
596             FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT, w) ;
597
598           FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
599           FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
600           if(fHaveTPCRP)
601             FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
602
603           FillHistogram(Form("hPi0Both_cen%d",fCentBin),m,pt, w) ;
604           FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),mcore,ptcore, w) ;
605
606           if(a<kAlphaCut){
607             snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
608             FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),m,pt, w) ;
609           }
610           if(ph1->Module()==1 && ph2->Module()==1)
611             FillHistogram("hPi0M11", m, pt, w) ;
612           else if(ph1->Module()==2 && ph2->Module()==2)
613             FillHistogram("hPi0M22", m, pt, w) ;
614           else if(ph1->Module()==3 && ph2->Module()==3)
615             FillHistogram("hPi0M33", m, pt, w) ;
616           else if(ph1->Module()==1 && ph2->Module()==2)
617             FillHistogram("hPi0M12", m, pt, w) ;
618           else if(ph1->Module()==1 && ph2->Module()==3)
619             FillHistogram("hPi0M13", m, pt, w) ;
620           else if(ph1->Module()==2 && ph2->Module()==3)
621             FillHistogram("hPi0M23", m, pt, w) ;
622
623         }
624         
625       }
626       
627       
628       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
629         FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),m,pt,w) ;
630         FillHistogram(Form("hPi0Disp2core_cen%d",fCentBin),mcore, ptcore, w) ;  
631
632         FillHistogram(Form("hMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA, w) ;
633         FillHistogram(Form("hMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC, w) ;
634         if(fHaveTPCRP)
635           FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT, w) ;
636
637         FillHistogram(Form("hMassPtV0ADisp2core_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
638         FillHistogram(Form("hMassPtV0CDisp2core_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
639         if(fHaveTPCRP)
640           FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
641           
642         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
643           FillHistogram(Form("hMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA, w) ;
644           FillHistogram(Form("hMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC, w) ;
645           if(fHaveTPCRP)
646             FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT, w) ;
647
648           FillHistogram(Form("hMassPtV0ABoth2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
649           FillHistogram(Form("hMassPtV0CBoth2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
650           if(fHaveTPCRP)
651             FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
652
653           FillHistogram(Form("hPi0Both2_cen%d",fCentBin),m,pt, w) ;
654           FillHistogram(Form("hPi0Both2core_cen%d",fCentBin),mcore,ptcore, w) ;
655         }
656
657       }
658     } // end of loop i2
659   } // end of loop i1
660 }
661
662 //________________________________________________________________________
663 void AliAnalysisTaskPi0FlowMCAOD::ConsiderPi0sMix()
664 {
665   char key[55];
666
667   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
668
669   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
670     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
671     const Double_t w1 = ph1->GetWeight();
672     for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
673       TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
674       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
675         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
676         TLorentzVector p12  = *ph1  + *ph2;
677         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
678         
679         const Double_t w2 = ph2->GetWeight();
680         Double_t w = TMath::Sqrt(w1*w2);
681
682         Double_t dphiA=p12.Phi()-fRPV0A ;
683         while(dphiA<0)dphiA+=TMath::Pi() ;
684         while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
685
686         Double_t dphiC=p12.Phi()-fRPV0C ;
687         while(dphiC<0)dphiC+=TMath::Pi() ;
688         while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
689
690         Double_t dphiT=p12.Phi()-fRP ;
691         while(dphiT<0)dphiT+=TMath::Pi() ;
692         while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
693
694
695         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
696         Double_t m=p12.M() ;
697         Double_t mcore=pv12.M() ;
698         Double_t pt=p12.Pt() ;
699         Double_t ptcore=pv12.Pt() ;
700         Double_t pt1=ph1->Pt() ;
701         Double_t pt2=ph2->Pt() ;
702         Double_t ptcore1=ph1->GetMomV2()->Pt() ;
703         Double_t ptcore2=ph2->GetMomV2()->Pt() ;
704
705
706         snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
707         FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA, w) ;
708         FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC, w) ;
709         if(fHaveTPCRP)
710           FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT, w) ;
711
712         FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),mcore, ptcore, dphiA, w) ;
713         FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),mcore, ptcore, dphiC, w) ;
714         if(fHaveTPCRP)
715           FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),mcore, ptcore, dphiT, w) ;
716
717         FillHistogram(Form("hMiPi0All_cen%d",fCentBin),m,pt, w) ;
718         FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),mcore,ptcore, w) ;
719         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
720           FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),m,pt, w) ;
721           FillHistogram(Form("hMiMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA, w) ;
722           FillHistogram(Form("hMiMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC, w) ;
723           if(fHaveTPCRP)
724             FillHistogram(Form("hMiMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT, w) ;
725         }
726
727         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt1, w) ;
728         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt2, w) ;
729         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore1, w) ;
730         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore2, w) ;
731         if(ph1->IsntUnfolded())
732           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt1, w) ;
733         if(ph2->IsntUnfolded())
734           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt2, w) ;
735         if(ph1->IsCPVOK()){
736           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt1, w) ;
737           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore1, w) ;
738         }
739         if(ph2->IsCPVOK()){
740           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt2, w) ;
741           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore2, w) ;
742         }
743         if(ph1->IsCPV2OK()){
744           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt1, w) ;
745           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore1, w) ;
746         }
747         if(ph2->IsCPV2OK()){
748           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt2, w) ;
749           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
750         }
751         if(ph1->IsDispOK()){
752           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt1, w) ;
753           if(ph1->IsntUnfolded()){
754             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt1, w) ;
755           }
756           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore1, w) ;
757         }
758         if(ph2->IsDispOK()){
759           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt2, w) ;
760           if(ph1->IsntUnfolded()){
761             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt2, w) ;
762           }
763           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore2, w) ;
764         }
765         if(ph1->IsDisp2OK()){
766           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt1, w) ;
767           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
768         }
769         if(ph2->IsDisp2OK()){
770           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt2, w) ;
771           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore2, w) ;
772         }
773         if(ph1->IsDispOK() && ph1->IsCPVOK()){
774           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
775           FillHistogram(key,m,pt1, w) ;
776           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin);
777           FillHistogram(key,mcore,ptcore1, w) ;
778         }
779         if(ph2->IsDispOK() && ph2->IsCPVOK()){
780           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin);
781           FillHistogram(key,m,pt2, w) ;
782           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin);
783           FillHistogram(key,mcore,ptcore2, w) ;
784         }
785         if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
786           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt1, w) ;
787           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore1, w) ;
788         }
789         if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
790           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt2, w) ;
791           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore2, w) ;
792         }
793
794
795
796         if(a<kAlphaCut){
797           FillHistogram(Form("hMiPi0All_a07_cen%d",fCentBin),m,pt, w) ;
798         }
799         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
800           FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA, w) ;
801           FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC, w) ;
802           if(fHaveTPCRP)
803             FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT, w) ;
804
805           FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
806           FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
807           if(fHaveTPCRP)
808             FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
809
810           FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),m,pt, w) ;
811           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),mcore, ptcore, w) ;
812
813           if(a<kAlphaCut){
814             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),m,pt, w) ;
815           }
816         }
817         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
818           FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),m,pt, w) ;
819           FillHistogram(Form("hMiPi0CPV2core_cen%d",fCentBin),mcore, ptcore, w) ;
820
821           FillHistogram(Form("hMiMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA, w) ;
822           FillHistogram(Form("hMiMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC, w) ;
823           if(fHaveTPCRP)
824             FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT, w) ;
825           FillHistogram(Form("hMiMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
826           FillHistogram(Form("hMiMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
827           if(fHaveTPCRP)
828             FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
829
830           if(a<kAlphaCut){
831             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),m,pt, w) ;
832           }
833         }
834         if(ph1->IsDispOK() && ph2->IsDispOK()){
835           FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA, w) ;
836           FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC, w) ;
837           if(fHaveTPCRP)
838             FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT, w) ;
839
840           FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
841           FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
842           if(fHaveTPCRP)
843             FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
844
845
846           FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),m,pt, w) ;
847           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
848           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
849             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),m,pt, w) ;
850             FillHistogram(Form("hMiMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA, w) ;
851             FillHistogram(Form("hMiMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC, w) ;
852             if(fHaveTPCRP)
853               FillHistogram(Form("hMiMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT, w) ;
854           }
855
856           if(a<kAlphaCut){
857             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),m,pt, w) ;
858           }
859           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
860             FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA, w) ;
861             FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC, w) ;
862             if(fHaveTPCRP)
863               FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT, w) ;
864
865             FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
866             FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
867             if(fHaveTPCRP)
868               FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
869
870             FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),m,pt, w) ;
871             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
872
873             if(a<kAlphaCut){
874               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),m,pt, w) ;
875             }
876           }
877         }
878         
879         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
880           FillHistogram(Form("hMiMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA, w) ;
881           FillHistogram(Form("hMiMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC, w) ;
882           if(fHaveTPCRP)
883             FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT, w) ;
884
885           FillHistogram(Form("hMiMassPtV0ADisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
886           FillHistogram(Form("hMiMassPtV0CDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
887           if(fHaveTPCRP)
888             FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
889
890
891           FillHistogram(Form("hMiPi0Disp2_cen%d",fCentBin),m,pt, w) ;
892           FillHistogram(Form("hMiPi0Disp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
893
894           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
895             FillHistogram(Form("hMiMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA, w) ;
896             FillHistogram(Form("hMiMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC, w) ;
897             if(fHaveTPCRP)
898               FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT, w) ;
899
900             FillHistogram(Form("hMiMassPtV0ABoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
901             FillHistogram(Form("hMiMassPtV0CBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
902             if(fHaveTPCRP)
903               FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
904
905             FillHistogram(Form("hMiPi0Both2_cen%d",fCentBin),m,pt, w) ;
906             FillHistogram(Form("hMiPi0Both2core_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
907
908           }
909         }
910       } // end of loop i2
911     }
912   } // end of loop i1
913 }
914
915
916 void AliAnalysisTaskPi0FlowMCAOD::ProcessMC()
917 {
918   FillMCHist();
919   FillSecondaries() ;
920 }
921
922
923
924 //___________________________________________________________________________
925 void AliAnalysisTaskPi0FlowMCAOD::FillMCHist(){
926   //fill histograms for efficiensy etc. calculation
927   
928   //---------First pi0/eta-----------------------------
929   char partName[10] ; char hkey[55] ;
930   Int_t ntrack = fMcArray->GetEntriesFast(); 
931   
932   for(Int_t i=0;i<ntrack;i++){
933     AliAODMCParticle* particle = GetParticle(i);
934       
935     if(particle->PdgCode() == kPi0)
936       snprintf(partName,10,"pi0") ;
937     else
938       if(particle->PdgCode() == kEta)
939         snprintf(partName,10,"eta") ;
940       else
941         if(particle->PdgCode() == kGamma)
942           snprintf(partName,10,"gamma") ;
943         else
944           if(particle->PdgCode() == 310)
945             snprintf(partName,10,"K0S") ;
946           else
947             continue ;
948     //Primary particle
949     Double_t r  = R(particle) ;
950     Double_t pt = particle->Pt() ;
951     //Distribution over vertex
952     FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
953     
954     if(r >kRCut)
955       continue ;
956
957     Double_t phi=particle->Phi() ;
958     while(phi<0.)phi+=TMath::TwoPi() ;
959     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
960     
961     Double_t phig = 180./TMath::Pi()*phi; // phi in deg
962     
963     //Total number of pi0 with creation radius <1 cm
964     Double_t weight = PrimaryParticleWeight(particle) ;  
965     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCentBin) ;
966     
967     FillHistogram(hkey,pt,weight) ;
968     
969     if(TMath::Abs(particle->Y())<0.135 && phig>260. && phig<320.){
970       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCentBin) ;
971       FillHistogram(hkey,pt,weight) ;
972       
973       snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCentBin) ;
974       FillHistogram(hkey,particle->Y(),weight) ;
975       
976       snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCentBin) ;
977       FillHistogram(hkey,phi,weight) ;
978     }  
979   }   
980   
981 }
982  
983
984
985 //________________________________________________________________________
986 void AliAnalysisTaskPi0FlowMCAOD::FillSecondaries()
987 {
988   //Sort secondaires
989   
990   //Fill spectra of primary particles 
991   //with proper weight
992   if( fDebug )
993     AliInfo("start");
994   
995   Int_t ntrack = fMcArray->GetEntriesFast(); 
996   
997   for(Int_t i=0; i<ntrack; i++){
998     AliAODMCParticle* p = GetParticle(i);
999     
1000     if(R(p)>kRCut)
1001       continue ;
1002     if( TMath::Abs(p->Pt())<1.e-6 )
1003       continue;
1004     if(TMath::Abs(p->Y())>0.5)
1005       continue ;
1006     Double_t w = PrimaryParticleWeight(p) ;  
1007     Int_t primPdgCode=p->PdgCode() ;
1008     switch(primPdgCode){
1009     case  kGamma: FillHistogram(Form("hPrimPhot_cen%d",fCentBin),p->Pt(),w); 
1010       break ;
1011     case  kElectron: 
1012     case -kElectron: 
1013       FillHistogram(Form("hPrimEl_cen%d",fCentBin),p->Pt(),w); 
1014       break ;
1015     case  kPi0: 
1016       FillHistogram(Form("hPrimPi0_cen%d",fCentBin),p->Pt(),w); 
1017       break ;
1018     case  kEta: 
1019       FillHistogram(Form("hPrimEta_cen%d",fCentBin),p->Pt(),w); 
1020       break ;
1021     case  kPiPlus: 
1022     case  kPiMinus: 
1023       FillHistogram(Form("hPrimPipm_cen%d",fCentBin),p->Pt(),w); 
1024       break ;             
1025     case  kProton:  //p 
1026       FillHistogram(Form("hPrimP_cen%d",fCentBin),p->Pt(),w); 
1027       break ;             
1028     case kProtonBar:  //pbar
1029       FillHistogram(Form("hPrimPbar_cen%d",fCentBin),p->Pt(),w); 
1030       break ;             
1031     case  kNeutron:  //n 
1032       FillHistogram(Form("hPrimN_cen%d",fCentBin),p->Pt(),w); 
1033       break ;             
1034     case  kNeutronBar:  //nbar
1035       FillHistogram(Form("hPrimNbar_cen%d",fCentBin),p->Pt(),w); 
1036       break ;
1037     case  310:  //nbar
1038       FillHistogram(Form("hPrimK0S_cen%d",fCentBin),p->Pt(),w); 
1039       break ;
1040     case  130:  //nbar
1041       FillHistogram(Form("hPrimK0L_cen%d",fCentBin),p->Pt(),w); 
1042       break ;
1043     case  321:  //K+
1044     case -321:  //K-
1045       FillHistogram(Form("hPrimKpm_cen%d",fCentBin),p->Pt(),w); 
1046       break ;
1047     default:       //other
1048       FillHistogram(Form("hPrimOther_cen%d",fCentBin),p->Pt(),w);    
1049     }
1050   }
1051   if(fDebug)
1052     AliInfo("Origins of secondary pi0s");
1053   //Origins of secondary pi0s
1054   for(Int_t i=0; i<ntrack; i++){
1055     AliAODMCParticle* p = GetParticle(i);
1056     if(p->PdgCode()!=111)
1057       continue ;
1058     FillHistogram("Vertex",p->Pt(),R(p));
1059     if(R(p)<kRCut)
1060       continue ;
1061     Double_t phi=p->Phi() ;
1062     while(phi<0.)phi+=TMath::TwoPi() ;
1063     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1064     FillHistogram("hSecondPi0RphiZ",R(p),phi,p->Zv()) ;   
1065     Double_t w = PrimaryParticleWeight(p) ;  
1066     FillHistogram("hSecondPi0RE",R(p),p->Pt(),w) ;   
1067   }
1068
1069   TLorentzVector p1;
1070
1071   Int_t inPHOS=fCaloPhotonsPHOS->GetEntries() ;
1072   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1073     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1074     Double_t w1=ph1->GetWeight() ;
1075     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1076       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
1077       TLorentzVector p12  = *ph1  + *ph2;
1078       Double_t w2=ph2->GetWeight() ;
1079       Double_t w = TMath::Sqrt(w1*w2) ;
1080       FillHistogram(Form("hParentAll_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1081       Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
1082       if(prim>-1){
1083         AliAODMCParticle * particle = GetParticle(prim);
1084         FillHistogram("hMass_R",p12.M(),p12.Pt(),R(particle)) ;
1085                 
1086         
1087         Int_t pdgCode=particle->PdgCode() ;
1088         if(pdgCode!=111){ //common parent - not pi0
1089           if(pdgCode==22)  
1090             FillHistogram(Form("hParentGamma_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1091           else{             
1092             if(pdgCode==11 || pdgCode==-11){   
1093               FillHistogram(Form("hParentEl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1094             }
1095             else{
1096               if(InPi0mass(p12.M() ,p12.Pt())){
1097                 if(fDebug >1) AliInfo(Form("Common parent: %d",pdgCode));
1098               }
1099               FillHistogram(Form("hParentOther_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1100             }
1101           }//Not photons
1102         }//Common parent not pi0
1103         else{ //common parent - pi0
1104           FillHistogram(Form("hParentPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; 
1105           FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),R(particle),w) ;     
1106           FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Zv(),w) ;  
1107           if(R(particle)<kRCut && TMath::Abs(particle->Zv())<fMaxAbsVertexZ){
1108             FillHistogram(Form("hParentDirPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1109             continue ;
1110           }
1111           //Common particle pi0, created off-vertex
1112           Int_t primPi0 = particle->GetMother();
1113           if(primPi0==-1){
1114             FillHistogram(Form("hParentPi0NoPrim_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1115           }
1116           else{
1117             Int_t primPdgCode=((AliAODMCParticle *)fMcArray->At(primPi0))->PdgCode();
1118             switch(primPdgCode){
1119             case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //eta
1120               break ;
1121             case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //omega
1122               break ;
1123             case  211:  //pi+-
1124             case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1125               break ;
1126             case  321:  //K+-
1127             case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1128               break ;
1129             case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0S
1130               break ;
1131             case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0L
1132               break ;
1133             case  2212:  //p 
1134             case  2112:  //n 
1135               FillHistogram(Form("hParentPi0pn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
1136               break ;
1137             case -2212:  //pbar
1138             case -2112:  //nbar
1139               FillHistogram(Form("hParentPi0antipn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
1140               break ;
1141             default:       //other
1142               FillHistogram(Form("hParentPi0Other_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1143             }//switch     
1144           }//pi0 with primary
1145         }//common parent - pi0
1146       }//there is common primary 
1147     }//seond photon loop
1148   }//first photon loop
1149   
1150   
1151   //Now look at photon contaiminations
1152   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1153   
1154     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1155     Int_t iprim = ph1->GetPrimary() ;
1156     if(iprim<0)
1157       FillAllHistograms(Form("hGammaNoPrim_cen%d",fCentBin),ph1) ; //
1158     else{
1159       //Find primary at vertex
1160       AliAODMCParticle * primPHOS = GetParticle(iprim);
1161       Int_t iprimV = primPHOS->GetMother();
1162       AliAODMCParticle * primVtx = primPHOS ;
1163       while((iprimV>-1) && R(primVtx)>kRCut){
1164         primVtx = GetParticle(iprimV);
1165         iprimV = primVtx->GetMother();
1166       }
1167     
1168       //photon
1169       Int_t primPdgCode=primVtx->PdgCode() ;
1170       switch(primPdgCode){
1171       case  22: FillAllHistograms("hGammaPhot",ph1); 
1172         break ;
1173       case  11: 
1174       case -11: 
1175         FillAllHistograms("hGammaEl",ph1); 
1176         break ;
1177       case  111: 
1178         FillAllHistograms("hGammaPi0",ph1); 
1179         break ;
1180       case  221: 
1181         FillAllHistograms("hGammaEta",ph1); 
1182         break ;
1183       case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
1184         break ;
1185       case  211: 
1186       case -211: 
1187         FillAllHistograms("hGammaPipm",ph1); 
1188         //Find particle entered PHOS
1189         if(primVtx == primPHOS)
1190           FillAllHistograms("hGammaPipmDirect",ph1); 
1191         else{
1192           Int_t primPdgPHOS=primPHOS->PdgCode() ;
1193           if(primPdgPHOS==22){
1194             FillAllHistograms("hGammaPipmGamma",ph1); 
1195             FillHistogram("hPipmGammaConvR",ph1->Pt(),R(primPHOS));
1196             FillHistogram("hPipmGammaConvRZ",primPHOS->Zv(),R(primPHOS));
1197             break ;               
1198           }
1199           if(TMath::Abs(primPdgPHOS)==11){
1200             FillAllHistograms("hGammaPipmEl",ph1); 
1201             FillHistogram("hPipmElConvR",ph1->Pt(),R(primPHOS));
1202             break ;               
1203           }
1204           if(TMath::Abs(primPdgPHOS)==2212){
1205             FillAllHistograms("hGammaPipmp",ph1); 
1206             FillHistogram("hPipmNConvR",ph1->Pt(),R(primPHOS));
1207             break ;               
1208           }
1209           if(TMath::Abs(primPdgPHOS)==2112){
1210             FillAllHistograms("hGammaPipmn",ph1); 
1211             FillHistogram("hPipmNConvR",ph1->Pt(),R(primPHOS));
1212             break ;               
1213           }
1214           FillAllHistograms("hGammaPipmOther",ph1); 
1215           FillHistogram("hPipmOtherConvR",ph1->Pt(),R(primPHOS));                   
1216         }
1217         break ;           
1218       case  2212:  //p 
1219         FillAllHistograms("hGammaP",ph1); 
1220         break ;           
1221       case -2212:  //pbar
1222         FillAllHistograms("hGammaPbar",ph1); 
1223         break ;           
1224       case  2112:  //n 
1225         FillAllHistograms("hGammaN",ph1); 
1226         break ;           
1227       case -2112:  //nbar
1228         FillAllHistograms("hGammaNbar",ph1) ; // pn
1229         break ;
1230       case  310:  //nbar
1231         FillAllHistograms("hGammaK0S",ph1) ; // pn
1232         break ;
1233       case  130:  //nbar
1234         FillAllHistograms("hGammaK0L",ph1) ; // pn
1235         break ;
1236       case  321:  //K+
1237       case -321:  //K-
1238         FillAllHistograms("hGammaKpm",ph1) ; // pn
1239         break ;
1240       case -323: 
1241       case  323: 
1242       case -313: 
1243       case  313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
1244         break ;
1245                   
1246       case -2224 : //Deltas
1247       case  2224 : //Deltas
1248       case -2214 : //Deltas
1249       case  2214 : //Deltas
1250       case -2114 : //Deltas
1251       case  2114 : //Deltas
1252       case -1114 : //Deltas
1253       case  1114 : //Deltas
1254         FillAllHistograms("hGammaDelta",ph1) ; // pn
1255         break ;           
1256       default:     //other
1257         if(primVtx->Charge())
1258           FillAllHistograms("hGammaOtherCharged",ph1) ; //
1259         else
1260           FillAllHistograms("hGammaOtherNeutral",ph1) ; //
1261       }
1262     }
1263   
1264   }//single photons
1265   
1266 }
1267
1268 //_____________________________________________________________________________
1269 void AliAnalysisTaskPi0FlowMCAOD::FillAllHistograms(const char * particleName,AliCaloPhoton * ph)
1270 {
1271   //Fill All PID histograms
1272         
1273   Double_t w=ph->GetWeight() ;
1274   Double_t pt = ph->Pt() ;
1275   Double_t ptC=ph->GetMomV2()->Pt() ;
1276   FillHistogram(Form("%s_All_cen%d",particleName,fCentBin),pt,w) ;
1277   FillHistogram(Form("%s_Allcore_cen%d",particleName,fCentBin),ptC,w) ;
1278   if(ph->IsCPVOK()){
1279     FillHistogram(Form("%s_CPV_cen%d",particleName,fCentBin),pt,w) ;
1280     FillHistogram(Form("%s_CPVcore_cen%d",particleName,fCentBin),ptC,w) ;
1281   }
1282   if(ph->IsCPV2OK()){
1283     FillHistogram(Form("%s_CPV2_cen%d",particleName,fCentBin),pt,w) ;
1284     FillHistogram(Form("%s_CPV2core_cen%d",particleName,fCentBin),ptC,w) ;
1285   }
1286   if(ph->IsDispOK()){     
1287     FillHistogram(Form("%s_Disp_cen%d",particleName,fCentBin),pt,w) ;
1288     FillHistogram(Form("%s_Dispcore_cen%d",particleName,fCentBin),ptC,w) ;
1289     if(ph->IsDisp2OK()){
1290       FillHistogram(Form("%s_Disp2_cen%d",particleName,fCentBin),pt,w) ;
1291       FillHistogram(Form("%s_Disp2core_cen%d",particleName,fCentBin),ptC,w) ;
1292     }
1293     if(ph->IsCPVOK()){
1294       FillHistogram(Form("%s_Both_cen%d",particleName,fCentBin),pt,w) ;
1295       FillHistogram(Form("%s_Bothcore_cen%d",particleName,fCentBin),ptC,w) ;
1296     }
1297   }  
1298 }
1299
1300
1301 //___________________________________________________________________________
1302 Double_t AliAnalysisTaskPi0FlowMCAOD::PrimaryWeight(Int_t /* primary */){
1303   return 1.;
1304 }
1305 //________________________________________________________________________
1306 Double_t AliAnalysisTaskPi0FlowMCAOD::PrimaryParticleWeight(AliAODMCParticle * /* particle */){
1307   return 1.;
1308 }
1309
1310 //________________________________________________________________________
1311 TClonesArray* AliAnalysisTaskPi0FlowMCAOD::GetMCArray()
1312 {
1313   fMcArray = 0;
1314   AliAODInputHandler* aodHandler=dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
1315   
1316   if (aodHandler){
1317     AliAODEvent *aod=aodHandler->GetEvent();
1318     if (aod) {
1319       fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1320       if (!fMcArray) AliError("Could not retrieve MC array!");
1321     }
1322     else AliError("Could not retrieve AOD event!");
1323   }
1324   
1325   return fMcArray;
1326 }
1327
1328 AliAODMCParticle* AliAnalysisTaskPi0FlowMCAOD::GetParticle(Int_t particlepos) {
1329   //Returns particle at given position.
1330   
1331   if(!fMcArray){
1332     AliError("MC array is not initialized, run GetMCArray() first!");
1333     return 0;
1334   }
1335
1336   return (AliAODMCParticle *) fMcArray->At(particlepos);
1337 }
1338
1339
1340 Int_t AliAnalysisTaskPi0FlowMCAOD::FindPrimary(AliVCluster*clu,  Bool_t&sure){
1341   //Finds primary and estimates if it unique one?
1342   //First check can it be photon/electron
1343
1344   const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
1345   Int_t n = clu->GetNLabels();
1346
1347   for(Int_t i=0;  i<n;  i++){
1348     Int_t label = clu->GetLabelAt(i);
1349     AliAODMCParticle*  p=  GetParticle(label) ;
1350     Int_t pdg = p->PdgCode() ;
1351     if(pdg==22  ||  pdg==11 || pdg == -11){
1352       if(p->E()>emFraction*clu->E()){
1353         sure=kTRUE ;
1354         return label;
1355       }
1356     }
1357   }
1358
1359   Double_t*  Ekin=  new  Double_t[n] ;
1360
1361   for(Int_t i=0;  i<n;  i++){
1362     AliAODMCParticle*  p = GetParticle(clu->GetLabelAt(i)) ;
1363     Ekin[i]=p->P() ;  // estimate of kinetic energy
1364     if(p->PdgCode()==-2212  ||  p->PdgCode()==-2112){
1365       Ekin[i]+=1.8  ;  //due to annihilation
1366     }
1367   }
1368   Int_t iMax=0;
1369   Double_t eMax=0.,eSubMax=0. ;
1370   for(Int_t i=0;  i<n;  i++){
1371     if(Ekin[i]>eMax){
1372       eSubMax=eMax;
1373       eMax=Ekin[i];
1374       iMax=i;
1375     }
1376   }
1377   if(eSubMax>0.8*eMax)//not obvious primary
1378     sure=kFALSE;
1379   else
1380     sure=kTRUE;
1381   delete[]  Ekin;
1382   return clu->GetLabelAt(iMax);
1383 }
1384
1385 //________________________________________________________________________
1386 Int_t AliAnalysisTaskPi0FlowMCAOD::FindCommonParent(Int_t iPart, Int_t jPart){
1387   //check if there is a common parent for particles i and j
1388   // -1: no common parent or wrong iPart/jPart
1389   
1390   Int_t ntrack = fMcArray->GetEntriesFast();
1391   if(iPart==-1 || iPart>=ntrack || jPart==-1 || jPart>=ntrack) return -1;
1392   
1393   Int_t iprim1 = iPart;
1394   
1395   while(iprim1>-1){  
1396     Int_t iprim2=jPart;
1397     
1398     while(iprim2>-1){
1399       if(iprim1==iprim2)
1400         return iprim1 ;
1401       iprim2 = GetParticle(iprim2)->GetMother();
1402     }
1403     
1404     iprim1 = GetParticle(iprim1)->GetMother();
1405   }
1406   
1407   return -1;
1408 }
1409 //________________________________________________________________________
1410 Bool_t AliAnalysisTaskPi0FlowMCAOD::HaveParent(Int_t iPart, Int_t pdgParent){
1411   //check if there is a common parent for particles i and j
1412   // -1: no common parent or wrong iPart/jPart
1413   
1414   Int_t ntrack = fMcArray->GetEntriesFast();
1415   if(iPart==-1 || iPart>=ntrack) return -1;
1416   
1417   Int_t iprim1=iPart;
1418   while(iprim1>-1){  
1419     AliAODMCParticle * tmp = GetParticle(iprim1) ;
1420     if(tmp->PdgCode()==pdgParent)
1421       return kTRUE ;
1422     iprim1 = tmp->GetMother();
1423   }
1424   return kFALSE;
1425 }
1426 //________________________________________________________________________
1427 Bool_t AliAnalysisTaskPi0FlowMCAOD::InPi0mass(Double_t m, Double_t /*pt*/){
1428
1429  return TMath::Abs(m-0.135)<0.007*2.5 ;
1430 }
1431
1432 //________________________________________________________________________
1433 Double32_t AliAnalysisTaskPi0FlowMCAOD::R(AliAODMCParticle* p) {
1434   //Radius of vertex.
1435
1436   Double32_t x = p->Xv();
1437   Double32_t y = p->Yv();
1438   Double32_t z = p->Zv();
1439
1440   return x*x + y*y + z*z;
1441 }