Bug in flattening correction is fixed
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0Flow.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 #include "TChain.h"
17 #include "TTree.h"
18 #include "TObjArray.h"
19 #include "TF1.h"
20 #include "TFile.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 #include "TH2I.h"
24 #include "TH3F.h"
25 #include "TParticle.h"
26 #include "TCanvas.h"
27 #include "TStyle.h"
28 #include "TRandom.h"
29
30 #include "AliAnalysisManager.h"
31 #include "AliMCEventHandler.h"
32 #include "AliMCEvent.h"
33 #include "AliStack.h"
34 #include "AliAnalysisTaskSE.h"
35 #include "AliAnalysisTaskPi0Flow.h"
36 #include "AliCaloPhoton.h"
37 #include "AliPHOSGeometry.h"
38 #include "TGeoManager.h"
39 #include "AliPHOSEsdCluster.h"
40 #include "AliPHOSCalibData.h"
41 #include "AliESDEvent.h"
42 #include "AliESDCaloCells.h"
43 #include "AliESDVertex.h"
44 #include "AliESDtrackCuts.h"
45 #include "AliAODEvent.h"
46 #include "AliLog.h"
47 #include "AliPID.h"
48 #include "AliCDBManager.h"
49 #include <AliAODCaloCluster.h>
50 #include "AliCentrality.h"
51 #include "AliESDtrackCuts.h"
52 #include "AliEventplane.h"
53 #include "TProfile.h"
54 #include "AliOADBContainer.h"
55 #include "AliPHOSEPFlattener.h"
56
57 // Analysis task to fill histograms with PHOS ESD or AOD clusters and cells
58 // Authors : Dmitri Peressounko
59 // Date    : 28.05.2011
60 // Modified: 03.08.2012 Henrik Qvigstad
61 /* $Id$ */
62
63 ClassImp(AliAnalysisTaskPi0Flow);
64
65 //________________________________________________________________________
66 Double_t rnlin(Double_t *x, Double_t * /*par*/)
67 {
68   //a = par[0], b = par[1].
69   //1+a*exp(-e/b)
70
71 // return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
72  return 1.015*(0.0241+1.0504*x[0]+0.000249*x[0]*x[0]) ;
73
74 }
75
76 //________________________________________________________________________
77 AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
78 : AliAnalysisTaskSE(name),
79   fCentEdges(10),
80   fCentNMixed(10),
81   fNEMRPBins(9),
82   fPeriod(period),
83   fOutputContainer(0x0),
84   fNonLinCorr(0),
85   fEvent(0x0),
86   fEventESD(0x0),
87   fEventAOD(0x0),
88   fMCStack(0x0),
89   fRunNumber(-999),
90   fInternalRunNumber(0),
91   fPHOSGeo(0),
92   fMultV0(0x0),
93   fV0Cpol(0.),fV0Apol(0.),
94   fESDtrackCuts(0x0),
95   fPHOSCalibData(0x0),
96   fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"), 
97   fTPCFlat(0x0),
98   fV0AFlat(0x0),
99   fV0CFlat(0x0),
100   fVertexVector(),
101   fVtxBin(0),
102   fCentralityV0M(0.),
103   fCentBin(0),
104   fHaveTPCRP(0),
105   fRP(0),
106   fRPV0A(0),
107   fRPV0C(0),
108   fEMRPBin(0),
109   fCaloPhotonsPHOS(0x0),
110   fCaloPhotonsPHOSLists(0x0)
111 {
112   const int nbins = 9;
113   Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
114   fCentEdges = TArrayD(nbins+1, edges);
115   Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
116   fCentNMixed = TArrayI(nbins, nMixed);
117   
118   for(Int_t i=0;i<kNCenBins;i++){
119     for(Int_t j=0;j<2; j++)
120       for(Int_t k=0; k<2; k++) {
121         fMeanQ[i][j][k]=0.;
122         fWidthQ[i][j][k]=0.;
123       }
124   }
125   
126   fVertex[0]=0; fVertex[1]=0; fVertex[2]=0; 
127
128   // Output slots #0 write into a TH1 container
129   DefineOutput(1,TList::Class());
130
131
132
133   // Set bad channel map
134   char key[55] ;
135   for(Int_t i=0; i<6; i++){
136     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
137     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
138   }
139
140
141   // Initialize non-linrarity correction
142   fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
143
144
145 }
146 //___________________________________________________________________________
147 AliAnalysisTaskPi0Flow::~AliAnalysisTaskPi0Flow()
148 {
149   delete fNonLinCorr;
150   delete fESDtrackCuts;
151   delete fPHOSCalibData;
152   delete fCaloPhotonsPHOSLists;
153   if(fTPCFlat)delete fTPCFlat;  fTPCFlat=0x0;
154   if(fV0AFlat)delete fV0AFlat;  fV0AFlat=0x0;
155   if(fV0CFlat)delete fV0CFlat;  fV0CFlat=0x0;
156   
157 }
158 //________________________________________________________________________
159 void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
160 {
161   // Create histograms
162   // Called once
163   const Int_t nRuns=200 ;
164
165   // histograms
166   if(fOutputContainer != NULL){
167     delete fOutputContainer;
168   }
169   fOutputContainer = new TList();
170   fOutputContainer->SetOwner(kTRUE);
171
172   //========QA histograms=======
173
174   //Event selection
175   fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 12,0.,13.,nRuns,0.,float(nRuns))) ;
176   fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 12,0.,12.)) ;
177
178   //vertex distribution
179   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
180
181   //Centrality
182   fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
183   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
184   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
185   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
186   fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,10.,100,0.,100.)) ;
187
188
189   //Reaction plane
190   fOutputContainer->Add(new TH3F("hPHOSphi","cos" ,10,0.,100.,20,0.,10.,100,-TMath::Pi(),TMath::Pi()));
191
192   fOutputContainer->Add(new TH2F("cos2AC","RP correlation between TPC subs", 100,-1.,1.,20,0.,100.)) ;
193   fOutputContainer->Add(new TH2F("cos2V0AC","RP correlation between VO A and C sides", 100,-1.,1.,20,0.,100.)) ;
194   fOutputContainer->Add(new TH2F("cos2V0ATPC","RP correlation between TPC and V0A", 100,-1.,1.,20,0.,100.)) ;
195   fOutputContainer->Add(new TH2F("cos2V0CTPC","RP correlation between TPC and V0C", 100,-1.,1.,20,0.,100.)) ;
196
197   fOutputContainer->Add(new TH2F("phiRP","RP distribution with TPC", 100,0.,TMath::Pi(),20,0.,100.)) ;
198   fOutputContainer->Add(new TH2F("phiRPflat","RP distribution with TPC flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
199   fOutputContainer->Add(new TH2F("phiRPV0A","RP distribution with V0A", 100,0.,TMath::Pi(),20,0.,100.)) ;
200   fOutputContainer->Add(new TH2F("phiRPV0C","RP distribution with V0C", 100,0.,TMath::Pi(),20,0.,100.)) ;
201   fOutputContainer->Add(new TH3F("phiRPV0AC","RP distribution with V0A and V0C", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
202   fOutputContainer->Add(new TH2F("phiRPV0Aflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
203   fOutputContainer->Add(new TH2F("phiRPV0Cflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
204   fOutputContainer->Add(new TH3F("phiRPV0ATPC","RP distribution with V0A + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
205   fOutputContainer->Add(new TH3F("phiRPV0CTPC","RP distribution with V0C + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
206
207
208   //PHOS QA
209   fOutputContainer->Add(new TH1I("hCellMultEvent"  ,"PHOS cell multiplicity per event"    ,2000,0,2000));
210   fOutputContainer->Add(new TH1I("hCellMultEventM1","PHOS cell multiplicity per event, M1",2000,0,2000));
211   fOutputContainer->Add(new TH1I("hCellMultEventM2","PHOS cell multiplicity per event, M2",2000,0,2000));
212   fOutputContainer->Add(new TH1I("hCellMultEventM3","PHOS cell multiplicity per event, M3",2000,0,2000));
213
214   fOutputContainer->Add(new TH1F("hCellEnergy"  ,"Cell energy"            ,3000,0.,30.));
215   fOutputContainer->Add(new TH1F("hCellEnergyM1","Cell energy in module 1",3000,0.,30.));
216   fOutputContainer->Add(new TH1F("hCellEnergyM2","Cell energy in module 2",3000,0.,30.));
217   fOutputContainer->Add(new TH1F("hCellEnergyM3","Cell energy in module 3",3000,0.,30.));
218
219   fOutputContainer->Add(new TH2F("hCellNXZM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
220   fOutputContainer->Add(new TH2F("hCellNXZM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
221   fOutputContainer->Add(new TH2F("hCellNXZM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
222   fOutputContainer->Add(new TH2F("hCellEXZM1","Cell E(X,Z), M1",64,0.5,64.5, 56,0.5,56.5));
223   fOutputContainer->Add(new TH2F("hCellEXZM2","Cell E(X,Z), M2",64,0.5,64.5, 56,0.5,56.5));
224   fOutputContainer->Add(new TH2F("hCellEXZM3","Cell E(X,Z), M3",64,0.5,64.5, 56,0.5,56.5));
225
226   //Bad Map
227   fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
228   fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
229   fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
230
231   fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
232   fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
233   fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
234
235   fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
236   fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
237   fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
238
239   fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
240   fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
241   fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
242
243
244   //Single photon and pi0 spectrum
245   const Int_t nPtPhot = 300 ;
246   const Double_t ptPhotMax = 30 ;
247   const Int_t nM       = 500;
248   const Double_t mMin  = 0.0;
249   const Double_t mMax  = 1.0;
250
251   //PHOS calibration QA
252   fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
253   fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
254   fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
255   fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
256   fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
257   fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
258
259   // Histograms for different centralities
260   char key[55] ;
261   for(Int_t cent=0; cent < fCentEdges.GetSize(); cent++){
262     snprintf(key,55,"hPhotAll_cen%d",cent) ;
263     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
264     snprintf(key,55,"hPhotAllcore_cen%d",cent) ;
265     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
266     snprintf(key,55,"hPhotAllwou_cen%d",cent) ;
267     fOutputContainer->Add(new TH1F(key,"All clusters",nPtPhot,0.,ptPhotMax));
268     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
269     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
270     snprintf(key,55,"hPhotDisp2_cen%d",cent) ;
271     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
272     snprintf(key,55,"hPhotDispcore_cen%d",cent) ;
273     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
274     snprintf(key,55,"hPhotDispwou_cen%d",cent) ;
275     fOutputContainer->Add(new TH1F(key,"Disp clusters",nPtPhot,0.,ptPhotMax));
276     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
277     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
278     snprintf(key,55,"hPhotCPVcore_cen%d",cent) ;
279     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
280     snprintf(key,55,"hPhotCPV2_cen%d",cent) ;
281     fOutputContainer->Add(new TH1F(key,"CPV clusters",nPtPhot,0.,ptPhotMax));
282     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
283     fOutputContainer->Add(new TH1F(key,"Both clusters",nPtPhot,0.,ptPhotMax));
284     snprintf(key,55,"hPhotBothcore_cen%d",cent) ;
285     fOutputContainer->Add(new TH1F(key,"Both clusters",nPtPhot,0.,ptPhotMax));
286
287     snprintf(key,55,"hPi0All_cen%d",cent) ;
288     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
289     snprintf(key,55,"hPi0Allcore_cen%d",cent) ;
290     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
291     snprintf(key,55,"hPi0Allwou_cen%d",cent) ;
292     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
293     snprintf(key,55,"hPi0Disp_cen%d",cent) ;
294     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
295     snprintf(key,55,"hPi0Disp2_cen%d",cent) ;
296     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
297     snprintf(key,55,"hPi0Dispcore_cen%d",cent) ;
298     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
299     snprintf(key,55,"hPi0Dispwou_cen%d",cent) ;
300     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
301     snprintf(key,55,"hPi0CPV_cen%d",cent) ;
302     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
303     snprintf(key,55,"hPi0CPVcore_cen%d",cent) ;
304     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
305     snprintf(key,55,"hPi0CPV2_cen%d",cent) ;
306     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
307     snprintf(key,55,"hPi0Both_cen%d",cent) ;
308     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
309     snprintf(key,55,"hPi0Bothcore_cen%d",cent) ;
310     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
311
312     snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
313     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
314     snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
315     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
316     snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
317     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
318     snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
319     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
320     snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
321     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
322
323     snprintf(key,55,"hSingleAll_cen%d",cent) ;
324     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
325     snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
326     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
327     snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
328     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
329     snprintf(key,55,"hSingleDisp_cen%d",cent) ;
330     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
331     snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
332     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
333     snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
334     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
335     snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
336     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
337     snprintf(key,55,"hSingleCPV_cen%d",cent) ;
338     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
339     snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
340     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
341     snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
342     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
343     snprintf(key,55,"hSingleBoth_cen%d",cent) ;
344     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
345     snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
346     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
347
348
349     snprintf(key,55,"hMiPi0All_cen%d",cent) ;
350     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
351     snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
352     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
353     snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
354     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
355     snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
356     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
357     snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
358     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
359     snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
360     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
361     snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
362     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
363     snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
364     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
365     snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
366     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
367     snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
368     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
369     snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
370     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
371     snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
372     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
373
374     snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
375     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
376     snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
377     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
378     snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
379     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
380     snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
381     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
382     snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
383     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
384
385     snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
386     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
387     snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
388     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
389     snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
390     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
391     snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
392     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
393     snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
394     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
395     snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
396     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
397     snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
398     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
399     snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
400     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
401     snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
402     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
403     snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
404     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
405     snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
406     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
407     snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
408     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
409   }
410
411
412
413   const Int_t nPt      = 20;
414   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.} ;
415   const Int_t nPhi=10 ;
416   Double_t xPhi[nPhi+1] ;
417   for(Int_t i=0;i<=nPhi;i++)
418     xPhi[i]=i*TMath::Pi() /nPhi ;
419   const Int_t nMm=200 ;
420   Double_t xM[nMm+1] ;
421   for(Int_t i=0;i<=nMm;i++)
422     xM[i]=i*0.5 /nMm;
423
424   char phiTitle[15] ;
425   for(Int_t iRP=0; iRP<3; iRP++){
426     if(iRP==0)
427       snprintf(phiTitle,15,"TPC") ;
428     if(iRP==1)
429       snprintf(phiTitle,15,"V0A") ;
430     if(iRP==2)
431       snprintf(phiTitle,15,"V0C") ;
432     for(Int_t cent=0; cent<fCentEdges.GetSize(); cent++){
433       snprintf(key,55,"hPhotPhi%sAll_cen%d",phiTitle,cent) ;
434       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
435       snprintf(key,55,"hPhotPhi%sAllcore_cen%d",phiTitle,cent) ;
436       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
437       snprintf(key,55,"hPhotPhi%sDisp_cen%d",phiTitle,cent) ;
438       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
439       snprintf(key,55,"hPhotPhi%sDispcore_cen%d",phiTitle,cent) ;
440       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
441       snprintf(key,55,"hPhotPhi%sCPV_cen%d",phiTitle,cent) ;
442       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
443       snprintf(key,55,"hPhotPhi%sCPVcore_cen%d",phiTitle,cent) ;
444       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
445       snprintf(key,55,"hPhotPhi%sBoth_cen%d",phiTitle,cent) ;
446       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
447       snprintf(key,55,"hPhotPhi%sBothcore_cen%d",phiTitle,cent) ;
448       fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPt,xPt,nPhi,xPhi));
449
450       //Pions
451       snprintf(key,55,"hMassPt%sAll_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,"hMassPt%sAllcore_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,"hMassPt%sCPV_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,"hMassPt%sCPVcore_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,"hMassPt%sDisp_cen%d",phiTitle,cent) ;
460       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
461       snprintf(key,55,"hMassPt%sDispcore_cen%d",phiTitle,cent) ;
462       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
463       snprintf(key,55,"hMassPt%sBoth_cen%d",phiTitle,cent) ;
464       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
465       snprintf(key,55,"hMassPt%sBothcore_cen%d",phiTitle,cent) ;
466       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
467
468       //Mixed
469       snprintf(key,55,"hMiMassPt%sAll_cen%d",phiTitle,cent) ;
470       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
471       snprintf(key,55,"hMiMassPt%sAllcore_cen%d",phiTitle,cent) ;
472       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
473       snprintf(key,55,"hMiMassPt%sCPV_cen%d",phiTitle,cent) ;
474       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
475       snprintf(key,55,"hMiMassPt%sCPVcore_cen%d",phiTitle,cent) ;
476       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
477       snprintf(key,55,"hMiMassPt%sDisp_cen%d",phiTitle,cent) ;
478       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
479       snprintf(key,55,"hMiMassPt%sDispcore_cen%d",phiTitle,cent) ;
480       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
481       snprintf(key,55,"hMiMassPt%sBoth_cen%d",phiTitle,cent) ;
482       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
483       snprintf(key,55,"hMiMassPt%sBothcore_cen%d",phiTitle,cent) ;
484       fOutputContainer->Add(new TH3F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nMm,xM,nPt,xPt,nPhi,xPhi));
485     }
486   }
487   
488   // Setup photon lists
489   Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
490   fCaloPhotonsPHOSLists = new TObjArray(kapacity);
491   fCaloPhotonsPHOSLists->SetOwner();
492   
493   PostData(1, fOutputContainer);
494 }
495
496 //________________________________________________________________________
497 void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
498 {
499   // Main loop, called for each event
500   // Analyze ESD/AOD
501
502
503   // Step 0: Event Objects
504   fEvent = GetEvent();
505   fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
506   fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
507   fMCStack = GetMCStack();
508   LogProgress(0);
509
510
511   // Step 1: Run Number, Misalignment Matrix, and Calibration
512   // fRunNumber, fInternalRunNumber, fMultV0, fV0Cpol, fV0Apol, fMeanQ, fWidthQ
513   if( fRunNumber != fEvent->GetRunNumber()) { // Check run number
514     // this should run only at first call of UserExec(),
515     // or if task runs over multiple runs, which should not occur in normal use.
516
517     // if run number has changed, set run variables
518     fRunNumber = fEvent->GetRunNumber();
519     fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
520     // then set misalignment and V0 calibration
521     SetGeometry();
522     SetMisalignment();
523     SetV0Calibration();
524     SetESDTrackCuts();
525     SetPHOSCalibData();
526     SetFlatteningData();
527   }
528   LogProgress(1);
529   LogSelection(0, fInternalRunNumber);
530
531
532   // Step 2: Vertex
533   // fVertex, fVertexVector, fVtxBin
534   SetVertex();
535   if( RejectEventVertex() ) {
536     PostData(1, fOutputContainer);
537     return; // Reject!
538   }
539   LogProgress(2);
540
541 // Step 3:
542 //   if(event->IsPileupFromSPD()){
543 //     PostData(1, fOutputContainer);
544 //     return; // Reject!
545 //   }
546   LogProgress(3);
547
548
549   // Step 4: Centrality
550   // fCentralityV0M, fCentBin
551   SetCentrality();
552   if( RejectCentrality() ){
553     PostData(1, fOutputContainer);
554     return; // Reject! 
555   }
556   LogProgress(4);
557
558
559   // Step 5: Reaction Plane
560   // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
561   EvalReactionPlane(); //TODO: uncomment this, or at least deal with it
562   EvalV0ReactionPlane(); //TODO: uncomment this, or at least deal with it
563   fEMRPBin = GetRPBin(); //TODO: uncomment this, or at least deal with it
564   LogProgress(5);
565
566
567   // Step 6: MC
568   //  ProcessMC() ;
569   LogProgress(6);
570
571
572   // Step 7: QA PHOS cells
573   FillPHOSCellQAHists();
574   LogProgress(7);
575
576
577   // Step 8: Event Photons (PHOS Clusters) selection
578   SelectPhotonClusters();
579   FillSelectedClusterHistograms();
580   LogProgress(8);
581
582   // Step 9: Consider pi0 (photon/cluster) pairs.
583   ConsiderPi0s();
584   LogProgress(9);
585
586   // Step 10; Mixing
587   ConsiderPi0sMix();
588   LogProgress(10);
589   
590   // Step 11: Update lists
591   UpdateLists();
592   LogProgress(11);
593
594   
595   // Post output data.
596   PostData(1, fOutputContainer);
597 }
598
599 //________________________________________________________________________
600 // void AliAnalysisTaskPi0Flow::Terminate(Option_t *)
601 // {
602 //   // Draw result to the screen
603 //   // Called once at the end of the query
604 //   // new TCanvas;
605 //   // TH1 * hTotSelEvents = dynamic_cast<TH1*>(fOutputContainer->FindObject("hTotSelEvents"));
606 //   // hTotSelEvents->Draw();
607 // }
608 //________________________________________________________________________
609 void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges)
610 {
611   // Define centrality bins by their edges
612
613   int last = edges.GetSize()-1;
614   if( edges.At(0) < 0.) 
615     AliError("lower edge less then 0");
616   if( 90. < edges.At(last)  )
617     AliError("upper edge larger then 90.");
618   for(int i=0; i<last-1; ++i)
619     if(edges.At(i) > edges.At(i+1))
620       AliError("edges are not sorted");
621   
622   fCentEdges = edges;
623 }
624 //________________________________________________________________________
625 void AliAnalysisTaskPi0Flow::SetNMixedPerCentrality(const TArrayI& nMixed)
626 {
627   // Set number of mixed events for all centrality bins
628
629   fCentNMixed = nMixed;
630 }
631
632 //________________________________________________________________________
633 void AliAnalysisTaskPi0Flow::SetPHOSBadMap(Int_t mod, TH2I* badMapHist)
634   {
635     if(fPHOSBadMap[mod])
636       delete fPHOSBadMap[mod];
637
638     fPHOSBadMap[mod]=new TH2I(*badMapHist);
639     if(fDebug)
640       AliInfo(Form("Setting Bad Map Histogram  %s",fPHOSBadMap[mod]->GetName()));
641   }
642
643 //________________________________________________________________________
644 Bool_t AliAnalysisTaskPi0Flow::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
645 {
646   //Check if this channel belogs to the good ones
647
648   if(strcmp(det,"PHOS")==0){
649     if(mod>5 || mod<1){
650       AliError(Form("No bad map for PHOS module %d",mod)) ;
651       return kTRUE ;
652     }
653     if(!fPHOSBadMap[mod]){
654       AliError(Form("No Bad map for PHOS module %d, !fPHOSBadMap[mod]",mod)) ;
655       return kTRUE ;
656     }
657     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
658       return kFALSE ;
659     else
660       return kTRUE ;
661   }
662   else{
663     AliError(Form("Can not find bad channels for detector %s ",det)) ;
664   }
665
666   //Remove 6 noisy channels in run 139036, LHC10h
667   if( 139036 == fRunNumber
668     && mod==1
669     && (ix==9||ix==10||ix==11)
670     && (iz==45 || iz==46))
671     return kFALSE;
672
673   return kTRUE ;
674 }
675 //_____________________________________________________________________________
676 void AliAnalysisTaskPi0Flow::FillPHOSCellQAHists()
677 {
678   // Fill cell occupancy per module
679
680   AliVCaloCells * cells = fEvent->GetPHOSCells();
681
682   FillHistogram("hCenPHOSCells",fCentralityV0M,cells->GetNumberOfCells()) ;
683   FillHistogram("hCenTrack",fCentralityV0M,fEvent->GetNumberOfTracks()) ;
684
685
686   Int_t nCellModule[3] = {0,0,0};
687   for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
688     Int_t cellAbsId = cells->GetCellNumber(iCell);
689     Int_t relId[4] = {0,0,0,0};
690     fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
691     Int_t mod1  = relId[0];
692     Int_t cellX = relId[2];
693     Int_t cellZ = relId[3] ;
694     Float_t energy = cells->GetAmplitude(iCell);
695     FillHistogram("hCellEnergy",energy);
696     if(mod1==1) {
697       nCellModule[0]++;
698       FillHistogram("hCellEnergyM1",cells->GetAmplitude(iCell));
699       FillHistogram("hCellNXZM1",cellX,cellZ,1.);
700       FillHistogram("hCellEXZM1",cellX,cellZ,energy);
701     }
702     else if (mod1==2) {
703       nCellModule[1]++;
704       FillHistogram("hCellEnergyM2",cells->GetAmplitude(iCell));
705       FillHistogram("hCellNXZM2",cellX,cellZ,1.);
706       FillHistogram("hCellEXZM2",cellX,cellZ,energy);
707     }
708     else if (mod1==3) {
709       nCellModule[2]++;
710       FillHistogram("hCellEnergyM3",cells->GetAmplitude(iCell));
711       FillHistogram("hCellNXZM3",cellX,cellZ,1.);
712       FillHistogram("hCellEXZM3",cellX,cellZ,energy);
713     }
714   }
715   FillHistogram("hCellMultEventM1",nCellModule[0]);
716   FillHistogram("hCellMultEventM2",nCellModule[1]);
717   FillHistogram("hCellMultEventM3",nCellModule[2]);
718
719 }
720 //_____________________________________________________________________________
721 void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
722 {
723   // clear (or create) array for holding events photons/clusters
724   if(fCaloPhotonsPHOS)
725     fCaloPhotonsPHOS->Clear();
726   else {
727     fCaloPhotonsPHOS = new TObjArray(200);
728     fCaloPhotonsPHOS->SetOwner();
729   }
730
731   
732   AliVCaloCells* cells = dynamic_cast<AliVCaloCells*> (fEvent->GetPHOSCells());
733   for (Int_t i=0; i<fEvent->GetNumberOfCaloClusters(); i++) {
734     AliVCluster *clu = fEvent->GetCaloCluster(i);
735
736     if ( !clu->IsPHOS() || clu->E()< kMinClusterEnergy) continue; // reject cluster
737
738     // check if cell/channel is good.
739     Float_t  position[3];
740     clu->GetPosition(position);
741     TVector3 global(position) ;
742     Int_t relId[4] ;
743     fPHOSGeo->GlobalPos2RelId(global,relId) ;
744     Int_t mod  = relId[0] ;
745     Int_t cellX = relId[2];
746     Int_t cellZ = relId[3] ;
747     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
748       continue ; // reject if not.
749
750
751     FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
752
753     if(clu->GetNCells() < kMinNCells) continue ;
754     if(clu->GetM02() < kMinM02)   continue ;
755
756     TLorentzVector lorentzMomentum;
757     Double_t ecore;
758     ecore = CoreEnergy(clu,cells);
759
760     //if ESD, Apply re-Calibreation
761     Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0)
762     if( fEventESD ) {
763       AliPHOSEsdCluster cluPHOS1( *(AliESDCaloCluster*) (clu) );
764       cluPHOS1.Recalibrate(fPHOSCalibData, dynamic_cast<AliESDCaloCells*> (cells)); // modify the cell energies
765       Reclusterize(&cluPHOS1) ;
766       cluPHOS1.EvalAll(kLogWeight, fVertexVector);         // recalculate the cluster parameters
767       cluPHOS1.SetE(fNonLinCorr->Eval(cluPHOS1.E()));// Users's nonlinearity
768
769       if(cluPHOS1.E()<0.3) continue; // check energy again
770
771       //correct misalignment
772       TVector3 localPos;
773       const Float_t shiftX[6]={0.,-2.3,-2.11,-1.53,0.,0.} ;
774       const Float_t shiftZ[6]={0.,-0.4, 0.52, 0.8,0.,0.} ;
775       fPHOSGeo->Global2Local(localPos,global,mod) ;
776       fPHOSGeo->Local2Global(mod,localPos.X()+shiftX[mod],localPos.Z()+shiftZ[mod],global);
777       position[0]=global.X() ;
778       position[1]=global.Y() ;
779       position[2]=global.Z() ;
780       cluPHOS1.SetPosition(position);
781
782       cluPHOS1.GetMomentum(lorentzMomentum ,origo);
783       
784       //TODO: Check, this may be LHC10h specific:
785       if(mod==2) lorentzMomentum*=135.5/134.0 ;
786       if(mod==3) lorentzMomentum*=135.5/137.2 ;
787       if(mod==2) ecore*=135.5/134.0 ;
788       if(mod==3) ecore*=135.5/137.2 ;
789         
790     }
791     else if (fEventAOD) { // is ! ESD, AOD.
792       AliESDCaloCluster* aodCluster = (AliESDCaloCluster*) (clu);
793       aodCluster->GetMomentum(lorentzMomentum ,origo);
794     }
795     else {
796       AliError("(Calo)Cluster is neither ESD nor AOD");
797       continue;
798     }
799
800
801     char skey[55];
802     snprintf(skey,55,"hCluLowM%d",mod) ;
803     FillHistogram(skey,cellX,cellZ,1.);
804     if(lorentzMomentum.E()>1.5){
805       sprintf(skey,"hCluHighM%d",mod) ;
806       FillHistogram(skey,cellX,cellZ,1.);
807     }
808
809     fCaloPhotonsPHOS->Add(new  AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
810     AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
811
812     ph->SetModule(mod) ;
813     lorentzMomentum*= ecore/lorentzMomentum.E() ;
814     ph->SetMomV2(&lorentzMomentum) ;
815     ph->SetNCells(clu->GetNCells());
816     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
817     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
818     if(ph->IsDispOK()){
819       sprintf(skey,"hCluDispM%d",mod) ;
820       FillHistogram(skey,cellX,cellZ,1.);
821     }
822
823     // Track Matching
824     Double_t dx=clu->GetTrackDx() ;
825     Double_t dz=clu->GetTrackDz() ;
826     Bool_t cpvBit=kTRUE ; //No track matched by default
827     Bool_t cpvBit2=kTRUE ; //More Strict criterion
828     if( fEventESD ) {
829       
830       TArrayI * itracks = static_cast<AliESDCaloCluster*> (clu)->GetTracksMatched() ;
831       if(itracks->GetSize()>0){
832         Int_t iTr = itracks->At(0);
833         if(iTr>=0 && iTr<fEvent->GetNumberOfTracks()){
834           AliVParticle* track = fEvent->GetTrack(iTr);
835           Double_t pt = track->Pt() ;
836           Short_t charge = track->Charge() ;
837           Double_t r=TestCPV(dx, dz, pt, charge) ;
838           cpvBit=(r>2.) ;
839           cpvBit2=(r>4.) ;
840         }
841       }
842     }
843     else if ( fEventAOD ) {
844       int nTracksMatched = clu->GetNTracksMatched();
845       if(nTracksMatched > 0) {
846         AliVTrack* track = dynamic_cast<AliVTrack*> (clu->GetTrackMatched(0));
847         if ( track ) {
848           Double_t pt = track->Pt();
849           Short_t charge = track->Charge();
850           Double_t r = TestCPV(dx, dz, pt, charge) ;
851           cpvBit=(r>2.) ;
852           cpvBit2=(r>4.) ;
853         }
854       }
855     }
856     ph->SetCPVBit(cpvBit) ;
857     ph->SetCPV2Bit(cpvBit2) ;
858     if(cpvBit){
859       sprintf(skey,"hCluVetoM%d",mod) ;
860       FillHistogram(skey,cellX,cellZ,1.);
861     }
862     ph->SetEMCx(float(cellX)) ;
863     ph->SetEMCz(float(cellZ)) ;
864     //    ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
865     ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolde
866
867   }
868   FillHistogram("hCenPHOS",fCentralityV0M, fCaloPhotonsPHOS->GetEntriesFast()) ;
869 }
870 //_____________________________________________________________________________
871 void AliAnalysisTaskPi0Flow::FillSelectedClusterHistograms()
872 {
873   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
874     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
875
876     Double_t dphiA=ph1->Phi()-fRPV0A ;
877     while(dphiA<0)dphiA+=TMath::Pi() ;
878     while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
879
880     Double_t dphiC=ph1->Phi()-fRPV0C ;
881     while(dphiC<0)dphiC+=TMath::Pi() ;
882     while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
883
884     Double_t dphiT=ph1->Phi()-fRP ;
885     while(dphiT<0)dphiT+=TMath::Pi() ;
886     while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
887
888     FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),ph1->Pt(),dphiA) ;
889     FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),ph1->Pt(),dphiC) ;
890     if(fHaveTPCRP)
891       FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),ph1->Pt(),dphiT) ;
892
893     FillHistogram(Form("hPhotAll_cen%d",fCentBin),ph1->Pt()) ;
894     FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
895     if(ph1->IsntUnfolded()){
896       FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),ph1->Pt()) ;
897     }
898     if(ph1->IsCPVOK()){
899       FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),ph1->Pt(),dphiA) ;
900       FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),ph1->Pt(),dphiC) ;
901       if(fHaveTPCRP)
902         FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),ph1->Pt(),dphiT) ;
903
904       FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
905       FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
906       if(fHaveTPCRP)
907         FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
908
909       FillHistogram(Form("hPhotCPV_cen%d",fCentBin),ph1->Pt()) ;
910       FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
911     }
912     if(ph1->IsCPV2OK()){
913       FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),ph1->Pt()) ;
914     }
915     if(ph1->IsDispOK()){
916       FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),ph1->Pt(),dphiA) ;
917       FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),ph1->Pt(),dphiC) ;
918       if(fHaveTPCRP)
919         FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),ph1->Pt(),dphiT) ;
920
921       FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
922       FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
923       if(fHaveTPCRP)
924         FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
925
926       FillHistogram(Form("hPhotDisp_cen%d",fCentBin),ph1->Pt()) ;
927       if(ph1->IsDisp2OK()){
928         FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),ph1->Pt()) ;
929       }
930       FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
931       if(ph1->IsntUnfolded()){
932         FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),ph1->Pt()) ;
933       }
934       if(ph1->IsCPVOK()){
935         FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),ph1->Pt(),dphiA) ;
936         FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),ph1->Pt(),dphiC) ;
937         if(fHaveTPCRP)
938           FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),ph1->Pt(),dphiT) ;
939
940         FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiA) ;
941         FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiC) ;
942         if(fHaveTPCRP)
943           FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt(),dphiT) ;
944
945         FillHistogram(Form("hPhotBoth_cen%d",fCentBin),ph1->Pt()) ;
946         FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ph1->GetMomV2()->Pt()) ;
947       }
948     }
949   }
950 }
951 //_____________________________________________________________________________
952 void AliAnalysisTaskPi0Flow::ConsiderPi0s()
953 {
954   char key[55];
955   for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
956     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
957     for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
958       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
959       TLorentzVector p12  = *ph1  + *ph2;
960       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
961       FillHistogram("hPHOSphi",fCentralityV0M,p12.Pt(),p12.Phi());
962       Double_t dphiA=p12.Phi()-fRPV0A ;
963       while(dphiA<0)dphiA+=TMath::Pi() ;
964       while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
965
966       Double_t dphiC=p12.Phi()-fRPV0C ;
967       while(dphiC<0)dphiC+=TMath::Pi() ;
968       while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
969
970       Double_t dphiT=p12.Phi()-fRP ;
971       while(dphiT<0)dphiT+=TMath::Pi() ;
972       while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
973
974       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
975
976       FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
977       FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
978       if(fHaveTPCRP)
979         FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
980
981       FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
982       FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
983       if(fHaveTPCRP)
984         FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
985
986
987       FillHistogram(Form("hPi0All_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
988       FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
989       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
990         FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
991       }
992
993       FillHistogram(Form("hSingleAll_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
994       FillHistogram(Form("hSingleAll_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
995       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
996       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
997       if(ph1->IsntUnfolded())
998         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
999       if(ph2->IsntUnfolded())
1000         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1001       if(ph1->IsCPVOK()){
1002         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1003         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1004       }
1005       if(ph2->IsCPVOK()){
1006         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1007         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1008       }
1009       if(ph1->IsCPV2OK()){
1010         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1011       }
1012       if(ph2->IsCPV2OK()){
1013         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1014       }
1015       if(ph1->IsDispOK()){
1016         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1017         if(ph1->IsntUnfolded()){
1018           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1019         }
1020         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1021       }
1022       if(ph2->IsDispOK()){
1023         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1024         if(ph1->IsntUnfolded()){
1025           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1026         }
1027         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1028       }
1029       if(ph1->IsDisp2OK()){
1030         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1031       }
1032       if(ph2->IsDisp2OK()){
1033         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1034       }
1035       if(ph1->IsDispOK() && ph1->IsCPVOK()){
1036         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1037         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1038       }
1039       if(ph2->IsDispOK() && ph2->IsCPVOK()){
1040         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1041         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1042       }
1043
1044
1045       if(a<kAlphaCut){
1046         FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1047       }
1048
1049       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1050         snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
1051         FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1052         FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1053         if(fHaveTPCRP)
1054           FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1055
1056         FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
1057         FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
1058         if(fHaveTPCRP)
1059           FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
1060
1061         FillHistogram(Form("hPi0CPV_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1062         FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1063
1064         if(a<kAlphaCut){
1065           FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1066         }
1067       }
1068       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1069         FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),p12.M(),p12.Pt()) ;
1070         if(a<kAlphaCut){
1071           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1072         }
1073       }
1074       if(ph1->IsDispOK() && ph2->IsDispOK()){
1075         snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
1076         FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1077         FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1078         if(fHaveTPCRP)
1079           FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1080
1081         FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiA) ;
1082         FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiC) ;
1083         if(fHaveTPCRP)
1084           FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiT) ;
1085
1086         FillHistogram(Form("hPi0Disp_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1087         FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1088         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1089           FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1090         }
1091         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1092           FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),p12.M(), p12.Pt()) ;
1093         }
1094
1095         if(a<kAlphaCut){
1096           FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1097         }
1098
1099         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1100           FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1101           FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1102           if(fHaveTPCRP)
1103             FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1104
1105           FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiA) ;
1106           FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiC) ;
1107           if(fHaveTPCRP)
1108             FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt(),dphiT) ;
1109
1110           FillHistogram(Form("hPi0Both_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1111           FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
1112
1113           if(a<kAlphaCut){
1114             snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
1115             FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1116           }
1117           if(fCentralityV0M>20.){
1118           if(ph1->Module()==1 && ph2->Module()==1)
1119             FillHistogram("hPi0M11",p12.M(),p12.Pt() );
1120           else if(ph1->Module()==2 && ph2->Module()==2)
1121             FillHistogram("hPi0M22",p12.M(),p12.Pt() );
1122           else if(ph1->Module()==3 && ph2->Module()==3)
1123             FillHistogram("hPi0M33",p12.M(),p12.Pt() );
1124           else if(ph1->Module()==1 && ph2->Module()==2)
1125             FillHistogram("hPi0M12",p12.M(),p12.Pt() );
1126           else if(ph1->Module()==1 && ph2->Module()==3)
1127             FillHistogram("hPi0M13",p12.M(),p12.Pt() );
1128           else if(ph1->Module()==2 && ph2->Module()==3)
1129             FillHistogram("hPi0M23",p12.M(),p12.Pt() );
1130           }
1131
1132         }
1133       }
1134     } // end of loop i2
1135   } // end of loop i1
1136 }
1137 //_____________________________________________________________________________
1138 void AliAnalysisTaskPi0Flow::ConsiderPi0sMix()
1139 {
1140   char key[55];
1141
1142   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1143
1144   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
1145     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1146     for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1147       TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
1148       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1149         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1150         TLorentzVector p12  = *ph1  + *ph2;
1151         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1152
1153         Double_t dphiA=p12.Phi()-fRPV0A ;
1154         while(dphiA<0)dphiA+=TMath::Pi() ;
1155         while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
1156
1157         Double_t dphiC=p12.Phi()-fRPV0C ;
1158         while(dphiC<0)dphiC+=TMath::Pi() ;
1159         while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
1160
1161         Double_t dphiT=p12.Phi()-fRP ;
1162         while(dphiT<0)dphiT+=TMath::Pi() ;
1163         while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
1164
1165
1166         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1167
1168         snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
1169         FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiA) ;
1170         FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiC) ;
1171         if(fHaveTPCRP)
1172           FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),p12.M() ,p12.Pt(),dphiT) ;
1173
1174         FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiA) ;
1175         FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiC) ;
1176         if(fHaveTPCRP)
1177           FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),pv12.M(), pv12.Pt(), dphiT) ;
1178
1179         FillHistogram(Form("hMiPi0All_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1180         FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),pv12.M() ,pv12.Pt()) ;
1181         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1182           FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1183         }
1184
1185         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1186         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1187         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1188         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1189         if(ph1->IsntUnfolded())
1190           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1191         if(ph2->IsntUnfolded())
1192           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1193         if(ph1->IsCPVOK()){
1194           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1195           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1196         }
1197         if(ph2->IsCPVOK()){
1198           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1199           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1200         }
1201         if(ph1->IsCPV2OK()){
1202           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1203         }
1204         if(ph2->IsCPV2OK()){
1205           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1206         }
1207         if(ph1->IsDispOK()){
1208           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1209           if(ph1->IsntUnfolded()){
1210             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1211           }
1212           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),pv12.M(),ph1->GetMomV2()->Pt()) ;
1213         }
1214         if(ph2->IsDispOK()){
1215           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1216           if(ph1->IsntUnfolded()){
1217             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1218           }
1219           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),pv12.M(),ph2->GetMomV2()->Pt()) ;
1220         }
1221         if(ph1->IsDisp2OK()){
1222           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),p12.M(),ph1->Pt()) ;
1223         }
1224         if(ph2->IsDisp2OK()){
1225           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),p12.M(),ph2->Pt()) ;
1226         }
1227         if(ph1->IsDispOK() && ph1->IsCPVOK()){
1228           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1229           FillHistogram(key,p12.M(),ph1->Pt()) ;
1230           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1231           FillHistogram(key,pv12.M(),ph1->GetMomV2()->Pt()) ;
1232         }
1233         if(ph2->IsDispOK() && ph2->IsCPVOK()){
1234           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1235           FillHistogram(key,p12.M(),ph2->Pt()) ;
1236           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1237           FillHistogram(key,pv12.M(),ph2->GetMomV2()->Pt()) ;
1238         }
1239
1240
1241
1242         if(a<kAlphaCut){
1243           snprintf(key,55,"hMiPi0All_a07_cen%d",fCentBin) ;
1244           FillHistogram(key,p12.M() ,p12.Pt()) ;
1245         }
1246         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1247           FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiA) ;
1248           FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiC) ;
1249           if(fHaveTPCRP)
1250             FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),p12.M(), p12.Pt(),dphiT) ;
1251
1252           FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiA) ;
1253           FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiC) ;
1254           if(fHaveTPCRP)
1255             FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt(),dphiT) ;
1256
1257           FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),p12.M(), p12.Pt()) ;
1258           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),pv12.M(), pv12.Pt()) ;
1259
1260           if(a<kAlphaCut){
1261             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1262           }
1263         }
1264         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1265           FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1266
1267           if(a<kAlphaCut){
1268             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1269           }
1270         }
1271         if(ph1->IsDispOK() && ph2->IsDispOK()){
1272           FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiA) ;
1273           FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiC) ;
1274           if(fHaveTPCRP)
1275             FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),p12.M(),p12.Pt(),dphiT) ;
1276
1277           FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1278           FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1279           if(fHaveTPCRP)
1280             FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1281
1282
1283           FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1284           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1285           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1286             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),p12.M(),p12.Pt()) ;
1287           }
1288
1289           if(a<kAlphaCut){
1290             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1291           }
1292           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1293             FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiA) ;
1294             FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiC) ;
1295             if(fHaveTPCRP)
1296               FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),p12.M(),p12.Pt(),dphiT) ;
1297
1298             FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1299             FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1300             if(fHaveTPCRP)
1301               FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1302
1303             FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1304             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1305
1306             if(a<kAlphaCut){
1307               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),p12.M() ,p12.Pt()) ;
1308             }
1309           }
1310         }
1311       } // end of loop i2
1312     }
1313   } // end of loop i1
1314 }
1315 //_____________________________________________________________________________
1316 void AliAnalysisTaskPi0Flow::UpdateLists()
1317 {
1318   //Now we either add current events to stack or remove
1319   //If no photons in current event - no need to add it to mixed
1320
1321   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1322   if( fDebug >= 2 )
1323     AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1324   if(fCaloPhotonsPHOS->GetEntriesFast()>0){
1325     arrayList->AddFirst(fCaloPhotonsPHOS) ;
1326     fCaloPhotonsPHOS=0;
1327     if(arrayList->GetEntries() > fCentNMixed[fCentBin]){ // Remove redundant events
1328       TObjArray * tmp = static_cast<TObjArray*>(arrayList->Last()) ;
1329       arrayList->RemoveLast() ;
1330       delete tmp ; // TODO: may conflict with delete done by list being owner.
1331     }
1332   }
1333   else
1334     fCaloPhotonsPHOS->Clear(); // TODO: redundant???
1335 }
1336 //_____________________________________________________________________________
1337 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x)const{
1338   //FillHistogram
1339   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
1340   if(tmpI){
1341     tmpI->Fill(x) ;
1342     return ;
1343   }
1344   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
1345   if(tmpF){
1346     tmpF->Fill(x) ;
1347     return ;
1348   }
1349   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
1350   if(tmpD){
1351     tmpD->Fill(x) ;
1352     return ;
1353   }
1354   AliInfo(Form("can not find histogram <%s> ",key)) ;
1355 }
1356 //_____________________________________________________________________________
1357 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y)const{
1358   //FillHistogram
1359   TObject * tmp = fOutputContainer->FindObject(key) ;
1360   if(!tmp){
1361     AliInfo(Form("can not find histogram <%s> ",key)) ;
1362     return ;
1363   }
1364   if(tmp->IsA() == TClass::GetClass("TH1F")){
1365     ((TH1F*)tmp)->Fill(x,y) ;
1366     return ;
1367   }
1368   if(tmp->IsA() == TClass::GetClass("TH2F")){
1369     ((TH2F*)tmp)->Fill(x,y) ;
1370     return ;
1371   }
1372   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1373 }
1374
1375 //_____________________________________________________________________________
1376 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1377   //Fills 1D histograms with key
1378   TObject * tmp = fOutputContainer->FindObject(key) ;
1379   if(!tmp){
1380     AliInfo(Form("can not find histogram <%s> ",key)) ;
1381     return ;
1382   }
1383   if(tmp->IsA() == TClass::GetClass("TH2F")){
1384     ((TH2F*)tmp)->Fill(x,y,z) ;
1385     return ;
1386   }
1387   if(tmp->IsA() == TClass::GetClass("TH3F")){
1388     ((TH3F*)tmp)->Fill(x,y,z) ;
1389     return ;
1390   }
1391 }
1392
1393 //_____________________________________________________________________________
1394 AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
1395 {
1396   fEvent = InputEvent();
1397   if( ! fEvent ) {
1398     AliError("Event could not be retrieved");
1399     PostData(1, fOutputContainer);
1400   }
1401   return fEvent;
1402 }
1403
1404
1405 //___________________________________________________________________________
1406 AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
1407 {
1408   fMCStack = 0;
1409   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
1410   if(eventHandler){
1411     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
1412     if( mcEventHandler)
1413       fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
1414   }
1415   return fMCStack;
1416 }
1417
1418 //___________________________________________________________________________
1419 Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
1420 {
1421  /* fCentBin=1+Int_t(centralityV0M/100. *kNCenBins) ;
1422   if(centralityV0M < 5. || fCentBin < 0)
1423    fCentBin=0 ;
1424   if(fCentBin > kNCenBins-1)
1425     fCentBin = kNCenBins-1 ;
1426  */
1427   int lastBinUpperIndex = fCentEdges.GetSize() -1;
1428   if( centralityV0M > fCentEdges[lastBinUpperIndex] ) {
1429     if( fDebug >= 1 )
1430       AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1431     return lastBinUpperIndex-1;
1432   }
1433   if( centralityV0M < fCentEdges[0] ) {
1434     if( fDebug >= 1 )
1435       AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1436     return 0;
1437   }
1438   
1439   fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1440   return fCentBin;
1441 }
1442
1443 //___________________________________________________________________________
1444 Int_t AliAnalysisTaskPi0Flow::GetRPBin()
1445 {
1446   Double_t averageRP;
1447   if(fHaveTPCRP)
1448     averageRP = fRPV0A+fRPV0C+fRP /3.;
1449   else
1450     averageRP = fRPV0A+fRPV0C /2.;
1451
1452   fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1453
1454   if(fEMRPBin> (Int_t) fNEMRPBins-1)
1455     fEMRPBin=fNEMRPBins-1 ;
1456   else if(fEMRPBin<0)
1457     fEMRPBin=0;
1458
1459   if ( fDebug >= 2 )
1460     AliInfo(Form("Event Mixing Reaction Plane bin is: %d", fEMRPBin));
1461
1462   return fEMRPBin;
1463 }
1464
1465
1466 //_____________________________________________________________________________
1467 void AliAnalysisTaskPi0Flow::LogProgress(int step)
1468 {
1469   if(fDebug >= 2) {
1470     AliInfo(Form("step %d completed", step));
1471   }
1472   // the +0.5 is not realy neccisarry, but oh well... -henrik
1473   //FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1474   //FillHistogram("hTotSelEvents", step+0.5);
1475 }
1476
1477 void AliAnalysisTaskPi0Flow::LogSelection(int step, int internalRunNumber)
1478 {
1479   // if(fDebug > 1) {
1480   //   AliInfo(Form("step %d completed", step));
1481   // }
1482   // the +0.5 is not realy neccisarry, but oh well... -henrik
1483   FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1484   FillHistogram("hTotSelEvents", step+0.5);
1485 }
1486
1487
1488 //___________________________________________________________________________
1489 Int_t AliAnalysisTaskPi0Flow::ConvertToInternalRunNumber(Int_t run){
1490   if( kLHC11h == fPeriod ) {
1491     switch(run){
1492     case  170593 : return 179 ;
1493     case  170572 : return 178 ;
1494     case  170556 : return 177 ;
1495     case  170552 : return 176 ;
1496     case  170546 : return 175 ;
1497     case  170390 : return 174 ;
1498     case  170389 : return 173 ;
1499     case  170388 : return 172 ;
1500     case  170387 : return 171 ;
1501     case  170315 : return 170 ;
1502     case  170313 : return 169 ;
1503     case  170312 : return 168 ;
1504     case  170311 : return 167 ;
1505     case  170309 : return 166 ;
1506     case  170308 : return 165 ;
1507     case  170306 : return 164 ;
1508     case  170270 : return 163 ;
1509     case  170269 : return 162 ;
1510     case  170268 : return 161 ;
1511     case  170267 : return 160 ;
1512     case  170264 : return 159 ;
1513     case  170230 : return 158 ;
1514     case  170228 : return 157 ;
1515     case  170208 : return 156 ;
1516     case  170207 : return 155 ;
1517     case  170205 : return 154 ;
1518     case  170204 : return 153 ;
1519     case  170203 : return 152 ;
1520     case  170195 : return 151 ;
1521     case  170193 : return 150 ;
1522     case  170163 : return 149 ;
1523     case  170162 : return 148 ;
1524     case  170159 : return 147 ;
1525     case  170155 : return 146 ;
1526     case  170152 : return 145 ;
1527     case  170091 : return 144 ;
1528     case  170089 : return 143 ;
1529     case  170088 : return 142 ;
1530     case  170085 : return 141 ;
1531     case  170084 : return 140 ;
1532     case  170083 : return 139 ;
1533     case  170081 : return 138 ;
1534     case  170040 : return 137 ;
1535     case  170038 : return 136 ;
1536     case  170036 : return 135 ;
1537     case  170027 : return 134 ;
1538     case  169981 : return 133 ;
1539     case  169975 : return 132 ;
1540     case  169969 : return 131 ;
1541     case  169965 : return 130 ;
1542     case  169961 : return 129 ;
1543     case  169956 : return 128 ;
1544     case  169926 : return 127 ;
1545     case  169924 : return 126 ;
1546     case  169923 : return 125 ;
1547     case  169922 : return 124 ;
1548     case  169919 : return 123 ;
1549     case  169918 : return 122 ;
1550     case  169914 : return 121 ;
1551     case  169859 : return 120 ;
1552     case  169858 : return 119 ;
1553     case  169855 : return 118 ;
1554     case  169846 : return 117 ;
1555     case  169838 : return 116 ;
1556     case  169837 : return 115 ;
1557     case  169835 : return 114 ;
1558     case  169683 : return 113 ;
1559     case  169628 : return 112 ;
1560     case  169591 : return 111 ;
1561     case  169590 : return 110 ;
1562     case  169588 : return 109 ;
1563     case  169587 : return 108 ;
1564     case  169586 : return 107 ;
1565     case  169584 : return 106 ;
1566     case  169557 : return 105 ;
1567     case  169555 : return 104 ;
1568     case  169554 : return 103 ;
1569     case  169553 : return 102 ;
1570     case  169550 : return 101 ;
1571     case  169515 : return 100 ;
1572     case  169512 : return 99 ;
1573     case  169506 : return 98 ;
1574     case  169504 : return 97 ;
1575     case  169498 : return 96 ;
1576     case  169475 : return 95 ;
1577     case  169420 : return 94 ;
1578     case  169419 : return 93 ;
1579     case  169418 : return 92 ;
1580     case  169417 : return 91 ;
1581     case  169415 : return 90 ;
1582     case  169411 : return 89 ;
1583     case  169238 : return 88 ;
1584     case  169236 : return 87 ;
1585     case  169167 : return 86 ;
1586     case  169160 : return 85 ;
1587     case  169156 : return 84 ;
1588     case  169148 : return 83 ;
1589     case  169145 : return 82 ;
1590     case  169144 : return 81 ;
1591     case  169143 : return 80 ;
1592     case  169138 : return 79 ;
1593     case  169099 : return 78 ;
1594     case  169094 : return 77 ;
1595     case  169091 : return 76 ;
1596     case  169045 : return 75 ;
1597     case  169044 : return 74 ;
1598     case  169040 : return 73 ;
1599     case  169035 : return 72 ;
1600     case  168992 : return 71 ;
1601     case  168988 : return 70 ;
1602     case  168984 : return 69 ;
1603     case  168826 : return 68 ;
1604     case  168777 : return 67 ;
1605     case  168514 : return 66 ;
1606     case  168512 : return 65 ;
1607     case  168511 : return 64 ;
1608     case  168467 : return 63 ;
1609     case  168464 : return 62 ;
1610     case  168461 : return 61 ;
1611     case  168460 : return 60 ;
1612     case  168458 : return 59 ;
1613     case  168362 : return 58 ;
1614     case  168361 : return 57 ;
1615     case  168356 : return 56 ;
1616     case  168342 : return 55 ;
1617     case  168341 : return 54 ;
1618     case  168325 : return 53 ;
1619     case  168322 : return 52 ;
1620     case  168318 : return 51 ;
1621     case  168311 : return 50 ;
1622     case  168310 : return 49 ;
1623     case  168213 : return 48 ;
1624     case  168212 : return 47 ;
1625     case  168208 : return 46 ;
1626     case  168207 : return 45 ;
1627     case  168206 : return 44 ;
1628     case  168205 : return 43 ;
1629     case  168204 : return 42 ;
1630     case  168203 : return 41 ;
1631     case  168181 : return 40 ;
1632     case  168177 : return 39 ;
1633     case  168175 : return 38 ;
1634     case  168173 : return 37 ;
1635     case  168172 : return 36 ;
1636     case  168171 : return 35 ;
1637     case  168115 : return 34 ;
1638     case  168108 : return 33 ;
1639     case  168107 : return 32 ;
1640     case  168105 : return 31 ;
1641     case  168104 : return 30 ;
1642     case  168103 : return 29 ;
1643     case  168076 : return 28 ;
1644     case  168069 : return 27 ;
1645     case  168068 : return 26 ;
1646     case  168066 : return 25 ;
1647     case  167988 : return 24 ;
1648     case  167987 : return 23 ;
1649     case  167986 : return 22 ;
1650     case  167985 : return 21 ;
1651     case  167921 : return 20 ;
1652     case  167920 : return 19 ;
1653     case  167915 : return 18 ;
1654     case  167909 : return 17 ;
1655     case  167903 : return 16 ;
1656     case  167902 : return 15 ;
1657     case  167818 : return 14 ;
1658     case  167814 : return 13 ;
1659     case  167813 : return 12 ;
1660     case  167808 : return 11 ;
1661     case  167807 : return 10 ;
1662     case  167806 : return 9 ;
1663     case  167713 : return 8 ;
1664     case  167712 : return 7 ;
1665     case  167711 : return 6 ;
1666     case  167706 : return 5 ;
1667     case  167693 : return 4 ;
1668     case  166532 : return 3 ;
1669     case  166530 : return 2 ;
1670     case  166529 : return 1 ;
1671
1672     default : return 199;
1673     }
1674   }
1675   if( kLHC10h == fPeriod ) {
1676     switch(run){
1677     case  139517 : return 137;
1678     case  139514 : return 136;
1679     case  139513 : return 135;
1680     case  139511 : return 134;
1681     case  139510 : return 133;
1682     case  139507 : return 132;
1683     case  139505 : return 131;
1684     case  139504 : return 130;
1685     case  139503 : return 129;
1686     case  139470 : return 128;
1687     case  139467 : return 127;
1688     case  139466 : return 126;
1689     case  139465 : return 125;
1690     case  139440 : return 124;
1691     case  139439 : return 123;
1692     case  139438 : return 122;
1693     case  139437 : return 121;
1694     case  139360 : return 120;
1695     case  139329 : return 119;
1696     case  139328 : return 118;
1697     case  139314 : return 117;
1698     case  139311 : return 116;
1699     case  139310 : return 115;
1700     case  139309 : return 114;
1701     case  139308 : return 113;
1702     case  139173 : return 112;
1703     case  139172 : return 111;
1704     case  139110 : return 110;
1705     case  139107 : return 109;
1706     case  139105 : return 108;
1707     case  139104 : return 107;
1708     case  139042 : return 106;
1709     case  139038 : return 105;
1710     case  139037 : return 104;
1711     case  139036 : return 103;
1712     case  139029 : return 102;
1713     case  139028 : return 101;
1714     case  138983 : return 100;
1715     case  138982 : return 99;
1716     case  138980 : return 98;
1717     case  138979 : return 97;
1718     case  138978 : return 96;
1719     case  138977 : return 95;
1720     case  138976 : return 94;
1721     case  138973 : return 93;
1722     case  138972 : return 92;
1723     case  138965 : return 91;
1724     case  138924 : return 90;
1725     case  138872 : return 89;
1726     case  138871 : return 88;
1727     case  138870 : return 87;
1728     case  138837 : return 86;
1729     case  138830 : return 85;
1730     case  138828 : return 84;
1731     case  138826 : return 83;
1732     case  138796 : return 82;
1733     case  138795 : return 81;
1734     case  138742 : return 80;
1735     case  138732 : return 79;
1736     case  138730 : return 78;
1737     case  138666 : return 77;
1738     case  138662 : return 76;
1739     case  138653 : return 75;
1740     case  138652 : return 74;
1741     case  138638 : return 73;
1742     case  138624 : return 72;
1743     case  138621 : return 71;
1744     case  138583 : return 70;
1745     case  138582 : return 69;
1746     case  138579 : return 68;
1747     case  138578 : return 67;
1748     case  138534 : return 66;
1749     case  138469 : return 65;
1750     case  138442 : return 64;
1751     case  138439 : return 63;
1752     case  138438 : return 62;
1753     case  138396 : return 61;
1754     case  138364 : return 60;
1755     case  138359 : return 59;
1756     case  138275 : return 58;
1757     case  138225 : return 57;
1758     case  138201 : return 56;
1759     case  138200 : return 55;
1760     case  138197 : return 54;
1761     case  138192 : return 53;
1762     case  138190 : return 52;
1763     case  138154 : return 51;
1764     case  138153 : return 50;
1765     case  138151 : return 49;
1766     case  138150 : return 48;
1767     case  138126 : return 47;
1768     case  138125 : return 46;
1769     case  137848 : return 45;
1770     case  137847 : return 44;
1771     case  137844 : return 43;
1772     case  137843 : return 42;
1773     case  137752 : return 41;
1774     case  137751 : return 40;
1775     case  137748 : return 39;
1776     case  137724 : return 38;
1777     case  137722 : return 37;
1778     case  137718 : return 36;
1779     case  137704 : return 35;
1780     case  137693 : return 34;
1781     case  137692 : return 33;
1782     case  137691 : return 32;
1783     case  137689 : return 31;
1784     case  137686 : return 30;
1785     case  137685 : return 29;
1786     case  137639 : return 28;
1787     case  137638 : return 27;
1788     case  137608 : return 26;
1789     case  137595 : return 25;
1790     case  137549 : return 24;
1791     case  137546 : return 23;
1792     case  137544 : return 22;
1793     case  137541 : return 21;
1794     case  137539 : return 20;
1795     case  137531 : return 19;
1796     case  137530 : return 18;
1797     case  137443 : return 17;
1798     case  137441 : return 16;
1799     case  137440 : return 15;
1800     case  137439 : return 14;
1801     case  137434 : return 13;
1802     case  137432 : return 12;
1803     case  137431 : return 11;
1804     case  137430 : return 10;
1805     case  137366 : return 9;
1806     case  137243 : return 8;
1807     case  137236 : return 7;
1808     case  137235 : return 6;
1809     case  137232 : return 5;
1810     case  137231 : return 4;
1811     case  137165 : return 3;
1812     case  137162 : return 2;
1813     case  137161 : return 1;
1814     default : return 199;
1815     }
1816   }
1817   if(kUndefinedPeriod && fDebug >= 1 ) {
1818     AliWarning("Period not defined");
1819   }
1820   return 1;
1821 }
1822 //_____________________________________________________________________________
1823 Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1824
1825   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1826   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1827   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1828   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1829   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1830   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1831     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1832     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1833   return (R2<2.5*2.5) ;
1834
1835 }
1836 //_____________________________________________________________________________
1837 Bool_t AliAnalysisTaskPi0Flow::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1838
1839   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1840   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1841   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1842   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1843   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1844   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1845     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1846     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1847   return (R2<1.5*1.5) ;
1848
1849 }
1850 //____________________________________________________________________________
1851 TList* AliAnalysisTaskPi0Flow::GetCaloPhotonsPHOSList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin)
1852 {
1853   int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins
1854               + centBin * fNEMRPBins
1855               + rpBin;
1856   if( fCaloPhotonsPHOSLists->At(offset) ) { // list exists
1857     TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1858     if( ! list )
1859       AliError("object in fCaloPhotonsPHOSLists at %i did not cast");
1860     return list;
1861   }
1862   else {// no list for this bin has been created, yet
1863     TList* list = new TList();
1864     list->SetOwner();
1865     fCaloPhotonsPHOSLists->AddAt(list, offset);
1866     return list;
1867   }
1868 }
1869
1870 //____________________________________________________________________________
1871 Double_t AliAnalysisTaskPi0Flow::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1872   //Parameterization of LHC10h period
1873   //_true if neutral_
1874
1875   Double_t meanX=0;
1876   Double_t meanZ=0.;
1877   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1878                          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);
1879   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) ;
1880   AliVEvent *event = InputEvent();
1881   Double_t mf = 0.; //
1882   if(event) mf = event->GetMagneticField(); //Positive for ++ and negative for --
1883
1884   if(mf<0.){ //field --
1885     meanZ = -0.468318 ;
1886     if(charge>0)
1887       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)) ;
1888     else
1889       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)) ;
1890   }
1891   else{ //Field ++
1892     meanZ= -0.468318;
1893     if(charge>0)
1894       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)) ;
1895     else
1896       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)) ;
1897   }
1898
1899   Double_t rz=(dz-meanZ)/sz ;
1900   Double_t rx=(dx-meanX)/sx ;
1901   return TMath::Sqrt(rx*rx+rz*rz) ;
1902 }
1903 //____________________________________________________________________________
1904 void AliAnalysisTaskPi0Flow::SetFlatteningData(){
1905   //Read objects with flattening parameters 
1906   AliOADBContainer flatContainer("phosFlat");
1907   flatContainer.InitFromFile(fEPcalibFileName.Data(),"phosFlat");
1908   TObjArray *maps = (TObjArray*)flatContainer.GetObject(fRunNumber,"phosFlat");
1909   if(!maps){
1910       AliError(Form("Can not read Flattening for run %d. \n From file >%s<\n",fRunNumber,fEPcalibFileName.Data())) ;    
1911   }
1912   else{
1913     AliInfo(Form("Setting PHOS flattening with name %s \n",maps->GetName())) ;
1914     AliPHOSEPFlattener * h = (AliPHOSEPFlattener*)maps->At(0) ;  
1915     if(fTPCFlat) delete fTPCFlat ;
1916     fTPCFlat = new AliPHOSEPFlattener() ;
1917     fTPCFlat = h ;
1918     h = (AliPHOSEPFlattener*)maps->At(1) ;  
1919     if(fV0AFlat) delete fV0AFlat ;
1920     fV0AFlat = new AliPHOSEPFlattener() ;
1921     fV0AFlat = h ;
1922     h = (AliPHOSEPFlattener*)maps->At(2) ;  
1923     if(fV0CFlat) delete fV0CFlat ;
1924     fV0CFlat = new AliPHOSEPFlattener() ;
1925     fV0CFlat = h ;
1926   }    
1927   
1928 }
1929  //____________________________________________________________________________
1930 Double_t  AliAnalysisTaskPi0Flow::ApplyFlattening(Double_t phi, Double_t c){
1931   
1932   if(fTPCFlat)
1933     return fTPCFlat->MakeFlat(phi,c);
1934   return phi ;
1935
1936 }
1937 //____________________________________________________________________________
1938 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0A(Double_t phi, Double_t c){
1939   
1940   if(fV0AFlat)
1941     return fV0AFlat->MakeFlat(phi,c);
1942   return phi ;
1943
1944 }
1945 //____________________________________________________________________________
1946 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0C(Double_t phi, Double_t c){
1947   
1948   if(fV0CFlat)
1949     return fV0CFlat->MakeFlat(phi,c);
1950   return phi ;
1951
1952 }
1953 //____________________________________________________________________________
1954 Double_t  AliAnalysisTaskPi0Flow::CoreEnergy(AliVCluster * clu, AliVCaloCells * cells)
1955 {
1956   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1957
1958   //Can not use already calculated coordinates?
1959   //They have incidence correction...
1960   const Double_t distanceCut =3.5 ;
1961   const Double_t logWeight=4.5 ;
1962
1963   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1964 // Calculates the center of gravity in the local PHOS-module coordinates
1965   Float_t wtot = 0;
1966   Double_t xc[100]={0} ;
1967   Double_t zc[100]={0} ;
1968   Double_t x = 0 ;
1969   Double_t z = 0 ;
1970   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1971   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1972     Int_t relid[4] ;
1973     Float_t xi ;
1974     Float_t zi ;
1975     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1976     fPHOSGeo->RelPosInModule(relid, xi, zi);
1977     xc[iDigit]=xi ;
1978     zc[iDigit]=zi ;
1979     elist[iDigit] *= cells->GetCellAmplitude(clu->GetCellsAbsId()[iDigit]);
1980     if( fDebug >= 3 )
1981       printf("%f ",elist[iDigit]);
1982     if (clu->E()>0 && elist[iDigit]>0) {
1983       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1984       x    += xc[iDigit] * w ;
1985       z    += zc[iDigit] * w ;
1986       wtot += w ;
1987     }
1988   }
1989   if (wtot>0) {
1990     x /= wtot ;
1991     z /= wtot ;
1992   }
1993   Double_t coreE=0. ;
1994   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1995     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1996     if(distance < distanceCut)
1997       coreE += elist[iDigit] ;
1998   }
1999   //Apply non-linearity correction
2000   return fNonLinCorr->Eval(coreE) ;
2001 }
2002 //____________________________________________________________________________
2003 Bool_t  AliAnalysisTaskPi0Flow::AreNeibors(Int_t id1,Int_t id2){
2004   // return true if absId are "Neighbors" (adjacent, including diagornaly,)
2005   // false if not.
2006
2007   Int_t relid1[4] ;
2008   fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
2009
2010   Int_t relid2[4] ;
2011   fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
2012
2013   // if inside the same PHOS module
2014   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) {
2015     const Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
2016     const Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
2017
2018     // and if diff in both direction is 1 or less
2019     if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
2020       return true; // are neighbors
2021   }
2022
2023   // else false
2024   return false;
2025 }
2026 //____________________________________________________________________________
2027 void  AliAnalysisTaskPi0Flow::Reclusterize(AliVCluster * clu){
2028   //Re-clusterize to make continues cluster
2029
2030   const Int_t oldMulDigit=clu->GetNCells() ;
2031   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2032   UShort_t * dlist = clu->GetCellsAbsId();
2033
2034   Int_t index[oldMulDigit] ;
2035   Bool_t used[oldMulDigit] ;
2036   for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
2037   Int_t inClu=0 ;
2038   Double_t eMax=0. ;
2039   //find maximum
2040   for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
2041     if(eMax<elist[iDigit]){
2042       eMax=elist[iDigit];
2043       index[0]=iDigit ;
2044       inClu=1 ;
2045     }
2046   }
2047   if(inClu==0){ //empty cluster
2048     return ;
2049   }
2050   used[index[0]]=kTRUE ; //mark as used
2051   for(Int_t i=0; i<inClu; i++){
2052     for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
2053        if(used[iDigit]) //already used
2054          continue ;
2055        if(AreNeibors(dlist[index[i]],dlist[iDigit])){
2056          index[inClu]= iDigit ;
2057          inClu++ ;
2058          used[iDigit]=kTRUE ;
2059        }
2060     }
2061   }
2062
2063   if(inClu==oldMulDigit) //no need to modify
2064     return ;
2065
2066   clu->SetNCells(inClu);
2067   //copy
2068   UShort_t tmpD[oldMulDigit] ;
2069   Double_t tmpE[oldMulDigit] ;
2070   for(Int_t i=0; i<oldMulDigit; i++){
2071     tmpD[i]=dlist[i] ;
2072     tmpE[i]=elist[i] ;
2073   }
2074   //change order of digits in list so that
2075   //first inClu cells were true ones
2076   for(Int_t i=0; i<inClu; i++){
2077     dlist[i]=tmpD[index[i]] ;
2078     elist[i]=tmpE[index[i]] ;
2079   }
2080
2081
2082 }
2083
2084 //_____________________________________________________________________________
2085 void AliAnalysisTaskPi0Flow::SetMisalignment(){
2086   // sets the misalignment vertex if ESD
2087   if( fEventESD ) {
2088     for(Int_t mod=0; mod<5; mod++) {
2089       const TGeoHMatrix* modMatrix = fEvent->GetPHOSMatrix(mod);
2090       if( ! modMatrix) {
2091         if( fDebug )
2092           AliInfo(Form("no PHOS Geometric Misalignment Matrix for module %d", mod));
2093         continue;
2094       }
2095       else {
2096         fPHOSGeo->SetMisalMatrix(modMatrix, mod);
2097         if( fDebug )
2098           AliInfo(Form("PHOS Geometric Misalignment Matrix set for module %d", mod));
2099       }
2100     }
2101   }
2102 }
2103
2104 //_____________________________________________________________________________
2105 void AliAnalysisTaskPi0Flow::SetV0Calibration(){
2106     // assigns: fMultV0, fV0Cpol, fV0Apol, fMeanQ, and fWidthQ
2107
2108     int runNumber = this->fRunNumber;
2109     
2110     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
2111     TFile *foadb = TFile::Open(oadbfilename.Data());
2112
2113     if(!foadb){
2114         AliError(Form("OADB file %s cannot be opened\n", oadbfilename.Data()));
2115         AliError("V0 Calibration not set !\n");
2116         return;
2117     }
2118
2119     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
2120     if(!cont){
2121         AliError("OADB object hMultV0BefCorr is not available in the file");
2122         AliError("V0 Calibration not set!\n");
2123         return;
2124     }
2125
2126     if(!(cont->GetObject(runNumber))){
2127         AliError(Form("OADB object hMultV0BefCorr is not available for run %i, trying 137366)",runNumber));
2128         runNumber = 137366;
2129     }
2130     if(!(cont->GetObject(runNumber))){
2131         AliError(Form("OADB object hMultV0BefCorr is not available for run %i ",runNumber));
2132         AliError("V0 Calibration not set!\n");
2133         return;
2134     }
2135
2136     if( fDebug )  AliInfo("Setting V0 calibration") ;
2137     fMultV0 = ((TH2F *) cont->GetObject(runNumber))->ProfileX();
2138
2139     TF1 *fpol0 = new TF1("fpol0","pol0");
2140     fMultV0->Fit(fpol0,"Q0","",0,31);
2141     fV0Cpol = fpol0->GetParameter(0);
2142     fMultV0->Fit(fpol0,"Q0","",32,64);
2143     fV0Apol = fpol0->GetParameter(0);
2144
2145     for(Int_t iside=0;iside<2;iside++){
2146         for(Int_t icoord=0;icoord<2;icoord++){
2147             for(Int_t i=0;i  < kNCenBins;i++){
2148                 char namecont[100];
2149                 if(iside==0 && icoord==0)
2150                     snprintf(namecont,100,"hQxc2_%i",i);
2151                 else if(iside==1 && icoord==0)
2152                     snprintf(namecont,100,"hQxa2_%i",i);
2153                 else if(iside==0 && icoord==1)
2154                     snprintf(namecont,100,"hQyc2_%i",i);
2155                 else if(iside==1 && icoord==1)
2156                     snprintf(namecont,100,"hQya2_%i",i);
2157
2158                 cont = (AliOADBContainer*) foadb->Get(namecont);
2159                 if(!cont){
2160                     AliError(Form("OADB object %s is not available in the file %s", namecont, oadbfilename.Data()));
2161                     AliError("V0 Calibration not fully set!\n");
2162                     return;
2163                 }
2164
2165                 if(!(cont->GetObject(runNumber))){
2166                     AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2167                     runNumber = 137366;
2168                 }
2169                 if(!(cont->GetObject(runNumber))){
2170                   AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2171                   AliError("V0 Calibration not fully set!\n");
2172                   return;
2173                 }
2174                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetMean();
2175                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetRMS();
2176
2177                 //for v3
2178 //              if(iside==0 && icoord==0)
2179 //                  snprintf(namecont,100,"hQxc3_%i",i);
2180 //              else if(iside==1 && icoord==0)
2181 //                  snprintf(namecont,100,"hQxa3_%i",i);
2182 //              else if(iside==0 && icoord==1)
2183 //                  snprintf(namecont,100,"hQyc3_%i",i);
2184 //              else if(iside==1 && icoord==1)
2185 //                  snprintf(namecont,100,"hQya3_%i",i);
2186 // 
2187 //              cont = (AliOADBContainer*) foadb->Get(namecont);
2188 //              if(!cont){
2189 //                  AliError(Form("OADB object %s is not available in the file",namecont));
2190 //                  AliError("V0 Calibration not fully set!\n");
2191 //                  return;
2192 //              }
2193 // 
2194 //              if(!(cont->GetObject(runNumber))){
2195 //                  AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2196 //                  runNumber = 137366;
2197 //              }
2198 //              if(!(cont->GetObject(runNumber))){
2199 //                AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2200 //                AliError("V0 Calibration not fully set!\n");
2201 //                return;
2202 //              }
2203 //              fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2204 //              fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2205
2206             }
2207         }
2208     }
2209
2210     delete fpol0; fpol0=0;
2211 }
2212
2213 //_____________________________________________________________________________
2214 void AliAnalysisTaskPi0Flow::SetESDTrackCuts()
2215 {
2216   if( fEventESD ) {
2217     // Create ESD track cut
2218     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
2219     //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
2220     fESDtrackCuts->SetRequireTPCRefit(kTRUE);
2221   }
2222 }
2223
2224 //_____________________________________________________________________________
2225 void AliAnalysisTaskPi0Flow::SetGeometry()
2226 {
2227   // Initialize the PHOS geometry
2228   if( kLHC10h == fPeriod && fEventESD ) {
2229     TGeoManager::Import("geometry.root"); //TODO: should perhaps not be done
2230     fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
2231     if( ! fPHOSGeo )
2232       AliError("geometry (fPHOSGeo) not initialised");
2233   }
2234
2235   //Init geometry
2236   if(!fPHOSGeo){
2237     AliOADBContainer geomContainer("phosGeo");
2238     geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
2239     TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
2240     fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
2241     for(Int_t mod=0; mod<5; mod++) {
2242       if(!matrixes->At(mod)) {
2243         continue;
2244         if( fDebug )
2245           AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2246       }
2247       else {
2248         fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
2249         if( fDebug >1 )
2250           AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2251       }
2252     }
2253   } 
2254 }
2255
2256 //_____________________________________________________________________________
2257 void AliAnalysisTaskPi0Flow::SetPHOSCalibData()
2258 {
2259   if( fPHOSCalibData )
2260     delete fPHOSCalibData; 
2261   fPHOSCalibData = 0;
2262   
2263   // Calibration only needed for ESD
2264   if( fEventESD /*&& */ ) {
2265     if( kLHC10h == fPeriod && fEventESD ) {
2266       //We have to apply re-calibration for pass1 LCH10h
2267       // Initialize decalibration factors in the form of the OCDB object
2268       AliCDBManager * man = AliCDBManager::Instance();
2269       man->SetRun(140000) ; //TODO; revise, this should probably not b done.
2270       man->SetDefaultStorage("local://OCDB");
2271     }
2272     fPHOSCalibData = new AliPHOSCalibData();
2273   }
2274 }
2275
2276 //_____________________________________________________________________________
2277 void AliAnalysisTaskPi0Flow::SetVertex()
2278 {
2279   const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
2280   if( primaryVertex ) {
2281     fVertex[0] = primaryVertex->GetX();
2282     fVertex[1] = primaryVertex->GetY();
2283     fVertex[2] = primaryVertex->GetZ();
2284   }
2285   else {
2286     AliError("Event has 0x0 Primary Vertex, defaulting to origo");
2287     fVertex[0] = 0;
2288     fVertex[1] = 0;
2289     fVertex[2] = 0;
2290   }
2291   fVertexVector = TVector3(fVertex);
2292   FillHistogram("hZvertex", fVertexVector.z(), fInternalRunNumber-0.5);
2293   
2294   if( fDebug >= 2 )
2295     AliInfo(Form("Vertex is set to (%.1f,%.1f,%.1f)", fVertex[0], fVertex[1], fVertex[2]));
2296
2297   fVtxBin=0 ;// No support for vtx binning implemented.
2298 }
2299
2300 //_____________________________________________________________________________
2301 Bool_t AliAnalysisTaskPi0Flow::RejectEventVertex()
2302 {
2303   if( ! fEvent->GetPrimaryVertex() )
2304     return true; // reject
2305   LogSelection(1, fInternalRunNumber);
2306
2307   if ( TMath::Abs(fVertexVector.z()) > 10. )
2308     return true; // reject
2309   LogSelection(2, fInternalRunNumber);
2310
2311   return false; // accept event.
2312 }
2313
2314 //_____________________________________________________________________________
2315 void AliAnalysisTaskPi0Flow::SetCentrality()
2316 {
2317   AliCentrality *centrality = fEvent->GetCentrality();
2318   if( centrality )
2319     fCentralityV0M=centrality->GetCentralityPercentile("V0M");
2320   else {
2321     AliError("Event has 0x0 centrality");
2322     fCentralityV0M = -1.;
2323   }
2324   FillHistogram("hCentrality",fCentralityV0M,fInternalRunNumber-0.5) ;
2325
2326   fCentBin = GetCentralityBin(fCentralityV0M);
2327
2328   if ( fDebug >= 2 )
2329     AliInfo(Form("Centrality (bin) is: %f (%d)", fCentralityV0M, fCentBin));
2330 }
2331
2332 //_____________________________________________________________________________
2333 Bool_t AliAnalysisTaskPi0Flow::RejectCentrality()
2334 {
2335   if( ! fEvent->GetCentrality() )
2336     return true; // reject
2337   LogSelection(3, fInternalRunNumber);
2338
2339 //   if( fCentralityV0M <= 0. || fCentralityV0M>80. )
2340 //     return true; // reject
2341     
2342   int lastBinUpperIndex = fCentEdges.GetSize() -1;
2343   if( fCentralityV0M > fCentEdges[lastBinUpperIndex] ) {
2344     if( fDebug )
2345       AliInfo("Rejecting due to centrality outside of binning.");
2346     return true; // reject
2347   }
2348   LogSelection(4, fInternalRunNumber);
2349
2350   if( fCentralityV0M < fCentEdges[0] ) {
2351     if( fDebug )
2352       AliInfo("Rejecting due to centrality outside of binning.");
2353     return true; // reject
2354   }
2355   LogSelection(5, fInternalRunNumber);
2356
2357   return false;
2358 }
2359
2360
2361 //_____________________________________________________________________________
2362 void AliAnalysisTaskPi0Flow::EvalReactionPlane()
2363 {
2364   // assigns: fHaveTPCRP and fRP
2365   // also does a few histogram fills
2366
2367   AliEventplane *eventPlane = fEvent->GetEventplane();
2368   if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2369   
2370   Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
2371   FillHistogram("phiRP",reactionPlaneQ,fCentralityV0M) ;
2372
2373   if(reactionPlaneQ==999 || reactionPlaneQ < 0.){ //reaction plain was not defined
2374     if( fDebug ) AliInfo(Form("No Q Reaction Plane, value is %f", reactionPlaneQ));
2375     fHaveTPCRP = kFALSE;
2376   }
2377   else{
2378     if( fDebug >= 2 ) AliInfo(Form("Q Reaction Plane is %f", reactionPlaneQ));
2379     fHaveTPCRP = kTRUE;
2380   }
2381
2382   if(fHaveTPCRP){
2383     fRP = ApplyFlattening(reactionPlaneQ, fCentralityV0M) ;
2384
2385     while(fRP<0)  fRP+=TMath::Pi();
2386     while(fRP>TMath::Pi())  fRP-=TMath::Pi();
2387     FillHistogram("phiRPflat",fRP,fCentralityV0M) ;
2388     Double_t dPsi = eventPlane->GetQsubRes() ;
2389     FillHistogram("cos2AC",TMath::Cos(2.*dPsi),fCentralityV0M) ;
2390   }
2391   else
2392     fRP=0.;
2393 }
2394
2395
2396 //____________________________________________________________________________
2397 void  AliAnalysisTaskPi0Flow::EvalV0ReactionPlane(){
2398   // set: fRPV0A and fRPV0C
2399
2400   //VZERO data
2401   AliVVZERO* v0 = fEvent->GetVZEROData();
2402
2403   //reset Q vector info
2404   Double_t Qxa2 = 0, Qya2 = 0;
2405   Double_t Qxc2 = 0, Qyc2 = 0;
2406
2407   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
2408     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
2409     Float_t multv0 = v0->GetMultiplicity(iv0);
2410     if (iv0 < 32){ // V0C
2411       Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2412       Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2413     } else {       // V0A
2414       Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2415       Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2416     }
2417   }
2418
2419   Int_t iC = -1;
2420   // centrality bins
2421   if(fCentralityV0M < 5) iC = 0;
2422   else if(fCentralityV0M < 10) iC = 1;
2423   else if(fCentralityV0M < 20) iC = 2;
2424   else if(fCentralityV0M < 30) iC = 3;
2425   else if(fCentralityV0M < 40) iC = 4;
2426   else if(fCentralityV0M < 50) iC = 5;
2427   else if(fCentralityV0M < 60) iC = 6;
2428   else if(fCentralityV0M < 70) iC = 7;
2429   else iC = 8;
2430
2431     //grab for each centrality the proper histo with the Qx and Qy to do the recentering
2432   Double_t Qxamean2 = fMeanQ[iC][1][0];
2433   Double_t Qxarms2  = fWidthQ[iC][1][0];
2434   Double_t Qyamean2 = fMeanQ[iC][1][1];
2435   Double_t Qyarms2  = fWidthQ[iC][1][1];
2436
2437   Double_t Qxcmean2 = fMeanQ[iC][0][0];
2438   Double_t Qxcrms2  = fWidthQ[iC][0][0];
2439   Double_t Qycmean2 = fMeanQ[iC][0][1];
2440   Double_t Qycrms2  = fWidthQ[iC][0][1];
2441
2442   Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
2443   Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
2444   Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
2445   Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
2446
2447   fRPV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
2448   fRPV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
2449
2450   while(fRPV0A<0)fRPV0A+=TMath::Pi() ;
2451   while(fRPV0A>TMath::Pi())fRPV0A-=TMath::Pi() ;
2452   while(fRPV0C<0)fRPV0C+=TMath::Pi() ;
2453   while(fRPV0C>TMath::Pi())fRPV0C-=TMath::Pi() ;
2454
2455   // Flattening
2456   fRPV0A=ApplyFlatteningV0A(fRPV0A,fCentralityV0M) ;
2457   while(fRPV0A<0)fRPV0A+=TMath::Pi() ;
2458   while(fRPV0A>TMath::Pi())fRPV0A-=TMath::Pi() ;
2459   fRPV0C=ApplyFlatteningV0C(fRPV0C,fCentralityV0M) ;
2460   while(fRPV0C<0)fRPV0C+=TMath::Pi() ;
2461   while(fRPV0C>TMath::Pi())fRPV0C-=TMath::Pi() ;
2462   
2463   if( fDebug >= 2 )
2464     AliInfo(Form("V0 Reaction Plane is: A side: %f, C side: %f", fRPV0A, fRPV0C));
2465
2466   FillHistogram("phiRPV0A",fRPV0A,fCentralityV0M);
2467   FillHistogram("phiRPV0C",fRPV0C,fCentralityV0M);
2468
2469   FillHistogram("phiRPV0Aflat",fRPV0A,fCentralityV0M) ;
2470   FillHistogram("phiRPV0AC",fRPV0A,fRPV0C,fCentralityV0M) ;
2471   FillHistogram("cos2V0AC",TMath::Cos(2.*(fRPV0A-fRPV0C)),fCentralityV0M) ;
2472   if(fHaveTPCRP){
2473     FillHistogram("phiRPV0ATPC",fRP,fRPV0A,fCentralityV0M) ;
2474     FillHistogram("cos2V0ATPC",TMath::Cos(2.*(fRP-fRPV0A)),fCentralityV0M) ;
2475   }
2476
2477   FillHistogram("phiRPV0Cflat",fRPV0C,fCentralityV0M) ;
2478   if(fHaveTPCRP){
2479     FillHistogram("phiRPV0CTPC",fRP,fRPV0C,fCentralityV0M) ;
2480     FillHistogram("cos2V0CTPC",TMath::Cos(2.*(fRP-fRPV0C)),fCentralityV0M) ;
2481   }
2482 }