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