PWGGA/GammaConv/AliAnalysisTaskPi0v2.cxx

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