PWGGA/GammaConv/AliAnalysisTaskPi0v2.cxx

   1 #include <exception>
   2 #include "TRandom3.h"
   3 #include "TChain.h"
   4 #include "AliLog.h"
   5 #include "TList.h"
   6 #include "TChain.h"
   7 #include "TDirectory.h"
   8 #include "TTree.h"
   9 #include "TH1.h"
  10 #include "TH1F.h"
  11 #include "THnSparse.h"
  12 #include "AliAnalysisTask.h"
  13 #include "AliAnalysisManager.h"
  14 #include "AliAnalysisTaskPi0v2.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliAODEvent.h"
  18 #include "AliCentrality.h"
  19 #include <iostream>
  20
  21 #include "TFile.h"
  22 #include "AliOADBContainer.h"
  23
  24 #include "AliEPSelectionTask.h"
  25
  26 // Author Daniel Lohner (Daniel.Lohner@cern.ch)
  27
  28 using namespace std;
  29
  30 ClassImp(AliAnalysisTaskPi0v2)
  31
  32 //________________________________________________________________________
  33     AliAnalysisTaskPi0v2::AliAnalysisTaskPi0v2(const char *name,Int_t harmonic) : AliAnalysisTaskSE(name),
  34     fV0Reader(NULL),
  35     fNCuts(0),
  36     fConversionSelection(NULL),
  37     fConversionGammas(NULL),
  38     fNCentralityBins(1),
  39     fCentralityBins(NULL),
  40     fCentrality(-1),
  41     fCentralityBin(0),
  42     fNBinsPhi(6),
  43     fEP(NULL),
  44     fUseTPCOnlyTracks(kTRUE),
  45     fEtaMax(0.75),
  46     fEtaGap(1),
  47     fRPTPCEtaA(0),
  48     fRPTPCEtaC(0),
  49     fRPV0A(0),
  50     fRPV0C(0),
  51     fRPTPC(0),
  52     fRPTPCEtaABF(0),
  53     fRPTPCEtaCBF(0),
  54     fRPV0ABF(0),
  55     fRPV0CBF(0),
  56     fRPTPCBF(0),
  57     fConversionCuts(NULL),
  58     fRandomizer(NULL),
  59     fOutputList(NULL),
  60     fMesonPDGCode(kPi0),
  61     fDeltaPsiRP(0),
  62     fRunNumber(0),
  63     fRunIndex(0),
  64     fNEPMethods(knEPMethod),
  65     fFillQA(kTRUE),
  66     fHarmonic(harmonic),
  67     fPsiMax(2*TMath::Pi()/Double_t(harmonic)),
  68     fPeriod("LHC10h"),
  69     fIsAOD(kFALSE),
  70     fSparseDist(NULL),
  71     fHruns(NULL),
  72     fDoEPFlattening(kTRUE),
  73
  74     hNEvents(NULL),
  75     hEventSelection(NULL),
  76     hRPTPC(NULL),
  77     hRPV0A(NULL),
  78     hRPV0C(NULL),
  79     hRPTPCAC(NULL),
  80     hRPV0ATPC(NULL),
  81     hRPV0CTPC(NULL),
  82     hRPV0AC(NULL),
  83     hCos2TPC(NULL),
  84     hCos2V0ATPC(NULL),
  85     hCos2V0CTPC(NULL),
  86     hCos2V0AC(NULL),
  87     hRPTPCEtaA(NULL),
  88     hRPTPCEtaC(NULL),
  89     hRPTPCEtaAC(NULL),
  90     hCos2TPCEta(NULL),
  91     hCos2V0ATPCEtaA(NULL),
  92     hCos2V0ATPCEtaC(NULL),
  93     hCos2V0CTPCEtaA(NULL),
  94     hCos2V0CTPCEtaC(NULL),
  95     hCos2SumWeights(NULL),
  96     hEtaTPCEP(NULL),
  97     hGammaMultCent(NULL),
  98     hGammaPhi(NULL),
  99     hMultChargedvsNGamma(NULL),
 100     hMultChargedvsVZERO(NULL),
 101     hMultChargedvsSPD(NULL),
 102     hGammadNdPhi(NULL),
 103     hGammaMultdPhiTRUE(NULL),
 104     hGammaMultdPhiRECOTRUE(NULL),
 105     hGammaMultTRUE(NULL),
 106     hGammaMultRECOTRUE(NULL),
 107     hGammaMultdPhi(NULL),
 108     hGammaMult(NULL),
 109     hGamma(NULL),
 110     hGammaFull(NULL),
 111     hCharged(NULL),
 112     hPi0(NULL),
 113     hPi0BG(NULL),
 114
 115     fMultV0(0x0),
 116     fV0Cpol(0.),
 117     fV0Apol(0.),
 118     //hEPVertex(NULL)
 119     hEPQA(NULL)
 120 {
 121     fInvMassRange[0]=0.05;
 122     fInvMassRange[1]=0.3;
 123
 124     for(Int_t ii=0;ii<knEPMethod;ii++)fEPSelectionMask[ii]=1;
 125
 126     fRandomizer=new TRandom3();
 127     fRandomizer->SetSeed(0);
 128
 129     for(Int_t i = 0; i < 4; ++i) {
 130         fPhiDist[i] = 0;
 131     }
 132
 133     // Define input and output slots here
 134     DefineInput(0, TChain::Class());
 135     DefineOutput(1, TList::Class());
 136 }
 137
 138 //________________________________________________________________________
 139 AliAnalysisTaskPi0v2::~AliAnalysisTaskPi0v2(){
 140
 141     if(fRandomizer){
 142         delete fRandomizer;
 143         fRandomizer=NULL;
 144     }
 145     if(fCentralityBins){
 146         delete fCentralityBins;
 147         fCentralityBins=NULL;
 148     }
 149
 150     if(fConversionSelection){
 151         for(Int_t ii=0;ii<fNCuts;ii++)delete fConversionSelection[ii];
 152         delete[] fConversionSelection;
 153         fConversionSelection=NULL;
 154     }
 155
 156     if (fPeriod.CompareTo("LHC11h")==0){
 157         for(Int_t i = 0; i < 4; i++) {
 158             if(fPhiDist[i]){
 159                 delete fPhiDist[i];
 160                 fPhiDist[i] = 0;
 161             }
 162         }
 163         if(fHruns) delete fHruns;
 164     }
 165 }
 166
 167 //________________________________________________________________________
 168 void AliAnalysisTaskPi0v2::UserCreateOutputObjects()
 169 {
 170     OpenFile(1);
 171
 172     // GetConversionCuts
 173     fConversionCuts=fV0Reader->GetConversionCuts();
 174
 175     // Flags
 176
 177     Bool_t IsMC=AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
 178     Bool_t IsHeavyIon=fConversionCuts->IsHeavyIon();
 179
 180     if(!IsHeavyIon||IsMC)fNEPMethods=1;
 181
 182     if(!fCentralityBins){
 183         fCentralityBins=new Double_t[fNCentralityBins+1];
 184         fCentralityBins[0]=-0.5;
 185         fCentralityBins[1]=0.5;
 186     }
 187
 188     // Create histograms
 189
 190     if(fOutputList != NULL){
 191         delete fOutputList;
 192         fOutputList = NULL;
 193     }
 194     if(fOutputList == NULL){
 195         fOutputList = new TList();
 196         fOutputList->SetOwner(kTRUE);
 197     }
 198
 199     Int_t kGCnXBinsSpectra = Int_t((fInvMassRange[1]-fInvMassRange[0])*500);  //500 for range 0 - 1
 200     Double_t kGCfirstXBinSpectra = fInvMassRange[0];
 201     Double_t kGClastXBinSpectra = fInvMassRange[1];
 202
 203     Int_t nbinspi0[knbinsPi0]={kGCnYBinsSpectra,kGCnXBinsSpectra,fNBinsPhi,fNCentralityBins,fNEPMethods};
 204     Double_t minpi0[knbinsPi0]={kGCfirstYBinSpectra,kGCfirstXBinSpectra,0,-0.5,-0.5};
 205     Double_t maxpi0[knbinsPi0]={kGClastYBinSpectra,kGClastXBinSpectra,fPsiMax/2.,fNCentralityBins-0.5,fNEPMethods-0.5};
 206     const char *binpi0[knbinsPi0]={"pt","mass","dPhi","centr","EPm"};
 207
 208     Int_t nbinsg[knbinsGamma]={kGCnYBinsSpectra,fNBinsPhi,fNCentralityBins,fNEPMethods};
 209     Double_t ming[knbinsGamma]={kGCfirstYBinSpectra,0,-0.5,-0.5};
 210     Double_t maxg[knbinsGamma]={kGClastYBinSpectra,fPsiMax/2.,fNCentralityBins-0.5,fNEPMethods-0.5};
 211     const char *bingamma[knbinsGamma]={"pt","dPhi","centr","EPm"};
 212
 213     // Define Binning
 214
 215     if(!IsMC){
 216
 217         hPi0=new THnSparseF*[fNCuts+1];
 218         hPi0BG=new THnSparseF*[fNCuts+1];
 219         hGamma=new THnSparseF*[fNCuts+1];
 220
 221         // Photon Cuts
 222         for(Int_t iCut=0;iCut<fNCuts;iCut++){
 223
 224             TList *fCutOutputList=new TList();
 225             fCutOutputList->SetName(fConversionSelection[iCut]->GetCutString().Data());
 226             fCutOutputList->SetOwner(kTRUE);
 227             fOutputList->Add(fCutOutputList);
 228
 229             hPi0[iCut]=new THnSparseF("Pi0_Sparse","Pi0_Sparse",knbinsPi0,nbinspi0,minpi0,maxpi0);
 230             for(Int_t i=0;i<knbinsPi0;i++) hPi0[iCut]->GetAxis(i)->SetName(binpi0[i]);
 231             hPi0[iCut]->Sumw2();
 232             fCutOutputList->Add(hPi0[iCut]);
 233
 234             hPi0BG[iCut]=new THnSparseF("Pi0BG_Sparse","Pi0BG_Sparse",knbinsPi0,nbinspi0,minpi0,maxpi0);
 235             for(Int_t i=0;i<knbinsPi0;i++) hPi0BG[iCut]->GetAxis(i)->SetName(binpi0[i]);
 236             hPi0BG[iCut]->Sumw2();
 237             fCutOutputList->Add(hPi0BG[iCut]);
 238
 239             // Gamma
 240
 241             hGamma[iCut]=new THnSparseF("Gamma_Sparse","Gamma_Sparse",knbinsGamma,nbinsg,ming,maxg);
 242             for(Int_t i=0;i<knbinsGamma;i++) hGamma[iCut]->GetAxis(i)->SetName(bingamma[i]);
 243             hGamma[iCut]->Sumw2();
 244             fCutOutputList->Add(hGamma[iCut]);
 245         }
 246
 247         // no EP Flattening
 248         Int_t iCut=0;
 249
 250         TList *fCutOutputList=new TList();
 251         fCutOutputList->SetName(Form("%s_BF",fConversionSelection[iCut]->GetCutString().Data()));
 252         fCutOutputList->SetOwner(kTRUE);
 253         fOutputList->Add(fCutOutputList);
 254
 255         iCut=fNCuts;
 256
 257         hPi0[iCut]=new THnSparseF("Pi0_Sparse","Pi0_Sparse",knbinsPi0,nbinspi0,minpi0,maxpi0);
 258         for(Int_t i=0;i<knbinsPi0;i++) hPi0[iCut]->GetAxis(i)->SetName(binpi0[i]);
 259         hPi0[iCut]->Sumw2();
 260         fCutOutputList->Add(hPi0[iCut]);
 261
 262         hPi0BG[iCut]=new THnSparseF("Pi0BG_Sparse","Pi0BG_Sparse",knbinsPi0,nbinspi0,minpi0,maxpi0);
 263         for(Int_t i=0;i<knbinsPi0;i++) hPi0BG[iCut]->GetAxis(i)->SetName(binpi0[i]);
 264         hPi0BG[iCut]->Sumw2();
 265         fCutOutputList->Add(hPi0BG[iCut]);
 266
 267         // Gamma
 268
 269         hGamma[iCut]=new THnSparseF("Gamma_Sparse","Gamma_Sparse",knbinsGamma,nbinsg,ming,maxg);
 270         for(Int_t i=0;i<knbinsGamma;i++) hGamma[iCut]->GetAxis(i)->SetName(bingamma[i]);
 271         hGamma[iCut]->Sumw2();
 272         fCutOutputList->Add(hGamma[iCut]);
 273
 274     }
 275
 276     if(IsHeavyIon&&!IsMC){
 277
 278         // RP Calculation
 279         TList *fRPList=new TList();
 280         fRPList->SetName("Event Plane");
 281         fRPList->SetOwner(kTRUE);
 282         fOutputList->Add(fRPList);
 283
 284         hRPTPC=new TH2F("TPCAC" ,"TPC_AC" , fNCentralityBins,fCentralityBins, 180, 0, fPsiMax);
 285         hRPTPC->Sumw2();
 286         fRPList->Add(hRPTPC);
 287         hRPTPCEtaA=new TH2F("TPCEtaA" ,"TPC_EtaA" , fNCentralityBins,fCentralityBins, 180, 0, fPsiMax);
 288         hRPTPCEtaA->Sumw2();
 289         fRPList->Add(hRPTPCEtaA);
 290         hRPTPCEtaC=new TH2F("TPCEtaC" ,"TPC_EtaC" , fNCentralityBins,fCentralityBins, 180, 0, fPsiMax);
 291         hRPTPCEtaC->Sumw2();
 292         fRPList->Add(hRPTPCEtaC);
 293         hRPV0A=new TH2F("V0A" ,"VZERO_A" , fNCentralityBins,fCentralityBins, 180, 0, fPsiMax);
 294         hRPV0A->Sumw2();
 295         fRPList->Add(hRPV0A);
 296         hRPV0C=new TH2F("V0C" ,"VZERO_C" , fNCentralityBins,fCentralityBins, 180, 0, fPsiMax);
 297         hRPV0C->Sumw2();
 298         fRPList->Add(hRPV0C);
 299
 300         hRPTPCAC=new TH2F("TPCA_TPCC" ,"TPCA_TPCC" , 180, 0, fPsiMax, 180, 0, fPsiMax);
 301         hRPTPCAC->Sumw2();
 302         fRPList->Add(hRPTPCAC);
 303         hRPV0ATPC=new TH2F("V0A_TPC" ,"V0A_TPC" , 180, 0, fPsiMax, 180, 0, fPsiMax);
 304         hRPV0ATPC->Sumw2();
 305         fRPList->Add(hRPV0ATPC);
 306         hRPV0CTPC=new TH2F("V0C_TPC" ,"V0C_TPC" , 180, 0, fPsiMax, 180, 0, fPsiMax);
 307         hRPV0CTPC->Sumw2();
 308         fRPList->Add(hRPV0CTPC);
 309         hRPV0AC=new TH2F("V0A_V0C" ,"V0A_V0C" , 180, 0, fPsiMax, 180, 0, fPsiMax);
 310         hRPV0AC->Sumw2();
 311         fRPList->Add(hRPV0AC);
 312         hRPTPCEtaAC=new TH2F("TPCEtaA_TPCEtaC" ,"TPCEtaA_TPCEtaC" , 180, 0, fPsiMax, 180, 0, fPsiMax);
 313         hRPTPCEtaAC->Sumw2();
 314         fRPList->Add(hRPTPCEtaAC);
 315
 316         Int_t nsystep=4;// 3 different weights + unflattened EP
 317
 318         hCos2TPC=new TH2F("Cos2_TPCAC" ,"Cos2_TPCAC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 319         hCos2TPC->Sumw2();
 320         fRPList->Add(hCos2TPC);
 321         hCos2TPCEta=new TH2F("Cos2_TPCEtaAC" ,"Cos2_TPCEtaAC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 322         hCos2TPCEta->Sumw2();
 323         fRPList->Add(hCos2TPCEta);
 324         hCos2V0ATPC=new TH2F("Cos2_V0ATPC" ,"Cos2_V0ATPC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 325         hCos2V0ATPC->Sumw2();
 326         fRPList->Add(hCos2V0ATPC);
 327         hCos2V0CTPC=new TH2F("Cos2_V0CTPC" ,"Cos2_V0CTPC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 328         hCos2V0CTPC->Sumw2();
 329         fRPList->Add(hCos2V0CTPC);
 330         hCos2V0AC=new TH2F("Cos2_V0AC" ,"Cos2_V0AC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 331         hCos2V0AC->Sumw2();
 332         fRPList->Add(hCos2V0AC);
 333         hCos2V0ATPCEtaA=new TH2F("Cos2_V0ATPCEtaA" ,"Cos2_V0ATPCEtaA" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 334         hCos2V0ATPCEtaA->Sumw2();
 335         fRPList->Add(hCos2V0ATPCEtaA);
 336         hCos2V0ATPCEtaC=new TH2F("Cos2_V0ATPCEtaC" ,"Cos2_V0ATPCEtaC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 337         hCos2V0ATPCEtaC->Sumw2();
 338         fRPList->Add(hCos2V0ATPCEtaC);
 339         hCos2V0CTPCEtaA=new TH2F("Cos2_V0CTPCEtaA" ,"Cos2_V0CTPCEtaA" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 340         hCos2V0CTPCEtaA->Sumw2();
 341         fRPList->Add(hCos2V0CTPCEtaA);
 342         hCos2V0CTPCEtaC=new TH2F("Cos2_V0CTPCEtaC" ,"Cos2_V0CTPCEtaC" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 343         hCos2V0CTPCEtaC->Sumw2();
 344         fRPList->Add(hCos2V0CTPCEtaC);
 345         hCos2SumWeights=new TH2F("Cos2_SumWeights" ,"Cos2_SumWeights" ,fNCentralityBins,fCentralityBins,nsystep+1,-0.5,nsystep+0.5);
 346         hCos2SumWeights->Sumw2();
 347         fRPList->Add(hCos2SumWeights);
 348
 349         hEtaTPCEP=new TH2F("Eta_TPCEP" ,"Eta_TPCEP" ,fNCentralityBins,fCentralityBins,100,-1,1);
 350         hEtaTPCEP->Sumw2();
 351         fRPList->Add(hEtaTPCEP);
 352
 353
 354         /*const Int_t nepbins=4;
 355          Int_t nbinsep[nepbins]={30,30,40,180};
 356          Double_t minep[nepbins]={-0.015,0.17,-10,0};
 357          Double_t maxep[nepbins]={0.015,0.20,10,fPsiMax};
 358
 359          hEPVertex=new THnSparseF("EP_Vertex","EP_Vertex",nepbins,nbinsep,minep,maxep);
 360          fRPList->Add(hEPVertex);
 361          */
 362
 363         const Int_t nRuns=275;
 364         const Int_t nepbins=4;
 365         Int_t nbinsep[nepbins]={fNCentralityBins,180,nRuns,5};
 366         Double_t minep[nepbins]={-0.5,0,0.5,-0.5};
 367         Double_t maxep[nepbins]={fNCentralityBins-0.5,fPsiMax,Double_t(nRuns)+0.5,4.5};
 368         hEPQA=new THnSparseF("EP_QA","EP_QA",nepbins,nbinsep,minep,maxep);
 369         fRPList->Add(hEPQA);
 370     }
 371
 372     TList *fPhotonQAList=new TList();
 373     fPhotonQAList->SetName("Gamma_QA");
 374     fPhotonQAList->SetOwner(kTRUE);
 375     fOutputList->Add(fPhotonQAList);
 376
 377     if(fFillQA){
 378         // Gamma QA
 379         hGammaPhi=new TH2F*[fNCentralityBins];
 380         for(Int_t m=0;m<fNCentralityBins;m++){
 381             hGammaPhi[m]=new TH2F(Form("%d_GammaPhi",m),"GammaPhi",kGCnYBinsSpectra,kGCfirstYBinSpectra,kGClastYBinSpectra,360,0,2*TMath::Pi());
 382             hGammaPhi[m]->Sumw2();
 383             fPhotonQAList->Add(hGammaPhi[m]);
 384         }
 385
 386         hGammaMultCent=new TH2F("GammaMultvsCent","GammaMultvsCent",fNCentralityBins,fCentralityBins, 60,-0.5,59.5);
 387         hGammaMultCent->Sumw2();
 388         fPhotonQAList->Add(hGammaMultCent);
 389
 390         hMultChargedvsSPD=new TH2F("Mult_ChargedvsSPD","Mult_ChargedvsSPD",250,0,2500, 250,0,5000);
 391         hMultChargedvsSPD->Sumw2();
 392         fPhotonQAList->Add(hMultChargedvsSPD);
 393         hMultChargedvsVZERO=new TH2F("Mult_ChargedvsVZERO","Mult_ChargedvsVZERO",250,0,2500, 200,0,20000);
 394         hMultChargedvsVZERO->Sumw2();
 395         fPhotonQAList->Add(hMultChargedvsVZERO);
 396         hMultChargedvsNGamma=new TH2F("Mult_ChargedvsNGamma","Mult_ChargedvsNGamma",250,0,2500,60,-0.5,59.5);
 397         hMultChargedvsNGamma->Sumw2();
 398         fPhotonQAList->Add(hMultChargedvsNGamma);
 399
 400         Int_t nbinsgmult[knbinsGammaMult]={kGCnYBinsSpectra,400,fNCentralityBins};
 401         Double_t mingmult[knbinsGammaMult]={kGCfirstYBinSpectra,0,-0.5};
 402         Double_t maxgmult[knbinsGammaMult]={kGClastYBinSpectra,8000,fNCentralityBins-0.5};
 403         Double_t maxgmultdPhi[knbinsGammaMult]={kGClastYBinSpectra,2000,fNCentralityBins-0.5};
 404         const char *bingammamult[knbinsGammaMult]={"pt","gmult","centr"};
 405
 406         if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 407
 408             hGammaMultdPhiTRUE=new THnSparseF("Gamma_MultdPhi_TRUE","Gamma_MultdPhi_TRUE",knbinsGammaMult,nbinsgmult,mingmult,maxgmultdPhi);
 409             for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMultdPhiTRUE->GetAxis(i)->SetName(bingammamult[i]);
 410             hGammaMultdPhiTRUE->Sumw2();
 411             fPhotonQAList->Add(hGammaMultdPhiTRUE);
 412
 413         hGammaMultTRUE=new THnSparseF("Gamma_Mult_TRUE","Gamma_Mult_TRUE",knbinsGammaMult,nbinsgmult,mingmult,maxgmult);
 414         for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMultTRUE->GetAxis(i)->SetName(bingammamult[i]);
 415         hGammaMultTRUE->Sumw2();
 416         fPhotonQAList->Add(hGammaMultTRUE);
 417
 418         hGammaMultdPhiRECOTRUE=new THnSparseF("Gamma_MultdPhi_RECOTRUE","Gamma_MultdPhi_RECOTRUE",knbinsGammaMult,nbinsgmult,mingmult,maxgmultdPhi);
 419         for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMultdPhiRECOTRUE->GetAxis(i)->SetName(bingammamult[i]);
 420         hGammaMultdPhiRECOTRUE->Sumw2();
 421         fPhotonQAList->Add(hGammaMultdPhiRECOTRUE);
 422
 423         hGammaMultRECOTRUE=new THnSparseF("Gamma_Mult_RECOTRUE","Gamma_Mult_RECOTRUE",knbinsGammaMult,nbinsgmult,mingmult,maxgmult);
 424         for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMultRECOTRUE->GetAxis(i)->SetName(bingammamult[i]);
 425         hGammaMultRECOTRUE->Sumw2();
 426         fPhotonQAList->Add(hGammaMultRECOTRUE);
 427         }
 428
 429         hGammaMultdPhi=new THnSparseF*[fNEPMethods];
 430         hGammaMult=new THnSparseF*[fNEPMethods];
 431
 432         hGammadNdPhi=new THnSparseF("Gamma_dNdPhi","Gamma_dNdPhi",knbinsGamma,nbinsg,ming,maxg);
 433         for(Int_t i=0;i<knbinsGamma;i++) hGammadNdPhi->GetAxis(i)->SetName(bingamma[i]);
 434         hGammadNdPhi->Sumw2();
 435         fPhotonQAList->Add(hGammadNdPhi);
 436
 437         for(Int_t iEP=0;iEP<fNEPMethods;iEP++){
 438             hGammaMultdPhi[iEP]=new THnSparseF(Form("Gamma_MultdPhi_%d",iEP),"Gamma_MultdPhi",knbinsGammaMult,nbinsgmult,mingmult,maxgmultdPhi);
 439             for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMultdPhi[iEP]->GetAxis(i)->SetName(bingammamult[i]);
 440             hGammaMultdPhi[iEP]->Sumw2();
 441             fPhotonQAList->Add(hGammaMultdPhi[iEP]);
 442
 443             hGammaMult[iEP]=new THnSparseF(Form("Gamma_Mult_%d",iEP),"Gamma_Mult",knbinsGammaMult,nbinsgmult,mingmult,maxgmult);
 444             for(Int_t i=0;i<knbinsGammaMult;i++) hGammaMult[iEP]->GetAxis(i)->SetName(bingammamult[i]);
 445             hGammaMult[iEP]->Sumw2();
 446             fPhotonQAList->Add(hGammaMult[iEP]);
 447         }
 448
 449         const Int_t knbinsCharged=3;
 450         Int_t nbinscharged[knbinsCharged]={fNBinsPhi,fNCentralityBins,fNEPMethods};
 451         Double_t mincharged[knbinsCharged]={0,-0.5,-0.5};
 452         Double_t maxcharged[knbinsCharged]={fPsiMax/2,fNCentralityBins-0.5,fNEPMethods-0.5};
 453         hCharged=new THnSparseF("Charged","Charged",knbinsCharged,nbinscharged,mincharged,maxcharged);
 454         hCharged->GetAxis(0)->SetName("dPhi");
 455         hCharged->GetAxis(1)->SetName("centr");
 456         hCharged->GetAxis(2)->SetName("EPm");
 457         hCharged->Sumw2();
 458         fPhotonQAList->Add(hCharged);
 459
 460         Int_t nbinsgfull[knbinsGamma]={kGCnYBinsSpectra,24,fNCentralityBins,fNEPMethods};
 461         Double_t mingfull[knbinsGamma]={kGCfirstYBinSpectra,0,-0.5,-0.5};
 462         Double_t maxgfull[knbinsGamma]={kGClastYBinSpectra,2*TMath::Pi(),fNCentralityBins-0.5,fNEPMethods-0.5};
 463         hGammaFull=new THnSparseF("Gamma_Sparse_Full","Gamma_Sparse_Full",knbinsGamma,nbinsgfull,mingfull,maxgfull);
 464         for(Int_t i=0;i<knbinsGamma;i++) hGammaFull->GetAxis(i)->SetName(bingamma[i]);
 465         hGammaFull->Sumw2();
 466         fPhotonQAList->Add(hGammaFull);
 467
 468     }
 469     hNEvents=new TH1F("NEvents","NEvents",fNCentralityBins,fCentralityBins);
 470     fPhotonQAList->Add(hNEvents);
 471     hEventSelection=new TH1F("EventSelection","EventSelection",kEventSelected,0.5,kEventSelected+0.5);
 472     hEventSelection->GetXaxis()->SetBinLabel(kEventIn,"in");
 473     hEventSelection->GetXaxis()->SetBinLabel(kEventSelV0Reader,"v0reader");
 474     hEventSelection->GetXaxis()->SetBinLabel(kEventCentrality,"centr");
 475     hEventSelection->GetXaxis()->SetBinLabel(kEventRun,"run");
 476     hEventSelection->GetXaxis()->SetBinLabel(kEventNoTPCEP,"no ep");
 477     hEventSelection->GetXaxis()->SetBinLabel(kEventProcessEvent,"proc ev");
 478     hEventSelection->GetXaxis()->SetBinLabel(kEventSelected,"selected");
 479     fPhotonQAList->Add(hEventSelection);
 480
 481     // V0 Reader Cuts
 482     TList *fV0ReaderCuts=fConversionCuts->GetCutHistograms();
 483     fV0ReaderCuts->SetOwner(kTRUE);
 484     fOutputList->Add(fV0ReaderCuts);
 485
 486     PostData(1, fOutputList);
 487
 488 }
 489
 490 //________________________________________________________________________
 491 Bool_t AliAnalysisTaskPi0v2::InitEvent(){
 492
 493     if(!fV0Reader){AliError("Error: No V0 Reader and Pi0 Reconstructor");return kFALSE;}
 494     if(!fV0Reader->IsEventSelected()){
 495         hEventSelection->Fill(kEventSelV0Reader);
 496         return kFALSE;
 497     }
 498
 499     fConversionGammas=fV0Reader->GetReconstructedGammas();
 500
 501     fIsAOD=(fInputEvent->IsA()==AliAODEvent::Class());
 502
 503     if(!fConversionSelection){
 504         AliError("No Cut Selection");
 505         return kFALSE;
 506     }
 507
 508     if(!SetCentrality()){
 509         hEventSelection->Fill(kEventCentrality);
 510         return kFALSE;
 511     }
 512
 513     if(fConversionCuts->IsHeavyIon()&&!fMCEvent){
 514
 515         if(fRunNumber!=fInputEvent->GetRunNumber()){
 516             fRunNumber=fInputEvent->GetRunNumber();
 517             if (fRunNumber >= 136851 && fRunNumber <= 139515){fPeriod = "LHC10h";}
 518             if (fRunNumber >= 166529 && fRunNumber <= 170593){fPeriod = "LHC11h";}
 519             fRunIndex=GetRunIndex(fRunNumber);
 520             LoadVZEROCalibration(fRunNumber); // Load Calibration for V0 Event Plane
 521             LoadTPCCalibration(fRunNumber); // Load Calibration for TPC Event Plane
 522         }
 523         if(fRunIndex<0){
 524             AliInfo("Run not selected");
 525             hEventSelection->Fill(kEventRun);
 526             return kFALSE;
 527         }
 528
 529         // TPC Event Plane
 530         if(!GetTPCEventPlane()){
 531             hEventSelection->Fill(kEventNoTPCEP);
 532             return kFALSE;
 533         }
 534         //fEP=fInputEvent->GetEventplane();
 535         //if(!fEP)return kFALSE;
 536
 537         fRPTPCBF=GetCorrectedTPCEPAngle(NULL,NULL,kFALSE);
 538         fRPTPC=GetEventPlaneAngle(kTPC);
 539
 540         // TPC Eta Sub Events
 541         fRPTPCEtaABF=GetTPCSubEPEta(kEPTPCEtaA);
 542         fRPTPCEtaA=ApplyFlattening(fRPTPCEtaABF,kEPTPCEtaA);
 543
 544         fRPTPCEtaCBF=GetTPCSubEPEta(kEPTPCEtaC);
 545         fRPTPCEtaC=ApplyFlattening(fRPTPCEtaCBF,kEPTPCEtaC);
 546
 547         // GetV0 Event Plane
 548         GetV0EP(fInputEvent,fRPV0ABF,fRPV0CBF);
 549         fRPV0A=ApplyFlattening(fRPV0ABF,kEPV0A);
 550         fRPV0C=ApplyFlattening(fRPV0CBF,kEPV0C);
 551
 552     }
 553     return kTRUE;
 554 }
 555
 556
 557 //________________________________________________________________________
 558 void AliAnalysisTaskPi0v2::UserExec(Option_t *)
 559 {
 560     cout<<hNEvents->GetEntries()<<endl;
 561     if(hNEvents->GetEntries()>5)return;
 562
 563     hEventSelection->Fill(kEventIn);
 564
 565     if(!InitEvent())return;
 566
 567     // Process Cuts
 568     for(Int_t iCut=0;iCut<fNCuts;iCut++){
 569
 570         if(fConversionSelection[iCut]->ProcessEvent(fConversionGammas,fInputEvent,fMCEvent)){    // May only be called once per event!
 571
 572             if(!fMCEvent){
 573
 574                 // Process EP methods
 575                 for(Int_t iEP=0;iEP<fNEPMethods;iEP++){
 576
 577                     if(!fEPSelectionMask[iEP])continue; // dont fill THnSparse if not needed-> Save Memory
 578
 579                     ProcessPi0s(iCut,EEventPlaneMethod(iEP));
 580
 581                     ProcessGammas(iCut,EEventPlaneMethod(iEP));
 582                 }
 583             }
 584
 585             // QA
 586             if(fFillQA&&iCut==0)ProcessQA();
 587         }
 588         else{
 589             hEventSelection->Fill(kEventProcessEvent);
 590         }
 591     }
 592
 593     // Fill N Events
 594     hEventSelection->Fill(kEventSelected);
 595     hNEvents->Fill(fCentrality);
 596
 597     // EventPlaneResolution
 598     ProcessEventPlane();
 599
 600     PostData(1, fOutputList);
 601 }
 602
 603 //________________________________________________________________________
 604 void AliAnalysisTaskPi0v2::ProcessPi0s(Int_t iCut,EEventPlaneMethod iEP){
 605
 606     if(!fConversionSelection[iCut])return;
 607
 608     if(fConversionSelection[iCut]->GetNumberOfPhotons()==0)return;
 609
 610     // Process Pi0s
 611
 612     for(Int_t ii=0;ii<fConversionSelection[iCut]->GetNumberOfPi0s();ii++){
 613
 614         AliAODConversionMother *pi0cand=fConversionSelection[iCut]->GetPi0(ii);
 615
 616         if(!pi0cand)continue;
 617
 618         Double_t val[knbinsPi0];
 619         val[kPi0Pt]=pi0cand->Pt();
 620         val[kPi0Mass]=pi0cand->M();
 621         val[kPi0dPhi]=GetPi0PhiwrtRP(pi0cand,iEP);
 622         val[kPi0Cent]=fCentralityBin;
 623         val[kPi0EPM]=Int_t(iEP);
 624
 625         hPi0[iCut]->Fill(val);
 626
 627         if(iCut==0){
 628             // no flattening
 629             val[kPi0dPhi]=GetPi0PhiwrtRP(pi0cand,iEP,kFALSE);
 630             hPi0[fNCuts]->Fill(val);
 631         }
 632     }
 633
 634     // Pi0 BG
 635     for(Int_t ii=0;ii<fConversionSelection[iCut]->GetNumberOfBGs();ii++){
 636
 637         AliAODConversionMother *pi0cand=fConversionSelection[iCut]->GetBG(ii);
 638
 639         if(!pi0cand)continue;
 640
 641         Double_t val[knbinsPi0];
 642         val[kPi0Pt]=pi0cand->Pt();
 643         val[kPi0Mass]=pi0cand->M();
 644         val[kPi0dPhi]=GetPi0PhiwrtRP(pi0cand,iEP);
 645         val[kPi0Cent]=fCentralityBin;
 646         val[kPi0EPM]=Int_t(iEP);
 647
 648         hPi0BG[iCut]->Fill(val,pi0cand->GetWeight());
 649
 650         if(iCut==0){
 651             // no flattening
 652             val[kPi0dPhi]=GetPi0PhiwrtRP(pi0cand,iEP,kFALSE);
 653             hPi0BG[fNCuts]->Fill(val);
 654         }
 655     }
 656 }
 657
 658 //________________________________________________________________________
 659 void AliAnalysisTaskPi0v2::ProcessGammas(Int_t iCut,EEventPlaneMethod iEP){
 660
 661     if(!fConversionSelection[iCut]){
 662         AliWarning("Conversion Selection does not exist");
 663         return;
 664     }
 665
 666     for(Int_t ii=0;ii<fConversionSelection[iCut]->GetNumberOfPhotons();ii++){
 667
 668         AliAODConversionPhoton *gamma=fConversionSelection[iCut]->GetPhoton(ii);
 669
 670         Double_t val[knbinsGamma];
 671         val[kGammaPt]=gamma->Pt();
 672         val[kGammadPhi]=GetPhotonPhiwrtRP(gamma,iEP);
 673         val[kGammaCent]=fCentralityBin;
 674         val[kGammaEPM]=Int_t(iEP);
 675
 676         hGamma[iCut]->Fill(val);
 677
 678         if(iCut==0){
 679             // no flattening
 680             val[kGammadPhi]=GetPhotonPhiwrtRP(gamma,iEP,kFALSE);
 681             hGamma[fNCuts]->Fill(val);
 682         }
 683
 684         if(iCut==0&&fFillQA){
 685             hGammadNdPhi->Fill(val);
 686
 687             Double_t EPAngle=GetEventPlaneAngle(iEP,gamma->Eta(),gamma,NULL);
 688             Double_t dPhi=gamma->Phi()-EPAngle;
 689             if(dPhi>=(2*TMath::Pi()))dPhi-=2*TMath::Pi();
 690             if(dPhi<0)dPhi+=2*TMath::Pi();
 691             val[kGammadPhi]=dPhi;
 692             hGammaFull->Fill(val);
 693         }
 694     }
 695 }
 696
 697 //________________________________________________________________________
 698 void AliAnalysisTaskPi0v2::ProcessQA(){
 699
 700     if(!fConversionSelection[0]){
 701         AliWarning("Conversion Selection does not exist");
 702         return;
 703     }
 704
 705
 706     AliStack *fMCStack=NULL;
 707     if(fMCEvent)fMCStack=fMCEvent->Stack();
 708
 709     // Multiplicity
 710
 711     Double_t multcharged=fConversionSelection[0]->GetNumberOfChargedTracks(fInputEvent);
 712     Double_t multVZERO=fConversionSelection[0]->GetVZEROMult(fInputEvent);
 713     Double_t multSPD=fConversionSelection[0]->GetSPDMult(fInputEvent);
 714
 715     hMultChargedvsNGamma->Fill(multcharged,fConversionSelection[0]->GetNumberOfPhotons());
 716     hMultChargedvsVZERO->Fill(multcharged,multVZERO);
 717     hMultChargedvsSPD->Fill(multcharged,multSPD);
 718
 719     // Efficiency Purity
 720
 721     Double_t valdPhi[knbinsGammaMult];
 722     Double_t val[knbinsGammaMult];
 723
 724     Int_t dNdPhi[fNBinsPhi];
 725     Int_t ncharged;
 726
 727     for(Int_t iEP=0;iEP<fNEPMethods;iEP++){
 728         GetChargeddNdPhi(&dNdPhi[0],ncharged,iEP);
 729
 730         // Reco
 731         for(Int_t ii=0;ii<fConversionSelection[0]->GetNumberOfPhotons();ii++){
 732             AliAODConversionPhoton *gamma=fConversionSelection[0]->GetPhoton(ii);
 733             val[0]=gamma->Pt();
 734             val[1]=ncharged;
 735             val[2]=fCentralityBin;
 736
 737             valdPhi[0]=gamma->Pt();
 738             valdPhi[1]=dNdPhi[GetPhotonPhiBin(gamma,iEP)];
 739             valdPhi[2]=fCentralityBin;
 740
 741             hGammaMult[iEP]->Fill(val);
 742             hGammaMultdPhi[iEP]->Fill(valdPhi);
 743
 744             // Gamma Phi
 745             hGammaPhi[fCentralityBin]->Fill(gamma->Pt(),gamma->Phi());
 746             hGammaMultCent->Fill(fCentrality,Float_t(fConversionSelection[0]->GetNumberOfPhotons()));
 747
 748             if(fMCStack){
 749                 if(gamma->IsTruePhoton(fMCStack)){
 750                     hGammaMultRECOTRUE->Fill(val);
 751                     hGammaMultdPhiRECOTRUE->Fill(valdPhi);
 752                 }
 753             }
 754         }
 755
 756         // MC Truth
 757         if(fMCEvent){
 758             for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
 759                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
 760                 if (!particle) continue;
 761                 if(fConversionCuts->PhotonIsSelectedMC(particle,fMCStack)){
 762                     TParticle *daughter=(TParticle *)fMCStack->Particle(particle->GetDaughter(0));
 763                     if(daughter){
 764                         val[0]=particle->Pt();
 765                         val[1]=ncharged;
 766                         val[2]=fCentralityBin;
 767
 768                         valdPhi[0]=particle->Pt();
 769                         valdPhi[1]=dNdPhi[GetPhiBin(GetMCPhotonPhiwrtRP(particle,EEventPlaneMethod(iEP)))];
 770                         valdPhi[2]=fCentralityBin;
 771
 772                         hGammaMultTRUE->Fill(val);
 773                         hGammaMultdPhiTRUE->Fill(valdPhi);
 774                     }
 775                 }
 776             }
 777         }
 778     }
 779 }
 780
 781 //________________________________________________________________________
 782 void AliAnalysisTaskPi0v2::GetPhotondNdPhi(Int_t *dNdPhi,Int_t iEP,Int_t iCut){
 783
 784     for(Int_t iPhi=0;iPhi<fNBinsPhi;iPhi++)dNdPhi[iPhi]=0;
 785
 786     for(Int_t ii=0;ii<fConversionSelection[iCut]->GetNumberOfPhotons();ii++){
 787         AliAODConversionPhoton *gamma=fConversionSelection[iCut]->GetPhoton(ii);
 788         Int_t phibin=GetPhotonPhiBin(gamma,iEP);
 789         dNdPhi[phibin]++;
 790     }
 791 }
 792
 793 //________________________________________________________________________
 794 void AliAnalysisTaskPi0v2::GetChargeddNdPhi(Int_t *dNdPhi,Int_t &ntot,Int_t iEP){
 795
 796     for(Int_t iPhi=0;iPhi<fNBinsPhi;iPhi++)dNdPhi[iPhi]=0;
 797     ntot=0;
 798
 799     Double_t val[3];
 800     val[1]=fCentralityBin;
 801     val[2]=Int_t(iEP);
 802
 803     for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 804         AliVTrack* currentTrack = (AliVTrack*)fInputEvent->GetTrack(iTracks);
 805         if(!currentTrack) continue;
 806         if(TMath::Abs(currentTrack->Eta())>fEtaMax)continue;
 807
 808         Double_t phiwrt=GetChargedPhiwrtRP(currentTrack,EEventPlaneMethod(iEP));
 809         Int_t phibin=GetPhiBin(phiwrt);
 810
 811         val[0]=phiwrt;
 812         hCharged->Fill(val);
 813
 814         dNdPhi[phibin]++;
 815         ntot++;
 816     }
 817 }
 818
 819 //________________________________________________________________________
 820 Int_t AliAnalysisTaskPi0v2::GetPhiBin(Double_t phiwrt){
 821     Double_t binrange=TMath::Pi()/(Double_t(fHarmonic*fNBinsPhi));
 822     for(Int_t iPhi=0;iPhi<fNBinsPhi;iPhi++){
 823         if(phiwrt>=(binrange*iPhi)&&phiwrt<(binrange*(iPhi+1)))return iPhi;
 824     }
 825     return -1;
 826 }
 827
 828 //________________________________________________________________________
 829 Int_t AliAnalysisTaskPi0v2::GetPhotonPhiBin(AliAODConversionPhoton *gamma,Int_t iEP){
 830     Double_t phiwrt=GetPhotonPhiwrtRP(gamma,EEventPlaneMethod(iEP));
 831     return GetPhiBin(phiwrt);
 832 }
 833
 834 //________________________________________________________________________
 835 Double_t AliAnalysisTaskPi0v2::GetPi0PhiwrtRP(AliAODConversionMother *pi0,EEventPlaneMethod iEP,Bool_t bDoFlattening){
 836
 837     AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(pi0->GetLabel1()));
 838     AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(pi0->GetLabel2()));
 839
 840     Double_t EPAngle=GetEventPlaneAngle(iEP,pi0->Eta(),gamma0,gamma1,bDoFlattening);
 841
 842     return GetPhiwrtRP(pi0->Phi()-EPAngle);
 843 }
 844
 845 //________________________________________________________________________
 846 Double_t AliAnalysisTaskPi0v2::GetPhotonPhiwrtRP(AliAODConversionPhoton *gamma,EEventPlaneMethod iEP,Bool_t bDoFlattening){
 847
 848     Double_t EPAngle=GetEventPlaneAngle(iEP,gamma->Eta(),gamma,NULL,bDoFlattening);
 849
 850     return GetPhiwrtRP(gamma->Phi()-EPAngle);
 851 }
 852
 853 //________________________________________________________________________
 854 Double_t AliAnalysisTaskPi0v2::GetMCPhotonPhiwrtRP(TParticle *gamma,EEventPlaneMethod iEP,Bool_t bDoFlattening){
 855
 856     Double_t EPAngle=GetEventPlaneAngle(iEP,gamma->Eta(),NULL,NULL,bDoFlattening);
 857
 858     return GetPhiwrtRP(gamma->Phi()-EPAngle);
 859 }
 860 //________________________________________________________________________
 861 Double_t AliAnalysisTaskPi0v2::GetChargedPhiwrtRP(AliVTrack *track,EEventPlaneMethod iEP,Bool_t bDoFlattening){
 862
 863     Double_t EPAngle=GetEventPlaneAngle(iEP,track->Eta(),NULL,NULL,bDoFlattening);
 864
 865     return GetPhiwrtRP(track->Phi()-EPAngle);
 866 }
 867 //________________________________________________________________________
 868 Double_t AliAnalysisTaskPi0v2::GetPhiwrtRP(Double_t dPhi){
 869     Double_t newdPhi=TMath::Abs(dPhi); // Cos is symmetric
 870     while(newdPhi>fPsiMax/2)newdPhi=TMath::Abs(newdPhi-fPsiMax);
 871     return newdPhi;
 872 }
 873
 874 //________________________________________________________________________
 875 void AliAnalysisTaskPi0v2::Terminate(Option_t *)
 876 {
 877
 878 }
 879
 880 //________________________________________________________________________
 881 void AliAnalysisTaskPi0v2::ProcessEventPlane()
 882 {
 883     if(!fMCEvent&&fConversionCuts->IsHeavyIon()){
 884
 885         /*  Double_t val[4];
 886         val[0]=fInputEvent->GetPrimaryVertex()->GetX();
 887         val[1]=fInputEvent->GetPrimaryVertex()->GetY();
 888         val[2]=fInputEvent->GetPrimaryVertex()->GetZ();
 889         val[3]=GetEventPlaneAngle(kTPC);
 890         hEPVertex->Fill(val);   */
 891
 892         // Run by run monitoring (before flattening)
 893
 894         Double_t val[4];
 895         val[0]=fCentralityBin;
 896         val[2]=fRunIndex;
 897
 898         val[1]=fRPTPCBF;
 899         val[3]=kEPTPC;
 900         hEPQA->Fill(val);
 901
 902         val[1]=fRPTPCEtaABF;
 903         val[3]=kEPTPCEtaA;
 904         hEPQA->Fill(val);
 905         val[1]=fRPTPCEtaCBF;
 906         val[3]=kEPTPCEtaC;
 907         hEPQA->Fill(val);
 908
 909         val[1]=fRPV0ABF;
 910         val[3]=kEPV0A;
 911         hEPQA->Fill(val);
 912         val[1]=fRPV0CBF;
 913         val[3]=kEPV0C;
 914         hEPQA->Fill(val);
 915
 916         // After Flattening
 917
 918         // TPC EP
 919         Double_t PsiRP1BF=fEP->GetQsub1()->Phi()/Double_t(fHarmonic);
 920         Double_t PsiRP2BF=fEP->GetQsub2()->Phi()/Double_t(fHarmonic);
 921         Double_t PsiRP1=ApplyFlattening(PsiRP1BF,kEPTPC);
 922         Double_t PsiRP2=ApplyFlattening(PsiRP2BF,kEPTPC);
 923
 924         hRPTPC->Fill(fCentrality,fRPTPC);
 925         hRPTPCAC->Fill(PsiRP1,PsiRP2);
 926
 927         // TPC Eta Gap
 928         hRPTPCEtaA->Fill(fCentrality,fRPTPCEtaA);
 929         hRPTPCEtaC->Fill(fCentrality,fRPTPCEtaC);
 930         hRPTPCEtaAC->Fill(fRPTPCEtaA,fRPTPCEtaC);
 931
 932         // V0
 933
 934         hRPV0A->Fill(fCentrality,fRPV0A);
 935         hRPV0C->Fill(fCentrality,fRPV0C);
 936
 937         hRPV0ATPC->Fill(fRPV0A,fRPTPC);
 938         hRPV0CTPC->Fill(fRPV0C,fRPTPC);
 939         hRPV0AC->Fill(fRPV0A,fRPV0C);
 940
 941         Double_t cos2V0ATPC=TMath::Cos(Double_t(fHarmonic)*(fRPTPC-fRPV0A));
 942         Double_t cos2V0CTPC=TMath::Cos(Double_t(fHarmonic)*(fRPTPC-fRPV0C));
 943         Double_t cos2V0AV0C=TMath::Cos(Double_t(fHarmonic)*(fRPV0C-fRPV0A));
 944         Double_t cos2TPCEta=TMath::Cos(Double_t(fHarmonic)*(fRPTPCEtaA-fRPTPCEtaC));
 945         Double_t cos2TPC=TMath::Cos(Double_t(fHarmonic)*(PsiRP1-PsiRP2));
 946         Double_t cos2V0ATPCEtaA=TMath::Cos(Double_t(fHarmonic)*(fRPV0A-fRPTPCEtaA));
 947         Double_t cos2V0CTPCEtaA=TMath::Cos(Double_t(fHarmonic)*(fRPV0C-fRPTPCEtaA));
 948         Double_t cos2V0ATPCEtaC=TMath::Cos(Double_t(fHarmonic)*(fRPV0A-fRPTPCEtaC));
 949         Double_t cos2V0CTPCEtaC=TMath::Cos(Double_t(fHarmonic)*(fRPV0C-fRPTPCEtaC));
 950
 951         const Int_t nfill=4;
 952         Double_t weight[nfill];
 953         weight[0]=1.;// Fill unweighted
 954         weight[1]=Float_t(fConversionSelection[0]->GetNumberOfPhotons());// Weight with Photon Mult
 955         weight[2]=Float_t(fConversionSelection[0]->GetNumberOfChargedTracks(fInputEvent)); // Weight with charged Track Mult
 956         weight[3]=Float_t(fConversionSelection[0]->GetVZEROMult(fInputEvent)); // Weight with V0 mult
 957
 958         for(Int_t i=0;i<nfill;i++){
 959
 960             hCos2V0ATPC->Fill(fCentrality,i,cos2V0ATPC*weight[i]);
 961             hCos2V0CTPC->Fill(fCentrality,i,cos2V0CTPC*weight[i]);
 962             hCos2V0AC->Fill(fCentrality,i,cos2V0AV0C*weight[i]);
 963             hCos2TPCEta->Fill(fCentrality,i,cos2TPCEta*weight[i]);
 964             hCos2TPC->Fill(fCentrality,i,cos2TPC*weight[i]);
 965             hCos2V0ATPCEtaA->Fill(fCentrality,i,cos2V0ATPCEtaA*weight[i]);
 966             hCos2V0ATPCEtaC->Fill(fCentrality,i,cos2V0ATPCEtaC*weight[i]);
 967             hCos2V0CTPCEtaA->Fill(fCentrality,i,cos2V0CTPCEtaA*weight[i]);
 968             hCos2V0CTPCEtaC->Fill(fCentrality,i,cos2V0CTPCEtaC*weight[i]);
 969
 970             hCos2SumWeights->Fill(fCentrality,i,weight[i]);
 971         }
 972
 973         // Fill Resolution before EP Flattening
 974         Double_t cos2V0ATPCBF=TMath::Cos(Double_t(fHarmonic)*(fRPTPCBF-fRPV0ABF));
 975         Double_t cos2V0CTPCBF=TMath::Cos(Double_t(fHarmonic)*(fRPTPCBF-fRPV0CBF));
 976         Double_t cos2V0AV0CBF=TMath::Cos(Double_t(fHarmonic)*(fRPV0CBF-fRPV0ABF));
 977         Double_t cos2TPCEtaBF=TMath::Cos(Double_t(fHarmonic)*(fRPTPCEtaABF-fRPTPCEtaCBF));
 978         Double_t cos2TPCBF=TMath::Cos(Double_t(fHarmonic)*(PsiRP1BF-PsiRP2BF));
 979         Double_t cos2V0ATPCEtaABF=TMath::Cos(Double_t(fHarmonic)*(fRPV0ABF-fRPTPCEtaABF));
 980         Double_t cos2V0CTPCEtaABF=TMath::Cos(Double_t(fHarmonic)*(fRPV0CBF-fRPTPCEtaABF));
 981         Double_t cos2V0ATPCEtaCBF=TMath::Cos(Double_t(fHarmonic)*(fRPV0ABF-fRPTPCEtaCBF));
 982         Double_t cos2V0CTPCEtaCBF=TMath::Cos(Double_t(fHarmonic)*(fRPV0CBF-fRPTPCEtaCBF));
 983
 984         hCos2V0ATPC->Fill(fCentrality,nfill,cos2V0ATPCBF);
 985         hCos2V0CTPC->Fill(fCentrality,nfill,cos2V0CTPCBF);
 986         hCos2V0AC->Fill(fCentrality,nfill,cos2V0AV0CBF);
 987         hCos2TPCEta->Fill(fCentrality,nfill,cos2TPCEtaBF);
 988         hCos2TPC->Fill(fCentrality,nfill,cos2TPCBF);
 989         hCos2V0ATPCEtaA->Fill(fCentrality,nfill,cos2V0ATPCEtaABF);
 990         hCos2V0ATPCEtaC->Fill(fCentrality,nfill,cos2V0ATPCEtaCBF);
 991         hCos2V0CTPCEtaA->Fill(fCentrality,nfill,cos2V0CTPCEtaABF);
 992         hCos2V0CTPCEtaC->Fill(fCentrality,nfill,cos2V0CTPCEtaCBF);
 993
 994         hCos2SumWeights->Fill(fCentrality,nfill);
 995
 996     }
 997 }
 998
 999 //________________________________________________________________________
1000 TVector2 AliAnalysisTaskPi0v2::GetEPContribution(AliAODConversionPhoton *gamma){
1001     TVector2 q;
1002     for(Int_t ii=0;ii<2;ii++){
1003         AliVTrack *fCurrentTrack=AliConversionCuts::GetTrack(fInputEvent,gamma->GetTrackLabel(ii));
1004         TVector2 qtrack=GetContributionEP(fCurrentTrack);
1005         q+=qtrack;
1006     }
1007     return q;
1008 }
1009
1010 //________________________________________________________________________
1011 Double_t AliAnalysisTaskPi0v2::GetCorrectedTPCEPAngle(AliAODConversionPhoton *gamma0,AliAODConversionPhoton *gamma1,Bool_t bDoFlattening){
1012     // Correct Event Plane for Dilepton Tracks
1013     TVector2 q0(*fEP->GetQVector());
1014     if(gamma0)q0-=GetEPContribution(gamma0);
1015     if(gamma1)q0-=GetEPContribution(gamma1);
1016     Double_t EPangle=GetPsiInRange(q0.Phi()/Double_t(fHarmonic));
1017     if(bDoFlattening)EPangle=ApplyFlattening(EPangle,kEPTPC);
1018
1019     return EPangle;
1020 }
1021
1022 //________________________________________________________________________
1023 Double_t AliAnalysisTaskPi0v2::GetTPCSubEPEta(EEventPlane ep){
1024
1025     Double_t etamin,etamax;
1026     switch(ep){
1027     case kEPTPCEtaA:
1028         etamin=fEtaGap/2;
1029         etamax=1;
1030         break;
1031     case kEPTPCEtaC:
1032         etamin=-1;
1033         etamax=-fEtaGap/2;
1034         break;
1035     default:
1036         return 0;
1037     }
1038
1039     TVector2 q;
1040     for(Int_t ii=0;ii<fInputEvent->GetNumberOfTracks();ii++){
1041         AliVTrack *fCurrentTrack=dynamic_cast<AliVTrack*>(fInputEvent->GetTrack(ii));
1042         if(!fCurrentTrack)continue;
1043         if(fCurrentTrack->Eta()>=etamin&&fCurrentTrack->Eta()<=etamax){
1044             TVector2 qtrack=GetContributionEP(fCurrentTrack);
1045             q+=qtrack;
1046         }
1047     }
1048
1049     Double_t phi=GetPsiInRange(q.Phi()/Double_t(fHarmonic));
1050
1051     return phi;
1052 }
1053
1054 //________________________________________________________________________
1055 TVector2 AliAnalysisTaskPi0v2::GetContributionEP(AliVTrack *track){
1056
1057     TVector2 q(0,0);
1058
1059     TArrayF *fQContributionX=fEP->GetQContributionXArray();
1060     TArrayF *fQContributionY=fEP->GetQContributionYArray();
1061
1062     Int_t trackID=track->GetID();
1063
1064     if(fIsAOD){
1065         if((trackID>-1)&&fUseTPCOnlyTracks)return q;
1066         if((trackID<0)&&!fUseTPCOnlyTracks)return q;
1067         if (fUseTPCOnlyTracks) trackID = trackID*(-1) - 1;
1068     }
1069
1070     q.Set(fQContributionX->GetAt(trackID),fQContributionY->GetAt(trackID));
1071
1072     return q;
1073 }
1074
1075 //________________________________________________________________________
1076 Double_t AliAnalysisTaskPi0v2::GetPsiInRange(Double_t phi){
1077
1078     Double_t newphi=phi;
1079     while(newphi<0)newphi+=fPsiMax;
1080     while(newphi>fPsiMax)newphi-=fPsiMax;
1081     return newphi;
1082 }
1083
1084 //________________________________________________________________________
1085 Double_t AliAnalysisTaskPi0v2::GetEventPlaneAngle(EEventPlaneMethod EPmethod,Double_t eta,AliAODConversionPhoton *gamma0,AliAODConversionPhoton *gamma1,Bool_t bDoFlattening)
1086 {
1087     // If arguments are not null, the contribution of these photons is subtracted from the TPC EP
1088
1089     if(fConversionCuts->IsHeavyIon()){
1090
1091         // For MC select random EP angle in order to avoid correlations due to azimuth dependent efficiencies (ITS holes)
1092         if(fMCEvent){
1093             return fRandomizer->Uniform(0,fPsiMax);
1094         }
1095         switch(EPmethod){
1096         case kTPC:
1097             return GetCorrectedTPCEPAngle(gamma0,gamma1,bDoFlattening);
1098         case kTPCEtaGap:
1099             // Use opposite EP
1100             if(bDoFlattening){
1101                 if(eta<0)return fRPTPCEtaA;
1102                 else return fRPTPCEtaC;
1103             }
1104             else{
1105                 if(eta<0)return fRPTPCEtaABF;
1106                 else return fRPTPCEtaCBF;
1107             }
1108         case kV0A:
1109             if(bDoFlattening)return fRPV0A;
1110             else return fRPV0ABF;
1111         case kV0C:
1112             if(bDoFlattening)return fRPV0C;
1113             else return fRPV0CBF;
1114         default:
1115             return 0;
1116         }
1117     }
1118
1119     // NO EP in pp mode
1120     return 0;
1121 }
1122
1123 ///________________________________________________________________________
1124 Bool_t AliAnalysisTaskPi0v2::SetCentrality(){
1125
1126     // Set centrality bin for current event
1127
1128     if(!fConversionCuts->IsHeavyIon()){
1129         fCentrality=0;
1130         fCentralityBin=0;
1131         return kTRUE;
1132     }
1133
1134     fCentrality=fConversionCuts->GetCentrality(fInputEvent);
1135
1136     if(fNCentralityBins>1){
1137         for(fCentralityBin=0;fCentralityBin<fNCentralityBins;fCentralityBin++){
1138             if(fCentrality>=fCentralityBins[fCentralityBin]&&fCentrality<fCentralityBins[fCentralityBin+1])return kTRUE;
1139         }
1140         return kFALSE;
1141     }
1142     fCentralityBin=0;
1143     return kTRUE;
1144 }
1145
1146 //________________________________________________________________________
1147 void AliAnalysisTaskPi0v2::SetCentralityBins(Double_t *bins,Int_t nbins)
1148 {
1149     // Set Centrality bins for analysis
1150
1151     fNCentralityBins=nbins;
1152
1153     if(fCentralityBins){
1154         delete[] fCentralityBins;
1155         fCentralityBins=NULL;
1156     }
1157
1158     fCentralityBins=new Double_t[fNCentralityBins+1];
1159     for(Int_t ii=0;ii<=fNCentralityBins;ii++){
1160         fCentralityBins[ii]=bins[ii];
1161     }
1162 }
1163
1164 //________________________________________________________________________
1165 void AliAnalysisTaskPi0v2::SetCuts(AliConversionSelection **conversionselection,Int_t ncuts){
1166
1167     if(fConversionSelection){
1168         for(Int_t ii=0;ii<fNCuts;ii++)delete fConversionSelection;
1169         delete[] fConversionSelection;
1170         fConversionSelection=NULL;
1171     }
1172     fNCuts=ncuts;
1173     fConversionSelection=new AliConversionSelection*[fNCuts];
1174     for(Int_t ii=0;ii<fNCuts;ii++)fConversionSelection[ii]=new AliConversionSelection(*conversionselection[ii]);
1175 }
1176
1177 //___________________________________________________________________________
1178 Int_t AliAnalysisTaskPi0v2::GetRunIndex(Int_t run){
1179
1180     switch(run){
1181
1182         //LHC11h (131 runs)
1183     case 167902 : return 140;
1184     case 167903 : return 141;
1185     case 167909 : return 142;
1186     case 167915 : return 143;
1187     case 167920 : return 144;
1188     case 167985 : return 145;
1189     case 167986 : return 146;
1190     case 167987 : return 147;
1191     case 167988 : return 148;
1192     case 168066 : return 149;
1193     case 168068 : return 150;
1194     case 168069 : return 151;
1195     case 168076 : return 152;
1196     case 168103 : return 153;
1197     case 168104 : return 154;
1198     case 168105 : return 155;
1199     case 168107 : return 156;
1200     case 168108 : return 157;
1201     case 168115 : return 158;
1202     case 168212 : return 159;
1203     case 168310 : return 160;
1204     case 168311 : return 161;
1205     case 168322 : return 162;
1206     case 168325 : return 163;
1207     case 168341 : return 164;
1208     case 168342 : return 165;
1209     case 168361 : return 166;
1210     case 168362 : return 167;
1211     case 168458 : return 168;
1212     case 168460 : return 169;
1213     case 168461 : return 170;
1214     case 168464 : return 171;
1215     case 168467 : return 172;
1216     case 168511 : return 173;
1217     case 168512 : return 174;
1218     case 168514 : return 175;
1219     case 168777 : return 176;
1220     case 168826 : return 177;
1221     case 168984 : return 178;
1222     case 168988 : return 179;
1223     case 168992 : return 180;
1224     case 169035 : return 181;
1225     case 169040 : return 182;
1226     case 169044 : return 183;
1227     case 169045 : return 184;
1228     case 169091 : return 185;
1229     case 169094 : return 186;
1230     case 169099 : return 187;
1231     case 169138 : return 188;
1232     case 169143 : return 189;
1233     case 169144 : return 190;
1234     case 169145 : return 191;
1235     case 169148 : return 192;
1236     case 169156 : return 193;
1237     case 169160 : return 194;
1238     case 169167 : return 195;
1239     case 169238 : return 196;
1240     case 169411 : return 197;
1241     case 169415 : return 198;
1242     case 169417 : return 199;
1243     case 169418 : return 200;
1244     case 169419 : return 201;
1245     case 169420 : return 202;
1246     case 169475 : return 203;
1247     case 169498 : return 204;
1248     case 169504 : return 205;
1249     case 169506 : return 206;
1250     case 169512 : return 207;
1251     case 169515 : return 208;
1252     case 169550 : return 209;
1253     case 169553 : return 210;
1254     case 169554 : return 211;
1255     case 169555 : return 212;
1256     case 169557 : return 213;
1257     case 169584 : return 214;
1258     case 169586 : return 215;
1259     case 169587 : return 216;
1260     case 169588 : return 217;
1261     case 169590 : return 218;
1262     case 169591 : return 219;
1263     case 169835 : return 220;
1264     case 169837 : return 221;
1265     case 169838 : return 222;
1266     case 169846 : return 223;
1267     case 169855 : return 224;
1268     case 169858 : return 225;
1269     case 169859 : return 226;
1270     case 169922 : return 227;
1271     case 169923 : return 228;
1272     case 169956 : return 229;
1273     case 169965 : return 230;
1274     case 169975 : return 231;
1275     case 169981 : return 232;
1276     case 170027 : return 233;
1277     case 170036 : return 234;
1278     case 170038 : return 235;
1279     case 170040 : return 236;
1280     case 170081 : return 237;
1281     case 170083 : return 238;
1282     case 170084 : return 239;
1283     case 170085 : return 240;
1284     case 170088 : return 241;
1285     case 170089 : return 242;
1286     case 170091 : return 243;
1287     case 170152 : return 244;
1288     case 170155 : return 245;
1289     case 170159 : return 246;
1290     case 170163 : return 247;
1291     case 170193 : return 248;
1292     case 170195 : return 249;
1293     case 170203 : return 250;
1294     case 170204 : return 251;
1295     case 170207 : return 252;
1296     case 170208 : return 253;
1297     case 170228 : return 254;
1298     case 170230 : return 255;
1299     case 170268 : return 256;
1300     case 170269 : return 257;
1301     case 170270 : return 258;
1302     case 170306 : return 259;
1303     case 170308 : return 260;
1304     case 170309 : return 261;
1305     case 170311 : return 262;
1306     case 170312 : return 263;
1307     case 170313 : return 264;
1308     case 170315 : return 265;
1309     case 170387 : return 266;
1310     case 170388 : return 267;
1311     case 170556 : return 268;
1312     case 170572 : return 269;
1313     case 170593 : return 270;
1314
1315     //LHC10h (137 runs)
1316     case  139517 : return 137;
1317     case  139514 : return 136;
1318     case  139513 : return 135;
1319     case  139511 : return 134;
1320     case  139510 : return 133;
1321     case  139507 : return 132;
1322     case  139505 : return 131;
1323     case  139504 : return 130;
1324     case  139503 : return 129;
1325     case  139470 : return 128;
1326     case  139467 : return 127;
1327     case  139466 : return 126;
1328     case  139465 : return 125;
1329     case  139440 : return 124;
1330     case  139439 : return 123;
1331     case  139438 : return 122;
1332     case  139437 : return 121;
1333     case  139360 : return 120;
1334     case  139329 : return 119;
1335     case  139328 : return 118;
1336     case  139314 : return 117;
1337     case  139311 : return 116;
1338     case  139310 : return 115;
1339     case  139309 : return 114;
1340     case  139308 : return 113;
1341     case  139173 : return 112;
1342     case  139172 : return 111;
1343     case  139110 : return 110;
1344     case  139107 : return 109;
1345     case  139105 : return 108;
1346     case  139104 : return 107;
1347     case  139042 : return 106;
1348     case  139038 : return 105;
1349     case  139037 : return 104;
1350     case  139036 : return 103;
1351     case  139029 : return 102;
1352     case  139028 : return 101;
1353     case  138983 : return 100;
1354     case  138982 : return 99;
1355     case  138980 : return 98;
1356     case  138979 : return 97;
1357     case  138978 : return 96;
1358     case  138977 : return 95;
1359     case  138976 : return 94;
1360     case  138973 : return 93;
1361     case  138972 : return 92;
1362     case  138965 : return 91;
1363     case  138924 : return 90;
1364     case  138872 : return 89;
1365     case  138871 : return 88;
1366     case  138870 : return 87;
1367     case  138837 : return 86;
1368     case  138830 : return 85;
1369     case  138828 : return 84;
1370     case  138826 : return 83;
1371     case  138796 : return 82;
1372     case  138795 : return 81;
1373     case  138742 : return 80;
1374     case  138732 : return 79;
1375     case  138730 : return 78;
1376     case  138666 : return 77;
1377     case  138662 : return 76;
1378     case  138653 : return 75;
1379     case  138652 : return 74;
1380     case  138638 : return 73;
1381     case  138624 : return 72;
1382     case  138621 : return 71;
1383     case  138583 : return 70;
1384     case  138582 : return 69;
1385     // case  138579 : return 68;
1386     // case  138578 : return 67;
1387     case  138534 : return 66;
1388     case  138469 : return 65;
1389     case  138442 : return 64;
1390     case  138439 : return 63;
1391     case  138438 : return 62;
1392     case  138396 : return 61;
1393     case  138364 : return 60;
1394     case  138359 : return 59;
1395     case  138275 : return 58;
1396     case  138225 : return 57;
1397     case  138201 : return 56;
1398     case  138200 : return 55;
1399     case  138197 : return 54;
1400     case  138192 : return 53;
1401     case  138190 : return 52;
1402     case  138154 : return 51;
1403     case  138153 : return 50;
1404     case  138151 : return 49;
1405     case  138150 : return 48;
1406     case  138126 : return 47;
1407     case  138125 : return 46;
1408     case  137848 : return 45;
1409     case  137847 : return 44;
1410     case  137844 : return 43;
1411     case  137843 : return 42;
1412     case  137752 : return 41;
1413     case  137751 : return 40;
1414     case  137748 : return 39;
1415     case  137724 : return 38;
1416     case  137722 : return 37;
1417     case  137718 : return 36;
1418     case  137704 : return 35;
1419     case  137693 : return 34;
1420     case  137692 : return 33;
1421     case  137691 : return 32;
1422     case  137689 : return 31;
1423     case  137686 : return 30;
1424     case  137685 : return 29;
1425     case  137639 : return 28;
1426     case  137638 : return 27;
1427     case  137608 : return 26;
1428     case  137595 : return 25;
1429     case  137549 : return 24;
1430     case  137546 : return 23;
1431     case  137544 : return 22;
1432     case  137541 : return 21;
1433     case  137539 : return 20;
1434     case  137531 : return 19;
1435     case  137530 : return 18;
1436     case  137443 : return 17;
1437     case  137441 : return 16;
1438     case  137440 : return 15;
1439     case  137439 : return 14;
1440     case  137434 : return 13;
1441     case  137432 : return 12;
1442     case  137431 : return 11;
1443     case  137430 : return 10;
1444     case  137366 : return 9;
1445     case  137243 : return 8;
1446     case  137236 : return 7;
1447     case  137235 : return 6;
1448     case  137232 : return 5;
1449     case  137231 : return 4;
1450     case  137165 : return 3;
1451     case  137162 : return 2;
1452     case  137161 : return 1;
1453
1454     // Default
1455     //default : return -1;
1456     default : return 275;
1457     }
1458 }
1459
1460 //____________________________________________________________________________
1461 void AliAnalysisTaskPi0v2::GetV0EP(AliVEvent * event,Double_t &rpv0a,Double_t &rpv0c){
1462
1463     if (fPeriod.CompareTo("LHC10h")==0){
1464         // Corrected VZERO EP (from AliAnalysisTaskPi0Flow)
1465         //VZERO data
1466         AliESDVZERO* esdV0 = (AliESDVZERO*)event->GetVZEROData();
1467
1468         //reset Q vector info
1469         Double_t Qxa2 = 0, Qya2 = 0;
1470         Double_t Qxc2 = 0, Qyc2 = 0;
1471
1472         for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1473             Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1474             Float_t multv0 = esdV0->GetMultiplicity(iv0);
1475             Double_t lqx=TMath::Cos(fHarmonic*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1476             Double_t lqy=TMath::Sin(fHarmonic*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1477             if (iv0 < 32){ // V0C
1478                 Qxc2 += lqx;
1479                 Qyc2 += lqy;
1480             } else {       // V0A
1481                 Qxa2 += lqx;
1482                 Qya2 += lqy;
1483             }
1484         }
1485
1486         Int_t iC = -1;
1487         // centrality bins
1488         if(fCentrality < 5) iC = 0;
1489         else if(fCentrality < 10) iC = 1;
1490         else if(fCentrality < 20) iC = 2;
1491         else if(fCentrality < 30) iC = 3;
1492         else if(fCentrality < 40) iC = 4;
1493         else if(fCentrality < 50) iC = 5;
1494         else if(fCentrality < 60) iC = 6;
1495         else if(fCentrality < 70) iC = 7;
1496         else iC = 8;
1497
1498         //grab for each centrality the proper histo with the Qx and Qy to do the recentering
1499         Double_t Qxamean2 = fMeanQ[iC][1][0];
1500         Double_t Qxarms2  = fWidthQ[iC][1][0];
1501         Double_t Qyamean2 = fMeanQ[iC][1][1];
1502         Double_t Qyarms2  = fWidthQ[iC][1][1];
1503
1504         Double_t Qxcmean2 = fMeanQ[iC][0][0];
1505         Double_t Qxcrms2  = fWidthQ[iC][0][0];
1506         Double_t Qycmean2 = fMeanQ[iC][0][1];
1507         Double_t Qycrms2  = fWidthQ[iC][0][1];
1508
1509         Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1510         Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1511         Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1512         Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1513         rpv0a = TMath::ATan2(QyaCor2, QxaCor2)/Double_t(fHarmonic);
1514         rpv0c = TMath::ATan2(QycCor2, QxcCor2)/Double_t(fHarmonic);
1515
1516         //rpv0a = TMath::ATan2(Qya2, Qxa2)/Double_t(fHarmonic);
1517         //rpv0c = TMath::ATan2(Qyc2, Qxc2)/Double_t(fHarmonic);
1518
1519         // cout<<"Compare v"<<fHarmonic<<" "<<rpv0a<<" "<<fInputEvent->GetEventplane()->GetEventplane("V0A",fInputEvent,fHarmonic)<<endl;
1520
1521     }
1522     if (fPeriod.CompareTo("LHC11h")==0){
1523
1524         AliEventplane *eventEP=fInputEvent->GetEventplane();
1525
1526         Double_t qx,qy;
1527         rpv0a=eventEP->CalculateVZEROEventPlane(fInputEvent,8,fHarmonic,qx,qy);
1528         rpv0c=eventEP->CalculateVZEROEventPlane(fInputEvent,9,fHarmonic,qx,qy);
1529     }
1530
1531     rpv0a=GetPsiInRange(rpv0a);
1532     rpv0c=GetPsiInRange(rpv0c);
1533 }
1534
1535 //_____________________________________________________________________________
1536 void AliAnalysisTaskPi0v2::LoadVZEROCalibration(Int_t run){
1537
1538     cout<<fPeriod.Data()<<endl;
1539
1540     // VZERO Phi Weights and Recentering
1541     if (fPeriod.CompareTo("LHC10h")==0){
1542         TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1543         TFile *foadb = TFile::Open(oadbfilename.Data());
1544
1545         if(!foadb){
1546             printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1547             return;
1548         }
1549
1550         AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1551         if(!cont){
1552             printf("OADB object hMultV0BefCorr is not available in the file\n");
1553             return;
1554         }
1555
1556         if(!(cont->GetObject(run))){
1557             printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1558             run = 137366;
1559         }
1560         printf("Setting V0 calibration \n") ;
1561         fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1562
1563         TF1 *fpol0 = new TF1("fpol0","pol0");
1564         fMultV0->Fit(fpol0,"0","",0,31);
1565         fV0Cpol = fpol0->GetParameter(0);
1566         fMultV0->Fit(fpol0,"0","",32,64);
1567         fV0Apol = fpol0->GetParameter(0);
1568
1569         for(Int_t iside=0;iside<2;iside++){
1570             for(Int_t icoord=0;icoord<2;icoord++){
1571                 for(Int_t i=0;i  < nCentrBinV0;i++){
1572                     char namecont[100];
1573
1574                     if(iside==0 && icoord==0)
1575                         snprintf(namecont,100,"hQxc%i_%i",fHarmonic,i);
1576                     else if(iside==1 && icoord==0)
1577                         snprintf(namecont,100,"hQxa%i_%i",fHarmonic,i);
1578                     else if(iside==0 && icoord==1)
1579                         snprintf(namecont,100,"hQyc%i_%i",fHarmonic,i);
1580                     else if(iside==1 && icoord==1)
1581                         snprintf(namecont,100,"hQya%i_%i",fHarmonic,i);
1582
1583                     cont = (AliOADBContainer*) foadb->Get(namecont);
1584                     if(!cont){
1585                         printf("OADB object %s is not available in the file\n",namecont);
1586                         return;
1587                     }
1588
1589                     if(!(cont->GetObject(run))){
1590                         printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1591                         run = 137366;
1592                     }
1593                     fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1594                     fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1595                 }
1596             }
1597         }
1598     }
1599 }
1600
1601 //_____________________________________________________________________________
1602 void AliAnalysisTaskPi0v2::LoadTPCCalibration(Int_t run){
1603
1604     // TPC Event Plane Weights
1605     AliOADBContainer *fEPContainer=NULL;
1606     TString oadbfilename="";
1607
1608     if (fPeriod.CompareTo("LHC10h")==0){
1609         // LHC10h
1610
1611         if(fIsAOD){
1612             oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
1613         }
1614         else{
1615             oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
1616         }
1617
1618         TFile foadb(oadbfilename);
1619         if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
1620
1621         AliInfo("Using Standard OADB");
1622         fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
1623         if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
1624         foadb.Close();
1625
1626         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
1627
1628         Bool_t emptybins;
1629
1630         int iter = 0;
1631         while (iter<3){
1632             emptybins = kFALSE;
1633
1634             for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
1635                 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
1636                     emptybins = kTRUE;
1637                 }
1638             }
1639             if (emptybins) {
1640                 cout << "empty bins - rebinning!" << endl;
1641                 fPhiDist[0]->Rebin();
1642                 iter++;
1643             }
1644             else iter = 3;
1645         }
1646
1647         if (emptybins) {
1648             AliError("After Maximum of rebinning still empty Phi-bins!!!");
1649         }
1650     }
1651
1652     if (fPeriod.CompareTo("LHC11h")==0){
1653         // LHC11h
1654
1655         oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
1656         TFile *foadb = TFile::Open(oadbfilename);
1657         if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
1658
1659         AliInfo("Using Standard OADB");
1660         fSparseDist = (THnSparse*) foadb->Get("Default");
1661         if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
1662         foadb->Close();
1663         if(!fHruns){
1664           fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
1665            fHruns->SetName("runsHisto");
1666         }
1667
1668         Int_t runbin=fHruns->FindBin(fRunNumber);
1669         if (fHruns->GetBinContent(runbin) > 1){
1670             fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
1671         }
1672         else if(fHruns->GetBinContent(runbin) < 2){
1673             fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
1674             AliInfo("Using integrated Phi-weights for this run");
1675         }
1676         for (Int_t i = 0; i<4 ;i++)
1677         {
1678             if(fPhiDist[i]){
1679                 delete fPhiDist[i];
1680                 fPhiDist[i] = 0x0;
1681             }
1682             if(i == 0){
1683                 fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
1684                 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
1685             if(i == 1){
1686                 fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
1687                 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
1688             if(i == 2){
1689                 fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
1690                 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
1691             if(i == 3){
1692                 fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
1693                 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
1694             fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
1695             fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
1696             fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
1697             fSparseDist->GetAxis(2)->SetRange(1,2);
1698         }
1699         fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
1700     }
1701     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", run));
1702
1703 }
1704
1705 //_________________________________________________________________________
1706 Int_t AliAnalysisTaskPi0v2::GetAODEPTrackFilterBit(){
1707
1708     if(fUseTPCOnlyTracks){
1709         return 128;// TPC only with vertex constraint
1710     }
1711     return 1;// Use Global Tracks
1712 }
1713
1714 //_________________________________________________________________________
1715 TObjArray* AliAnalysisTaskPi0v2::GetEventPlaneTracks(Int_t &maxID)
1716 {
1717     TObjArray *tracklist=NULL;
1718
1719     AliESDtrackCuts *fEPESDtrackCuts=AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
1720     fEPESDtrackCuts->SetPtRange(0.15,20.);
1721     fEPESDtrackCuts->SetEtaRange(-0.8,0.8);
1722
1723     Int_t fAODfilterbit=GetAODEPTrackFilterBit();
1724
1725     if(fInputEvent->IsA()==AliESDEvent::Class()){
1726         maxID=fInputEvent->GetNumberOfTracks();
1727         tracklist=fEPESDtrackCuts->GetAcceptedTracks((AliESDEvent*)fInputEvent,fUseTPCOnlyTracks);
1728     }
1729     if(fInputEvent->IsA()==AliAODEvent::Class()){
1730
1731         // From AliEPSelectionTask
1732         tracklist = new TObjArray();
1733
1734         AliAODTrack *tr = 0;
1735         Int_t maxid1 = 0;
1736         Int_t maxidtemp = -1;
1737         Float_t ptlow = 0;
1738         Float_t ptup = 0;
1739         Float_t etalow = 0;
1740         Float_t etaup = 0;
1741         fEPESDtrackCuts->GetPtRange(ptlow,ptup);
1742         fEPESDtrackCuts->GetEtaRange(etalow,etaup);
1743         Int_t ntpc = fEPESDtrackCuts->GetMinNClusterTPC();
1744
1745         for (Int_t i = 0; i < fInputEvent->GetNumberOfTracks() ; i++){
1746             tr = (AliAODTrack*)fInputEvent->GetTrack(i);
1747             maxidtemp = tr->GetID();
1748             if(maxidtemp < 0 && fAODfilterbit != 128) continue;// id<0 means filter bit 128
1749             if(maxidtemp > -1 && fAODfilterbit == 128) continue;// id>01 means filter bit 1
1750             if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
1751             if (maxidtemp > maxid1) maxid1 = maxidtemp;
1752             if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
1753                 tracklist->Add(tr);
1754             }
1755         }
1756         maxID = maxid1;
1757     }
1758     delete fEPESDtrackCuts;
1759     if(!tracklist)AliError("No tracklist");
1760     return tracklist;
1761 }
1762
1763 //____________________________________________________________________________
1764 Bool_t AliAnalysisTaskPi0v2::GetTPCEventPlane(){
1765
1766     if(fEP){
1767         delete fEP;
1768         fEP=NULL;
1769     }
1770     fEP=new AliEventplane();
1771
1772     float mQx=0, mQy=0;
1773     float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
1774     AliVTrack* track;
1775     Double_t weight;
1776     Int_t idtemp = -1;
1777     int trackcounter1=0, trackcounter2=0;
1778
1779     Int_t maxID=0;
1780
1781     TObjArray *tracklist=GetEventPlaneTracks(maxID);
1782
1783     fEP->GetQContributionXArray()->Set(maxID);
1784     fEP->GetQContributionYArray()->Set(maxID);
1785     fEP->GetQContributionXArraysub1()->Set(maxID);
1786     fEP->GetQContributionYArraysub1()->Set(maxID);
1787     fEP->GetQContributionXArraysub2()->Set(maxID);
1788     fEP->GetQContributionYArraysub2()->Set(maxID);
1789
1790     int nt = tracklist->GetEntries();
1791     for (int i=0; i<nt; i++){
1792         weight = 1;
1793         track = dynamic_cast<AliVTrack*> (tracklist->At(i));
1794         if (track) {
1795             // Fill Eta Distribution
1796             hEtaTPCEP->Fill(fCentrality,track->Eta());
1797
1798             weight=GetWeight(track);
1799             idtemp = track->GetID();
1800             // TPC only tracks have negative id ((-1)*IDESD - 1) in AOD
1801             if (fIsAOD && (fUseTPCOnlyTracks)) idtemp = idtemp*(-1) - 1;
1802
1803             Double_t qx=weight*cos(Double_t(fHarmonic)*track->Phi());
1804             Double_t qy=weight*sin(Double_t(fHarmonic)*track->Phi());
1805             fEP->GetQContributionXArray()->AddAt(qx,idtemp);
1806             fEP->GetQContributionYArray()->AddAt(qy,idtemp);
1807
1808             mQx += (qx);
1809             mQy += (qy);
1810
1811             // This loop splits the track set into 2 random subsets
1812             if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
1813                 float random = fRandomizer->Rndm();
1814                 if(random < .5){
1815                     mQx1 += (qx);
1816                     mQy1 += (qy);
1817                     fEP->GetQContributionXArraysub1()->AddAt(qx,idtemp);
1818                     fEP->GetQContributionYArraysub1()->AddAt(qy,idtemp);
1819                     trackcounter1++;
1820                     }
1821                 else {
1822                     mQx2 += (qx);
1823                     mQy2 += (qy);
1824                     fEP->GetQContributionXArraysub2()->AddAt(qx,idtemp);
1825                     fEP->GetQContributionYArraysub2()->AddAt(qy,idtemp);
1826                     trackcounter2++;
1827                 }
1828             }
1829             else{
1830                 if( trackcounter1 >= int(nt/2.)){
1831                     mQx2 += (qx);
1832                     mQy2 += (qy);
1833                     fEP->GetQContributionXArraysub2()->AddAt(qx,idtemp);
1834                     fEP->GetQContributionYArraysub2()->AddAt(qy,idtemp);
1835                     trackcounter2++;
1836                 }
1837                 else {
1838                     mQx1 += (qx);
1839                     mQy1 += (qy);
1840                     fEP->GetQContributionXArraysub1()->AddAt(qx,idtemp);
1841                     fEP->GetQContributionYArraysub1()->AddAt(qy,idtemp);
1842                     trackcounter1++;
1843                 }
1844             }
1845         }
1846     }
1847
1848     tracklist->Clear();
1849     delete tracklist;
1850     tracklist = NULL;
1851
1852     TVector2 *mQ=new TVector2();
1853     mQ->Set(mQx,mQy);
1854     Double_t EPAngle=mQ->Phi()/Double_t(fHarmonic);
1855
1856     TVector2 *fQsub1=new TVector2();
1857     TVector2 *fQsub2=new TVector2();
1858     fQsub1->Set(mQx1,mQy1);
1859     fQsub2->Set(mQx2,mQy2);
1860
1861     fEP->SetQVector(mQ);
1862     fEP->SetEventplaneQ(EPAngle);
1863     fEP->SetQsub(fQsub1,fQsub2);
1864     fEP->SetQsubRes(fQsub1->Phi()/Double_t(fHarmonic) - fQsub2->Phi()/Double_t(fHarmonic));
1865
1866     Int_t ntracks=trackcounter1+trackcounter2;
1867
1868     if(ntracks<3)return kFALSE;// <3 -> no subevents
1869     return kTRUE;
1870 }
1871
1872
1873 //____________________________________________________________________________
1874 Double_t AliAnalysisTaskPi0v2::ApplyFlattening(Double_t phi,EEventPlane ep){
1875
1876     Double_t c2=0;
1877     Double_t s2=0;
1878     Double_t c4=0;
1879     Double_t s4=0;
1880
1881     // GetCorrection Parameter
1882     if(fHarmonic==2){
1883
1884         if(ep==kEPTPC){
1885
1886             Double_t cc2[5]={0.00904396,0.00472483,0.00306154,0.00218462,0.00167447};
1887             Double_t cs2[5]={0.00885519,0.00516223,0.00411065,0.00380145,0.00324424};
1888             Double_t cc4[5]={-0.00110933,-0.00110521,-0.00124342,0.00104131,0.000651779};
1889             Double_t cs4[5]={0.00163869,-0.00053565,0.000878745,-0.000563657,-0.000604021};
1890
1891             c2=cc2[fCentralityBin];
1892             s2=cs2[fCentralityBin];
1893             c4=cc4[fCentralityBin];
1894             s4=cs4[fCentralityBin];
1895         }
1896
1897         if(ep==kEPTPCEtaA){
1898
1899             Double_t cc2[5]={0.00529447,0.00278029,0.00315325,0.00173634,0.000763168};
1900             Double_t cs2[5]={0.00314285,0.00170173,0.00263333,0.0018509,0.00223784};
1901             Double_t cc4[5]={-0.000737254,-0.00037845,-0.000492715,0.000775897,0.000768656};
1902             Double_t cs4[5]={0.000347583,3.79872e-05,0.000387037,-0.000186129,0.000432698};
1903
1904             c2=cc2[fCentralityBin];
1905             s2=cs2[fCentralityBin];
1906             c4=cc4[fCentralityBin];
1907             s4=cs4[fCentralityBin];
1908         }
1909
1910         if(ep==kEPTPCEtaC){
1911
1912             Double_t cc2[5]={-0.00562282,-0.00456735,-0.00306068,-0.0027173,-0.00172432};
1913             Double_t cs2[5]={0.0101804,0.00430782,0.00394715,0.00350156,0.00302749};
1914             Double_t cc4[5]={0.00150831,-0.00159271,-0.000964157,0.000525894,9.93172e-05};
1915             Double_t cs4[5]={0.00119279,-4.74629e-05,0.000118845,0.000278554,3.20868e-05};
1916
1917             c2=cc2[fCentralityBin];
1918             s2=cs2[fCentralityBin];
1919             c4=cc4[fCentralityBin];
1920             s4=cs4[fCentralityBin];
1921         }
1922
1923         if(ep==kEPV0A){
1924
1925             Double_t cc2[5]={0.046427,0.0105401,-0.000152992,-0.00578274,-0.0108038};
1926             Double_t cs2[5]={0.00551503,0.0158159,0.00965148,0.00135414,-0.00548846};
1927             Double_t cc4[5]={0.00362833,0.00170777,0.000152998,0.00223823,0.00215164};
1928             Double_t cs4[5]={0.00349056,0.00142802,0.00123298,0.00207995,0.00145625};
1929
1930             c2=cc2[fCentralityBin];
1931             s2=cs2[fCentralityBin];
1932             c4=cc4[fCentralityBin];
1933             s4=cs4[fCentralityBin];
1934         }
1935
1936         if(ep==kEPV0C){
1937
1938             Double_t cc2[5]={-0.00473277,-0.000371313,0.000857122,-1.54263e-05,-0.000686139};
1939             Double_t cs2[5]={0.00408304,-0.00208615,-0.00149018,-0.000853616,-2.78855e-05};
1940             Double_t cc4[5]={-0.00451741,-0.00399036,-0.00318784,-0.00186472,-0.00106299};
1941             Double_t cs4[5]={0.00188045,-0.00713956,-0.00484254,-0.00448149,-0.00482164};
1942
1943             c2=cc2[fCentralityBin];
1944             s2=cs2[fCentralityBin];
1945             c4=cc4[fCentralityBin];
1946             s4=cs4[fCentralityBin];
1947         }
1948     }
1949
1950     if(fHarmonic==3){
1951
1952         if(ep==kEPTPC){
1953
1954             Double_t cc2[5]={0.0116542,0.0103631,0.00897965,0.00707409,0.00605151};
1955             Double_t cs2[5]={-0.0171191,-0.013024,-0.0114752,-0.0086613,-0.00706863};
1956             Double_t cc4[5]={-0.000602948,0.00144836,-0.000193641,0.000108773,-0.000518333};
1957             Double_t cs4[5]={-0.00164769,0.00134327,-0.00106369,7.96546e-06,-0.000261517};
1958
1959             c2=cc2[fCentralityBin];
1960             s2=cs2[fCentralityBin];
1961             c4=cc4[fCentralityBin];
1962             s4=cs4[fCentralityBin];
1963         }
1964
1965         if(ep==kEPTPCEtaA){
1966
1967             Double_t cc2[5]={0.000386277,0.000119225,0.00111969,0.000534801,0.000642703};
1968             Double_t cs2[5]={-0.00581604,-0.00607255,-0.00443819,-0.00268834,-0.00299961};
1969             Double_t cc4[5]={0.00051635,0.00036326,-0.000221272,4.66775e-05,-3.05784e-06};
1970             Double_t cs4[5]={1.43285e-05,0.000514099,0.000619339,0.00106466,0.000344196};
1971
1972             c2=cc2[fCentralityBin];
1973             s2=cs2[fCentralityBin];
1974             c4=cc4[fCentralityBin];
1975             s4=cs4[fCentralityBin];
1976         }
1977
1978         if(ep==kEPTPCEtaC){
1979
1980             Double_t cc2[5]={0.0116475,0.0102385,0.00801121,0.00552336,0.00423273};
1981             Double_t cs2[5]={-0.0112722,-0.00796059,-0.00683678,-0.00531097,-0.00430716};
1982             Double_t cc4[5]={-0.000609051,1.36573e-08,-0.000464961,-0.000387943,-2.28363e-05};
1983             Double_t cs4[5]={0.00125449,0.00168484,-0.000390491,-0.000219447,8.11997e-07};
1984
1985             c2=cc2[fCentralityBin];
1986             s2=cs2[fCentralityBin];
1987             c4=cc4[fCentralityBin];
1988             s4=cs4[fCentralityBin];
1989         }
1990
1991         if(ep==kEPV0A){
1992
1993             Double_t cc2[5]={-0.0057427,-0.00482728,-0.00565919,-0.000717094,-0.00933233};
1994             Double_t cs2[5]={0.0306554,-0.0144675,-0.0159243,-0.0120465,-0.00814124};
1995             Double_t cc4[5]={-0.002868,0.00159533,0.00754171,0.00683898,0.00689441};
1996             Double_t cs4[5]={0.00083196,0.00198133,4.68307e-05,-0.00018187,-0.0014258};
1997
1998             c2=cc2[fCentralityBin];
1999             s2=cs2[fCentralityBin];
2000             c4=cc4[fCentralityBin];
2001             s4=cs4[fCentralityBin];
2002         }
2003
2004         if(ep==kEPV0C){
2005
2006             Double_t cc2[5]={-0.00259141,-0.00115826,-0.000738658,-4.96667e-05,-0.000346694};
2007             Double_t cs2[5]={-0.0111001,0.00258109,0.00110959,-0.000147296,-0.000199817};
2008             Double_t cc4[5]={0.000968742,0.00157903,0.000206157,0.000444206,-0.00046573};
2009             Double_t cs4[5]={-0.00307319,-0.0047952,-0.00412117,-0.00320344,-0.00386629};
2010
2011             c2=cc2[fCentralityBin];
2012             s2=cs2[fCentralityBin];
2013             c4=cc4[fCentralityBin];
2014             s4=cs4[fCentralityBin];
2015         }
2016     }
2017
2018     // Do correction
2019     Double_t newphi=phi;
2020     newphi+=1/Double_t(fHarmonic)*(2*c2*sin(Double_t(fHarmonic)*phi)-2*s2*cos(Double_t(fHarmonic)*phi)+c4*sin(2*Double_t(fHarmonic)*phi)-s4*cos(2*Double_t(fHarmonic)*phi));
2021     newphi=GetPsiInRange(newphi);
2022
2023     return newphi;
2024 }
2025
2026 //________________________________________________________________________
2027 Double_t AliAnalysisTaskPi0v2::GetWeight(TObject* track1)
2028 {
2029     Double_t ptweight=1;
2030     AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
2031     if (track) {
2032         if (track->Pt()<2) ptweight=track->Pt();
2033         else ptweight=2;
2034     }
2035     return ptweight*GetPhiWeight(track);
2036 }
2037
2038 //________________________________________________________________________
2039 Double_t AliAnalysisTaskPi0v2::GetPhiWeight(TObject* track1)
2040 {
2041   Double_t phiweight=1;
2042   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
2043
2044   TH1F *phiDist = 0x0;
2045   if(track) phiDist = SelectPhiDist(track);
2046
2047   if (phiDist && track) {
2048       Double_t nParticles = phiDist->Integral();
2049       Double_t nPhibins = phiDist->GetNbinsX();
2050
2051       Double_t Phi = track->Phi();
2052
2053       while (Phi<0) Phi += TMath::TwoPi();
2054       while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
2055
2056       Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
2057
2058       if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
2059   }
2060   return phiweight;
2061 }
2062
2063 //_________________________________________________________________________
2064 TH1F* AliAnalysisTaskPi0v2::SelectPhiDist(AliVTrack *track)
2065 {
2066   if (fPeriod.CompareTo("LHC10h")==0) return fPhiDist[0];
2067   else if(fPeriod.CompareTo("LHC11h")==0){
2068       if (track->Charge() < 0){
2069           if(track->Eta() < 0.)       return fPhiDist[0];
2070           else if (track->Eta() > 0.) return fPhiDist[2];
2071       }
2072       else if (track->Charge() > 0){
2073           if(track->Eta() < 0.)       return fPhiDist[1];
2074           else if (track->Eta() > 0.) return fPhiDist[3];
2075       }
2076   }
2077   return 0;
2078 }