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