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