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