PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0FlowMC.cxx

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