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