]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0Flow.cxx
In AliAnalysisTaskPi0Flow, changed logic for different length of past event buffer...
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0Flow.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TObjArray.h"
4 #include "TF1.h"
5 #include "TFile.h"
6 #include "TH1F.h"
7 #include "TH2F.h"
8 #include "TH2I.h"
9 #include "TH3F.h"
10 #include "TParticle.h"
11 #include "TCanvas.h"
12 #include "TStyle.h"
13 #include "TRandom.h"
14
15 #include "AliAnalysisManager.h"
16 #include "AliMCEventHandler.h"
17 #include "AliMCEvent.h"
18 #include "AliStack.h"
19 #include "AliAnalysisTaskSE.h"
20 #include "AliAnalysisTaskPi0Flow.h"
21 #include "AliCaloPhoton.h"
22 #include "AliPHOSGeometry.h"
23 #include "TGeoManager.h"
24 #include "AliPHOSEsdCluster.h"
25 #include "AliPHOSCalibData.h"
26 #include "AliESDEvent.h"
27 #include "AliESDCaloCells.h"
28 #include "AliESDVertex.h"
29 #include "AliESDtrackCuts.h"
30 #include "AliAODEvent.h"
31 #include "AliLog.h"
32 #include "AliPID.h"
33 #include "AliCDBManager.h"
34 #include <AliAODCaloCluster.h>
35 #include "AliCentrality.h"
36 #include "AliESDtrackCuts.h"
37 #include "AliEventplane.h"
38 #include "TProfile.h"
39 #include "AliOADBContainer.h"
40
41 // Analysis task to fill histograms with PHOS ESD or AOD clusters and cells
42 // Authors : Dmitri Peressounko
43 // Date    : 28.05.2011
44 // Modified: 03.08.2012 Henrik Qvigstad
45 /* $Id: AliAnalysisTaskPi0Flow.cxx 58214 2012-08-17 06:53:05Z kharlov $ */
46
47 ClassImp(AliAnalysisTaskPi0Flow);
48
49 //________________________________________________________________________
50 Double_t rnlin(Double_t *x, Double_t * /*par*/)
51 {
52   //a = par[0], b = par[1].
53   //1+a*exp(-e/b)
54
55 // return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
56  return 1.015*(0.0241+1.0504*x[0]+0.000249*x[0]*x[0]) ;
57
58 }
59
60 //________________________________________________________________________
61 AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
62 : AliAnalysisTaskSE(name),
63   fCentEdges(10),
64   fCentNMixed(10),
65   fNEMRPBins(9),
66   fPeriod(period),
67   fOutputContainer(0x0),
68   fNonLinCorr(0),
69   fEvent(0x0),
70   fEventESD(0x0),
71   fEventAOD(0x0),
72   fMCStack(0x0),
73   fRunNumber(-999),
74   fInternalRunNumber(0),
75   fPHOSGeo(0),
76   fMultV0(0x0),
77   fV0Cpol(0.),fV0Apol(0.),
78   fESDtrackCuts(0x0),
79   fPHOSCalibData(0x0),
80   fVertexVector(),
81   fVtxBin(0),
82   fCentralityV0M(0.),
83   fCentBin(0),
84   fHaveTPCRP(0),
85   fRP(0),
86   fRPV0A(0),
87   fRPV0C(0),
88   fEMRPBin(0),
89   fCaloPhotonsPHOS(0x0),
90   fCaloPhotonsPHOSLists(0x0)
91 {
92   const int elength = 10;
93   Double_t edges[elength] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
94   fCentEdges = TArrayD(elength, edges);
95   Int_t nMixed[elength] = {4,4,6,10,20,30,50,100,100,100};
96   fCentNMixed = TArrayI(elength, nMixed);
97   
98   for(Int_t i=0;i<kNCenBins;i++){
99     for(Int_t j=0;j<2; j++)
100       for(Int_t k=0; k<2; k++) {
101         fMeanQ[i][j][k]=0.;
102         fWidthQ[i][j][k]=0.;
103       }
104   }
105   
106   fVertex[0]=0; fVertex[1]=0; fVertex[2]=0; 
107
108   // Output slots #0 write into a TH1 container
109   DefineOutput(1,TList::Class());
110
111
112
113   // Set bad channel map
114   char key[55] ;
115   for(Int_t i=0; i<6; i++){
116     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
117     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
118   }
119
120
121   // Initialize non-linrarity correction
122   fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
123
124
125 }
126 //___________________________________________________________________________
127 AliAnalysisTaskPi0Flow::~AliAnalysisTaskPi0Flow()
128 {
129   delete fNonLinCorr;
130   delete fESDtrackCuts;
131   delete fPHOSCalibData;
132   delete fCaloPhotonsPHOSLists;
133 }
134 //________________________________________________________________________
135 void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
136 {
137   // Create histograms
138   // Called once
139   const Int_t nRuns=200 ;
140
141   // histograms
142   if(fOutputContainer != NULL){
143     delete fOutputContainer;
144   }
145   fOutputContainer = new TList();
146   fOutputContainer->SetOwner(kTRUE);
147
148   //========QA histograms=======
149
150   //Event selection
151   fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 12,0.,13.,nRuns,0.,float(nRuns))) ;
152   fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 12,0.,12.)) ;
153
154   //vertex distribution
155   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
156
157   //Centrality
158   fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
159   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
160   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
161   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
162   fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,10.,100,0.,100.)) ;
163
164
165   //Reaction plane
166   fOutputContainer->Add(new TH3F("hPHOSphi","cos" ,10,0.,100.,20,0.,10.,100,-TMath::Pi(),TMath::Pi()));
167
168   fOutputContainer->Add(new TH2F("cos2AC","RP correlation between TPC subs", 100,-1.,1.,20,0.,100.)) ;
169   fOutputContainer->Add(new TH2F("cos2V0AC","RP correlation between VO A and C sides", 100,-1.,1.,20,0.,100.)) ;
170   fOutputContainer->Add(new TH2F("cos2V0ATPC","RP correlation between TPC and V0A", 100,-1.,1.,20,0.,100.)) ;
171   fOutputContainer->Add(new TH2F("cos2V0CTPC","RP correlation between TPC and V0C", 100,-1.,1.,20,0.,100.)) ;
172
173   fOutputContainer->Add(new TH2F("phiRP","RP distribution with TPC", 100,0.,TMath::Pi(),20,0.,100.)) ;
174   fOutputContainer->Add(new TH2F("phiRPflat","RP distribution with TPC flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
175   fOutputContainer->Add(new TH2F("phiRPV0A","RP distribution with V0A", 100,0.,TMath::Pi(),20,0.,100.)) ;
176   fOutputContainer->Add(new TH2F("phiRPV0C","RP distribution with V0C", 100,0.,TMath::Pi(),20,0.,100.)) ;
177   fOutputContainer->Add(new TH3F("phiRPV0AC","RP distribution with V0A and V0C", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
178   fOutputContainer->Add(new TH2F("phiRPV0Aflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
179   fOutputContainer->Add(new TH2F("phiRPV0Cflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
180   fOutputContainer->Add(new TH3F("phiRPV0ATPC","RP distribution with V0A + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
181   fOutputContainer->Add(new TH3F("phiRPV0CTPC","RP distribution with V0C + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
182
183
184   //PHOS QA
185   fOutputContainer->Add(new TH1I("hCellMultEvent"  ,"PHOS cell multiplicity per event"    ,2000,0,2000));
186   fOutputContainer->Add(new TH1I("hCellMultEventM1","PHOS cell multiplicity per event, M1",2000,0,2000));
187   fOutputContainer->Add(new TH1I("hCellMultEventM2","PHOS cell multiplicity per event, M2",2000,0,2000));
188   fOutputContainer->Add(new TH1I("hCellMultEventM3","PHOS cell multiplicity per event, M3",2000,0,2000));
189
190   fOutputContainer->Add(new TH1F("hCellEnergy"  ,"Cell energy"            ,3000,0.,30.));
191   fOutputContainer->Add(new TH1F("hCellEnergyM1","Cell energy in module 1",3000,0.,30.));
192   fOutputContainer->Add(new TH1F("hCellEnergyM2","Cell energy in module 2",3000,0.,30.));
193   fOutputContainer->Add(new TH1F("hCellEnergyM3","Cell energy in module 3",3000,0.,30.));
194
195   fOutputContainer->Add(new TH2F("hCellNXZM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
196   fOutputContainer->Add(new TH2F("hCellNXZM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
197   fOutputContainer->Add(new TH2F("hCellNXZM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
198   fOutputContainer->Add(new TH2F("hCellEXZM1","Cell E(X,Z), M1",64,0.5,64.5, 56,0.5,56.5));
199   fOutputContainer->Add(new TH2F("hCellEXZM2","Cell E(X,Z), M2",64,0.5,64.5, 56,0.5,56.5));
200   fOutputContainer->Add(new TH2F("hCellEXZM3","Cell E(X,Z), M3",64,0.5,64.5, 56,0.5,56.5));
201
202   //Bad Map
203   fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
204   fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
205   fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
206
207   fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
208   fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
209   fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
210
211   fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
212   fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
213   fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
214
215   fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
216   fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
217   fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
218
219
220   //Single photon and pi0 spectrum
221   const Int_t nPtPhot = 300 ;
222   const Double_t ptPhotMax = 30 ;
223   const Int_t nM       = 500;
224   const Double_t mMin  = 0.0;
225   const Double_t mMax  = 1.0;
226
227   //PHOS calibration QA
228   fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
229   fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
230   fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
231   fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
232   fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
233   fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
234
235   // Histograms for different centralities
236   char key[55] ;
237   for(Int_t cent=0; cent < fCentEdges.GetSize(); cent++){
238     snprintf(key,55,"hPhotAll_cen%d",cent) ;
239     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
240     snprintf(key,55,"hPhotAllcore_cen%d",cent) ;
241     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
242     snprintf(key,55,"hPhotAllwou_cen%d",cent) ;
243     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
244     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
245     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
246     snprintf(key,55,"hPhotDisp2_cen%d",cent) ;
247     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
248     snprintf(key,55,"hPhotDispcore_cen%d",cent) ;
249     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
250     snprintf(key,55,"hPhotDispwou_cen%d",cent) ;
251     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
252     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
253     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
254     snprintf(key,55,"hPhotCPVcore_cen%d",cent) ;
255     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
256     snprintf(key,55,"hPhotCPV2_cen%d",cent) ;
257     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
258     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
259     fOutputContainer->Add(new TH1F(key,"Both clusters",nPtPhot,0.,ptPhotMax));
260     snprintf(key,55,"hPhotBothcore_cen%d",cent) ;
261     fOutputContainer->Add(new TH1F(key,"Both clusters",nPtPhot,0.,ptPhotMax));
262
263     snprintf(key,55,"hPi0All_cen%d",cent) ;
264     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
265     snprintf(key,55,"hPi0Allcore_cen%d",cent) ;
266     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
267     snprintf(key,55,"hPi0Allwou_cen%d",cent) ;
268     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
269     snprintf(key,55,"hPi0Disp_cen%d",cent) ;
270     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
271     snprintf(key,55,"hPi0Disp2_cen%d",cent) ;
272     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
273     snprintf(key,55,"hPi0Dispcore_cen%d",cent) ;
274     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
275     snprintf(key,55,"hPi0Dispwou_cen%d",cent) ;
276     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
277     snprintf(key,55,"hPi0CPV_cen%d",cent) ;
278     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
279     snprintf(key,55,"hPi0CPVcore_cen%d",cent) ;
280     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
281     snprintf(key,55,"hPi0CPV2_cen%d",cent) ;
282     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
283     snprintf(key,55,"hPi0Both_cen%d",cent) ;
284     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
285     snprintf(key,55,"hPi0Bothcore_cen%d",cent) ;
286     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
287
288     snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
289     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
290     snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
291     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
292     snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
293     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
294     snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
295     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
296     snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
297     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
298
299     snprintf(key,55,"hSingleAll_cen%d",cent) ;
300     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
301     snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
302     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
303     snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
304     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
305     snprintf(key,55,"hSingleDisp_cen%d",cent) ;
306     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
307     snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
308     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
309     snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
310     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
311     snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
312     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
313     snprintf(key,55,"hSingleCPV_cen%d",cent) ;
314     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
315     snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
316     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
317     snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
318     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
319     snprintf(key,55,"hSingleBoth_cen%d",cent) ;
320     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
321     snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
322     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
323
324
325     snprintf(key,55,"hMiPi0All_cen%d",cent) ;
326     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
327     snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
328     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
329     snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
330     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
331     snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
332     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
333     snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
334     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
335     snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
336     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
337     snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
338     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
339     snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
340     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
341     snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
342     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
343     snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
344     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
345     snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
346     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
347     snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
348     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
349
350     snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
351     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
352     snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
353     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
354     snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
355     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
356     snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
357     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
358     snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
359     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
360
361     snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
362     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
363     snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
364     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
365     snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
366     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
367     snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
368     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
369     snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
370     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
371     snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
372     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
373     snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
374     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
375     snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
376     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
377     snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
378     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
379     snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
380     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
381     snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
382     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
383     snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
384     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
385   }
386
387
388
389   const Int_t nPt      = 20;
390   const Double_t xPt[21]={0.6,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.} ;
391   const Int_t nPhi=10 ;
392   Double_t xPhi[nPhi+1] ;
393   for(Int_t i=0;i<=nPhi;i++)
394     xPhi[i]=i*TMath::Pi() /nPhi ;
395   const Int_t nMm=200 ;
396   Double_t xM[nMm+1] ;
397   for(Int_t i=0;i<=nMm;i++)
398     xM[i]=i*0.5 /nMm;
399
400   char phiTitle[15] ;
401   for(Int_t iRP=0; iRP<3; iRP++){
402     if(iRP==0)
403       snprintf(phiTitle,15,"TPC") ;
404     if(iRP==1)
405       snprintf(phiTitle,15,"V0A") ;
406     if(iRP==2)
407       snprintf(phiTitle,15,"V0C") ;
408     for(Int_t cent=0; cent<fCentEdges.GetSize(); cent++){
409       snprintf(key,55,"hPhotPhi%sAll_cen%d",phiTitle,cent) ;
410       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
411       snprintf(key,55,"hPhotPhi%sAllcore_cen%d",phiTitle,cent) ;
412       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
413       snprintf(key,55,"hPhotPhi%sDisp_cen%d",phiTitle,cent) ;
414       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
415       snprintf(key,55,"hPhotPhi%sDispcore_cen%d",phiTitle,cent) ;
416       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
417       snprintf(key,55,"hPhotPhi%sCPV_cen%d",phiTitle,cent) ;
418       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
419       snprintf(key,55,"hPhotPhi%sCPVcore_cen%d",phiTitle,cent) ;
420       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
421       snprintf(key,55,"hPhotPhi%sBoth_cen%d",phiTitle,cent) ;
422       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
423       snprintf(key,55,"hPhotPhi%sBothcore_cen%d",phiTitle,cent) ;
424       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
425
426       //Pions
427       snprintf(key,55,"hMassPt%sAll_cen%d",phiTitle,cent) ;
428       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
429       snprintf(key,55,"hMassPt%sAllcore_cen%d",phiTitle,cent) ;
430       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
431       snprintf(key,55,"hMassPt%sCPV_cen%d",phiTitle,cent) ;
432       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
433       snprintf(key,55,"hMassPt%sCPVcore_cen%d",phiTitle,cent) ;
434       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
435       snprintf(key,55,"hMassPt%sDisp_cen%d",phiTitle,cent) ;
436       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
437       snprintf(key,55,"hMassPt%sDispcore_cen%d",phiTitle,cent) ;
438       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
439       snprintf(key,55,"hMassPt%sBoth_cen%d",phiTitle,cent) ;
440       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
441       snprintf(key,55,"hMassPt%sBothcore_cen%d",phiTitle,cent) ;
442       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
443
444       //Mixed
445       snprintf(key,55,"hMiMassPt%sAll_cen%d",phiTitle,cent) ;
446       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
447       snprintf(key,55,"hMiMassPt%sAllcore_cen%d",phiTitle,cent) ;
448       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
449       snprintf(key,55,"hMiMassPt%sCPV_cen%d",phiTitle,cent) ;
450       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
451       snprintf(key,55,"hMiMassPt%sCPVcore_cen%d",phiTitle,cent) ;
452       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
453       snprintf(key,55,"hMiMassPt%sDisp_cen%d",phiTitle,cent) ;
454       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
455       snprintf(key,55,"hMiMassPt%sDispcore_cen%d",phiTitle,cent) ;
456       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
457       snprintf(key,55,"hMiMassPt%sBoth_cen%d",phiTitle,cent) ;
458       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
459       snprintf(key,55,"hMiMassPt%sBothcore_cen%d",phiTitle,cent) ;
460       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
461     }
462   }
463   
464   // Setup photon lists
465   Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
466   fCaloPhotonsPHOSLists = new TObjArray(kapacity);
467   fCaloPhotonsPHOSLists->SetOwner();
468   
469   PostData(1, fOutputContainer);
470 }
471
472 //________________________________________________________________________
473 void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
474 {
475   // Main loop, called for each event
476   // Analyze ESD/AOD
477
478
479   // Step 0: Event Objects
480   fEvent = GetEvent();
481   fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
482   fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
483   fMCStack = GetMCStack();
484   LogProgress(0, ConvertToInternalRunNumber(fEvent->GetRunNumber()));
485
486
487   // Step 1: Run Number, Misalignment Matrix, and Calibration
488   // fRunNumber, fInternalRunNumber, fMultV0, fV0Cpol, fV0Apol, fMeanQ, fWidthQ
489   if( fRunNumber != fEvent->GetRunNumber()) { // Check run number
490     // this should run only at first call of UserExec(),
491     // or if task runs over multiple runs, which should not occur in normal use.
492
493     // if run number has changed, set run variables
494     fRunNumber = fEvent->GetRunNumber();
495     fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
496     // then set misalignment and V0 calibration
497     SetGeometry();
498     SetMisalignment();
499     SetV0Calibration();
500     SetESDTrackCuts();
501     SetPHOSCalibData();
502   }
503   LogProgress(1, fInternalRunNumber);
504
505
506   // Step 2: Vertex
507   // fVertex, fVertexVector, fVtxBin
508   SetVertex();
509   if( RejectEventVertex() ) {
510     PostData(1, fOutputContainer);
511     return; // Reject!
512   }
513   LogProgress(2, fInternalRunNumber);
514
515 // Step 3:
516 //   if(event->IsPileupFromSPD()){
517 //     PostData(1, fOutputContainer);
518 //     return; // Reject!
519 //   }
520   LogProgress(3, fInternalRunNumber);
521
522
523   // Step 4: Centrality
524   // fCentralityV0M, fCentBin
525   SetCentrality();
526   if( RejectCentrality() ){
527     PostData(1, fOutputContainer);
528     return; // Reject! 
529   }
530   LogProgress(4, fInternalRunNumber);
531
532
533   // Step 5: Reaction Plane
534   // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
535   EvalReactionPlane(); //TODO: uncomment this, or at least deal with it
536   EvalV0ReactionPlane(); //TODO: uncomment this, or at least deal with it
537   fEMRPBin = GetRPBin(); //TODO: uncomment this, or at least deal with it
538   LogProgress(5, fInternalRunNumber);
539
540
541   // Step 6: MC
542   //  ProcessMC() ;
543   LogProgress(6, fInternalRunNumber);
544
545
546   // Step 7: QA PHOS cells
547   FillPHOSCellQAHists();
548   LogProgress(7, fInternalRunNumber);
549
550
551   // Step 8: Event Photons (PHOS Clusters) selection
552   SelectPhotonClusters();
553   FillSelectedClusterHistograms();
554   LogProgress(8, fInternalRunNumber);
555
556   // Step 9: Consider pi0 (photon/cluster) pairs.
557   ConsiderPi0s();
558   LogProgress(9, fInternalRunNumber);
559
560   // Step 10; Mixing
561   ConsiderPi0sMix();
562   LogProgress(10, fInternalRunNumber);
563   
564   // Step 11: Update lists
565   UpdateLists();
566   LogProgress(11,fInternalRunNumber);
567
568   
569   // Post output data.
570   PostData(1, fOutputContainer);
571 }
572
573 //________________________________________________________________________
574 // void AliAnalysisTaskPi0Flow::Terminate(Option_t *)
575 // {
576 //   // Draw result to the screen
577 //   // Called once at the end of the query
578 //   // new TCanvas;
579 //   // TH1 * hTotSelEvents = dynamic_cast<TH1*>(fOutputContainer->FindObject("hTotSelEvents"));
580 //   // hTotSelEvents->Draw();
581 // }
582 //________________________________________________________________________
583 void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges)
584 {
585   // Define centrality bins by their edges
586
587   int last = edges.GetSize()-1;
588   if( edges.At(0) < 0.) 
589     AliError("lower edge less then 0");
590   if( 90. < edges.At(last)  )
591     AliError("upper edge larger then 90.");
592   for(int i=0; i<last-1; ++i)
593     if(edges.At(i) > edges.At(i+1))
594       AliError("edges are not sorted");
595   
596   fCentEdges = edges;
597 }
598 //________________________________________________________________________
599 void AliAnalysisTaskPi0Flow::SetNMixedPerCentrality(const TArrayI& nMixed)
600 {
601   // Set number of mixed events for all centrality bins
602
603   fCentNMixed = nMixed;
604 }
605
606 //________________________________________________________________________
607 void AliAnalysisTaskPi0Flow::SetPHOSBadMap(Int_t mod, TH2I* badMapHist)
608   {
609     if(fPHOSBadMap[mod])
610       delete fPHOSBadMap[mod];
611
612     fPHOSBadMap[mod]=new TH2I(*badMapHist);
613     if(fDebug)
614       AliInfo(Form("Setting Bad Map Histogram  %s",fPHOSBadMap[mod]->GetName()));
615   }
616
617 //________________________________________________________________________
618 Bool_t AliAnalysisTaskPi0Flow::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
619 {
620   //Check if this channel belogs to the good ones
621
622   if(strcmp(det,"PHOS")==0){
623     if(mod>5 || mod<1){
624       AliError(Form("No bad map for PHOS module %d",mod)) ;
625       return kTRUE ;
626     }
627     if(!fPHOSBadMap[mod]){
628       AliError(Form("No Bad map for PHOS module %d, !fPHOSBadMap[mod]",mod)) ;
629       return kTRUE ;
630     }
631     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
632       return kFALSE ;
633     else
634       return kTRUE ;
635   }
636   else{
637     AliError(Form("Can not find bad channels for detector %s ",det)) ;
638   }
639
640   //Remove 6 noisy channels in run 139036, LHC10h
641   if( 139036 == fRunNumber
642     && mod==1
643     && (ix==9||ix==10||ix==11)
644     && (iz==45 || iz==46))
645     return kFALSE;
646
647   return kTRUE ;
648 }
649 //_____________________________________________________________________________
650 void AliAnalysisTaskPi0Flow::FillPHOSCellQAHists()
651 {
652   AliVCaloCells * cells = fEvent->GetPHOSCells();
653
654   FillHistogram("hCenPHOSCells",fCentralityV0M,cells->GetNumberOfCells()) ;
655   FillHistogram("hCenTrack",fCentralityV0M,fEvent->GetNumberOfTracks()) ;
656
657
658   Int_t nCellModule[3] = {0,0,0};
659   for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
660     Int_t cellAbsId = cells->GetCellNumber(iCell);
661     Int_t relId[4] = {0,0,0,0};
662     fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
663     Int_t mod1  = relId[0];
664     Int_t cellX = relId[2];
665     Int_t cellZ = relId[3] ;
666     Float_t energy = cells->GetAmplitude(iCell);
667     FillHistogram("hCellEnergy",energy);
668     if(mod1==1) {
669       nCellModule[0]++;
670       FillHistogram("hCellEnergyM1",cells->GetAmplitude(iCell));
671       FillHistogram("hCellNXZM1",cellX,cellZ,1.);
672       FillHistogram("hCellEXZM1",cellX,cellZ,energy);
673     }
674     else if (mod1==2) {
675       nCellModule[1]++;
676       FillHistogram("hCellEnergyM2",cells->GetAmplitude(iCell));
677       FillHistogram("hCellNXZM2",cellX,cellZ,1.);
678       FillHistogram("hCellEXZM2",cellX,cellZ,energy);
679     }
680     else if (mod1==3) {
681       nCellModule[2]++;
682       FillHistogram("hCellEnergyM3",cells->GetAmplitude(iCell));
683       FillHistogram("hCellNXZM3",cellX,cellZ,1.);
684       FillHistogram("hCellEXZM3",cellX,cellZ,energy);
685     }
686   }
687   FillHistogram("hCellMultEventM1",nCellModule[0]);
688   FillHistogram("hCellMultEventM2",nCellModule[1]);
689   FillHistogram("hCellMultEventM3",nCellModule[2]);
690
691 }
692 //_____________________________________________________________________________
693 void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
694 {
695   // clear (or create) array for holding events photons/clusters
696   if(fCaloPhotonsPHOS)
697     fCaloPhotonsPHOS->Clear();
698   else {
699     fCaloPhotonsPHOS = new TObjArray(200);
700     fCaloPhotonsPHOS->SetOwner();
701   }
702
703   
704   AliESDCaloCells* esdCells = dynamic_cast<AliESDCaloCells*> (fEvent->GetPHOSCells());
705   for (Int_t i=0; i<fEvent->GetNumberOfCaloClusters(); i++) {
706     AliVCluster *clu = fEvent->GetCaloCluster(i);
707
708     if ( !clu->IsPHOS() || clu->E()< kMinClusterEnergy) continue; // reject cluster
709
710     // check if cell/channel is good.
711     Float_t  position[3];
712     clu->GetPosition(position);
713     TVector3 global(position) ;
714     Int_t relId[4] ;
715     fPHOSGeo->GlobalPos2RelId(global,relId) ;
716     Int_t mod  = relId[0] ;
717     Int_t cellX = relId[2];
718     Int_t cellZ = relId[3] ;
719     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
720       continue ; // reject if not.
721
722
723     FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
724
725     if(clu->GetNCells() < kMinNCells) continue ;
726     if(clu->GetM02() < kMinM02)   continue ;
727
728     TLorentzVector lorentzMomentum;
729     Double_t ecore;
730
731     //if ESD, Apply re-Calibreation
732     Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0)
733     if( fEventESD ) {
734       AliPHOSEsdCluster cluPHOS1( *(AliESDCaloCluster*) (clu) );
735       cluPHOS1.Recalibrate(fPHOSCalibData, esdCells); // modify the cell energies
736       Reclusterize(&cluPHOS1) ;
737       cluPHOS1.EvalAll(kLogWeight, fVertexVector);         // recalculate the cluster parameters
738       cluPHOS1.SetE(fNonLinCorr->Eval(cluPHOS1.E()));// Users's nonlinearity
739
740       if(cluPHOS1.E()<0.3) continue; // check energy again
741
742       //correct misalignment
743       TVector3 localPos;
744       const Float_t shiftX[6]={0.,-2.3,-2.11,-1.53,0.,0.} ;
745       const Float_t shiftZ[6]={0.,-0.4, 0.52, 0.8,0.,0.} ;
746       fPHOSGeo->Global2Local(localPos,global,mod) ;
747       fPHOSGeo->Local2Global(mod,localPos.X()+shiftX[mod],localPos.Z()+shiftZ[mod],global);
748       position[0]=global.X() ;
749       position[1]=global.Y() ;
750       position[2]=global.Z() ;
751       cluPHOS1.SetPosition(position);
752
753       cluPHOS1.GetMomentum(lorentzMomentum ,origo);
754       ecore = CoreEnergy(&cluPHOS1);
755       
756       //TODO: Check, this may be LHC10h specific:
757       if(mod==2) lorentzMomentum*=135.5/134.0 ;
758       if(mod==3) lorentzMomentum*=135.5/137.2 ;
759       if(mod==2) ecore*=135.5/134.0 ;
760       if(mod==3) ecore*=135.5/137.2 ;
761         
762     }
763     else if (fEventAOD) { // is ! ESD, AOD.
764       AliESDCaloCluster* aodCluster = (AliESDCaloCluster*) (clu);
765       aodCluster->GetMomentum(lorentzMomentum ,origo);
766       ecore = CoreEnergy(clu);
767     }
768     else {
769       AliError("(Calo)Cluster is neither ESD nor AOD");
770       continue;
771     }
772
773
774     char skey[55];
775     snprintf(skey,55,"hCluLowM%d",mod) ;
776     FillHistogram(skey,cellX,cellZ,1.);
777     if(lorentzMomentum.E()>1.5){
778       sprintf(skey,"hCluHighM%d",mod) ;
779       FillHistogram(skey,cellX,cellZ,1.);
780     }
781
782     fCaloPhotonsPHOS->Add(new  AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
783     AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
784
785     ph->SetModule(mod) ;
786     lorentzMomentum*= ecore/lorentzMomentum.E() ;
787     ph->SetMomV2(&lorentzMomentum) ;
788     ph->SetNCells(clu->GetNCells());
789     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
790     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
791     if(ph->IsDispOK()){
792       sprintf(skey,"hCluDispM%d",mod) ;
793       FillHistogram(skey,cellX,cellZ,1.);
794     }
795
796     // Track Matching
797     Double_t dx=clu->GetTrackDx() ;
798     Double_t dz=clu->GetTrackDz() ;
799     Bool_t cpvBit=kTRUE ; //No track matched by default
800     Bool_t cpvBit2=kTRUE ; //More Strict criterion
801     if( fEventESD ) {
802       
803       TArrayI * itracks = static_cast<AliESDCaloCluster*> (clu)->GetTracksMatched() ;
804       if(itracks->GetSize()>0){
805         Int_t iTr = itracks->At(0);
806         if(iTr>=0 && iTr<fEvent->GetNumberOfTracks()){
807           AliVParticle* track = fEvent->GetTrack(iTr);
808           Double_t pt = track->Pt() ;
809           Short_t charge = track->Charge() ;
810           Double_t r=TestCPV(dx, dz, pt, charge) ;
811           cpvBit=(r>2.) ;
812           cpvBit2=(r>4.) ;
813         }
814       }
815     }
816     else if ( fEventAOD ) {
817       int nTracksMatched = clu->GetNTracksMatched();
818       if(nTracksMatched > 0) {
819         AliVTrack* track = dynamic_cast<AliVTrack*> (clu->GetTrackMatched(0));
820         if ( track ) {
821           Double_t pt = track->Pt();
822           Short_t charge = track->Charge();
823           Double_t r = TestCPV(dx, dz, pt, charge) ;
824           cpvBit=(r>2.) ;
825           cpvBit2=(r>4.) ;
826         }
827       }
828     }
829     ph->SetCPVBit(cpvBit) ;
830     ph->SetCPV2Bit(cpvBit2) ;
831     if(cpvBit){
832       sprintf(skey,"hCluVetoM%d",mod) ;
833       FillHistogram(skey,cellX,cellZ,1.);
834     }
835     ph->SetEMCx(float(cellX)) ;
836     ph->SetEMCz(float(cellZ)) ;
837     //    ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
838     ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolde
839
840   }
841   FillHistogram("hCenPHOS",fCentralityV0M, fCaloPhotonsPHOS->GetEntriesFast()) ;
842 }
843 //_____________________________________________________________________________
844 void AliAnalysisTaskPi0Flow::FillSelectedClusterHistograms()
845 {
846   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
847     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
848
849     Double_t dphiA=ph1->Phi()-fRPV0A ;
850     while(dphiA<0)dphiA+=TMath::Pi() ;
851     while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
852
853     Double_t dphiC=ph1->Phi()-fRPV0C ;
854     while(dphiC<0)dphiC+=TMath::Pi() ;
855     while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
856
857     Double_t dphiT=ph1->Phi()-fRP ;
858     while(dphiT<0)dphiT+=TMath::Pi() ;
859     while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
860
861     FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),ph1->Pt(),dphiA) ;
862     FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),ph1->Pt(),dphiC) ;
863     if(fHaveTPCRP)
864       FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),ph1->Pt(),dphiT) ;
865
866     FillHistogram(Form("hPhotAll_cen%d",fCentBin),ph1->Pt()) ;
867     FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
868     if(ph1->IsntUnfolded()){
869       FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),ph1->Pt()) ;
870     }
871     if(ph1->IsCPVOK()){
872       FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),ph1->Pt(),dphiA) ;
873       FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),ph1->Pt(),dphiC) ;
874       if(fHaveTPCRP)
875         FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),ph1->Pt(),dphiT) ;
876
877       FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
878       FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
879       if(fHaveTPCRP)
880         FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
881
882       FillHistogram(Form("hPhotCPV_cen%d",fCentBin),ph1->Pt()) ;
883       FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
884     }
885     if(ph1->IsCPV2OK()){
886       FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),ph1->Pt()) ;
887     }
888     if(ph1->IsDispOK()){
889       FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),ph1->Pt(),dphiA) ;
890       FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),ph1->Pt(),dphiC) ;
891       if(fHaveTPCRP)
892         FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),ph1->Pt(),dphiT) ;
893
894       FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
895       FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
896       if(fHaveTPCRP)
897         FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
898
899       FillHistogram(Form("hPhotDisp_cen%d",fCentBin),ph1->Pt()) ;
900       if(ph1->IsDisp2OK()){
901         FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),ph1->Pt()) ;
902       }
903       FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
904       if(ph1->IsntUnfolded()){
905         FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),ph1->Pt()) ;
906       }
907       if(ph1->IsCPVOK()){
908         FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),ph1->Pt(),dphiA) ;
909         FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),ph1->Pt(),dphiC) ;
910         if(fHaveTPCRP)
911           FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),ph1->Pt(),dphiT) ;
912
913         FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
914         FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
915         if(fHaveTPCRP)
916           FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
917
918         FillHistogram(Form("hPhotBoth_cen%d",fCentBin),ph1->Pt()) ;
919         FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
920       }
921     }
922   }
923 }
924 //_____________________________________________________________________________
925 void AliAnalysisTaskPi0Flow::ConsiderPi0s()
926 {
927   char key[55];
928   for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
929     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
930     for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
931       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
932       TLorentzVector p12  = *ph1  + *ph2;
933       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
934       FillHistogram("hPHOSphi",fCentralityV0M,p12.Pt(),p12.Phi());
935       Double_t dphiA=p12.Phi()-fRPV0A ;
936       while(dphiA<0)dphiA+=TMath::Pi() ;
937       while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
938
939       Double_t dphiC=p12.Phi()-fRPV0C ;
940       while(dphiC<0)dphiC+=TMath::Pi() ;
941       while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
942
943       Double_t dphiT=p12.Phi()-fRP ;
944       while(dphiT<0)dphiT+=TMath::Pi() ;
945       while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
946
947       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
948
949       FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
950       FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
951       if(fHaveTPCRP)
952         FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
953
954       FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
955       FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
956       if(fHaveTPCRP)
957         FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
958
959
960       FillHistogram(Form("hPi0All_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
961       FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
962       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
963         FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
964       }
965
966       FillHistogram(Form("hSingleAll_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
967       FillHistogram(Form("hSingleAll_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
968       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
969       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
970       if(ph1->IsntUnfolded())
971         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
972       if(ph2->IsntUnfolded())
973         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
974       if(ph1->IsCPVOK()){
975         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
976         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
977       }
978       if(ph2->IsCPVOK()){
979         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
980         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
981       }
982       if(ph1->IsCPV2OK()){
983         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
984       }
985       if(ph2->IsCPV2OK()){
986         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
987       }
988       if(ph1->IsDispOK()){
989         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
990         if(ph1->IsntUnfolded()){
991           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
992         }
993         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
994       }
995       if(ph2->IsDispOK()){
996         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
997         if(ph1->IsntUnfolded()){
998           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
999         }
1000         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1001       }
1002       if(ph1->IsDisp2OK()){
1003         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1004       }
1005       if(ph2->IsDisp2OK()){
1006         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1007       }
1008       if(ph1->IsDispOK() && ph1->IsCPVOK()){
1009         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1010         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1011       }
1012       if(ph2->IsDispOK() && ph2->IsCPVOK()){
1013         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1014         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1015       }
1016
1017
1018       if(a<kAlphaCut){
1019         FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1020       }
1021
1022       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1023         snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
1024         FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1025         FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1026         if(fHaveTPCRP)
1027           FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1028
1029         FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
1030         FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
1031         if(fHaveTPCRP)
1032           FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
1033
1034         FillHistogram(Form("hPi0CPV_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1035         FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1036
1037         if(a<kAlphaCut){
1038           FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1039         }
1040       }
1041       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1042         FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),p12.M(),p12.Pt()) ;
1043         if(a<kAlphaCut){
1044           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1045         }
1046       }
1047       if(ph1->IsDispOK() && ph2->IsDispOK()){
1048         snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
1049         FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1050         FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1051         if(fHaveTPCRP)
1052           FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1053
1054         FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiA) ;
1055         FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiC) ;
1056         if(fHaveTPCRP)
1057           FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiT) ;
1058
1059         FillHistogram(Form("hPi0Disp_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1060         FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1061         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1062           FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1063         }
1064         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1065           FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),p12.M(), p12.Pt()) ;
1066         }
1067
1068         if(a<kAlphaCut){
1069           FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1070         }
1071
1072         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1073           FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1074           FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1075           if(fHaveTPCRP)
1076             FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1077
1078           FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
1079           FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
1080           if(fHaveTPCRP)
1081             FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
1082
1083           FillHistogram(Form("hPi0Both_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1084           FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
1085
1086           if(a<kAlphaCut){
1087             snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
1088             FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1089           }
1090           if(fCentralityV0M>20.){
1091           if(ph1->Module()==1 && ph2->Module()==1)
1092             FillHistogram("hPi0M11",p12.M(),p12.Pt() );
1093           else if(ph1->Module()==2 && ph2->Module()==2)
1094             FillHistogram("hPi0M22",p12.M(),p12.Pt() );
1095           else if(ph1->Module()==3 && ph2->Module()==3)
1096             FillHistogram("hPi0M33",p12.M(),p12.Pt() );
1097           else if(ph1->Module()==1 && ph2->Module()==2)
1098             FillHistogram("hPi0M12",p12.M(),p12.Pt() );
1099           else if(ph1->Module()==1 && ph2->Module()==3)
1100             FillHistogram("hPi0M13",p12.M(),p12.Pt() );
1101           else if(ph1->Module()==2 && ph2->Module()==3)
1102             FillHistogram("hPi0M23",p12.M(),p12.Pt() );
1103           }
1104
1105         }
1106       }
1107     } // end of loop i2
1108   } // end of loop i1
1109 }
1110 //_____________________________________________________________________________
1111 void AliAnalysisTaskPi0Flow::ConsiderPi0sMix()
1112 {
1113   char key[55];
1114
1115   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1116
1117   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
1118     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1119     for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1120       TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
1121       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1122         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1123         TLorentzVector p12  = *ph1  + *ph2;
1124         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1125
1126         Double_t dphiA=p12.Phi()-fRPV0A ;
1127         while(dphiA<0)dphiA+=TMath::Pi() ;
1128         while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
1129
1130         Double_t dphiC=p12.Phi()-fRPV0C ;
1131         while(dphiC<0)dphiC+=TMath::Pi() ;
1132         while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
1133
1134         Double_t dphiT=p12.Phi()-fRP ;
1135         while(dphiT<0)dphiT+=TMath::Pi() ;
1136         while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
1137
1138
1139         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1140
1141         snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
1142         FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1143         FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1144         if(fHaveTPCRP)
1145           FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1146
1147         FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiA) ;
1148         FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiC) ;
1149         if(fHaveTPCRP)
1150           FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiT) ;
1151
1152         FillHistogram(Form("hMiPi0All_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1153         FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
1154         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1155           FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1156         }
1157
1158         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1159         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1160         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1161         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1162         if(ph1->IsntUnfolded())
1163           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1164         if(ph2->IsntUnfolded())
1165           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1166         if(ph1->IsCPVOK()){
1167           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1168           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1169         }
1170         if(ph2->IsCPVOK()){
1171           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1172           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1173         }
1174         if(ph1->IsCPV2OK()){
1175           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1176         }
1177         if(ph2->IsCPV2OK()){
1178           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1179         }
1180         if(ph1->IsDispOK()){
1181           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1182           if(ph1->IsntUnfolded()){
1183             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1184           }
1185           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1186         }
1187         if(ph2->IsDispOK()){
1188           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1189           if(ph1->IsntUnfolded()){
1190             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1191           }
1192           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1193         }
1194         if(ph1->IsDisp2OK()){
1195           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1196         }
1197         if(ph2->IsDisp2OK()){
1198           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1199         }
1200         if(ph1->IsDispOK() && ph1->IsCPVOK()){
1201           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1202           FillHistogram(key,p12.M(),ph1->Pt()) ;
1203           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1204           FillHistogram(key,pv12.M(),ph1->GetMomV2()->Pt()) ;
1205         }
1206         if(ph2->IsDispOK() && ph2->IsCPVOK()){
1207           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1208           FillHistogram(key,p12.M(),ph2->Pt()) ;
1209           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1210           FillHistogram(key,pv12.M(),ph2->GetMomV2()->Pt()) ;
1211         }
1212
1213
1214
1215         if(a<kAlphaCut){
1216           snprintf(key,55,"hMiPi0All_a07_cen%d",fCentBin) ;
1217           FillHistogram(key,p12.M() ,p12.Pt()) ;
1218         }
1219         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1220           FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiA) ;
1221           FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiC) ;
1222           if(fHaveTPCRP)
1223             FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiT) ;
1224
1225           FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiA) ;
1226           FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiC) ;
1227           if(fHaveTPCRP)
1228             FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiT) ;
1229
1230           FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),p12.M(), p12.Pt()) ;
1231           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1232
1233           if(a<kAlphaCut){
1234             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1235           }
1236         }
1237         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1238           FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1239
1240           if(a<kAlphaCut){
1241             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1242           }
1243         }
1244         if(ph1->IsDispOK() && ph2->IsDispOK()){
1245           FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiA) ;
1246           FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiC) ;
1247           if(fHaveTPCRP)
1248             FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiT) ;
1249
1250           FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1251           FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1252           if(fHaveTPCRP)
1253             FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1254
1255
1256           FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1257           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1258           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1259             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),p12.M(),p12.Pt()) ;
1260           }
1261
1262           if(a<kAlphaCut){
1263             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1264           }
1265           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1266             FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiA) ;
1267             FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiC) ;
1268             if(fHaveTPCRP)
1269               FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiT) ;
1270
1271             FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1272             FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1273             if(fHaveTPCRP)
1274               FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1275
1276             FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1277             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1278
1279             if(a<kAlphaCut){
1280               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1281             }
1282           }
1283         }
1284       } // end of loop i2
1285     }
1286   } // end of loop i1
1287 }
1288 //_____________________________________________________________________________
1289 void AliAnalysisTaskPi0Flow::UpdateLists()
1290 {
1291   //Now we either add current events to stack or remove
1292   //If no photons in current event - no need to add it to mixed
1293
1294   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1295   if( fDebug )
1296     AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1297   if(fCaloPhotonsPHOS->GetEntriesFast()>0){
1298     arrayList->AddFirst(fCaloPhotonsPHOS) ;
1299     fCaloPhotonsPHOS=0;
1300     if(arrayList->GetEntries() > fCentNMixed[fCentBin]){ // Remove redundant events
1301       TObjArray * tmp = static_cast<TObjArray*>(arrayList->Last()) ;
1302       arrayList->RemoveLast() ;
1303       delete tmp ; // TODO: may conflict with delete done by list being owner.
1304     }
1305   }
1306   else
1307     fCaloPhotonsPHOS->Clear(); // TODO: redundant???
1308 }
1309 //_____________________________________________________________________________
1310 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x)const{
1311   //FillHistogram
1312   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
1313   if(tmpI){
1314     tmpI->Fill(x) ;
1315     return ;
1316   }
1317   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
1318   if(tmpF){
1319     tmpF->Fill(x) ;
1320     return ;
1321   }
1322   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
1323   if(tmpD){
1324     tmpD->Fill(x) ;
1325     return ;
1326   }
1327   AliInfo(Form("can not find histogram <%s> ",key)) ;
1328 }
1329 //_____________________________________________________________________________
1330 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y)const{
1331   //FillHistogram
1332   TObject * tmp = fOutputContainer->FindObject(key) ;
1333   if(!tmp){
1334     AliInfo(Form("can not find histogram <%s> ",key)) ;
1335     return ;
1336   }
1337   if(tmp->IsA() == TClass::GetClass("TH1F")){
1338     ((TH1F*)tmp)->Fill(x,y) ;
1339     return ;
1340   }
1341   if(tmp->IsA() == TClass::GetClass("TH2F")){
1342     ((TH2F*)tmp)->Fill(x,y) ;
1343     return ;
1344   }
1345   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1346 }
1347
1348 //_____________________________________________________________________________
1349 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1350   //Fills 1D histograms with key
1351   TObject * tmp = fOutputContainer->FindObject(key) ;
1352   if(!tmp){
1353     AliInfo(Form("can not find histogram <%s> ",key)) ;
1354     return ;
1355   }
1356   if(tmp->IsA() == TClass::GetClass("TH2F")){
1357     ((TH2F*)tmp)->Fill(x,y,z) ;
1358     return ;
1359   }
1360   if(tmp->IsA() == TClass::GetClass("TH3F")){
1361     ((TH3F*)tmp)->Fill(x,y,z) ;
1362     return ;
1363   }
1364 }
1365
1366 //_____________________________________________________________________________
1367 AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
1368 {
1369   fEvent = InputEvent();
1370   if( ! fEvent ) {
1371     AliError("Event could not be retrieved");
1372     PostData(1, fOutputContainer);
1373   }
1374   return fEvent;
1375 }
1376
1377
1378 //___________________________________________________________________________
1379 AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
1380 {
1381   fMCStack = 0;
1382   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
1383   if(eventHandler){
1384     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
1385     if( mcEventHandler)
1386       fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
1387   }
1388   return fMCStack;
1389 }
1390
1391 //___________________________________________________________________________
1392 Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
1393 {
1394  /* fCentBin=1+Int_t(centralityV0M/100. *kNCenBins) ;
1395   if(centralityV0M < 5. || fCentBin < 0)
1396    fCentBin=0 ;
1397   if(fCentBin > kNCenBins-1)
1398     fCentBin = kNCenBins-1 ;
1399  */
1400   int lastBinUpperIndex = fCentEdges.GetSize() -1;
1401   if( centralityV0M > fCentEdges[lastBinUpperIndex] ) {
1402     if(fDebug)
1403       AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1404     return lastBinUpperIndex-1;
1405   }
1406   if( centralityV0M < fCentEdges[0] ) {
1407     if( fDebug )
1408       AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1409     return 0;
1410   }
1411   
1412   fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1413   return fCentBin;
1414 }
1415
1416 //___________________________________________________________________________
1417 Int_t AliAnalysisTaskPi0Flow::GetRPBin()
1418 {
1419   Double_t averageRP;
1420   if(fHaveTPCRP)
1421     averageRP = fRPV0A+fRPV0C+fRP /3.;
1422   else
1423     averageRP = fRPV0A+fRPV0C /2.;
1424
1425   fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1426
1427   if(fEMRPBin> (Int_t) fNEMRPBins-1)
1428     fEMRPBin=fNEMRPBins-1 ;
1429   else if(fEMRPBin<0)
1430     fEMRPBin=0;
1431
1432   if ( fDebug )
1433     AliInfo(Form("Event Mixing Reaction Plane bin is: %d", fEMRPBin));
1434
1435   return fEMRPBin;
1436 }
1437
1438
1439 //_____________________________________________________________________________
1440 void AliAnalysisTaskPi0Flow::LogProgress(int step, int internalRunNumber)
1441 {
1442   if(fDebug > 1) {
1443     AliInfo(Form("step %d completed", step));
1444   }
1445   // the +0.5 is not realy neccisarry, but oh well... -henrik
1446   FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1447   FillHistogram("hTotSelEvents", step+0.5);
1448 }
1449
1450 //___________________________________________________________________________
1451 Int_t AliAnalysisTaskPi0Flow::ConvertToInternalRunNumber(Int_t run){
1452   if( kLHC11h == fPeriod ) {
1453     switch(run){
1454     case  170593 : return 179 ;
1455     case  170572 : return 178 ;
1456     case  170556 : return 177 ;
1457     case  170552 : return 176 ;
1458     case  170546 : return 175 ;
1459     case  170390 : return 174 ;
1460     case  170389 : return 173 ;
1461     case  170388 : return 172 ;
1462     case  170387 : return 171 ;
1463     case  170315 : return 170 ;
1464     case  170313 : return 169 ;
1465     case  170312 : return 168 ;
1466     case  170311 : return 167 ;
1467     case  170309 : return 166 ;
1468     case  170308 : return 165 ;
1469     case  170306 : return 164 ;
1470     case  170270 : return 163 ;
1471     case  170269 : return 162 ;
1472     case  170268 : return 161 ;
1473     case  170267 : return 160 ;
1474     case  170264 : return 159 ;
1475     case  170230 : return 158 ;
1476     case  170228 : return 157 ;
1477     case  170208 : return 156 ;
1478     case  170207 : return 155 ;
1479     case  170205 : return 154 ;
1480     case  170204 : return 153 ;
1481     case  170203 : return 152 ;
1482     case  170195 : return 151 ;
1483     case  170193 : return 150 ;
1484     case  170163 : return 149 ;
1485     case  170162 : return 148 ;
1486     case  170159 : return 147 ;
1487     case  170155 : return 146 ;
1488     case  170152 : return 145 ;
1489     case  170091 : return 144 ;
1490     case  170089 : return 143 ;
1491     case  170088 : return 142 ;
1492     case  170085 : return 141 ;
1493     case  170084 : return 140 ;
1494     case  170083 : return 139 ;
1495     case  170081 : return 138 ;
1496     case  170040 : return 137 ;
1497     case  170038 : return 136 ;
1498     case  170036 : return 135 ;
1499     case  170027 : return 134 ;
1500     case  169981 : return 133 ;
1501     case  169975 : return 132 ;
1502     case  169969 : return 131 ;
1503     case  169965 : return 130 ;
1504     case  169961 : return 129 ;
1505     case  169956 : return 128 ;
1506     case  169926 : return 127 ;
1507     case  169924 : return 126 ;
1508     case  169923 : return 125 ;
1509     case  169922 : return 124 ;
1510     case  169919 : return 123 ;
1511     case  169918 : return 122 ;
1512     case  169914 : return 121 ;
1513     case  169859 : return 120 ;
1514     case  169858 : return 119 ;
1515     case  169855 : return 118 ;
1516     case  169846 : return 117 ;
1517     case  169838 : return 116 ;
1518     case  169837 : return 115 ;
1519     case  169835 : return 114 ;
1520     case  169683 : return 113 ;
1521     case  169628 : return 112 ;
1522     case  169591 : return 111 ;
1523     case  169590 : return 110 ;
1524     case  169588 : return 109 ;
1525     case  169587 : return 108 ;
1526     case  169586 : return 107 ;
1527     case  169584 : return 106 ;
1528     case  169557 : return 105 ;
1529     case  169555 : return 104 ;
1530     case  169554 : return 103 ;
1531     case  169553 : return 102 ;
1532     case  169550 : return 101 ;
1533     case  169515 : return 100 ;
1534     case  169512 : return 99 ;
1535     case  169506 : return 98 ;
1536     case  169504 : return 97 ;
1537     case  169498 : return 96 ;
1538     case  169475 : return 95 ;
1539     case  169420 : return 94 ;
1540     case  169419 : return 93 ;
1541     case  169418 : return 92 ;
1542     case  169417 : return 91 ;
1543     case  169415 : return 90 ;
1544     case  169411 : return 89 ;
1545     case  169238 : return 88 ;
1546     case  169236 : return 87 ;
1547     case  169167 : return 86 ;
1548     case  169160 : return 85 ;
1549     case  169156 : return 84 ;
1550     case  169148 : return 83 ;
1551     case  169145 : return 82 ;
1552     case  169144 : return 81 ;
1553     case  169143 : return 80 ;
1554     case  169138 : return 79 ;
1555     case  169099 : return 78 ;
1556     case  169094 : return 77 ;
1557     case  169091 : return 76 ;
1558     case  169045 : return 75 ;
1559     case  169044 : return 74 ;
1560     case  169040 : return 73 ;
1561     case  169035 : return 72 ;
1562     case  168992 : return 71 ;
1563     case  168988 : return 70 ;
1564     case  168984 : return 69 ;
1565     case  168826 : return 68 ;
1566     case  168777 : return 67 ;
1567     case  168514 : return 66 ;
1568     case  168512 : return 65 ;
1569     case  168511 : return 64 ;
1570     case  168467 : return 63 ;
1571     case  168464 : return 62 ;
1572     case  168461 : return 61 ;
1573     case  168460 : return 60 ;
1574     case  168458 : return 59 ;
1575     case  168362 : return 58 ;
1576     case  168361 : return 57 ;
1577     case  168356 : return 56 ;
1578     case  168342 : return 55 ;
1579     case  168341 : return 54 ;
1580     case  168325 : return 53 ;
1581     case  168322 : return 52 ;
1582     case  168318 : return 51 ;
1583     case  168311 : return 50 ;
1584     case  168310 : return 49 ;
1585     case  168213 : return 48 ;
1586     case  168212 : return 47 ;
1587     case  168208 : return 46 ;
1588     case  168207 : return 45 ;
1589     case  168206 : return 44 ;
1590     case  168205 : return 43 ;
1591     case  168204 : return 42 ;
1592     case  168203 : return 41 ;
1593     case  168181 : return 40 ;
1594     case  168177 : return 39 ;
1595     case  168175 : return 38 ;
1596     case  168173 : return 37 ;
1597     case  168172 : return 36 ;
1598     case  168171 : return 35 ;
1599     case  168115 : return 34 ;
1600     case  168108 : return 33 ;
1601     case  168107 : return 32 ;
1602     case  168105 : return 31 ;
1603     case  168104 : return 30 ;
1604     case  168103 : return 29 ;
1605     case  168076 : return 28 ;
1606     case  168069 : return 27 ;
1607     case  168068 : return 26 ;
1608     case  168066 : return 25 ;
1609     case  167988 : return 24 ;
1610     case  167987 : return 23 ;
1611     case  167986 : return 22 ;
1612     case  167985 : return 21 ;
1613     case  167921 : return 20 ;
1614     case  167920 : return 19 ;
1615     case  167915 : return 18 ;
1616     case  167909 : return 17 ;
1617     case  167903 : return 16 ;
1618     case  167902 : return 15 ;
1619     case  167818 : return 14 ;
1620     case  167814 : return 13 ;
1621     case  167813 : return 12 ;
1622     case  167808 : return 11 ;
1623     case  167807 : return 10 ;
1624     case  167806 : return 9 ;
1625     case  167713 : return 8 ;
1626     case  167712 : return 7 ;
1627     case  167711 : return 6 ;
1628     case  167706 : return 5 ;
1629     case  167693 : return 4 ;
1630     case  166532 : return 3 ;
1631     case  166530 : return 2 ;
1632     case  166529 : return 1 ;
1633
1634     default : return 199;
1635     }
1636   }
1637   if( kLHC10h == fPeriod ) {
1638     switch(run){
1639     case  139517 : return 137;
1640     case  139514 : return 136;
1641     case  139513 : return 135;
1642     case  139511 : return 134;
1643     case  139510 : return 133;
1644     case  139507 : return 132;
1645     case  139505 : return 131;
1646     case  139504 : return 130;
1647     case  139503 : return 129;
1648     case  139470 : return 128;
1649     case  139467 : return 127;
1650     case  139466 : return 126;
1651     case  139465 : return 125;
1652     case  139440 : return 124;
1653     case  139439 : return 123;
1654     case  139438 : return 122;
1655     case  139437 : return 121;
1656     case  139360 : return 120;
1657     case  139329 : return 119;
1658     case  139328 : return 118;
1659     case  139314 : return 117;
1660     case  139311 : return 116;
1661     case  139310 : return 115;
1662     case  139309 : return 114;
1663     case  139308 : return 113;
1664     case  139173 : return 112;
1665     case  139172 : return 111;
1666     case  139110 : return 110;
1667     case  139107 : return 109;
1668     case  139105 : return 108;
1669     case  139104 : return 107;
1670     case  139042 : return 106;
1671     case  139038 : return 105;
1672     case  139037 : return 104;
1673     case  139036 : return 103;
1674     case  139029 : return 102;
1675     case  139028 : return 101;
1676     case  138983 : return 100;
1677     case  138982 : return 99;
1678     case  138980 : return 98;
1679     case  138979 : return 97;
1680     case  138978 : return 96;
1681     case  138977 : return 95;
1682     case  138976 : return 94;
1683     case  138973 : return 93;
1684     case  138972 : return 92;
1685     case  138965 : return 91;
1686     case  138924 : return 90;
1687     case  138872 : return 89;
1688     case  138871 : return 88;
1689     case  138870 : return 87;
1690     case  138837 : return 86;
1691     case  138830 : return 85;
1692     case  138828 : return 84;
1693     case  138826 : return 83;
1694     case  138796 : return 82;
1695     case  138795 : return 81;
1696     case  138742 : return 80;
1697     case  138732 : return 79;
1698     case  138730 : return 78;
1699     case  138666 : return 77;
1700     case  138662 : return 76;
1701     case  138653 : return 75;
1702     case  138652 : return 74;
1703     case  138638 : return 73;
1704     case  138624 : return 72;
1705     case  138621 : return 71;
1706     case  138583 : return 70;
1707     case  138582 : return 69;
1708     case  138579 : return 68;
1709     case  138578 : return 67;
1710     case  138534 : return 66;
1711     case  138469 : return 65;
1712     case  138442 : return 64;
1713     case  138439 : return 63;
1714     case  138438 : return 62;
1715     case  138396 : return 61;
1716     case  138364 : return 60;
1717     case  138359 : return 59;
1718     case  138275 : return 58;
1719     case  138225 : return 57;
1720     case  138201 : return 56;
1721     case  138200 : return 55;
1722     case  138197 : return 54;
1723     case  138192 : return 53;
1724     case  138190 : return 52;
1725     case  138154 : return 51;
1726     case  138153 : return 50;
1727     case  138151 : return 49;
1728     case  138150 : return 48;
1729     case  138126 : return 47;
1730     case  138125 : return 46;
1731     case  137848 : return 45;
1732     case  137847 : return 44;
1733     case  137844 : return 43;
1734     case  137843 : return 42;
1735     case  137752 : return 41;
1736     case  137751 : return 40;
1737     case  137748 : return 39;
1738     case  137724 : return 38;
1739     case  137722 : return 37;
1740     case  137718 : return 36;
1741     case  137704 : return 35;
1742     case  137693 : return 34;
1743     case  137692 : return 33;
1744     case  137691 : return 32;
1745     case  137689 : return 31;
1746     case  137686 : return 30;
1747     case  137685 : return 29;
1748     case  137639 : return 28;
1749     case  137638 : return 27;
1750     case  137608 : return 26;
1751     case  137595 : return 25;
1752     case  137549 : return 24;
1753     case  137546 : return 23;
1754     case  137544 : return 22;
1755     case  137541 : return 21;
1756     case  137539 : return 20;
1757     case  137531 : return 19;
1758     case  137530 : return 18;
1759     case  137443 : return 17;
1760     case  137441 : return 16;
1761     case  137440 : return 15;
1762     case  137439 : return 14;
1763     case  137434 : return 13;
1764     case  137432 : return 12;
1765     case  137431 : return 11;
1766     case  137430 : return 10;
1767     case  137366 : return 9;
1768     case  137243 : return 8;
1769     case  137236 : return 7;
1770     case  137235 : return 6;
1771     case  137232 : return 5;
1772     case  137231 : return 4;
1773     case  137165 : return 3;
1774     case  137162 : return 2;
1775     case  137161 : return 1;
1776     default : return 199;
1777     }
1778   }
1779   if(kUndefinedPeriod && fDebug > 0) {
1780     AliWarning("Period not defined");
1781   }
1782   return 1;
1783 }
1784 //_____________________________________________________________________________
1785 Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1786
1787   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1788   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1789   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1790   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1791   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1792   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1793     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1794     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1795   return (R2<2.5*2.5) ;
1796
1797 }
1798 //_____________________________________________________________________________
1799 Bool_t AliAnalysisTaskPi0Flow::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1800
1801   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1802   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1803   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1804   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1805   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1806   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1807     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1808     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1809   return (R2<1.5*1.5) ;
1810
1811 }
1812 //____________________________________________________________________________
1813 TList* AliAnalysisTaskPi0Flow::GetCaloPhotonsPHOSList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin)
1814 {
1815   int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins
1816               + centBin * fNEMRPBins
1817               + rpBin;
1818   if( fCaloPhotonsPHOSLists->At(offset) ) { // list exists
1819     TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1820     if( ! list )
1821       AliError("object in fCaloPhotonsPHOSLists at %i did not cast");
1822     return list;
1823   }
1824   else {// no list for this bin has been created, yet
1825     TList* list = new TList();
1826     list->SetOwner();
1827     fCaloPhotonsPHOSLists->AddAt(list, offset);
1828     return list;
1829   }
1830 }
1831
1832 //____________________________________________________________________________
1833 Double_t AliAnalysisTaskPi0Flow::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1834   //Parameterization of LHC10h period
1835   //_true if neutral_
1836
1837   Double_t meanX=0;
1838   Double_t meanZ=0.;
1839   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1840                          6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1841   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1842   AliVEvent *event = InputEvent();
1843   Double_t mf = 0.; //
1844   if(event) mf = event->GetMagneticField(); //Positive for ++ and negative for --
1845
1846   if(mf<0.){ //field --
1847     meanZ = -0.468318 ;
1848     if(charge>0)
1849       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1850     else
1851       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1852   }
1853   else{ //Field ++
1854     meanZ= -0.468318;
1855     if(charge>0)
1856       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1857     else
1858       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1859   }
1860
1861   Double_t rz=(dz-meanZ)/sz ;
1862   Double_t rx=(dx-meanX)/sx ;
1863   return TMath::Sqrt(rx*rx+rz*rz) ;
1864 }
1865 //____________________________________________________________________________
1866 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0A(Double_t phi, Double_t c){
1867   //LHC10h
1868
1869   Double_t v2c =0.;
1870   Double_t v2s =0;
1871
1872 if(fInternalRunNumber==1){
1873   v2c=0.00271301 + (0.00340046)*c + (-0.000376369)*c*c + (8.56362e-06)*c*c*c + (-5.70955e-08)*c*c*c*c ;
1874   v2s=0.024278 + (0.000586593)*c + (-1.07748e-06)*c*c + (-5.01567e-07)*c*c*c + (3.65956e-09)*c*c*c*c ;
1875 }
1876 if(fInternalRunNumber==2){
1877   v2c=0.111374 + (-0.0111331)*c + (0.000212342)*c*c + (-9.6774e-07)*c*c*c + (-2.54199e-09)*c*c*c*c ;
1878   v2s=-0.0252636 + (0.0152193)*c + (-0.000830793)*c*c + (1.5099e-05)*c*c*c + (-9.04687e-08)*c*c*c*c ;
1879 }
1880 if(fInternalRunNumber==4){
1881   v2c=-0.0323376 + (0.0058672)*c + (-0.000257135)*c*c + (3.54353e-06)*c*c*c + (-1.50431e-08)*c*c*c*c ;
1882   v2s=0.0596758 + (-0.00520203)*c + (0.000178495)*c*c + (-2.94143e-06)*c*c*c + (1.55424e-08)*c*c*c*c ;
1883 }
1884 if(fInternalRunNumber==5){
1885   v2c=-0.0415369 + (0.0098273)*c + (-0.000486907)*c*c + (8.00188e-06)*c*c*c + (-4.29107e-08)*c*c*c*c ;
1886   v2s=-0.0146049 + (0.00693627)*c + (-0.000348663)*c*c + (5.79331e-06)*c*c*c + (-3.33179e-08)*c*c*c*c ;
1887 }
1888 if(fInternalRunNumber==6){
1889   v2c=0.251902 + (-0.0415419)*c + (0.0018553)*c*c + (-3.15181e-05)*c*c*c + (1.80073e-07)*c*c*c*c ;
1890   v2s=0.0771839 + (-0.00962468)*c + (0.000394453)*c*c + (-7.22485e-06)*c*c*c + (4.53562e-08)*c*c*c*c ;
1891 }
1892 if(fInternalRunNumber==7){
1893   v2c=0.157038 + (-0.0204377)*c + (0.000847477)*c*c + (-1.4232e-05)*c*c*c + (8.04965e-08)*c*c*c*c ;
1894   v2s=0.0256337 + (-0.00242333)*c + (0.000140455)*c*c + (-3.5149e-06)*c*c*c + (2.44854e-08)*c*c*c*c ;
1895 }
1896 if(fInternalRunNumber==8){
1897   v2c=-0.0916108 + (0.0146996)*c + (-0.000624971)*c*c + (9.63331e-06)*c*c*c + (-4.8574e-08)*c*c*c*c ;
1898   v2s=0.0332113 + (-0.00332475)*c + (0.000151857)*c*c + (-2.86187e-06)*c*c*c + (1.54796e-08)*c*c*c*c ;
1899 }
1900 if(fInternalRunNumber==9){
1901   v2c=0.00925348 + (1.58854e-05)*c + (-3.85178e-05)*c*c + (6.75167e-07)*c*c*c + (-3.1502e-09)*c*c*c*c ;
1902   v2s=0.0249253 + (-0.00300902)*c + (0.000152653)*c*c + (-3.20923e-06)*c*c*c + (2.08818e-08)*c*c*c*c ;
1903 }
1904 if(fInternalRunNumber==10){
1905   v2c=0.040903 + (0.00142711)*c + (-0.000264675)*c*c + (6.50624e-06)*c*c*c + (-4.46628e-08)*c*c*c*c ;
1906   v2s=0.0914196 + (-0.00957901)*c + (0.000408998)*c*c + (-7.49412e-06)*c*c*c + (4.63915e-08)*c*c*c*c ;
1907 }
1908 if(fInternalRunNumber==11){
1909   v2c=-0.0130133 + (0.00296474)*c + (-0.000183518)*c*c + (3.65755e-06)*c*c*c + (-2.24484e-08)*c*c*c*c ;
1910   v2s=-0.00449031 + (0.00365474)*c + (-0.000226455)*c*c + (4.49792e-06)*c*c*c + (-2.96018e-08)*c*c*c*c ;
1911 }
1912 if(fInternalRunNumber==12){
1913   v2c=-0.021987 + (0.00506675)*c + (-0.000301768)*c*c + (5.66088e-06)*c*c*c + (-3.21447e-08)*c*c*c*c ;
1914   v2s=0.0830858 + (-0.00871637)*c + (0.000299212)*c*c + (-4.10705e-06)*c*c*c + (1.76888e-08)*c*c*c*c ;
1915 }
1916 if(fInternalRunNumber==13){
1917   v2c=0.0502952 + (-0.00152478)*c + (-0.000106289)*c*c + (3.39398e-06)*c*c*c + (-2.45015e-08)*c*c*c*c ;
1918   v2s=0.0639936 + (-0.00775191)*c + (0.000347449)*c*c + (-6.46454e-06)*c*c*c + (4.00439e-08)*c*c*c*c ;
1919 }
1920 if(fInternalRunNumber==14){
1921   v2c=-0.116922 + (0.0122105)*c + (-0.000373923)*c*c + (4.08464e-06)*c*c*c + (-1.05659e-08)*c*c*c*c ;
1922   v2s=0.128821 + (-0.0204399)*c + (0.000831009)*c*c + (-1.33147e-05)*c*c*c + (7.24876e-08)*c*c*c*c ;
1923 }
1924 if(fInternalRunNumber==15){
1925   v2c=0.0853741 + (-0.00924528)*c + (0.00034257)*c*c + (-5.69551e-06)*c*c*c + (3.44998e-08)*c*c*c*c ;
1926   v2s=0.0205492 + (-0.00052015)*c + (-5.42051e-05)*c*c + (1.82754e-06)*c*c*c + (-1.58393e-08)*c*c*c*c ;
1927 }
1928 if(fInternalRunNumber==16){
1929   v2c=-0.0233639 + (0.00703612)*c + (-0.000403314)*c*c + (7.45693e-06)*c*c*c + (-4.4522e-08)*c*c*c*c ;
1930   v2s=-0.0712605 + (0.0124469)*c + (-0.000543644)*c*c + (8.22618e-06)*c*c*c + (-4.1073e-08)*c*c*c*c ;
1931 }
1932 if(fInternalRunNumber==17){
1933   v2c=0.0113215 + (-0.0019788)*c + (9.81834e-05)*c*c + (-1.30009e-06)*c*c*c + (5.03998e-09)*c*c*c*c ;
1934   v2s=0.0459851 + (-0.00170929)*c + (-5.96459e-05)*c*c + (2.90146e-06)*c*c*c + (-2.85325e-08)*c*c*c*c ;
1935 }
1936 if(fInternalRunNumber==18){
1937   v2c=-0.0924521 + (0.0162203)*c + (-0.000702717)*c*c + (1.0269e-05)*c*c*c + (-4.632e-08)*c*c*c*c ;
1938   v2s=-0.0456605 + (0.00706138)*c + (-0.000243957)*c*c + (2.7335e-06)*c*c*c + (-8.83296e-09)*c*c*c*c ;
1939 }
1940 if(fInternalRunNumber==19){
1941   v2c=0.00150974 + (0.0113356)*c + (-0.000923111)*c*c + (2.0452e-05)*c*c*c + (-1.36949e-07)*c*c*c*c ;
1942   v2s=0.0673845 + (-0.0187379)*c + (0.00112044)*c*c + (-2.19997e-05)*c*c*c + (1.33115e-07)*c*c*c*c ;
1943 }
1944 if(fInternalRunNumber==20){
1945   v2c=0.0223891 + (-0.000921021)*c + (-3.56276e-05)*c*c + (9.95622e-07)*c*c*c + (-5.54486e-09)*c*c*c*c ;
1946   v2s=0.0657305 + (-0.00354234)*c + (2.90428e-05)*c*c + (5.82994e-07)*c*c*c + (-7.53656e-09)*c*c*c*c ;
1947 }
1948 if(fInternalRunNumber==21){
1949   v2c=0.00969434 + (-0.00195313)*c + (3.56409e-05)*c*c + (2.15514e-07)*c*c*c + (-5.5161e-09)*c*c*c*c ;
1950   v2s=0.000236584 + (0.00458234)*c + (-0.000306408)*c*c + (6.04476e-06)*c*c*c + (-3.79354e-08)*c*c*c*c ;
1951 }
1952 if(fInternalRunNumber==22){
1953   v2c=0.0973536 + (-0.00907174)*c + (0.000239041)*c*c + (-2.64399e-06)*c*c*c + (1.12038e-08)*c*c*c*c ;
1954   v2s=0.0506578 + (-0.00256521)*c + (1.53982e-05)*c*c + (4.18193e-07)*c*c*c + (-4.99622e-09)*c*c*c*c ;
1955 }
1956 if(fInternalRunNumber==23){
1957   v2c=0.120714 + (-0.0190965)*c + (0.000887955)*c*c + (-1.57144e-05)*c*c*c + (9.22823e-08)*c*c*c*c ;
1958   v2s=-0.00178837 + (0.00377706)*c + (-0.000358619)*c*c + (9.54949e-06)*c*c*c + (-7.39141e-08)*c*c*c*c ;
1959 }
1960 if(fInternalRunNumber==24){
1961   v2c=0.0593053 + (-0.00687397)*c + (0.000200607)*c*c + (-2.34483e-06)*c*c*c + (9.60642e-09)*c*c*c*c ;
1962   v2s=0.0235174 + (0.00219788)*c + (-0.000225065)*c*c + (4.9594e-06)*c*c*c + (-3.2695e-08)*c*c*c*c ;
1963 }
1964 if(fInternalRunNumber==25){
1965   v2c=0.0551353 + (-0.00340874)*c + (2.26718e-05)*c*c + (4.1675e-07)*c*c*c + (-3.37457e-09)*c*c*c*c ;
1966   v2s=0.0110038 + (-0.000180543)*c + (3.28645e-05)*c*c + (-1.08119e-06)*c*c*c + (6.7157e-09)*c*c*c*c ;
1967 }
1968 if(fInternalRunNumber==26){
1969   v2c=-0.0138673 + (0.00475006)*c + (-0.00030584)*c*c + (6.0294e-06)*c*c*c + (-3.69906e-08)*c*c*c*c ;
1970   v2s=0.0395065 + (-0.00331262)*c + (9.40592e-05)*c*c + (-1.1541e-06)*c*c*c + (4.38329e-09)*c*c*c*c ;
1971 }
1972 if(fInternalRunNumber==27){
1973   v2c=0.0547504 + (-0.00669067)*c + (0.000269785)*c*c + (-4.75217e-06)*c*c*c + (2.96254e-08)*c*c*c*c ;
1974   v2s=0.101508 + (-0.0107211)*c + (0.000385539)*c*c + (-6.28713e-06)*c*c*c + (3.75932e-08)*c*c*c*c ;
1975 }
1976 if(fInternalRunNumber==28){
1977   v2c=-0.102173 + (0.0158262)*c + (-0.000706394)*c*c + (1.14263e-05)*c*c*c + (-6.14646e-08)*c*c*c*c ;
1978   v2s=0.02013 + (9.82197e-05)*c + (-1.71572e-05)*c*c + (-7.69646e-08)*c*c*c + (1.45113e-09)*c*c*c*c ;
1979 }
1980 if(fInternalRunNumber==29){
1981   v2c=0.100539 + (-0.0115012)*c + (0.000295367)*c*c + (-2.57489e-06)*c*c*c + (5.41526e-09)*c*c*c*c ;
1982   v2s=0.138027 + (-0.00342397)*c + (-0.000392569)*c*c + (1.34626e-05)*c*c*c + (-1.06683e-07)*c*c*c*c ;
1983 }
1984 if(fInternalRunNumber==30){
1985   v2c=0.0695326 + (-0.00752755)*c + (0.000237872)*c*c + (-3.31513e-06)*c*c*c + (1.69653e-08)*c*c*c*c ;
1986   v2s=0.0290133 + (-0.00073999)*c + (-3.35653e-05)*c*c + (7.61029e-07)*c*c*c + (-2.95792e-09)*c*c*c*c ;
1987 }
1988 if(fInternalRunNumber==31){
1989   v2c=-0.103171 + (0.0134058)*c + (-0.000631602)*c*c + (1.10298e-05)*c*c*c + (-6.41177e-08)*c*c*c*c ;
1990   v2s=-0.137648 + (0.0111259)*c + (-0.000400035)*c*c + (5.90535e-06)*c*c*c + (-3.24676e-08)*c*c*c*c ;
1991 }
1992 if(fInternalRunNumber==32){
1993   v2c=0.0522957 + (-0.00718492)*c + (0.000288586)*c*c + (-5.14534e-06)*c*c*c + (3.31303e-08)*c*c*c*c ;
1994   v2s=0.0731685 + (-0.00469449)*c + (9.57641e-05)*c*c + (-1.11891e-06)*c*c*c + (6.13984e-09)*c*c*c*c ;
1995 }
1996 if(fInternalRunNumber==33){
1997   v2c=0.0882828 + (-0.00916109)*c + (0.00027821)*c*c + (-3.61312e-06)*c*c*c + (1.76418e-08)*c*c*c*c ;
1998   v2s=0.123302 + (-0.0156303)*c + (0.000623468)*c*c + (-9.83474e-06)*c*c*c + (5.24687e-08)*c*c*c*c ;
1999 }
2000 if(fInternalRunNumber==34){
2001   v2c=0.0848636 + (-0.0142808)*c + (0.000723319)*c*c + (-1.32627e-05)*c*c*c + (7.77064e-08)*c*c*c*c ;
2002   v2s=-0.138483 + (0.0268553)*c + (-0.0013585)*c*c + (2.56034e-05)*c*c*c + (-1.63065e-07)*c*c*c*c ;
2003 }
2004 if(fInternalRunNumber==35){
2005   v2c=-0.0134575 + (0.00267433)*c + (-0.000168784)*c*c + (3.24012e-06)*c*c*c + (-1.92766e-08)*c*c*c*c ;
2006   v2s=0.0596795 + (-0.00350808)*c + (2.97282e-05)*c*c + (3.52115e-07)*c*c*c + (-4.17462e-09)*c*c*c*c ;
2007 }
2008 if(fInternalRunNumber==36){
2009   v2c=0.0784813 + (-0.00718389)*c + (0.000145303)*c*c + (-1.40121e-06)*c*c*c + (7.52775e-09)*c*c*c*c ;
2010   v2s=-0.0439181 + (0.0101363)*c + (-0.000494075)*c*c + (8.43356e-06)*c*c*c + (-4.68549e-08)*c*c*c*c ;
2011 }
2012 if(fInternalRunNumber==37){
2013   v2c=-0.0118549 + (0.00328751)*c + (-0.000207156)*c*c + (3.92678e-06)*c*c*c + (-2.30059e-08)*c*c*c*c ;
2014   v2s=0.0344878 + (-0.000363669)*c + (-7.62481e-05)*c*c + (1.67135e-06)*c*c*c + (-9.64004e-09)*c*c*c*c ;
2015 }
2016 if(fInternalRunNumber==38){
2017   v2c=-0.0254171 + (0.00297494)*c + (-0.000122361)*c*c + (1.78861e-06)*c*c*c + (-7.87609e-09)*c*c*c*c ;
2018   v2s=-0.0321399 + (0.00705636)*c + (-0.000325906)*c*c + (5.17669e-06)*c*c*c + (-2.86143e-08)*c*c*c*c ;
2019 }
2020 if(fInternalRunNumber==39){
2021   v2c=-0.0368721 + (0.00253133)*c + (-0.000136793)*c*c + (2.51874e-06)*c*c*c + (-1.46055e-08)*c*c*c*c ;
2022   v2s=-0.0757763 + (0.005728)*c + (-0.000217489)*c*c + (3.11936e-06)*c*c*c + (-1.56919e-08)*c*c*c*c ;
2023 }
2024 if(fInternalRunNumber==40){
2025   v2c=0.0166884 + (0.000509256)*c + (-0.000117911)*c*c + (2.65839e-06)*c*c*c + (-1.60875e-08)*c*c*c*c ;
2026   v2s=0.0336629 + (-0.00123023)*c + (-4.90123e-05)*c*c + (2.17862e-06)*c*c*c + (-1.98445e-08)*c*c*c*c ;
2027 }
2028 if(fInternalRunNumber==41){
2029   v2c=-0.00887514 + (0.00491344)*c + (-0.000325782)*c*c + (6.54782e-06)*c*c*c + (-4.14143e-08)*c*c*c*c ;
2030   v2s=0.0720999 + (-0.00650524)*c + (0.000214752)*c*c + (-3.19681e-06)*c*c*c + (1.64573e-08)*c*c*c*c ;
2031 }
2032 if(fInternalRunNumber==42){
2033   v2c=-0.0495903 + (0.00428229)*c + (-0.000181618)*c*c + (2.3913e-06)*c*c*c + (-8.55469e-09)*c*c*c*c ;
2034   v2s=0.0127982 + (-0.00690065)*c + (0.000342707)*c*c + (-6.40515e-06)*c*c*c + (3.83973e-08)*c*c*c*c ;
2035 }
2036 if(fInternalRunNumber==43){
2037   v2c=0.0196783 + (-0.000221524)*c + (-7.09556e-05)*c*c + (1.76517e-06)*c*c*c + (-1.19612e-08)*c*c*c*c ;
2038   v2s=0.0350359 + (-0.0020458)*c + (5.18727e-05)*c*c + (-6.58401e-07)*c*c*c + (1.10936e-09)*c*c*c*c ;
2039 }
2040 if(fInternalRunNumber==44){
2041   v2c=-0.0196041 + (0.00257682)*c + (-0.000174854)*c*c + (3.2168e-06)*c*c*c + (-1.85569e-08)*c*c*c*c ;
2042   v2s=-0.0870278 + (0.00634528)*c + (-0.000273588)*c*c + (4.77944e-06)*c*c*c + (-2.98422e-08)*c*c*c*c ;
2043 }
2044 if(fInternalRunNumber==45){
2045   v2c=-0.0105827 + (0.00271714)*c + (-0.000152357)*c*c + (2.71899e-06)*c*c*c + (-1.5435e-08)*c*c*c*c ;
2046   v2s=0.135579 + (-0.0135952)*c + (0.000424184)*c*c + (-5.12316e-06)*c*c*c + (1.84497e-08)*c*c*c*c ;
2047 }
2048 if(fInternalRunNumber==46){
2049   v2c=-0.017755 + (0.00411137)*c + (-0.000231724)*c*c + (4.17641e-06)*c*c*c + (-2.39164e-08)*c*c*c*c ;
2050   v2s=-0.0198628 + (0.0010032)*c + (-4.9708e-05)*c*c + (7.12455e-07)*c*c*c + (-4.22231e-09)*c*c*c*c ;
2051 }
2052 if(fInternalRunNumber==47){
2053   v2c=-0.0111387 + (0.00258591)*c + (-0.000194287)*c*c + (3.84062e-06)*c*c*c + (-2.37192e-08)*c*c*c*c ;
2054   v2s=-0.0478616 + (0.00365349)*c + (-0.000152037)*c*c + (2.24323e-06)*c*c*c + (-1.23308e-08)*c*c*c*c ;
2055 }
2056 if(fInternalRunNumber==48){
2057   v2c=0.070344 + (-0.0137332)*c + (0.000638536)*c*c + (-1.188e-05)*c*c*c + (7.4801e-08)*c*c*c*c ;
2058   v2s=-0.00208823 + (-0.00314157)*c + (9.9553e-05)*c*c + (-8.73891e-07)*c*c*c + (-1.72944e-09)*c*c*c*c ;
2059 }
2060 if(fInternalRunNumber==51){
2061   v2c=0.017789 + (-0.00598224)*c + (0.000328477)*c*c + (-6.80787e-06)*c*c*c + (4.46905e-08)*c*c*c*c ;
2062   v2s=-0.0894734 + (0.007351)*c + (-0.000291085)*c*c + (4.38654e-06)*c*c*c + (-2.32261e-08)*c*c*c*c ;
2063 }
2064 if(fInternalRunNumber==52){
2065   v2c=0.0555614 + (-0.00430576)*c + (0.00011519)*c*c + (-1.81584e-06)*c*c*c + (1.19214e-08)*c*c*c*c ;
2066   v2s=0.0901479 + (-0.0091505)*c + (0.000343757)*c*c + (-5.71204e-06)*c*c*c + (3.31217e-08)*c*c*c*c ;
2067 }
2068 if(fInternalRunNumber==53){
2069   v2c=0.0181761 + (-0.00186378)*c + (4.60531e-05)*c*c + (-7.59581e-07)*c*c*c + (5.52465e-09)*c*c*c*c ;
2070   v2s=0.0456713 + (-0.00526813)*c + (0.000232978)*c*c + (-3.98724e-06)*c*c*c + (2.10384e-08)*c*c*c*c ;
2071 }
2072 if(fInternalRunNumber==54){
2073   v2c=0.0456225 + (-0.00611173)*c + (0.000228414)*c*c + (-3.70527e-06)*c*c*c + (2.10774e-08)*c*c*c*c ;
2074   v2s=-0.00553874 + (0.00501799)*c + (-0.000331323)*c*c + (6.96978e-06)*c*c*c + (-4.72472e-08)*c*c*c*c ;
2075 }
2076 if(fInternalRunNumber==55){
2077   v2c=-0.0578879 + (0.00664986)*c + (-0.000335443)*c*c + (6.0086e-06)*c*c*c + (-3.54293e-08)*c*c*c*c ;
2078   v2s=-0.0675269 + (0.00407581)*c + (-0.000148063)*c*c + (1.922e-06)*c*c*c + (-8.8635e-09)*c*c*c*c ;
2079 }
2080 if(fInternalRunNumber==56){
2081   v2c=-0.000898903 + (0.00199222)*c + (-0.000207257)*c*c + (4.46148e-06)*c*c*c + (-2.71361e-08)*c*c*c*c ;
2082   v2s=-0.103794 + (0.0198667)*c + (-0.000879042)*c*c + (1.43177e-05)*c*c*c + (-8.09672e-08)*c*c*c*c ;
2083 }
2084 if(fInternalRunNumber==57){
2085   v2c=0.00953087 + (0.000947327)*c + (-3.45822e-05)*c*c + (-4.22294e-07)*c*c*c + (9.08727e-09)*c*c*c*c ;
2086   v2s=0.0662314 + (-0.0076136)*c + (0.000288482)*c*c + (-4.82931e-06)*c*c*c + (2.67751e-08)*c*c*c*c ;
2087 }
2088 if(fInternalRunNumber==58){
2089   v2c=0.00682728 + (-0.000293407)*c + (-2.88037e-05)*c*c + (6.63629e-07)*c*c*c + (-3.93121e-09)*c*c*c*c ;
2090   v2s=0.0085365 + (0.00138919)*c + (-8.15068e-05)*c*c + (1.26437e-06)*c*c*c + (-7.50322e-09)*c*c*c*c ;
2091 }
2092 if(fInternalRunNumber==59){
2093   v2c=-0.0185226 + (0.00377881)*c + (-0.000226383)*c*c + (4.22897e-06)*c*c*c + (-2.51383e-08)*c*c*c*c ;
2094   v2s=-0.222605 + (0.0100699)*c + (-0.00029256)*c*c + (3.3855e-06)*c*c*c + (-1.33648e-08)*c*c*c*c ;
2095 }
2096 if(fInternalRunNumber==60){
2097   v2c=-0.00111322 + (0.00291259)*c + (-0.000215769)*c*c + (4.13355e-06)*c*c*c + (-2.36883e-08)*c*c*c*c ;
2098   v2s=0.033686 + (-0.00124534)*c + (1.52616e-05)*c*c + (-2.99302e-07)*c*c*c + (1.05198e-09)*c*c*c*c ;
2099 }
2100 if(fInternalRunNumber==61){
2101   v2c=0.00041873 + (-0.000231092)*c + (1.19226e-05)*c*c + (-6.29604e-07)*c*c*c + (7.26582e-09)*c*c*c*c ;
2102   v2s=0.00435438 + (0.0025265)*c + (-0.000153028)*c*c + (2.71169e-06)*c*c*c + (-1.67127e-08)*c*c*c*c ;
2103 }
2104 if(fInternalRunNumber==62){
2105   v2c=-0.0573899 + (0.00972389)*c + (-0.000495772)*c*c + (9.0103e-06)*c*c*c + (-5.45756e-08)*c*c*c*c ;
2106   v2s=0.00702691 + (0.00698504)*c + (-0.000541929)*c*c + (1.19971e-05)*c*c*c + (-8.26942e-08)*c*c*c*c ;
2107 }
2108 if(fInternalRunNumber==63){
2109   v2c=0.047222 + (-0.00604037)*c + (0.0002092)*c*c + (-3.07144e-06)*c*c*c + (1.68579e-08)*c*c*c*c ;
2110   v2s=0.0568815 + (-0.00607571)*c + (0.000244243)*c*c + (-4.22237e-06)*c*c*c + (2.38928e-08)*c*c*c*c ;
2111 }
2112 if(fInternalRunNumber==64){
2113   v2c=0.0464686 + (-0.00463482)*c + (0.000146379)*c*c + (-2.23008e-06)*c*c*c + (1.29549e-08)*c*c*c*c ;
2114   v2s=-0.00571995 + (0.00273485)*c + (-0.000131945)*c*c + (2.08986e-06)*c*c*c + (-1.26717e-08)*c*c*c*c ;
2115 }
2116 if(fInternalRunNumber==65){
2117   v2c=0.0544794 + (-0.00616918)*c + (0.000173808)*c*c + (-2.01338e-06)*c*c*c + (8.78123e-09)*c*c*c*c ;
2118   v2s=0.0557002 + (-0.00701287)*c + (0.000340418)*c*c + (-6.12754e-06)*c*c*c + (3.3384e-08)*c*c*c*c ;
2119 }
2120 if(fInternalRunNumber==66){
2121   v2c=0.0111176 + (0.000262254)*c + (-7.20349e-05)*c*c + (1.38755e-06)*c*c*c + (-6.65016e-09)*c*c*c*c ;
2122   v2s=0.0644184 + (-0.00518515)*c + (0.000174227)*c*c + (-2.7347e-06)*c*c*c + (1.29547e-08)*c*c*c*c ;
2123 }
2124 if(fInternalRunNumber==67){
2125   v2c=-0.0481037 + (0.00823686)*c + (-0.000418053)*c*c + (7.40186e-06)*c*c*c + (-4.19609e-08)*c*c*c*c ;
2126   v2s=0.0560171 + (-0.002556)*c + (2.54529e-05)*c*c + (-5.22269e-09)*c*c*c + (-2.19983e-09)*c*c*c*c ;
2127 }
2128 if(fInternalRunNumber==68){
2129   v2c=-0.0139809 + (0.00274305)*c + (-0.00017744)*c*c + (3.20286e-06)*c*c*c + (-1.74137e-08)*c*c*c*c ;
2130   v2s=-0.218673 + (0.00764806)*c + (-0.000222739)*c*c + (3.14656e-06)*c*c*c + (-1.73548e-08)*c*c*c*c ;
2131 }
2132 if(fInternalRunNumber==69){
2133   v2c=0.0736938 + (-0.0122004)*c + (0.000498227)*c*c + (-8.08779e-06)*c*c*c + (4.52749e-08)*c*c*c*c ;
2134   v2s=-0.236254 + (0.00979816)*c + (-0.000280698)*c*c + (3.4103e-06)*c*c*c + (-1.47299e-08)*c*c*c*c ;
2135 }
2136 if(fInternalRunNumber==70){
2137   v2c=-0.0366545 + (0.00399604)*c + (-0.000196876)*c*c + (3.49686e-06)*c*c*c + (-2.01142e-08)*c*c*c*c ;
2138   v2s=0.0220911 + (-0.000529459)*c + (9.10408e-06)*c*c + (-3.57334e-07)*c*c*c + (1.87909e-09)*c*c*c*c ;
2139 }
2140 if(fInternalRunNumber==71){
2141   v2c=-0.039698 + (0.00654615)*c + (-0.000323865)*c*c + (5.50998e-06)*c*c*c + (-3.00866e-08)*c*c*c*c ;
2142   v2s=0.0650078 + (-0.00502636)*c + (0.000128459)*c*c + (-1.3877e-06)*c*c*c + (3.21917e-09)*c*c*c*c ;
2143 }
2144 if(fInternalRunNumber==72){
2145   v2c=0.0594227 + (-0.00872513)*c + (0.000368684)*c*c + (-6.34986e-06)*c*c*c + (3.79523e-08)*c*c*c*c ;
2146   v2s=0.050348 + (-0.00512087)*c + (0.000228643)*c*c + (-4.40797e-06)*c*c*c + (2.62482e-08)*c*c*c*c ;
2147 }
2148 if(fInternalRunNumber==73){
2149   v2c=0.131281 + (-0.0162697)*c + (0.000584063)*c*c + (-8.64349e-06)*c*c*c + (4.61009e-08)*c*c*c*c ;
2150   v2s=0.0271721 + (-0.000487039)*c + (1.60624e-05)*c*c + (-5.38029e-07)*c*c*c + (2.08589e-09)*c*c*c*c ;
2151 }
2152 if(fInternalRunNumber==74){
2153   v2c=-0.0128651 + (0.00946426)*c + (-0.00055547)*c*c + (9.37245e-06)*c*c*c + (-4.78224e-08)*c*c*c*c ;
2154   v2s=0.0903476 + (-0.0164972)*c + (0.000927481)*c*c + (-1.88459e-05)*c*c*c + (1.20249e-07)*c*c*c*c ;
2155 }
2156 if(fInternalRunNumber==75){
2157   v2c=0.0204147 + (-0.00139139)*c + (-2.47732e-05)*c*c + (1.1214e-06)*c*c*c + (-7.22533e-09)*c*c*c*c ;
2158   v2s=0.0214247 + (-0.000147929)*c + (-7.58241e-06)*c*c + (-7.05565e-08)*c*c*c + (-4.86528e-10)*c*c*c*c ;
2159 }
2160 if(fInternalRunNumber==76){
2161   v2c=0.0469095 + (-0.00541793)*c + (0.000133238)*c*c + (-1.14334e-06)*c*c*c + (3.4826e-09)*c*c*c*c ;
2162   v2s=0.0304763 + (-0.00213957)*c + (8.27588e-05)*c*c + (-1.34495e-06)*c*c*c + (4.81134e-09)*c*c*c*c ;
2163 }
2164 if(fInternalRunNumber==77){
2165   v2c=0.113822 + (-0.0154286)*c + (0.000574533)*c*c + (-8.67874e-06)*c*c*c + (4.76727e-08)*c*c*c*c ;
2166   v2s=0.010623 + (0.00513467)*c + (-0.000372733)*c*c + (7.40692e-06)*c*c*c + (-4.64093e-08)*c*c*c*c ;
2167 }
2168 if(fInternalRunNumber==78){
2169   v2c=0.020658 + (-0.00274725)*c + (5.8872e-05)*c*c + (-3.43999e-07)*c*c*c + (-1.04273e-10)*c*c*c*c ;
2170   v2s=0.00219379 + (0.000972033)*c + (-1.41816e-05)*c*c + (-4.41149e-07)*c*c*c + (3.81037e-09)*c*c*c*c ;
2171 }
2172 if(fInternalRunNumber==79){
2173   v2c=0.227324 + (-0.0264915)*c + (0.000915298)*c*c + (-1.3343e-05)*c*c*c + (7.1044e-08)*c*c*c*c ;
2174   v2s=0.0474072 + (-0.00939515)*c + (0.000730108)*c*c + (-1.80811e-05)*c*c*c + (1.3003e-07)*c*c*c*c ;
2175 }
2176 if(fInternalRunNumber==81){
2177   v2c=-0.0368028 + (0.00368927)*c + (-0.000228321)*c*c + (4.48174e-06)*c*c*c + (-2.75298e-08)*c*c*c*c ;
2178   v2s=-0.238807 + (0.0108667)*c + (-0.000350649)*c*c + (4.95568e-06)*c*c*c + (-2.53022e-08)*c*c*c*c ;
2179 }
2180 if(fInternalRunNumber==82){
2181   v2c=-0.0214072 + (0.00119258)*c + (-7.80215e-05)*c*c + (1.37682e-06)*c*c*c + (-7.31103e-09)*c*c*c*c ;
2182   v2s=-0.233825 + (0.00735109)*c + (-0.000142737)*c*c + (9.15013e-07)*c*c*c + (-1.98463e-10)*c*c*c*c ;
2183 }
2184 if(fInternalRunNumber==83){
2185   v2c=-0.0716708 + (0.00855289)*c + (-0.00037249)*c*c + (6.41638e-06)*c*c*c + (-3.66318e-08)*c*c*c*c ;
2186   v2s=-0.329702 + (0.033231)*c + (-0.00131012)*c*c + (2.1101e-05)*c*c*c + (-1.18905e-07)*c*c*c*c ;
2187 }
2188 if(fInternalRunNumber==84){
2189   v2c=-0.0259022 + (0.00711773)*c + (-0.000434778)*c*c + (8.94774e-06)*c*c*c + (-5.84408e-08)*c*c*c*c ;
2190   v2s=-0.0936287 + (-0.00455282)*c + (0.000390982)*c*c + (-7.95224e-06)*c*c*c + (4.95612e-08)*c*c*c*c ;
2191 }
2192 if(fInternalRunNumber==85){
2193   v2c=0.0263223 + (-0.00475593)*c + (9.63236e-05)*c*c + (-4.95311e-07)*c*c*c + (-5.47428e-10)*c*c*c*c ;
2194   v2s=-0.220593 + (0.00823892)*c + (-0.000248477)*c*c + (3.28082e-06)*c*c*c + (-1.45098e-08)*c*c*c*c ;
2195 }
2196 if(fInternalRunNumber==86){
2197   v2c=-0.0157962 + (0.00408922)*c + (-0.000282087)*c*c + (5.81961e-06)*c*c*c + (-3.6481e-08)*c*c*c*c ;
2198   v2s=0.0340979 + (-0.00139471)*c + (6.48561e-06)*c*c + (3.02184e-08)*c*c*c + (-1.05879e-09)*c*c*c*c ;
2199 }
2200 if(fInternalRunNumber==87){
2201   v2c=0.0577366 + (-0.00471394)*c + (0.000107329)*c*c + (-7.28927e-07)*c*c*c + (-1.51199e-09)*c*c*c*c ;
2202   v2s=0.0350171 + (-0.00402687)*c + (0.00016046)*c*c + (-2.39625e-06)*c*c*c + (9.98146e-09)*c*c*c*c ;
2203 }
2204 if(fInternalRunNumber==88){
2205   v2c=-0.00547971 + (0.00161905)*c + (-0.00013831)*c*c + (2.94662e-06)*c*c*c + (-1.92504e-08)*c*c*c*c ;
2206   v2s=-0.0722251 + (0.0127634)*c + (-0.000612284)*c*c + (1.07553e-05)*c*c*c + (-6.34054e-08)*c*c*c*c ;
2207 }
2208 if(fInternalRunNumber==89){
2209   v2c=0.165574 + (-0.0241913)*c + (0.00101221)*c*c + (-1.62279e-05)*c*c*c + (8.77252e-08)*c*c*c*c ;
2210   v2s=0.0692954 + (-0.0147219)*c + (0.00068782)*c*c + (-1.11091e-05)*c*c*c + (5.37295e-08)*c*c*c*c ;
2211 }
2212 if(fInternalRunNumber==91){
2213   v2c=-0.0334666 + (0.00583866)*c + (-0.000329883)*c*c + (6.08932e-06)*c*c*c + (-3.59853e-08)*c*c*c*c ;
2214   v2s=-0.229775 + (0.010061)*c + (-0.00028144)*c*c + (3.3231e-06)*c*c*c + (-1.40799e-08)*c*c*c*c ;
2215 }
2216 if(fInternalRunNumber==95){
2217   v2c=-0.0993203 + (0.0150547)*c + (-0.000713531)*c*c + (1.22733e-05)*c*c*c + (-6.98669e-08)*c*c*c*c ;
2218   v2s=-0.217746 + (0.0068375)*c + (-0.000171962)*c*c + (2.05506e-06)*c*c*c + (-8.78703e-09)*c*c*c*c ;
2219 }
2220 if(fInternalRunNumber==96){
2221   v2c=-0.0398801 + (0.00406028)*c + (-0.000214066)*c*c + (3.69056e-06)*c*c*c + (-1.99339e-08)*c*c*c*c ;
2222   v2s=-0.240195 + (0.00782806)*c + (-0.000131365)*c*c + (3.18193e-07)*c*c*c + (4.88009e-09)*c*c*c*c ;
2223 }
2224 if(fInternalRunNumber==97){
2225   v2c=-0.0325538 + (0.00250986)*c + (-0.000136762)*c*c + (2.35572e-06)*c*c*c + (-1.23413e-08)*c*c*c*c ;
2226   v2s=-0.240415 + (0.010435)*c + (-0.000291848)*c*c + (3.35357e-06)*c*c*c + (-1.35833e-08)*c*c*c*c ;
2227 }
2228 if(fInternalRunNumber==98){
2229   v2c=-0.00692338 + (0.00379335)*c + (-0.000332833)*c*c + (7.35052e-06)*c*c*c + (-4.89956e-08)*c*c*c*c ;
2230   v2s=-0.263042 + (0.0134291)*c + (-0.000465334)*c*c + (7.09595e-06)*c*c*c + (-3.90332e-08)*c*c*c*c ;
2231 }
2232 if(fInternalRunNumber==101){
2233   v2c=-0.049409 + (0.0073433)*c + (-0.000376933)*c*c + (7.39945e-06)*c*c*c + (-4.967e-08)*c*c*c*c ;
2234   v2s=0.00564511 + (-0.00101958)*c + (0.000121077)*c*c + (-3.44889e-06)*c*c*c + (2.42733e-08)*c*c*c*c ;
2235 }
2236 if(fInternalRunNumber==102){
2237   v2c=-0.0393185 + (0.00610331)*c + (-0.000286099)*c*c + (4.61693e-06)*c*c*c + (-2.44104e-08)*c*c*c*c ;
2238   v2s=-0.00898987 + (0.00463766)*c + (-0.000269122)*c*c + (5.46456e-06)*c*c*c + (-3.85677e-08)*c*c*c*c ;
2239 }
2240 if(fInternalRunNumber==103){
2241   v2c=-0.0171815 + (0.00271616)*c + (-9.49063e-05)*c*c + (5.68888e-07)*c*c*c + (2.77777e-09)*c*c*c*c ;
2242   v2s=-0.0348362 + (0.00661647)*c + (-0.000196625)*c*c + (1.60914e-06)*c*c*c + (-4.34769e-09)*c*c*c*c ;
2243 }
2244 if(fInternalRunNumber==104){
2245   v2c=0.0453893 + (-0.00220981)*c + (-6.34642e-05)*c*c + (2.5268e-06)*c*c*c + (-1.81319e-08)*c*c*c*c ;
2246   v2s=0.0131555 + (0.00115036)*c + (-5.70975e-05)*c*c + (5.20224e-07)*c*c*c + (-3.11332e-09)*c*c*c*c ;
2247 }
2248 if(fInternalRunNumber==105){
2249   v2c=0.0052505 + (0.000251094)*c + (-3.95143e-05)*c*c + (5.6354e-07)*c*c*c + (-1.86603e-09)*c*c*c*c ;
2250   v2s=0.079882 + (-0.00841369)*c + (0.000342924)*c*c + (-5.64554e-06)*c*c*c + (2.89511e-08)*c*c*c*c ;
2251 }
2252 if(fInternalRunNumber==106){
2253   v2c=0.0561693 + (-0.00693389)*c + (0.000243013)*c*c + (-3.7816e-06)*c*c*c + (2.21666e-08)*c*c*c*c ;
2254   v2s=0.115847 + (-0.0145063)*c + (0.000597908)*c*c + (-9.37418e-06)*c*c*c + (4.57406e-08)*c*c*c*c ;
2255 }
2256 if(fInternalRunNumber==107){
2257   v2c=-0.160892 + (0.0259075)*c + (-0.00120084)*c*c + (1.99814e-05)*c*c*c + (-1.09708e-07)*c*c*c*c ;
2258   v2s=0.0922001 + (-0.00861465)*c + (0.000196846)*c*c + (-5.34939e-07)*c*c*c + (-1.20442e-08)*c*c*c*c ;
2259 }
2260 if(fInternalRunNumber==108){
2261   v2c=0.0325605 + (-0.00795036)*c + (0.000436987)*c*c + (-9.37794e-06)*c*c*c + (6.59912e-08)*c*c*c*c ;
2262   v2s=0.0894763 + (-0.00694006)*c + (0.000214928)*c*c + (-3.05745e-06)*c*c*c + (1.33679e-08)*c*c*c*c ;
2263 }
2264 if(fInternalRunNumber==109){
2265   v2c=-8.30505e-05 + (0.00304264)*c + (-0.000228365)*c*c + (4.43511e-06)*c*c*c + (-2.5788e-08)*c*c*c*c ;
2266   v2s=0.0276356 + (-0.000486726)*c + (-1.12152e-05)*c*c + (1.56258e-07)*c*c*c + (-2.17826e-09)*c*c*c*c ;
2267 }
2268 if(fInternalRunNumber==112){
2269   v2c=0.0273919 + (-0.00243837)*c + (2.41667e-05)*c*c + (3.06367e-07)*c*c*c + (-3.9028e-09)*c*c*c*c ;
2270   v2s=0.0119428 + (0.00161447)*c + (-0.00010526)*c*c + (1.76239e-06)*c*c*c + (-1.12844e-08)*c*c*c*c ;
2271 }
2272 if(fInternalRunNumber==114){
2273   v2c=0.114287 + (-0.016125)*c + (0.000673002)*c*c + (-1.11635e-05)*c*c*c + (6.39533e-08)*c*c*c*c ;
2274   v2s=0.0353709 + (-0.00169691)*c + (-3.93183e-05)*c*c + (2.08815e-06)*c*c*c + (-2.03159e-08)*c*c*c*c ;
2275 }
2276 if(fInternalRunNumber==115){
2277   v2c=-0.0109276 + (0.00378322)*c + (-0.000252651)*c*c + (4.599e-06)*c*c*c + (-2.36568e-08)*c*c*c*c ;
2278   v2s=0.0185448 + (0.00102852)*c + (-6.04966e-05)*c*c + (5.89933e-07)*c*c*c + (-2.5921e-09)*c*c*c*c ;
2279 }
2280 if(fInternalRunNumber==116){
2281   v2c=-0.117563 + (0.00405449)*c + (-0.000112806)*c*c + (4.01903e-06)*c*c*c + (-4.25358e-08)*c*c*c*c ;
2282   v2s=0.0280797 + (-0.00528592)*c + (0.000375784)*c*c + (-9.96825e-06)*c*c*c + (7.72147e-08)*c*c*c*c ;
2283 }
2284 if(fInternalRunNumber==117){
2285   v2c=0.0301869 + (-0.000751512)*c + (-9.43748e-05)*c*c + (2.60697e-06)*c*c*c + (-1.73948e-08)*c*c*c*c ;
2286   v2s=0.0594442 + (-0.00499442)*c + (0.000170617)*c*c + (-2.80903e-06)*c*c*c + (1.51075e-08)*c*c*c*c ;
2287 }
2288 if(fInternalRunNumber==118){
2289   v2c=0.0118942 + (-0.000669595)*c + (-3.6634e-05)*c*c + (1.40794e-06)*c*c*c + (-1.15253e-08)*c*c*c*c ;
2290   v2s=-0.0174395 + (0.00460355)*c + (-0.00023192)*c*c + (3.83941e-06)*c*c*c + (-2.30761e-08)*c*c*c*c ;
2291 }
2292 if(fInternalRunNumber==119){
2293   v2c=0.018229 + (-0.00144698)*c + (-1.07295e-06)*c*c + (6.22528e-07)*c*c*c + (-4.90504e-09)*c*c*c*c ;
2294   v2s=0.0122761 + (0.000985351)*c + (-6.02936e-05)*c*c + (9.41867e-07)*c*c*c + (-6.59756e-09)*c*c*c*c ;
2295 }
2296 if(fInternalRunNumber==120){
2297   v2c=-0.0968065 + (0.0131861)*c + (-0.000585224)*c*c + (9.82175e-06)*c*c*c + (-5.51313e-08)*c*c*c*c ;
2298   v2s=0.0719334 + (-0.00752912)*c + (0.000362464)*c*c + (-6.65114e-06)*c*c*c + (3.70295e-08)*c*c*c*c ;
2299 }
2300 if(fInternalRunNumber==121){
2301   v2c=0.0296329 + (-0.00448004)*c + (0.00010239)*c*c + (-2.40647e-07)*c*c*c + (-5.4555e-09)*c*c*c*c ;
2302   v2s=0.0774925 + (-0.00635491)*c + (0.0001939)*c*c + (-3.04755e-06)*c*c*c + (1.66055e-08)*c*c*c*c ;
2303 }
2304 if(fInternalRunNumber==122){
2305   v2c=0.00956939 + (0.000460793)*c + (-0.000105912)*c*c + (2.38286e-06)*c*c*c + (-1.446e-08)*c*c*c*c ;
2306   v2s=0.0152677 + (0.00345263)*c + (-0.000262399)*c*c + (5.5895e-06)*c*c*c + (-3.84974e-08)*c*c*c*c ;
2307 }
2308 if(fInternalRunNumber==123){
2309   v2c=0.0495618 + (-0.00059293)*c + (-0.00025227)*c*c + (6.61088e-06)*c*c*c + (-4.24063e-08)*c*c*c*c ;
2310   v2s=0.0862227 + (-0.0190161)*c + (0.0010647)*c*c + (-1.94957e-05)*c*c*c + (1.08665e-07)*c*c*c*c ;
2311 }
2312 if(fInternalRunNumber==124){
2313   v2c=-0.0408078 + (0.00565108)*c + (-0.00027528)*c*c + (5.2165e-06)*c*c*c + (-3.3244e-08)*c*c*c*c ;
2314   v2s=0.00169862 + (0.00371982)*c + (-0.000169757)*c*c + (2.07347e-06)*c*c*c + (-8.30142e-09)*c*c*c*c ;
2315 }
2316 if(fInternalRunNumber==125){
2317   v2c=0.0285065 + (-0.00348583)*c + (0.000106189)*c*c + (-1.56726e-06)*c*c*c + (9.40552e-09)*c*c*c*c ;
2318   v2s=0.0363645 + (-0.00252206)*c + (0.000122553)*c*c + (-2.74066e-06)*c*c*c + (1.66856e-08)*c*c*c*c ;
2319 }
2320 if(fInternalRunNumber==126){
2321   v2c=-0.00864424 + (0.00157555)*c + (-0.000127367)*c*c + (2.47447e-06)*c*c*c + (-1.42302e-08)*c*c*c*c ;
2322   v2s=-0.213003 + (0.00643156)*c + (-0.000121219)*c*c + (6.06838e-07)*c*c*c + (1.48381e-09)*c*c*c*c ;
2323 }
2324 if(fInternalRunNumber==127){
2325   v2c=-0.0336884 + (0.00299236)*c + (-0.000147557)*c*c + (2.47908e-06)*c*c*c + (-1.33495e-08)*c*c*c*c ;
2326   v2s=-0.271925 + (0.0158755)*c + (-0.000583802)*c*c + (9.14376e-06)*c*c*c + (-5.05302e-08)*c*c*c*c ;
2327 }
2328 if(fInternalRunNumber==128){
2329   v2c=-0.0639711 + (0.00822881)*c + (-0.000428416)*c*c + (7.98075e-06)*c*c*c + (-4.81995e-08)*c*c*c*c ;
2330   v2s=-0.22634 + (0.00908632)*c + (-0.00029)*c*c + (4.14095e-06)*c*c*c + (-2.15823e-08)*c*c*c*c ;
2331 }
2332 if(fInternalRunNumber==129){
2333   v2c=0.0652753 + (-0.00529737)*c + (0.000133883)*c*c + (-1.63339e-06)*c*c*c + (7.88294e-09)*c*c*c*c ;
2334   v2s=0.00357003 + (0.00503303)*c + (-0.000316531)*c*c + (6.45793e-06)*c*c*c + (-4.44467e-08)*c*c*c*c ;
2335 }
2336 if(fInternalRunNumber==131){
2337   v2c=0.117255 + (-0.0149009)*c + (0.000560178)*c*c + (-8.75754e-06)*c*c*c + (4.85024e-08)*c*c*c*c ;
2338   v2s=-0.00180786 + (0.00224258)*c + (-9.84056e-05)*c*c + (1.39525e-06)*c*c*c + (-7.63512e-09)*c*c*c*c ;
2339 }
2340 if(fInternalRunNumber==132){
2341   v2c=0.0100287 + (-0.00051338)*c + (-2.7702e-05)*c*c + (5.0205e-07)*c*c*c + (-1.12682e-09)*c*c*c*c ;
2342   v2s=0.00556895 + (0.00211082)*c + (-0.000125448)*c*c + (2.30409e-06)*c*c*c + (-1.63703e-08)*c*c*c*c ;
2343 }
2344 if(fInternalRunNumber==133){
2345   v2c=-0.00755913 + (0.00327793)*c + (-0.00024545)*c*c + (5.03417e-06)*c*c*c + (-3.11617e-08)*c*c*c*c ;
2346   v2s=0.0398477 + (-0.000728372)*c + (-2.63935e-05)*c*c + (7.85276e-07)*c*c*c + (-7.72703e-09)*c*c*c*c ;
2347 }
2348 if(fInternalRunNumber==134){
2349   v2c=0.00344056 + (0.000171769)*c + (-8.00887e-05)*c*c + (1.76696e-06)*c*c*c + (-1.08652e-08)*c*c*c*c ;
2350   v2s=-0.200474 + (0.00736532)*c + (-0.000221809)*c*c + (2.91763e-06)*c*c*c + (-1.40827e-08)*c*c*c*c ;
2351 }
2352 if(fInternalRunNumber==135){
2353   v2c=-0.0251554 + (0.00227996)*c + (-0.000133757)*c*c + (2.50604e-06)*c*c*c + (-1.48528e-08)*c*c*c*c ;
2354   v2s=-0.229223 + (0.0091122)*c + (-0.000245035)*c*c + (2.64826e-06)*c*c*c + (-9.57268e-09)*c*c*c*c ;
2355 }
2356 if(fInternalRunNumber==136){
2357   v2c=-0.0308547 + (0.00209237)*c + (-0.000183034)*c*c + (4.54164e-06)*c*c*c + (-3.26557e-08)*c*c*c*c ;
2358   v2s=-0.241702 + (0.00988355)*c + (-0.000263014)*c*c + (2.73376e-06)*c*c*c + (-9.79971e-09)*c*c*c*c ;
2359 }
2360
2361
2362
2363
2364
2365 /*
2366
2367   if(fRunNumber>=1 && fRunNumber<=45){ // 137161 -137848 //period 1
2368 //    v2c = 1.133676e-02-1.759454e-04*x-6.217033e-05*x*x+1.532507e-06*x*x*x-9.665265e-09*x*x*x*x;
2369 //    v2s = 2.173690e-02-8.410340e-04*x-1.391841e-05*x*x+5.003301e-07*x*x*x-4.175974e-09*x*x*x*x;
2370     v2c=2.295070e-02-1.050171e-03*x-3.432843e-05*x*x+1.156549e-06*x*x*x-7.781377e-09*x*x*x*x ;
2371     v2s=3.994886e-02-1.961957e-03*x+1.641958e-05*x*x+1.234025e-07*x*x*x-2.175698e-09*x*x*x*x ;
2372   }
2373   if(fRunNumber>=46 && fRunNumber<=58){ //138125-138275//period 2
2374     v2c =-1.394612e-03+1.333726e-03*x-1.115323e-04*x*x+2.103762e-06*x*x*x-1.186760e-08*x*x*x*x ;
2375     v2s =-7.788201e-03+7.178984e-04*x-4.553597e-05*x*x+6.941817e-07*x*x*x-4.465249e-09*x*x*x*x ;
2376   }
2377   if(fRunNumber>=59 && fRunNumber<=78){ //138359-138730//period 3
2378     v2c = 1.448809e-02+-7.894029e-04*x-2.693871e-05*x*x+7.944017e-07*x*x*x-4.272135e-09*x*x*x*x ;
2379     v2s = 3.442216e-03+-5.574560e-05*x-1.758343e-05*x*x+2.386547e-07*x*x*x-2.404210e-09*x*x*x*x ;
2380   }
2381   if(fRunNumber>=79){//period 4
2382     v2c =-1.553536e-03+7.761106e-04*x-9.529723e-05*x*x+2.026750e-06*x*x*x-1.217009e-08*x*x*x*x ;
2383     v2s =-7.787317e-02+3.843455e-03*x-1.221144e-04*x*x+1.478448e-06*x*x*x-7.442702e-09*x*x*x*x ;
2384   }
2385 */
2386
2387   return phi+v2c*TMath::Sin(2.*phi)-v2s*TMath::Cos(2.*phi) ;
2388
2389 }
2390 //____________________________________________________________________________
2391 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0C(Double_t phi, Double_t c){
2392   //LHC10h
2393
2394   Double_t v2c =0.;
2395   Double_t v2s =0;
2396
2397 if(fInternalRunNumber==1){
2398   v2c=-0.038901 + (0.00539521)*c + (-0.000213271)*c*c + (3.01147e-06)*c*c*c + (-1.35977e-08)*c*c*c*c ;
2399   v2s=0.0693294 + (-0.00851319)*c + (0.000337408)*c*c + (-5.33097e-06)*c*c*c + (2.90883e-08)*c*c*c*c ;
2400 }
2401 if(fInternalRunNumber==2){
2402   v2c=-0.0943695 + (0.0137829)*c + (-0.000569334)*c*c + (8.95522e-06)*c*c*c + (-4.78396e-08)*c*c*c*c ;
2403   v2s=0.00971525 + (-0.00181026)*c + (6.15798e-05)*c*c + (-8.51046e-07)*c*c*c + (4.57403e-09)*c*c*c*c ;
2404 }
2405 if(fInternalRunNumber==4){
2406   v2c=-0.00585475 + (0.000307793)*c + (1.89109e-06)*c*c + (-3.14804e-07)*c*c*c + (3.54341e-09)*c*c*c*c ;
2407   v2s=0.00374568 + (-0.00123435)*c + (8.01941e-05)*c*c + (-1.70807e-06)*c*c*c + (1.14576e-08)*c*c*c*c ;
2408 }
2409 if(fInternalRunNumber==5){
2410   v2c=-0.0301267 + (0.00206698)*c + (-4.82282e-05)*c*c + (4.39355e-07)*c*c*c + (-1.3747e-09)*c*c*c*c ;
2411   v2s=-0.0441947 + (0.00501774)*c + (-0.000181736)*c*c + (2.69235e-06)*c*c*c + (-1.37411e-08)*c*c*c*c ;
2412 }
2413 if(fInternalRunNumber==6){
2414   v2c=0.190652 + (-0.0254912)*c + (0.00105707)*c*c + (-1.63872e-05)*c*c*c + (8.51113e-08)*c*c*c*c ;
2415   v2s=0.0803124 + (-0.0119724)*c + (0.000575997)*c*c + (-1.07892e-05)*c*c*c + (6.67214e-08)*c*c*c*c ;
2416 }
2417 if(fInternalRunNumber==7){
2418   v2c=0.00875044 + (-0.00343147)*c + (0.0002198)*c*c + (-4.76292e-06)*c*c*c + (3.24844e-08)*c*c*c*c ;
2419   v2s=-0.081288 + (0.0146429)*c + (-0.000728327)*c*c + (1.32462e-05)*c*c*c + (-7.95362e-08)*c*c*c*c ;
2420 }
2421 if(fInternalRunNumber==8){
2422   v2c=0.0381367 + (-0.0082048)*c + (0.000410914)*c*c + (-7.28851e-06)*c*c*c + (4.22008e-08)*c*c*c*c ;
2423   v2s=-0.023338 + (0.003381)*c + (-0.000167791)*c*c + (3.13998e-06)*c*c*c + (-1.97233e-08)*c*c*c*c ;
2424 }
2425 if(fInternalRunNumber==9){
2426   v2c=-0.00734261 + (0.00121272)*c + (-5.31181e-05)*c*c + (7.88974e-07)*c*c*c + (-3.86356e-09)*c*c*c*c ;
2427   v2s=-0.0228333 + (0.00266052)*c + (-0.000112824)*c*c + (1.89515e-06)*c*c*c + (-1.07446e-08)*c*c*c*c ;
2428 }
2429 if(fInternalRunNumber==10){
2430   v2c=-0.0456385 + (0.0102956)*c + (-0.000529689)*c*c + (9.62809e-06)*c*c*c + (-5.7301e-08)*c*c*c*c ;
2431   v2s=0.0142197 + (-0.0031677)*c + (0.000190389)*c*c + (-4.10911e-06)*c*c*c + (2.81791e-08)*c*c*c*c ;
2432 }
2433 if(fInternalRunNumber==11){
2434   v2c=-0.00999249 + (0.00123939)*c + (-4.23895e-05)*c*c + (6.55415e-07)*c*c*c + (-3.98256e-09)*c*c*c*c ;
2435   v2s=-0.00809817 + (0.00199727)*c + (-0.000116718)*c*c + (2.33449e-06)*c*c*c + (-1.46407e-08)*c*c*c*c ;
2436 }
2437 if(fInternalRunNumber==12){
2438   v2c=0.0183566 + (-0.0020805)*c + (6.79726e-05)*c*c + (-7.96708e-07)*c*c*c + (3.1214e-09)*c*c*c*c ;
2439   v2s=-0.00714715 + (0.000721418)*c + (-3.50355e-05)*c*c + (9.09996e-07)*c*c*c + (-7.92198e-09)*c*c*c*c ;
2440 }
2441 if(fInternalRunNumber==13){
2442   v2c=0.0056915 + (0.000167211)*c + (-4.68018e-05)*c*c + (1.21891e-06)*c*c*c + (-8.70124e-09)*c*c*c*c ;
2443   v2s=0.00986223 + (-0.00247957)*c + (0.000133539)*c*c + (-2.51362e-06)*c*c*c + (1.5397e-08)*c*c*c*c ;
2444 }
2445 if(fInternalRunNumber==14){
2446   v2c=0.0479623 + (-0.00948196)*c + (0.000551795)*c*c + (-1.05692e-05)*c*c*c + (6.38754e-08)*c*c*c*c ;
2447   v2s=-0.00438456 + (-0.00152976)*c + (2.00645e-05)*c*c + (9.91017e-07)*c*c*c + (-1.35776e-08)*c*c*c*c ;
2448 }
2449 if(fInternalRunNumber==15){
2450   v2c=-0.00384088 + (0.00121719)*c + (-6.88665e-05)*c*c + (1.26138e-06)*c*c*c + (-7.13987e-09)*c*c*c*c ;
2451   v2s=0.0273017 + (-0.00314818)*c + (7.85031e-05)*c*c + (-9.07336e-08)*c*c*c + (-6.42421e-09)*c*c*c*c ;
2452 }
2453 if(fInternalRunNumber==16){
2454   v2c=-0.0313175 + (0.00368025)*c + (-0.000126895)*c*c + (1.7443e-06)*c*c*c + (-9.12442e-09)*c*c*c*c ;
2455   v2s=-0.0812178 + (0.0117402)*c + (-0.000506287)*c*c + (8.36766e-06)*c*c*c + (-4.69866e-08)*c*c*c*c ;
2456 }
2457 if(fInternalRunNumber==17){
2458   v2c=0.062893 + (-0.0128204)*c + (0.000591931)*c*c + (-1.00835e-05)*c*c*c + (5.72782e-08)*c*c*c*c ;
2459   v2s=-0.143105 + (0.0183921)*c + (-0.000787289)*c*c + (1.33237e-05)*c*c*c + (-7.71707e-08)*c*c*c*c ;
2460 }
2461 if(fInternalRunNumber==18){
2462   v2c=-0.0029353 + (-0.0016097)*c + (0.000199618)*c*c + (-5.42008e-06)*c*c*c + (4.08897e-08)*c*c*c*c ;
2463   v2s=0.0482015 + (-0.00115233)*c + (3.96019e-05)*c*c + (-7.7193e-07)*c*c*c + (3.90004e-09)*c*c*c*c ;
2464 }
2465 if(fInternalRunNumber==19){
2466   v2c=-0.0397728 + (0.0129154)*c + (-0.000703903)*c*c + (1.30972e-05)*c*c*c + (-7.93547e-08)*c*c*c*c ;
2467   v2s=0.181945 + (-0.0289807)*c + (0.00133685)*c*c + (-2.30866e-05)*c*c*c + (1.32746e-07)*c*c*c*c ;
2468 }
2469 if(fInternalRunNumber==20){
2470   v2c=-0.00927858 + (0.0021033)*c + (-0.000113971)*c*c + (2.17508e-06)*c*c*c + (-1.36812e-08)*c*c*c*c ;
2471   v2s=0.0201645 + (-0.00331311)*c + (0.000164916)*c*c + (-3.08376e-06)*c*c*c + (1.92565e-08)*c*c*c*c ;
2472 }
2473 if(fInternalRunNumber==21){
2474   v2c=-0.0254895 + (0.00388685)*c + (-0.000171429)*c*c + (2.87209e-06)*c*c*c + (-1.60625e-08)*c*c*c*c ;
2475   v2s=0.00237609 + (0.000664255)*c + (-4.93436e-05)*c*c + (8.86415e-07)*c*c*c + (-4.51913e-09)*c*c*c*c ;
2476 }
2477 if(fInternalRunNumber==22){
2478   v2c=0.0308671 + (-0.00365824)*c + (0.000148483)*c*c + (-2.52946e-06)*c*c*c + (1.49859e-08)*c*c*c*c ;
2479   v2s=-0.00805023 + (0.00154452)*c + (-6.89678e-05)*c*c + (1.14317e-06)*c*c*c + (-5.98991e-09)*c*c*c*c ;
2480 }
2481 if(fInternalRunNumber==23){
2482   v2c=0.122018 + (-0.0233973)*c + (0.00116385)*c*c + (-2.11232e-05)*c*c*c + (1.25836e-07)*c*c*c*c ;
2483   v2s=-0.0511879 + (0.00708815)*c + (-0.000227282)*c*c + (2.27542e-06)*c*c*c + (-4.83128e-09)*c*c*c*c ;
2484 }
2485 if(fInternalRunNumber==24){
2486   v2c=0.00888562 + (-0.000839847)*c + (2.52839e-05)*c*c + (-2.28921e-07)*c*c*c + (1.01057e-11)*c*c*c*c ;
2487   v2s=0.00814357 + (-0.000788997)*c + (2.43948e-05)*c*c + (-2.99558e-07)*c*c*c + (1.21675e-09)*c*c*c*c ;
2488 }
2489 if(fInternalRunNumber==25){
2490   v2c=0.00360833 + (0.000379977)*c + (-3.68193e-05)*c*c + (5.86861e-07)*c*c*c + (-2.01826e-09)*c*c*c*c ;
2491   v2s=-0.00642047 + (0.00060076)*c + (-2.42907e-05)*c*c + (4.12933e-07)*c*c*c + (-2.3557e-09)*c*c*c*c ;
2492 }
2493 if(fInternalRunNumber==26){
2494   v2c=-0.0131971 + (0.00189761)*c + (-7.7551e-05)*c*c + (1.25884e-06)*c*c*c + (-7.38989e-09)*c*c*c*c ;
2495   v2s=0.00679891 + (-0.000614692)*c + (1.04559e-06)*c*c + (4.34463e-07)*c*c*c + (-4.75821e-09)*c*c*c*c ;
2496 }
2497 if(fInternalRunNumber==27){
2498   v2c=-0.00887089 + (0.00126865)*c + (-7.15766e-05)*c*c + (1.43811e-06)*c*c*c + (-9.408e-09)*c*c*c*c ;
2499   v2s=-0.0302105 + (0.00349565)*c + (-0.000126108)*c*c + (1.835e-06)*c*c*c + (-9.27845e-09)*c*c*c*c ;
2500 }
2501 if(fInternalRunNumber==28){
2502   v2c=-0.0187616 + (0.0021974)*c + (-5.11412e-05)*c*c + (2.94628e-07)*c*c*c + (7.4194e-10)*c*c*c*c ;
2503   v2s=-0.0119204 + (-0.00156879)*c + (0.000132746)*c*c + (-2.67371e-06)*c*c*c + (1.61097e-08)*c*c*c*c ;
2504 }
2505 if(fInternalRunNumber==29){
2506   v2c=-0.121018 + (0.0228924)*c + (-0.00127093)*c*c + (2.48724e-05)*c*c*c + (-1.56121e-07)*c*c*c*c ;
2507   v2s=-0.0687197 + (0.015702)*c + (-0.000764074)*c*c + (1.38991e-05)*c*c*c + (-8.48812e-08)*c*c*c*c ;
2508 }
2509 if(fInternalRunNumber==30){
2510   v2c=0.01582 + (-0.0012293)*c + (1.62768e-05)*c*c + (2.36172e-07)*c*c*c + (-3.67153e-09)*c*c*c*c ;
2511   v2s=-0.00875651 + (0.00135376)*c + (-5.16763e-05)*c*c + (6.17614e-07)*c*c*c + (-1.64099e-09)*c*c*c*c ;
2512 }
2513 if(fInternalRunNumber==31){
2514   v2c=-0.044721 + (0.00756113)*c + (-0.000365759)*c*c + (6.15804e-06)*c*c*c + (-3.40075e-08)*c*c*c*c ;
2515   v2s=-0.0500828 + (-0.000695112)*c + (0.000130909)*c*c + (-3.38476e-06)*c*c*c + (2.41179e-08)*c*c*c*c ;
2516 }
2517 if(fInternalRunNumber==32){
2518   v2c=0.0306236 + (-0.00443307)*c + (0.000196517)*c*c + (-3.48059e-06)*c*c*c + (2.10266e-08)*c*c*c*c ;
2519   v2s=0.00819214 + (-0.000560558)*c + (-5.03783e-06)*c*c + (4.37862e-07)*c*c*c + (-3.93834e-09)*c*c*c*c ;
2520 }
2521 if(fInternalRunNumber==33){
2522   v2c=0.0120077 + (-0.00139976)*c + (5.85677e-05)*c*c + (-9.06209e-07)*c*c*c + (4.61403e-09)*c*c*c*c ;
2523   v2s=-0.0197677 + (0.00253589)*c + (-0.000101374)*c*c + (1.59347e-06)*c*c*c + (-8.65961e-09)*c*c*c*c ;
2524 }
2525 if(fInternalRunNumber==34){
2526   v2c=0.0759546 + (-0.0131538)*c + (0.000577001)*c*c + (-9.39636e-06)*c*c*c + (5.16814e-08)*c*c*c*c ;
2527   v2s=-0.0727864 + (0.00956655)*c + (-0.000406191)*c*c + (7.50453e-06)*c*c*c + (-4.93819e-08)*c*c*c*c ;
2528 }
2529 if(fInternalRunNumber==35){
2530   v2c=-0.00622564 + (0.00107647)*c + (-6.25644e-05)*c*c + (1.28553e-06)*c*c*c + (-8.46286e-09)*c*c*c*c ;
2531   v2s=0.0051704 + (-0.00165348)*c + (8.69531e-05)*c*c + (-1.69895e-06)*c*c*c + (1.09762e-08)*c*c*c*c ;
2532 }
2533 if(fInternalRunNumber==36){
2534   v2c=0.0117872 + (-0.000807391)*c + (-1.70999e-05)*c*c + (1.19089e-06)*c*c*c + (-1.13278e-08)*c*c*c*c ;
2535   v2s=-0.000162495 + (1.04465e-05)*c + (-1.07156e-05)*c*c + (2.54308e-07)*c*c*c + (-7.04948e-10)*c*c*c*c ;
2536 }
2537 if(fInternalRunNumber==37){
2538   v2c=-0.0119236 + (0.00290105)*c + (-0.000162086)*c*c + (3.22645e-06)*c*c*c + (-2.07075e-08)*c*c*c*c ;
2539   v2s=-0.0178535 + (0.00283024)*c + (-0.00012812)*c*c + (2.02404e-06)*c*c*c + (-1.0248e-08)*c*c*c*c ;
2540 }
2541 if(fInternalRunNumber==38){
2542   v2c=0.023002 + (-0.00433112)*c + (0.000229879)*c*c + (-4.52792e-06)*c*c*c + (2.92715e-08)*c*c*c*c ;
2543   v2s=-0.0217933 + (0.00235291)*c + (-0.000108054)*c*c + (2.08937e-06)*c*c*c + (-1.3524e-08)*c*c*c*c ;
2544 }
2545 if(fInternalRunNumber==39){
2546   v2c=0.00795585 + (-0.000372683)*c + (1.85243e-05)*c*c + (-4.28133e-07)*c*c*c + (3.02452e-09)*c*c*c*c ;
2547   v2s=-0.0849826 + (0.00529466)*c + (-0.000163148)*c*c + (2.1186e-06)*c*c*c + (-9.48208e-09)*c*c*c*c ;
2548 }
2549 if(fInternalRunNumber==40){
2550   v2c=0.00653019 + (-0.00127651)*c + (7.59264e-05)*c*c + (-1.6363e-06)*c*c*c + (1.11375e-08)*c*c*c*c ;
2551   v2s=-0.0107686 + (0.000982587)*c + (-2.97061e-05)*c*c + (4.04795e-07)*c*c*c + (-1.91289e-09)*c*c*c*c ;
2552 }
2553 if(fInternalRunNumber==41){
2554   v2c=0.012007 + (-0.00102703)*c + (3.48714e-05)*c*c + (-5.13011e-07)*c*c*c + (2.57827e-09)*c*c*c*c ;
2555   v2s=-0.00713083 + (0.00101952)*c + (-4.7675e-05)*c*c + (9.12723e-07)*c*c*c + (-5.99694e-09)*c*c*c*c ;
2556 }
2557 if(fInternalRunNumber==42){
2558   v2c=-0.0486696 + (0.00723345)*c + (-0.000324232)*c*c + (5.25317e-06)*c*c*c + (-2.78981e-08)*c*c*c*c ;
2559   v2s=-0.0485877 + (-0.00064506)*c + (9.15858e-05)*c*c + (-1.97486e-06)*c*c*c + (1.28831e-08)*c*c*c*c ;
2560 }
2561 if(fInternalRunNumber==43){
2562   v2c=0.0003541 + (-0.000236832)*c + (1.82448e-05)*c*c + (-4.20764e-07)*c*c*c + (2.96882e-09)*c*c*c*c ;
2563   v2s=0.00107616 + (-0.000462058)*c + (2.46761e-05)*c*c + (-4.33199e-07)*c*c*c + (2.47324e-09)*c*c*c*c ;
2564 }
2565 if(fInternalRunNumber==44){
2566   v2c=-0.0778259 + (0.0128797)*c + (-0.000514371)*c*c + (7.70095e-06)*c*c*c + (-3.89026e-08)*c*c*c*c ;
2567   v2s=-0.106839 + (0.00670369)*c + (-0.000184721)*c*c + (1.64581e-06)*c*c*c + (-1.47049e-09)*c*c*c*c ;
2568 }
2569 if(fInternalRunNumber==45){
2570   v2c=0.0513049 + (-0.00715356)*c + (0.000289106)*c*c + (-4.46938e-06)*c*c*c + (2.36599e-08)*c*c*c*c ;
2571   v2s=-0.0252464 + (0.00170124)*c + (-3.04376e-05)*c*c + (9.4935e-08)*c*c*c + (9.41807e-10)*c*c*c*c ;
2572 }
2573 if(fInternalRunNumber==46){
2574   v2c=0.0317697 + (-0.00121243)*c + (5.94543e-05)*c*c + (-1.26341e-06)*c*c*c + (8.49329e-09)*c*c*c*c ;
2575   v2s=-0.092276 + (0.00299941)*c + (-5.24713e-05)*c*c + (3.07331e-07)*c*c*c + (7.44139e-10)*c*c*c*c ;
2576 }
2577 if(fInternalRunNumber==47){
2578   v2c=0.0429796 + (-0.000861588)*c + (6.82123e-07)*c*c + (1.96388e-07)*c*c*c + (-1.56984e-09)*c*c*c*c ;
2579   v2s=-0.133047 + (0.00853153)*c + (-0.000287983)*c*c + (4.1155e-06)*c*c*c + (-2.00239e-08)*c*c*c*c ;
2580 }
2581 if(fInternalRunNumber==48){
2582   v2c=0.0480979 + (-0.000354969)*c + (-6.45176e-05)*c*c + (2.34807e-06)*c*c*c + (-1.9069e-08)*c*c*c*c ;
2583   v2s=-0.138315 + (0.0126874)*c + (-0.000672052)*c*c + (1.34449e-05)*c*c*c + (-8.64756e-08)*c*c*c*c ;
2584 }
2585 if(fInternalRunNumber==51){
2586   v2c=-0.0114422 + (0.00506157)*c + (-0.000227041)*c*c + (3.97457e-06)*c*c*c + (-2.35322e-08)*c*c*c*c ;
2587   v2s=-0.208135 + (0.0160488)*c + (-0.000536132)*c*c + (7.2714e-06)*c*c*c + (-3.28328e-08)*c*c*c*c ;
2588 }
2589 if(fInternalRunNumber==52){
2590   v2c=0.00113321 + (0.000922798)*c + (-6.5971e-05)*c*c + (1.40416e-06)*c*c*c + (-9.29905e-09)*c*c*c*c ;
2591   v2s=0.0550392 + (-0.00725661)*c + (0.000298984)*c*c + (-4.85473e-06)*c*c*c + (2.72063e-08)*c*c*c*c ;
2592 }
2593 if(fInternalRunNumber==53){
2594   v2c=-0.003508 + (6.54451e-06)*c + (1.28603e-05)*c*c + (-3.83674e-07)*c*c*c + (3.1053e-09)*c*c*c*c ;
2595   v2s=0.00714978 + (-0.00198147)*c + (0.000113909)*c*c + (-2.1949e-06)*c*c*c + (1.34933e-08)*c*c*c*c ;
2596 }
2597 if(fInternalRunNumber==54){
2598   v2c=-0.0024598 + (5.87979e-05)*c + (9.36516e-06)*c*c + (-2.66045e-07)*c*c*c + (1.80016e-09)*c*c*c*c ;
2599   v2s=-0.00651135 + (0.00109883)*c + (-4.41501e-05)*c*c + (6.19398e-07)*c*c*c + (-2.61977e-09)*c*c*c*c ;
2600 }
2601 if(fInternalRunNumber==55){
2602   v2c=0.0472936 + (-0.00156718)*c + (3.03713e-05)*c*c + (-1.97139e-07)*c*c*c + (-3.55205e-10)*c*c*c*c ;
2603   v2s=-0.155014 + (0.00717506)*c + (-0.000168137)*c*c + (1.55419e-06)*c*c*c + (-3.64579e-09)*c*c*c*c ;
2604 }
2605 if(fInternalRunNumber==56){
2606   v2c=0.011328 + (0.000924574)*c + (-6.803e-05)*c*c + (9.76956e-07)*c*c*c + (-2.90494e-09)*c*c*c*c ;
2607   v2s=-0.0412661 + (0.0107218)*c + (-0.000533108)*c*c + (9.75563e-06)*c*c*c + (-5.95608e-08)*c*c*c*c ;
2608 }
2609 if(fInternalRunNumber==57){
2610   v2c=0.00900452 + (0.000588343)*c + (-3.02208e-05)*c*c + (2.13239e-07)*c*c*c + (9.25727e-10)*c*c*c*c ;
2611   v2s=0.13086 + (-0.0204768)*c + (0.000928956)*c*c + (-1.63917e-05)*c*c*c + (9.8273e-08)*c*c*c*c ;
2612 }
2613 if(fInternalRunNumber==58){
2614   v2c=-0.0143951 + (0.00247373)*c + (-0.000107909)*c*c + (1.69603e-06)*c*c*c + (-9.05555e-09)*c*c*c*c ;
2615   v2s=-0.0209397 + (0.00180446)*c + (-5.97034e-05)*c*c + (8.17388e-07)*c*c*c + (-3.77989e-09)*c*c*c*c ;
2616 }
2617 if(fInternalRunNumber==59){
2618   v2c=-0.0619279 + (0.00393575)*c + (-0.00011867)*c*c + (1.45001e-06)*c*c*c + (-5.98792e-09)*c*c*c*c ;
2619   v2s=-0.248103 + (0.0158553)*c + (-0.000529631)*c*c + (7.46158e-06)*c*c*c + (-3.60091e-08)*c*c*c*c ;
2620 }
2621 if(fInternalRunNumber==60){
2622   v2c=-0.00886115 + (0.00152741)*c + (-6.94062e-05)*c*c + (1.11727e-06)*c*c*c + (-5.81351e-09)*c*c*c*c ;
2623   v2s=-0.000911417 + (0.000122344)*c + (-1.24605e-06)*c*c + (-6.29249e-08)*c*c*c + (8.81971e-10)*c*c*c*c ;
2624 }
2625 if(fInternalRunNumber==61){
2626   v2c=-0.00510166 + (0.00040596)*c + (-4.37341e-06)*c*c + (-1.15302e-07)*c*c*c + (1.53428e-09)*c*c*c*c ;
2627   v2s=-0.0170123 + (0.00220818)*c + (-9.50905e-05)*c*c + (1.59566e-06)*c*c*c + (-8.92053e-09)*c*c*c*c ;
2628 }
2629 if(fInternalRunNumber==62){
2630   v2c=-0.00332885 + (-0.000430498)*c + (5.60341e-05)*c*c + (-1.47922e-06)*c*c*c + (1.09723e-08)*c*c*c*c ;
2631   v2s=-0.0519534 + (0.00892818)*c + (-0.000434243)*c*c + (7.8407e-06)*c*c*c + (-4.69363e-08)*c*c*c*c ;
2632 }
2633 if(fInternalRunNumber==63){
2634   v2c=-0.00176529 + (0.00123097)*c + (-6.59677e-05)*c*c + (1.24591e-06)*c*c*c + (-7.90285e-09)*c*c*c*c ;
2635   v2s=-0.0181792 + (0.00222004)*c + (-7.96902e-05)*c*c + (1.13047e-06)*c*c*c + (-5.72272e-09)*c*c*c*c ;
2636 }
2637 if(fInternalRunNumber==64){
2638   v2c=0.00546184 + (-0.000449463)*c + (1.50833e-05)*c*c + (-1.72501e-07)*c*c*c + (3.58774e-10)*c*c*c*c ;
2639   v2s=0.013284 + (-0.00260671)*c + (0.000134403)*c*c + (-2.461e-06)*c*c*c + (1.47694e-08)*c*c*c*c ;
2640 }
2641 if(fInternalRunNumber==65){
2642   v2c=0.0108447 + (-0.0017699)*c + (4.62716e-05)*c*c + (-3.4647e-07)*c*c*c + (5.24536e-10)*c*c*c*c ;
2643   v2s=0.0302247 + (-0.00605136)*c + (0.000290676)*c*c + (-4.87706e-06)*c*c*c + (2.65104e-08)*c*c*c*c ;
2644 }
2645 if(fInternalRunNumber==66){
2646   v2c=-0.0100092 + (0.00144457)*c + (-6.512e-05)*c*c + (1.05356e-06)*c*c*c + (-5.52231e-09)*c*c*c*c ;
2647   v2s=-0.00380228 + (0.000610003)*c + (-3.47943e-05)*c*c + (7.21188e-07)*c*c*c + (-4.76727e-09)*c*c*c*c ;
2648 }
2649 if(fInternalRunNumber==67){
2650   v2c=-0.0229901 + (0.00400737)*c + (-0.000206024)*c*c + (3.93038e-06)*c*c*c + (-2.46388e-08)*c*c*c*c ;
2651   v2s=-0.00286041 + (0.000648607)*c + (-2.37397e-05)*c*c + (2.0372e-07)*c*c*c + (2.97679e-10)*c*c*c*c ;
2652 }
2653 if(fInternalRunNumber==68){
2654   v2c=-0.0603647 + (0.0045685)*c + (-0.000171555)*c*c + (2.66377e-06)*c*c*c + (-1.45892e-08)*c*c*c*c ;
2655   v2s=-0.21447 + (0.00909524)*c + (-0.000265298)*c*c + (3.73262e-06)*c*c*c + (-1.84484e-08)*c*c*c*c ;
2656 }
2657 if(fInternalRunNumber==69){
2658   v2c=-0.0207879 + (0.000139613)*c + (1.27802e-05)*c*c + (-6.40328e-07)*c*c*c + (6.35949e-09)*c*c*c*c ;
2659   v2s=-0.266459 + (0.0156802)*c + (-0.000481647)*c*c + (6.4597e-06)*c*c*c + (-3.047e-08)*c*c*c*c ;
2660 }
2661 if(fInternalRunNumber==70){
2662   v2c=-0.0324843 + (0.00376852)*c + (-0.000136873)*c*c + (2.01389e-06)*c*c*c + (-1.03877e-08)*c*c*c*c ;
2663   v2s=-0.0178378 + (0.00145491)*c + (-3.55421e-05)*c*c + (2.25993e-07)*c*c*c + (7.06258e-10)*c*c*c*c ;
2664 }
2665 if(fInternalRunNumber==71){
2666   v2c=-0.0138645 + (0.00167315)*c + (-6.68544e-05)*c*c + (1.07133e-06)*c*c*c + (-5.92659e-09)*c*c*c*c ;
2667   v2s=0.0359716 + (-0.00461076)*c + (0.000192868)*c*c + (-3.21409e-06)*c*c*c + (1.84468e-08)*c*c*c*c ;
2668 }
2669 if(fInternalRunNumber==72){
2670   v2c=0.0111712 + (-0.00175546)*c + (6.87862e-05)*c*c + (-8.99029e-07)*c*c*c + (3.44492e-09)*c*c*c*c ;
2671   v2s=-0.00352823 + (0.000480673)*c + (-2.66091e-05)*c*c + (5.54922e-07)*c*c*c + (-3.84534e-09)*c*c*c*c ;
2672 }
2673 if(fInternalRunNumber==73){
2674   v2c=-0.000939145 + (-0.000235945)*c + (2.51908e-05)*c*c + (-6.52435e-07)*c*c*c + (4.8218e-09)*c*c*c*c ;
2675   v2s=0.00666852 + (-0.000373105)*c + (-1.22838e-05)*c*c + (7.3621e-07)*c*c*c + (-7.2867e-09)*c*c*c*c ;
2676 }
2677 if(fInternalRunNumber==74){
2678   v2c=0.0474235 + (-0.00890967)*c + (0.000427741)*c*c + (-7.38723e-06)*c*c*c + (4.17242e-08)*c*c*c*c ;
2679   v2s=0.0142517 + (-0.00326338)*c + (0.000142729)*c*c + (-2.27831e-06)*c*c*c + (1.14084e-08)*c*c*c*c ;
2680 }
2681 if(fInternalRunNumber==75){
2682   v2c=0.0237808 + (-0.00310491)*c + (0.000134288)*c*c + (-2.27429e-06)*c*c*c + (1.32749e-08)*c*c*c*c ;
2683   v2s=-0.0258905 + (0.00417719)*c + (-0.000193901)*c*c + (3.41383e-06)*c*c*c + (-2.01695e-08)*c*c*c*c ;
2684 }
2685 if(fInternalRunNumber==76){
2686   v2c=0.0149865 + (-0.00140739)*c + (4.24537e-05)*c*c + (-4.06273e-07)*c*c*c + (5.44496e-10)*c*c*c*c ;
2687   v2s=0.00820134 + (-0.000519432)*c + (6.82211e-06)*c*c + (1.06291e-07)*c*c*c + (-1.65344e-09)*c*c*c*c ;
2688 }
2689 if(fInternalRunNumber==77){
2690   v2c=0.048123 + (-0.00608197)*c + (0.000246409)*c*c + (-4.05028e-06)*c*c*c + (2.33536e-08)*c*c*c*c ;
2691   v2s=-0.0214729 + (0.0029389)*c + (-0.000133826)*c*c + (2.29845e-06)*c*c*c + (-1.29253e-08)*c*c*c*c ;
2692 }
2693 if(fInternalRunNumber==78){
2694   v2c=-0.0172924 + (0.00202216)*c + (-6.478e-05)*c*c + (8.08739e-07)*c*c*c + (-3.73211e-09)*c*c*c*c ;
2695   v2s=-0.0355638 + (0.00121518)*c + (3.97839e-05)*c*c + (-1.62346e-06)*c*c*c + (1.27105e-08)*c*c*c*c ;
2696 }
2697 if(fInternalRunNumber==79){
2698   v2c=0.040762 + (-0.00311921)*c + (6.66868e-05)*c*c + (-1.18554e-06)*c*c*c + (1.07519e-08)*c*c*c*c ;
2699   v2s=-0.00558438 + (-0.00154821)*c + (0.000165732)*c*c + (-4.56377e-06)*c*c*c + (3.56988e-08)*c*c*c*c ;
2700 }
2701 if(fInternalRunNumber==81){
2702   v2c=-0.0247995 + (0.00179444)*c + (-0.000102085)*c*c + (2.31093e-06)*c*c*c + (-1.61924e-08)*c*c*c*c ;
2703   v2s=-0.249009 + (0.0133297)*c + (-0.000415717)*c*c + (5.84494e-06)*c*c*c + (-2.88179e-08)*c*c*c*c ;
2704 }
2705 if(fInternalRunNumber==82){
2706   v2c=-0.0249921 + (0.00199208)*c + (-7.86121e-05)*c*c + (1.40547e-06)*c*c*c + (-9.39651e-09)*c*c*c*c ;
2707   v2s=-0.254469 + (0.0150694)*c + (-0.000499848)*c*c + (7.08318e-06)*c*c*c + (-3.40806e-08)*c*c*c*c ;
2708 }
2709 if(fInternalRunNumber==83){
2710   v2c=-0.0721 + (0.00951559)*c + (-0.000431762)*c*c + (8.4632e-06)*c*c*c + (-5.79485e-08)*c*c*c*c ;
2711   v2s=-0.00864625 + (-0.00706196)*c + (0.000313773)*c*c + (-5.10066e-06)*c*c*c + (2.9662e-08)*c*c*c*c ;
2712 }
2713 if(fInternalRunNumber==84){
2714   v2c=-0.0198112 + (-0.0062788)*c + (0.000353589)*c*c + (-5.70525e-06)*c*c*c + (2.82617e-08)*c*c*c*c ;
2715   v2s=-0.0704647 + (0.0031015)*c + (-0.000244985)*c*c + (6.17258e-06)*c*c*c + (-4.43331e-08)*c*c*c*c ;
2716 }
2717 if(fInternalRunNumber==85){
2718   v2c=-0.0161955 + (-0.0017361)*c + (0.000130612)*c*c + (-2.9548e-06)*c*c*c + (2.05684e-08)*c*c*c*c ;
2719   v2s=-0.204119 + (0.00714989)*c + (-0.000162549)*c*c + (1.87459e-06)*c*c*c + (-7.79143e-09)*c*c*c*c ;
2720 }
2721 if(fInternalRunNumber==86){
2722   v2c=0.00995155 + (-0.00196822)*c + (0.00010617)*c*c + (-2.01071e-06)*c*c*c + (1.22728e-08)*c*c*c*c ;
2723   v2s=0.0148901 + (-0.002581)*c + (0.000132002)*c*c + (-2.44601e-06)*c*c*c + (1.48144e-08)*c*c*c*c ;
2724 }
2725 if(fInternalRunNumber==87){
2726   v2c=-0.00411031 + (0.000510639)*c + (-2.27814e-05)*c*c + (6.41192e-07)*c*c*c + (-5.85467e-09)*c*c*c*c ;
2727   v2s=0.0121118 + (-0.00281)*c + (0.000143754)*c*c + (-2.62557e-06)*c*c*c + (1.61207e-08)*c*c*c*c ;
2728 }
2729 if(fInternalRunNumber==88){
2730   v2c=-0.00570308 + (-0.000945116)*c + (7.42606e-05)*c*c + (-1.51212e-06)*c*c*c + (9.57993e-09)*c*c*c*c ;
2731   v2s=-0.0282666 + (0.00398475)*c + (-0.000165678)*c*c + (2.68542e-06)*c*c*c + (-1.4688e-08)*c*c*c*c ;
2732 }
2733 if(fInternalRunNumber==89){
2734   v2c=0.0845785 + (-0.012611)*c + (0.000561197)*c*c + (-9.39333e-06)*c*c*c + (5.16139e-08)*c*c*c*c ;
2735   v2s=0.149689 + (-0.0193648)*c + (0.000771629)*c*c + (-1.21901e-05)*c*c*c + (6.6443e-08)*c*c*c*c ;
2736 }
2737 if(fInternalRunNumber==91){
2738   v2c=-0.0289605 + (0.00127922)*c + (-4.47445e-06)*c*c + (-4.00293e-07)*c*c*c + (3.82966e-09)*c*c*c*c ;
2739   v2s=-0.241098 + (0.00906132)*c + (-0.000159552)*c*c + (6.47243e-07)*c*c*c + (4.61048e-09)*c*c*c*c ;
2740 }
2741 if(fInternalRunNumber==95){
2742   v2c=0.0144072 + (-0.000966056)*c + (6.47254e-05)*c*c + (-1.84164e-06)*c*c*c + (1.51588e-08)*c*c*c*c ;
2743   v2s=-0.306678 + (0.0184746)*c + (-0.000540218)*c*c + (6.36681e-06)*c*c*c + (-2.37512e-08)*c*c*c*c ;
2744 }
2745 if(fInternalRunNumber==96){
2746   v2c=0.00506182 + (-0.000632672)*c + (1.89752e-05)*c*c + (-2.14999e-07)*c*c*c + (3.12179e-10)*c*c*c*c ;
2747   v2s=-0.218693 + (0.00676786)*c + (-0.000115485)*c*c + (5.56737e-07)*c*c*c + (2.68891e-09)*c*c*c*c ;
2748 }
2749 if(fInternalRunNumber==97){
2750   v2c=-0.0225062 + (0.00204445)*c + (-6.76028e-05)*c*c + (8.92976e-07)*c*c*c + (-4.00576e-09)*c*c*c*c ;
2751   v2s=-0.260286 + (0.0139325)*c + (-0.000425631)*c*c + (5.64038e-06)*c*c*c + (-2.51554e-08)*c*c*c*c ;
2752 }
2753 if(fInternalRunNumber==98){
2754   v2c=-0.0261849 + (0.0063862)*c + (-0.000293076)*c*c + (4.37234e-06)*c*c*c + (-2.04899e-08)*c*c*c*c ;
2755   v2s=-0.17744 + (0.00236157)*c + (9.2782e-06)*c*c + (-5.63126e-07)*c*c*c + (4.9714e-09)*c*c*c*c ;
2756 }
2757 if(fInternalRunNumber==101){
2758   v2c=0.101445 + (-0.015664)*c + (0.000747255)*c*c + (-1.36979e-05)*c*c*c + (8.36903e-08)*c*c*c*c ;
2759   v2s=0.0152975 + (-0.00329946)*c + (0.00018627)*c*c + (-3.88572e-06)*c*c*c + (2.64895e-08)*c*c*c*c ;
2760 }
2761 if(fInternalRunNumber==102){
2762   v2c=-0.0344927 + (0.00470964)*c + (-0.000188941)*c*c + (3.0206e-06)*c*c*c + (-1.69735e-08)*c*c*c*c ;
2763   v2s=0.0142624 + (-0.00240623)*c + (9.16657e-05)*c*c + (-1.14663e-06)*c*c*c + (4.14195e-09)*c*c*c*c ;
2764 }
2765 if(fInternalRunNumber==103){
2766   v2c=-0.0238316 + (0.00322675)*c + (-0.000123462)*c*c + (1.76428e-06)*c*c*c + (-8.66176e-09)*c*c*c*c ;
2767   v2s=0.0227184 + (-0.00299954)*c + (0.00014422)*c*c + (-2.81018e-06)*c*c*c + (1.82546e-08)*c*c*c*c ;
2768 }
2769 if(fInternalRunNumber==104){
2770   v2c=0.0118014 + (-0.000998833)*c + (3.77271e-05)*c*c + (-6.19183e-07)*c*c*c + (3.49065e-09)*c*c*c*c ;
2771   v2s=0.000977264 + (0.00012542)*c + (1.40219e-08)*c*c + (-2.24472e-09)*c*c*c + (-3.39763e-10)*c*c*c*c ;
2772 }
2773 if(fInternalRunNumber==105){
2774   v2c=0.0312332 + (-0.00348485)*c + (0.000108654)*c*c + (-1.28443e-06)*c*c*c + (4.80188e-09)*c*c*c*c ;
2775   v2s=-0.0344234 + (0.00443313)*c + (-0.000173309)*c*c + (2.7124e-06)*c*c*c + (-1.42518e-08)*c*c*c*c ;
2776 }
2777 if(fInternalRunNumber==106){
2778   v2c=0.0369386 + (-0.00434063)*c + (0.00018894)*c*c + (-3.45688e-06)*c*c*c + (2.16617e-08)*c*c*c*c ;
2779   v2s=-0.00675251 + (0.000691284)*c + (-1.08033e-06)*c*c + (-5.57352e-07)*c*c*c + (6.28148e-09)*c*c*c*c ;
2780 }
2781 if(fInternalRunNumber==107){
2782   v2c=-0.0441483 + (0.00803486)*c + (-0.00034871)*c*c + (4.76349e-06)*c*c*c + (-1.81393e-08)*c*c*c*c ;
2783   v2s=-0.141046 + (0.0165183)*c + (-0.000592546)*c*c + (8.956e-06)*c*c*c + (-4.86662e-08)*c*c*c*c ;
2784 }
2785 if(fInternalRunNumber==108){
2786   v2c=-0.0531169 + (0.00637348)*c + (-0.000218334)*c*c + (2.66873e-06)*c*c*c + (-9.86859e-09)*c*c*c*c ;
2787   v2s=-0.00809622 + (-0.000552159)*c + (6.27984e-05)*c*c + (-1.33136e-06)*c*c*c + (7.45069e-09)*c*c*c*c ;
2788 }
2789 if(fInternalRunNumber==109){
2790   v2c=-0.00189593 + (0.00153994)*c + (-9.92105e-05)*c*c + (1.98756e-06)*c*c*c + (-1.24677e-08)*c*c*c*c ;
2791   v2s=0.0252107 + (-0.00330651)*c + (0.000128362)*c*c + (-1.89951e-06)*c*c*c + (9.45477e-09)*c*c*c*c ;
2792 }
2793 if(fInternalRunNumber==112){
2794   v2c=-0.00939028 + (0.00128226)*c + (-5.24125e-05)*c*c + (8.50382e-07)*c*c*c + (-4.65785e-09)*c*c*c*c ;
2795   v2s=0.00201633 + (-0.00116151)*c + (7.59727e-05)*c*c + (-1.59962e-06)*c*c*c + (1.04645e-08)*c*c*c*c ;
2796 }
2797 if(fInternalRunNumber==114){
2798   v2c=-0.00783351 + (0.00122277)*c + (-3.37736e-05)*c*c + (3.12136e-07)*c*c*c + (-9.19359e-10)*c*c*c*c ;
2799   v2s=0.0103833 + (-0.00215968)*c + (9.96942e-05)*c*c + (-1.71974e-06)*c*c*c + (1.01511e-08)*c*c*c*c ;
2800 }
2801 if(fInternalRunNumber==115){
2802   v2c=0.00588164 + (-0.00083683)*c + (5.13684e-05)*c*c + (-1.14177e-06)*c*c*c + (8.11899e-09)*c*c*c*c ;
2803   v2s=0.0037561 + (-0.000786163)*c + (3.88084e-05)*c*c + (-7.73678e-07)*c*c*c + (5.27396e-09)*c*c*c*c ;
2804 }
2805 if(fInternalRunNumber==116){
2806   v2c=0.0610685 + (-0.00392239)*c + (6.23272e-05)*c*c + (9.84282e-07)*c*c*c + (-1.75418e-08)*c*c*c*c ;
2807   v2s=-0.106295 + (0.0124333)*c + (-0.000510299)*c*c + (8.51235e-06)*c*c*c + (-4.9787e-08)*c*c*c*c ;
2808 }
2809 if(fInternalRunNumber==117){
2810   v2c=-0.0132264 + (0.0024143)*c + (-0.000106307)*c*c + (1.74632e-06)*c*c*c + (-9.79752e-09)*c*c*c*c ;
2811   v2s=-0.00788943 + (0.000443907)*c + (-1.83984e-05)*c*c + (3.81018e-07)*c*c*c + (-2.52604e-09)*c*c*c*c ;
2812 }
2813 if(fInternalRunNumber==118){
2814   v2c=-0.045835 + (0.00849686)*c + (-0.000408675)*c*c + (7.3031e-06)*c*c*c + (-4.34655e-08)*c*c*c*c ;
2815   v2s=-0.0455199 + (0.00718851)*c + (-0.000340502)*c*c + (6.00457e-06)*c*c*c + (-3.51883e-08)*c*c*c*c ;
2816 }
2817 if(fInternalRunNumber==119){
2818   v2c=-0.00627236 + (0.00146712)*c + (-9.05375e-05)*c*c + (1.96962e-06)*c*c*c + (-1.35867e-08)*c*c*c*c ;
2819   v2s=-0.00544096 + (0.000598445)*c + (-7.61032e-06)*c*c + (-2.09579e-07)*c*c*c + (3.00091e-09)*c*c*c*c ;
2820 }
2821 if(fInternalRunNumber==120){
2822   v2c=-0.0376823 + (0.00376965)*c + (-8.01573e-05)*c*c + (-2.80792e-09)*c*c*c + (6.79508e-09)*c*c*c*c ;
2823   v2s=-0.0675803 + (0.00849852)*c + (-0.000298266)*c*c + (3.93178e-06)*c*c*c + (-1.68582e-08)*c*c*c*c ;
2824 }
2825 if(fInternalRunNumber==121){
2826   v2c=0.0285153 + (-0.00340253)*c + (0.000130854)*c*c + (-2.02449e-06)*c*c*c + (1.09594e-08)*c*c*c*c ;
2827   v2s=0.0224991 + (-0.00524188)*c + (0.000269226)*c*c + (-4.93926e-06)*c*c*c + (2.97604e-08)*c*c*c*c ;
2828 }
2829 if(fInternalRunNumber==122){
2830   v2c=0.00379321 + (0.000532446)*c + (-6.16649e-05)*c*c + (1.67786e-06)*c*c*c + (-1.29277e-08)*c*c*c*c ;
2831   v2s=-0.011158 + (0.00256124)*c + (-0.000150921)*c*c + (3.22488e-06)*c*c*c + (-2.16809e-08)*c*c*c*c ;
2832 }
2833 if(fInternalRunNumber==123){
2834   v2c=-0.0607006 + (0.00831906)*c + (-0.000398705)*c*c + (6.66385e-06)*c*c*c + (-3.5222e-08)*c*c*c*c ;
2835   v2s=0.0389204 + (-0.00798361)*c + (0.000297135)*c*c + (-3.79751e-06)*c*c*c + (1.47914e-08)*c*c*c*c ;
2836 }
2837 if(fInternalRunNumber==124){
2838   v2c=-0.0288952 + (0.00325728)*c + (-0.000104126)*c*c + (1.29062e-06)*c*c*c + (-6.03539e-09)*c*c*c*c ;
2839   v2s=0.0318884 + (-0.00301585)*c + (0.000131913)*c*c + (-2.5702e-06)*c*c*c + (1.71313e-08)*c*c*c*c ;
2840 }
2841 if(fInternalRunNumber==125){
2842   v2c=-0.00413859 + (0.00024613)*c + (5.33902e-06)*c*c + (-3.86963e-07)*c*c*c + (3.89277e-09)*c*c*c*c ;
2843   v2s=-0.0492298 + (0.00624793)*c + (-0.000243844)*c*c + (3.69001e-06)*c*c*c + (-1.9156e-08)*c*c*c*c ;
2844 }
2845 if(fInternalRunNumber==126){
2846   v2c=-0.018816 + (0.00134076)*c + (-4.50486e-05)*c*c + (5.84259e-07)*c*c*c + (-2.34235e-09)*c*c*c*c ;
2847   v2s=-0.295912 + (0.0159608)*c + (-0.000497569)*c*c + (6.92442e-06)*c*c*c + (-3.35159e-08)*c*c*c*c ;
2848 }
2849 if(fInternalRunNumber==127){
2850   v2c=-0.039222 + (0.00400933)*c + (-0.0001341)*c*c + (1.65214e-06)*c*c*c + (-6.33829e-09)*c*c*c*c ;
2851   v2s=-0.290884 + (0.0147027)*c + (-0.000429021)*c*c + (5.64116e-06)*c*c*c + (-2.58784e-08)*c*c*c*c ;
2852 }
2853 if(fInternalRunNumber==128){
2854   v2c=0.0168685 + (-0.00285775)*c + (0.000141099)*c*c + (-2.48846e-06)*c*c*c + (1.41845e-08)*c*c*c*c ;
2855   v2s=-0.278776 + (0.0113771)*c + (-0.000279052)*c*c + (3.16144e-06)*c*c*c + (-1.18545e-08)*c*c*c*c ;
2856 }
2857 if(fInternalRunNumber==129){
2858   v2c=0.0490322 + (-0.00736827)*c + (0.000304529)*c*c + (-4.89472e-06)*c*c*c + (2.60746e-08)*c*c*c*c ;
2859   v2s=-0.0961198 + (0.0179264)*c + (-0.000914589)*c*c + (1.70145e-05)*c*c*c + (-1.03423e-07)*c*c*c*c ;
2860 }
2861 if(fInternalRunNumber==131){
2862   v2c=0.075924 + (-0.0108106)*c + (0.000493471)*c*c + (-8.39074e-06)*c*c*c + (4.67109e-08)*c*c*c*c ;
2863   v2s=0.0700446 + (-0.0105027)*c + (0.00048267)*c*c + (-8.4422e-06)*c*c*c + (4.91292e-08)*c*c*c*c ;
2864 }
2865 if(fInternalRunNumber==132){
2866   v2c=0.00752887 + (-0.00142672)*c + (7.37857e-05)*c*c + (-1.36894e-06)*c*c*c + (8.28054e-09)*c*c*c*c ;
2867   v2s=-0.0121719 + (0.00188256)*c + (-8.67104e-05)*c*c + (1.50514e-06)*c*c*c + (-8.79759e-09)*c*c*c*c ;
2868 }
2869 if(fInternalRunNumber==133){
2870   v2c=0.00242068 + (-0.000274576)*c + (1.06631e-05)*c*c + (-1.45829e-07)*c*c*c + (7.18001e-10)*c*c*c*c ;
2871   v2s=0.0113182 + (-0.000529549)*c + (-7.75839e-06)*c*c + (4.86382e-07)*c*c*c + (-3.98634e-09)*c*c*c*c ;
2872 }
2873 if(fInternalRunNumber==134){
2874   v2c=-0.0438721 + (0.00906564)*c + (-0.000407471)*c*c + (6.6254e-06)*c*c*c + (-3.61035e-08)*c*c*c*c ;
2875   v2s=-0.276264 + (0.0107322)*c + (-0.000281971)*c*c + (3.27634e-06)*c*c*c + (-1.17426e-08)*c*c*c*c ;
2876 }
2877 if(fInternalRunNumber==135){
2878   v2c=0.0224469 + (-0.00324762)*c + (0.000165861)*c*c + (-3.16496e-06)*c*c*c + (2.011e-08)*c*c*c*c ;
2879   v2s=-0.308891 + (0.0136721)*c + (-0.00032343)*c*c + (2.93575e-06)*c*c*c + (-5.7637e-09)*c*c*c*c ;
2880 }
2881 if(fInternalRunNumber==136){
2882   v2c=0.0279685 + (-0.00441)*c + (0.000294342)*c*c + (-6.56635e-06)*c*c*c + (4.53669e-08)*c*c*c*c ;
2883   v2s=-0.273287 + (0.00917778)*c + (-0.000144155)*c*c + (1.59881e-07)*c*c*c + (8.63874e-09)*c*c*c*c ;
2884 }
2885
2886 /*
2887   if(fRunNumber>=1 && fRunNumber<=45){ // 137161 -137848 //period 1
2888 //    v2c =-3.523975e-03+8.132703e-04*x+-4.107765e-05*x*x+7.308763e-07*x*x*x+-4.285144e-09*x*x*x*x;
2889 //    v2s =-1.555676e-02+1.075371e-03*x+-3.597818e-05*x*x+4.927188e-07*x*x*x+-2.217080e-09*x*x*x*x;
2890      v2c=-2.007559e-05+3.143290e-04*x-2.065995e-05*x*x+4.071671e-07*x*x*x-2.562114e-09*x*x*x*x ;
2891      v2s=-4.920033e-03+6.807606e-04*x-3.065278e-05*x*x+5.132261e-07*x*x*x-2.767644e-09*x*x*x*x ;
2892
2893   }
2894   if(fRunNumber>=46 && fRunNumber<=58){ //138125-138275//period 2
2895     v2c = 1.591471e-02+2.772445e-04*x-2.229627e-05*x*x+3.803341e-07*x*x*x-2.056551e-09*x*x*x*x;
2896     v2s =-5.800393e-02+2.007102e-03*x-3.292094e-05*x*x+7.266378e-08*x*x*x+1.828472e-09*x*x*x*x ;
2897   }
2898   if(fRunNumber>=59 && fRunNumber<=78){ //138359-138730//period 3
2899     v2c =-8.010228e-03+8.947184e-04*x-4.192609e-05*x*x+7.388210e-07*x*x*x-4.280541e-09*x*x*x*x ;
2900     v2s =-3.203888e-02+2.122650e-03*x-7.549389e-05*x*x+1.161380e-06*x*x*x-6.156621e-09*x*x*x*x ;
2901   }
2902   if(fRunNumber>=79){//period 4
2903     v2c =-1.025964e-02+1.128454e-03*x-4.603917e-05*x*x+7.386637e-07*x*x*x-4.035342e-09*x*x*x*x ;
2904     v2s =-9.746961e-02+4.672675e-03*x-1.311440e-04*x*x+1.650922e-06*x*x*x-7.083187e-09*x*x*x*x ;
2905   }
2906 */
2907   return phi+v2c*TMath::Sin(2.*phi)-v2s*TMath::Cos(2.*phi) ;
2908
2909 }
2910 //____________________________________________________________________________
2911 Double_t  AliAnalysisTaskPi0Flow::ApplyFlattening(Double_t phi, Double_t c){
2912
2913
2914   if(fInternalRunNumber==39){ //137748
2915     Double_t v2c= 0.0222974-0.00132297*c+0.000204021*c*c-4.49827e-06*c*c*c+2.76986e-08*c*c*c*c ;
2916     Double_t v2s=0.00615401-7.76462e-05*c ;
2917     return phi+v2c*TMath::Sin(2.*phi)-v2s*TMath::Cos(2.*phi) ;
2918   }
2919
2920
2921   //Periods 1,2,3
2922   //fRunNumber - run index
2923   if(fInternalRunNumber>=1 && fInternalRunNumber<79){
2924     Double_t v2c=4.40516e-04+TMath::Exp(-4.71923-7.62089e-02*c) ;
2925     Double_t v2s=1.79859e-03+TMath::Exp(-4.99649-7.87523e-02*c) ;
2926
2927     return phi+v2c*TMath::Sin(2.*phi)-v2s*TMath::Cos(2.*phi) ;
2928   }
2929
2930   //period4
2931    //138826, 138828, 138830
2932   if(fInternalRunNumber==83 ||fInternalRunNumber==84 ||fInternalRunNumber==85 ){
2933      Double_t v2c = -1.4e-03 ;
2934      Double_t v2s =  -5.60117e-02-TMath::Exp( -1.62827e+00-9.96071e-02*c);
2935     return phi + v2c*TMath::Sin(2.*phi) - v2s*TMath::Cos(2.*phi) ;
2936   }
2937   //Runs 138871,138872
2938   if(fInternalRunNumber==87 ||fInternalRunNumber==88 ){
2939     Double_t v2c=-0.00518371 ;
2940     Double_t v2s= 0.00633324 ;
2941     return phi + v2c*TMath::Sin(2.*phi) - v2s*TMath::Cos(2.*phi) ;
2942   }
2943   //run 139029
2944   if(fInternalRunNumber==102 ){
2945      Double_t v2c = -0.00354633 ;
2946      Double_t v2s =  0.00418512 ;
2947      return phi+v2c*TMath::Sin(2.*phi)-v2s*TMath::Cos(2.*phi) ;
2948   }
2949   //run 139110
2950 /*
2951   if(fRunNumber==110 ){
2952     Double_t v2c=4.78327e-04+TMath::Exp(-3.21625    -6.74306e-02*c) ;
2953     Double_t v2s=3.50176e-02+TMath::Exp(-9.73541e-01-9.19214e-02*c) ;
2954     return phi + v2c*TMath::Sin(2.*phi) - v2s*TMath::Cos(2.*phi) ;
2955   }
2956 */
2957   if(fInternalRunNumber>=79 && fInternalRunNumber<=139){
2958     Double_t v2c= 4.78327e-04 + TMath::Exp(-6.70587-2.33120e-02*c) ;
2959     Double_t v2s=-2.57731e-03 + TMath::Exp(-2.75493-1.05166e-01*c) ;
2960
2961     return phi + v2c*TMath::Sin(2.*phi) - v2s*TMath::Cos(2.*phi) ;
2962   }
2963   return phi ;
2964
2965 }
2966 //____________________________________________________________________________
2967 Double_t  AliAnalysisTaskPi0Flow::CoreEnergy(AliVCluster * clu){
2968   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
2969
2970   //Can not use already calculated coordinates?
2971   //They have incidence correction...
2972   const Double_t distanceCut =3.5 ;
2973   const Double_t logWeight=4.5 ;
2974
2975   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2976 // Calculates the center of gravity in the local PHOS-module coordinates
2977   Float_t wtot = 0;
2978   Double_t xc[100]={0} ;
2979   Double_t zc[100]={0} ;
2980   Double_t x = 0 ;
2981   Double_t z = 0 ;
2982   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
2983   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2984     Int_t relid[4] ;
2985     Float_t xi ;
2986     Float_t zi ;
2987     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
2988     fPHOSGeo->RelPosInModule(relid, xi, zi);
2989     xc[iDigit]=xi ;
2990     zc[iDigit]=zi ;
2991     if (clu->E()>0 && elist[iDigit]>0) {
2992       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
2993       x    += xc[iDigit] * w ;
2994       z    += zc[iDigit] * w ;
2995       wtot += w ;
2996     }
2997   }
2998   if (wtot>0) {
2999     x /= wtot ;
3000     z /= wtot ;
3001   }
3002   Double_t coreE=0. ;
3003   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
3004     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
3005     if(distance < distanceCut)
3006       coreE += elist[iDigit] ;
3007   }
3008   //Apply non-linearity correction
3009   return fNonLinCorr->Eval(coreE) ;
3010 }
3011 //____________________________________________________________________________
3012 Bool_t  AliAnalysisTaskPi0Flow::AreNeibors(Int_t id1,Int_t id2){
3013   // return true if absId are "Neighbors" (adjacent, including diagornaly,)
3014   // false if not.
3015
3016   Int_t relid1[4] ;
3017   fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
3018
3019   Int_t relid2[4] ;
3020   fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
3021
3022   // if inside the same PHOS module
3023   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) {
3024     const Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
3025     const Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
3026
3027     // and if diff in both direction is 1 or less
3028     if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
3029       return true; // are neighbors
3030   }
3031
3032   // else false
3033   return false;
3034 }
3035 //____________________________________________________________________________
3036 void  AliAnalysisTaskPi0Flow::Reclusterize(AliVCluster * clu){
3037   //Re-clusterize to make continues cluster
3038
3039   const Int_t oldMulDigit=clu->GetNCells() ;
3040   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
3041   UShort_t * dlist = clu->GetCellsAbsId();
3042
3043   Int_t index[oldMulDigit] ;
3044   Bool_t used[oldMulDigit] ;
3045   for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
3046   Int_t inClu=0 ;
3047   Double_t eMax=0. ;
3048   //find maximum
3049   for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
3050     if(eMax<elist[iDigit]){
3051       eMax=elist[iDigit];
3052       index[0]=iDigit ;
3053       inClu=1 ;
3054     }
3055   }
3056   if(inClu==0){ //empty cluster
3057     return ;
3058   }
3059   used[index[0]]=kTRUE ; //mark as used
3060   for(Int_t i=0; i<inClu; i++){
3061     for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
3062        if(used[iDigit]) //already used
3063          continue ;
3064        if(AreNeibors(dlist[index[i]],dlist[iDigit])){
3065          index[inClu]= iDigit ;
3066          inClu++ ;
3067          used[iDigit]=kTRUE ;
3068        }
3069     }
3070   }
3071
3072   if(inClu==oldMulDigit) //no need to modify
3073     return ;
3074
3075   clu->SetNCells(inClu);
3076   //copy
3077   UShort_t tmpD[oldMulDigit] ;
3078   Double_t tmpE[oldMulDigit] ;
3079   for(Int_t i=0; i<oldMulDigit; i++){
3080     tmpD[i]=dlist[i] ;
3081     tmpE[i]=elist[i] ;
3082   }
3083   //change order of digits in list so that
3084   //first inClu cells were true ones
3085   for(Int_t i=0; i<inClu; i++){
3086     dlist[i]=tmpD[index[i]] ;
3087     elist[i]=tmpE[index[i]] ;
3088   }
3089
3090
3091 }
3092
3093 //_____________________________________________________________________________
3094 void AliAnalysisTaskPi0Flow::SetMisalignment(){
3095   // sets the misalignment vertex if ESD
3096   if( fEventESD ) {
3097     for(Int_t mod=0; mod<5; mod++) {
3098       const TGeoHMatrix* modMatrix = fEvent->GetPHOSMatrix(mod);
3099       if( ! modMatrix) {
3100         if( fDebug )
3101           AliInfo(Form("no PHOS Geometric Misalignment Matrix for module %d", mod));
3102         continue;
3103       }
3104       else {
3105         fPHOSGeo->SetMisalMatrix(modMatrix, mod);
3106         if( fDebug )
3107           AliInfo("PHOS Geometric Misalignment Matrix");
3108       }
3109     }
3110   }
3111 }
3112
3113 //_____________________________________________________________________________
3114 void AliAnalysisTaskPi0Flow::SetV0Calibration(){
3115     // assigns: fMultV0, fV0Cpol, fV0Apol, fMeanQ, and fWidthQ
3116
3117     int runNumber = this->fRunNumber;
3118     
3119     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
3120     TFile *foadb = TFile::Open(oadbfilename.Data());
3121
3122     if(!foadb){
3123         AliError(Form("OADB file %s cannot be opened\n", oadbfilename.Data()));
3124         AliError("V0 Calibration not set !\n");
3125         return;
3126     }
3127
3128     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
3129     if(!cont){
3130         AliError("OADB object hMultV0BefCorr is not available in the file");
3131         AliError("V0 Calibration not set!\n");
3132         return;
3133     }
3134
3135     if(!(cont->GetObject(runNumber))){
3136         AliError(Form("OADB object hMultV0BefCorr is not available for run %i, trying 137366)",runNumber));
3137         runNumber = 137366;
3138     }
3139     if(!(cont->GetObject(runNumber))){
3140         AliError(Form("OADB object hMultV0BefCorr is not available for run %i ",runNumber));
3141         AliError("V0 Calibration not set!\n");
3142         return;
3143     }
3144
3145     if( fDebug )  AliInfo("Setting V0 calibration") ;
3146     fMultV0 = ((TH2F *) cont->GetObject(runNumber))->ProfileX();
3147
3148     TF1 *fpol0 = new TF1("fpol0","pol0");
3149     fMultV0->Fit(fpol0,"","",0,31);
3150     fV0Cpol = fpol0->GetParameter(0);
3151     fMultV0->Fit(fpol0,"","",32,64);
3152     fV0Apol = fpol0->GetParameter(0);
3153
3154     for(Int_t iside=0;iside<2;iside++){
3155         for(Int_t icoord=0;icoord<2;icoord++){
3156             for(Int_t i=0;i  < kNCenBins;i++){
3157                 char namecont[100];
3158                 if(iside==0 && icoord==0)
3159                     snprintf(namecont,100,"hQxc2_%i",i);
3160                 else if(iside==1 && icoord==0)
3161                     snprintf(namecont,100,"hQxa2_%i",i);
3162                 else if(iside==0 && icoord==1)
3163                     snprintf(namecont,100,"hQyc2_%i",i);
3164                 else if(iside==1 && icoord==1)
3165                     snprintf(namecont,100,"hQya2_%i",i);
3166
3167                 cont = (AliOADBContainer*) foadb->Get(namecont);
3168                 if(!cont){
3169                     AliError(Form("OADB object %s is not available in the file %s", namecont, oadbfilename.Data()));
3170                     AliError("V0 Calibration not fully set!\n");
3171                     return;
3172                 }
3173
3174                 if(!(cont->GetObject(runNumber))){
3175                     AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
3176                     runNumber = 137366;
3177                 }
3178                 if(!(cont->GetObject(runNumber))){
3179                   AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
3180                   AliError("V0 Calibration not fully set!\n");
3181                   return;
3182                 }
3183                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetMean();
3184                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetRMS();
3185
3186                 //for v3
3187 //              if(iside==0 && icoord==0)
3188 //                  snprintf(namecont,100,"hQxc3_%i",i);
3189 //              else if(iside==1 && icoord==0)
3190 //                  snprintf(namecont,100,"hQxa3_%i",i);
3191 //              else if(iside==0 && icoord==1)
3192 //                  snprintf(namecont,100,"hQyc3_%i",i);
3193 //              else if(iside==1 && icoord==1)
3194 //                  snprintf(namecont,100,"hQya3_%i",i);
3195 // 
3196 //              cont = (AliOADBContainer*) foadb->Get(namecont);
3197 //              if(!cont){
3198 //                  AliError(Form("OADB object %s is not available in the file",namecont));
3199 //                  AliError("V0 Calibration not fully set!\n");
3200 //                  return;
3201 //              }
3202 // 
3203 //              if(!(cont->GetObject(runNumber))){
3204 //                  AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
3205 //                  runNumber = 137366;
3206 //              }
3207 //              if(!(cont->GetObject(runNumber))){
3208 //                AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
3209 //                AliError("V0 Calibration not fully set!\n");
3210 //                return;
3211 //              }
3212 //              fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
3213 //              fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
3214
3215             }
3216         }
3217     }
3218
3219     delete fpol0; fpol0=0;
3220 }
3221
3222 //_____________________________________________________________________________
3223 void AliAnalysisTaskPi0Flow::SetESDTrackCuts()
3224 {
3225   if( fEventESD ) {
3226     // Create ESD track cut
3227     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
3228     //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
3229     fESDtrackCuts->SetRequireTPCRefit(kTRUE);
3230   }
3231 }
3232
3233 //_____________________________________________________________________________
3234 void AliAnalysisTaskPi0Flow::SetGeometry()
3235 {
3236   // Initialize the PHOS geometry
3237   if( kLHC10h == fPeriod && fEventESD ) {
3238     TGeoManager::Import("geometry.root"); //TODO: should perhaps not be done
3239     fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
3240     if( ! fPHOSGeo )
3241       AliError("geometry (fPHOSGeo) not initialised");
3242   }
3243
3244   //Init geometry
3245   if(!fPHOSGeo){
3246     AliOADBContainer geomContainer("phosGeo");
3247     geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
3248     TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
3249     fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
3250     for(Int_t mod=0; mod<5; mod++) {
3251       if(!matrixes->At(mod)) {
3252         continue;
3253         if( fDebug > 0 )
3254           AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
3255       }
3256       else {
3257         fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
3258         if( fDebug >1 )
3259           AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
3260       }
3261     }
3262   } 
3263 }
3264
3265 //_____________________________________________________________________________
3266 void AliAnalysisTaskPi0Flow::SetPHOSCalibData()
3267 {
3268   if( fPHOSCalibData )
3269     delete fPHOSCalibData; 
3270   fPHOSCalibData = 0;
3271   
3272   // Calibration only needed for ESD
3273   if( fEventESD /*&& */ ) {
3274     if( kLHC10h == fPeriod && fEventESD ) {
3275       //We have to apply re-calibration for pass1 LCH10h
3276       // Initialize decalibration factors in the form of the OCDB object
3277       AliCDBManager * man = AliCDBManager::Instance();
3278       man->SetRun(140000) ; //TODO; revise, this should probably not b done.
3279       man->SetDefaultStorage("local://OCDB");
3280     }
3281     fPHOSCalibData = new AliPHOSCalibData();
3282   }
3283 }
3284
3285 //_____________________________________________________________________________
3286 void AliAnalysisTaskPi0Flow::SetVertex()
3287 {
3288   const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
3289   if( primaryVertex ) {
3290     fVertex[0] = primaryVertex->GetX();
3291     fVertex[1] = primaryVertex->GetY();
3292     fVertex[2] = primaryVertex->GetZ();
3293   }
3294   else {
3295     AliError("Event has 0x0 Primary Vertex, defaulting to origo");
3296     fVertex[0] = 0;
3297     fVertex[1] = 0;
3298     fVertex[2] = 0;
3299   }
3300   fVertexVector = TVector3(fVertex);
3301   FillHistogram("hZvertex", fVertexVector.z(), fInternalRunNumber-0.5);
3302   
3303   if( fDebug > 1 )
3304     AliInfo(Form("Vertex is set to (%.1f,%.1f,%.1f)", fVertex[0], fVertex[1], fVertex[2]));
3305
3306   fVtxBin=0 ;// No support for vtx binning implemented.
3307 }
3308
3309 //_____________________________________________________________________________
3310 Bool_t AliAnalysisTaskPi0Flow::RejectEventVertex()
3311 {
3312   if( ! fEvent->GetPrimaryVertex() )
3313     return true; // reject
3314   if ( TMath::Abs(fVertexVector.z()) > 10. )
3315     return true; // reject
3316
3317   return false; // accept event.
3318 }
3319
3320 //_____________________________________________________________________________
3321 void AliAnalysisTaskPi0Flow::SetCentrality()
3322 {
3323   AliCentrality *centrality = fEvent->GetCentrality();
3324   if( centrality )
3325     fCentralityV0M=centrality->GetCentralityPercentile("V0M");
3326   else {
3327     AliError("Event has 0x0 centrality");
3328     fCentralityV0M = -1.;
3329   }
3330   FillHistogram("hCentrality",fCentralityV0M,fInternalRunNumber-0.5) ;
3331
3332   fCentBin = GetCentralityBin(fCentralityV0M);
3333
3334   if ( fDebug )
3335     AliInfo(Form("Centrality (bin) is: %f (%d)", fCentralityV0M, fCentBin));
3336 }
3337
3338 //_____________________________________________________________________________
3339 Bool_t AliAnalysisTaskPi0Flow::RejectCentrality()
3340 {
3341   if( ! fEvent->GetCentrality() )
3342     return true; // reject
3343 //   if( fCentralityV0M <= 0. || fCentralityV0M>80. )
3344 //     return true; // reject
3345     
3346   int lastBinUpperIndex = fCentEdges.GetSize() -1;
3347   if( fCentralityV0M > fCentEdges[lastBinUpperIndex] ) {
3348     if( fDebug )
3349       AliInfo("Rejecting due to centrality outside of binning.");
3350     return true; // reject
3351   }
3352   if( fCentralityV0M < fCentEdges[0] ) {
3353     if( fDebug )
3354       AliInfo("Rejecting due to centrality outside of binning.");
3355     return true; // reject
3356   }
3357
3358   return false;
3359 }
3360
3361
3362 //_____________________________________________________________________________
3363 void AliAnalysisTaskPi0Flow::EvalReactionPlane()
3364 {
3365   // assigns: fHaveTPCRP and fRP
3366   // also does a few histogram fills
3367
3368   AliEventplane *eventPlane = fEvent->GetEventplane();
3369   if( ! eventPlane ) { AliError("Event has no event plane"); return; }
3370   
3371   Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
3372   FillHistogram("phiRP",reactionPlaneQ,fCentralityV0M) ;
3373
3374   if(reactionPlaneQ==999 || reactionPlaneQ < 0.){ //reaction plain was not defined
3375     if( fDebug ) AliInfo(Form("No Q Reaction Plane, value is %f", reactionPlaneQ));
3376     fHaveTPCRP = kFALSE;
3377   }
3378   else{
3379     if( fDebug ) AliInfo(Form("Q Reaction Plane is %f", reactionPlaneQ));
3380     fHaveTPCRP = kTRUE;
3381   }
3382
3383   if(fHaveTPCRP){
3384     if( kLHC10h == fPeriod && fEventESD )
3385       fRP = ApplyFlattening(reactionPlaneQ, fCentralityV0M) ;
3386     else if ( fEventAOD )
3387       fRP = reactionPlaneQ;
3388     else
3389       AliError("Don't know how to set fRP");
3390
3391     while(fRP<0)  fRP+=TMath::Pi();
3392     while(fRP>TMath::Pi())  fRP-=TMath::Pi();
3393     FillHistogram("phiRPflat",fRP,fCentralityV0M) ;
3394     Double_t dPsi = eventPlane->GetQsubRes() ;
3395     FillHistogram("cos2AC",TMath::Cos(2.*dPsi),fCentralityV0M) ;
3396   }
3397   else
3398     fRP=0.;
3399 }
3400
3401
3402 //____________________________________________________________________________
3403 void  AliAnalysisTaskPi0Flow::EvalV0ReactionPlane(){
3404   // set: fRPV0A and fRPV0C
3405
3406   //VZERO data
3407   AliVVZERO* v0 = fEvent->GetVZEROData();
3408
3409   //reset Q vector info
3410   Double_t Qxa2 = 0, Qya2 = 0;
3411   Double_t Qxc2 = 0, Qyc2 = 0;
3412
3413   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
3414     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
3415     Float_t multv0 = v0->GetMultiplicity(iv0);
3416     if (iv0 < 32){ // V0C
3417       Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
3418       Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
3419     } else {       // V0A
3420       Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
3421       Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
3422     }
3423   }
3424
3425   Int_t iC = -1;
3426   // centrality bins
3427   if(fCentralityV0M < 5) iC = 0;
3428   else if(fCentralityV0M < 10) iC = 1;
3429   else if(fCentralityV0M < 20) iC = 2;
3430   else if(fCentralityV0M < 30) iC = 3;
3431   else if(fCentralityV0M < 40) iC = 4;
3432   else if(fCentralityV0M < 50) iC = 5;
3433   else if(fCentralityV0M < 60) iC = 6;
3434   else if(fCentralityV0M < 70) iC = 7;
3435   else iC = 8;
3436
3437     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
3438   Double_t Qxamean2 = fMeanQ[iC][1][0];
3439   Double_t Qxarms2  = fWidthQ[iC][1][0];
3440   Double_t Qyamean2 = fMeanQ[iC][1][1];
3441   Double_t Qyarms2  = fWidthQ[iC][1][1];
3442
3443   Double_t Qxcmean2 = fMeanQ[iC][0][0];
3444   Double_t Qxcrms2  = fWidthQ[iC][0][0];
3445   Double_t Qycmean2 = fMeanQ[iC][0][1];
3446   Double_t Qycrms2  = fWidthQ[iC][0][1];
3447
3448   Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
3449   Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
3450   Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
3451   Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
3452
3453   fRPV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
3454   fRPV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
3455
3456   while(fRPV0A<0)fRPV0A+=TMath::Pi() ;
3457   while(fRPV0A>TMath::Pi())fRPV0A-=TMath::Pi() ;
3458   while(fRPV0C<0)fRPV0C+=TMath::Pi() ;
3459   while(fRPV0C>TMath::Pi())fRPV0C-=TMath::Pi() ;
3460
3461   // Flattening
3462   if( kLHC10h == fPeriod && fEventESD ) {
3463     fRPV0A=ApplyFlatteningV0A(fRPV0A,fCentralityV0M) ;
3464     while(fRPV0A<0)fRPV0A+=TMath::Pi() ;
3465     while(fRPV0A>TMath::Pi())fRPV0A-=TMath::Pi() ;
3466     fRPV0C=ApplyFlatteningV0C(fRPV0C,fCentralityV0M) ;
3467     while(fRPV0C<0)fRPV0C+=TMath::Pi() ;
3468     while(fRPV0C>TMath::Pi())fRPV0C-=TMath::Pi() ;
3469   }
3470   else if ( fEventESD ) {
3471     AliError("Don't know how to, or if to, apply flattening for ESD");
3472   }
3473     
3474   if( fDebug )
3475     AliInfo(Form("V0 Reaction Plane is: A side: %f, C side: %f", fRPV0A, fRPV0C));
3476
3477   FillHistogram("phiRPV0A",fRPV0A,fCentralityV0M);
3478   FillHistogram("phiRPV0C",fRPV0C,fCentralityV0M);
3479
3480   FillHistogram("phiRPV0Aflat",fRPV0A,fCentralityV0M) ;
3481   FillHistogram("phiRPV0AC",fRPV0A,fRPV0C,fCentralityV0M) ;
3482   FillHistogram("cos2V0AC",TMath::Cos(2.*(fRPV0A-fRPV0C)),fCentralityV0M) ;
3483   if(fHaveTPCRP){
3484     FillHistogram("phiRPV0ATPC",fRP,fRPV0A,fCentralityV0M) ;
3485     FillHistogram("cos2V0ATPC",TMath::Cos(2.*(fRP-fRPV0A)),fCentralityV0M) ;
3486   }
3487
3488   FillHistogram("phiRPV0Cflat",fRPV0C,fCentralityV0M) ;
3489   if(fHaveTPCRP){
3490     FillHistogram("phiRPV0CTPC",fRP,fRPV0C,fCentralityV0M) ;
3491     FillHistogram("cos2V0CTPC",TMath::Cos(2.*(fRP-fRPV0C)),fCentralityV0M) ;
3492   }
3493 }