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