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