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