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