]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0Flow.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added FillHistogram for 3d Hists with weight
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0Flow.cxx
CommitLineData
b8bea077 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
2cde7444 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"
9ed0694f 29#include "THashList.h"
2cde7444 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"
ada6b882 39#include "TGeoManager.h"
2cde7444 40#include "AliPHOSEsdCluster.h"
41#include "AliPHOSCalibData.h"
42#include "AliESDEvent.h"
43#include "AliESDCaloCells.h"
44#include "AliESDVertex.h"
45#include "AliESDtrackCuts.h"
ada6b882 46#include "AliAODEvent.h"
2cde7444 47#include "AliLog.h"
48#include "AliPID.h"
49#include "AliCDBManager.h"
ada6b882 50#include <AliAODCaloCluster.h>
51#include "AliCentrality.h"
2cde7444 52#include "AliESDtrackCuts.h"
53#include "AliEventplane.h"
54#include "TProfile.h"
55#include "AliOADBContainer.h"
9d29cbcb 56#include "AliEPFlattener.h"
2cde7444 57
ada6b882 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
91b112bd 62/* $Id$ */
2cde7444 63
ada6b882 64ClassImp(AliAnalysisTaskPi0Flow);
2cde7444 65
66//________________________________________________________________________
67Double_t rnlin(Double_t *x, Double_t * /*par*/)
68{
69 //a = par[0], b = par[1].
70 //1+a*exp(-e/b)
71
ada6b882 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]) ;
2cde7444 74
75}
76
77//________________________________________________________________________
ada6b882 78AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
2cde7444 79: AliAnalysisTaskSE(name),
ada6b882 80 fCentEdges(10),
107849ac 81 fCentNMixed(10),
ada6b882 82 fNEMRPBins(9),
ada6b882 83 fPeriod(period),
10a684bf 84 fMaxAbsVertexZ(10.),
6a1a4e92 85 fManualV0EPCalc(false),
2cde7444 86 fOutputContainer(0x0),
2cde7444 87 fNonLinCorr(0),
ada6b882 88 fEvent(0x0),
89 fEventESD(0x0),
90 fEventAOD(0x0),
ada6b882 91 fRunNumber(-999),
92 fInternalRunNumber(0),
93 fPHOSGeo(0),
94 fMultV0(0x0),
95 fV0Cpol(0.),fV0Apol(0.),
96 fESDtrackCuts(0x0),
97 fPHOSCalibData(0x0),
b509e3a5 98 fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"),
99 fTPCFlat(0x0),
100 fV0AFlat(0x0),
101 fV0CFlat(0x0),
ada6b882 102 fVertexVector(),
103 fVtxBin(0),
104 fCentralityV0M(0.),
105 fCentBin(0),
106 fHaveTPCRP(0),
2cde7444 107 fRP(0),
108 fRPV0A(0),
109 fRPV0C(0),
ada6b882 110 fEMRPBin(0),
111 fCaloPhotonsPHOS(0x0),
112 fCaloPhotonsPHOSLists(0x0)
2cde7444 113{
bba0c9ef 114 const int nbins = 9;
115 Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
85394e40 116 TArrayD centEdges(nbins+1, edges);
bba0c9ef 117 Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
85394e40 118 TArrayI centNMixed(nbins, nMixed);
119 SetCentralityBinning(centEdges, centNMixed);
ada6b882 120
052d1b8e 121 for(int mod=1; mod <= kNMod; ++mod)
bb73874e 122 fModuleEnabled[mod-1] = kTRUE;
4c10872b 123
ada6b882 124 for(Int_t i=0;i<kNCenBins;i++){
2cde7444 125 for(Int_t j=0;j<2; j++)
ada6b882 126 for(Int_t k=0; k<2; k++) {
127 fMeanQ[i][j][k]=0.;
128 fWidthQ[i][j][k]=0.;
129 }
2cde7444 130 }
131
ada6b882 132 fVertex[0]=0; fVertex[1]=0; fVertex[2]=0;
133
2cde7444 134 // Output slots #0 write into a TH1 container
135 DefineOutput(1,TList::Class());
136
ada6b882 137
138
2cde7444 139 // Set bad channel map
140 char key[55] ;
141 for(Int_t i=0; i<6; i++){
53c2a02a 142 snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
2cde7444 143 fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
144 }
2cde7444 145
146
2cde7444 147 // Initialize non-linrarity correction
148 fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
149
150
151}
ada6b882 152//___________________________________________________________________________
153AliAnalysisTaskPi0Flow::~AliAnalysisTaskPi0Flow()
154{
155 delete fNonLinCorr;
156 delete fESDtrackCuts;
157 delete fPHOSCalibData;
158 delete fCaloPhotonsPHOSLists;
b509e3a5 159 if(fTPCFlat)delete fTPCFlat; fTPCFlat=0x0;
160 if(fV0AFlat)delete fV0AFlat; fV0AFlat=0x0;
161 if(fV0CFlat)delete fV0CFlat; fV0CFlat=0x0;
162
ada6b882 163}
2cde7444 164//________________________________________________________________________
165void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
166{
167 // Create histograms
168 // Called once
169 const Int_t nRuns=200 ;
ada6b882 170
171 // histograms
2cde7444 172 if(fOutputContainer != NULL){
173 delete fOutputContainer;
174 }
9ed0694f 175 fOutputContainer = new THashList();
2cde7444 176 fOutputContainer->SetOwner(kTRUE);
ada6b882 177
2cde7444 178 //========QA histograms=======
179
180 //Event selection
ada6b882 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
2cde7444 184 //vertex distribution
185 fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
ada6b882 186
2cde7444 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.)) ;
ada6b882 191 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
2cde7444 192 fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,10.,100,0.,100.)) ;
ada6b882 193
194
2cde7444 195 //Reaction plane
196 fOutputContainer->Add(new TH3F("hPHOSphi","cos" ,10,0.,100.,20,0.,10.,100,-TMath::Pi(),TMath::Pi()));
ada6b882 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.)) ;
f5e08ac9 207 fOutputContainer->Add(new TH3F("phiRPV0AC","RP distribution with V0A and V0C", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
ada6b882 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
2cde7444 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));
ada6b882 231
2cde7444 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));
ada6b882 240
2cde7444 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
ada6b882 249
2cde7444 250 //Single photon and pi0 spectrum
56f54833 251 const Int_t nPtPhot = 400 ;
252 const Double_t ptPhotMax = 40 ;
ada6b882 253 const Int_t nM = 500;
254 const Double_t mMin = 0.0;
255 const Double_t mMax = 1.0;
256
2cde7444 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));
ada6b882 264
265 // Histograms for different centralities
f0532ae0 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"};
f402cec6 268 char key[55];
269 TString name, title;
0374f417 270 for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
f402cec6 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
2cde7444 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) ;
ada6b882 330 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
2cde7444 331
2cde7444 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) ;
ada6b882 341 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
ada6b882 342 }
343
393743be 344 this->MakeMCHistograms();
345
ada6b882 346 // Setup photon lists
347 Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
348 fCaloPhotonsPHOSLists = new TObjArray(kapacity);
349 fCaloPhotonsPHOSLists->SetOwner();
2cde7444 350
351 PostData(1, fOutputContainer);
2cde7444 352}
353
393743be 354void AliAnalysisTaskPi0Flow::MakeMCHistograms()
355{
356 //empty, MC extensions should override
357}
358
359void AliAnalysisTaskPi0Flow::DoMC()
360{
361 //empty, MC extensions should override
362}
363
364
2cde7444 365//________________________________________________________________________
ada6b882 366void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
2cde7444 367{
368 // Main loop, called for each event
369 // Analyze ESD/AOD
ada6b882 370
371
372 // Step 0: Event Objects
373 fEvent = GetEvent();
374 fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
375 fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
c9fb2807 376 LogProgress(0);
ada6b882 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();
b509e3a5 394 SetFlatteningData();
2cde7444 395 }
c9fb2807 396 LogProgress(1);
397 LogSelection(0, fInternalRunNumber);
ada6b882 398
399
400 // Step 2: Vertex
401 // fVertex, fVertexVector, fVtxBin
402 SetVertex();
403 if( RejectEventVertex() ) {
2cde7444 404 PostData(1, fOutputContainer);
ada6b882 405 return; // Reject!
2cde7444 406 }
c9fb2807 407 LogProgress(2);
2cde7444 408
ada6b882 409// Step 3:
410// if(event->IsPileupFromSPD()){
411// PostData(1, fOutputContainer);
412// return; // Reject!
413// }
c9fb2807 414 LogProgress(3);
2cde7444 415
2cde7444 416
ada6b882 417 // Step 4: Centrality
418 // fCentralityV0M, fCentBin
419 SetCentrality();
420 if( RejectCentrality() ){
2cde7444 421 PostData(1, fOutputContainer);
ada6b882 422 return; // Reject!
2cde7444 423 }
c9fb2807 424 LogProgress(4);
2cde7444 425
ada6b882 426
427 // Step 5: Reaction Plane
428 // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
f5e08ac9 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
c9fb2807 432 LogProgress(5);
ada6b882 433
434
435 // Step 6: MC
436 // ProcessMC() ;
c9fb2807 437 LogProgress(6);
ada6b882 438
439
440 // Step 7: QA PHOS cells
441 FillPHOSCellQAHists();
c9fb2807 442 LogProgress(7);
ada6b882 443
444
445 // Step 8: Event Photons (PHOS Clusters) selection
446 SelectPhotonClusters();
447 FillSelectedClusterHistograms();
c9fb2807 448 LogProgress(8);
ada6b882 449
3d7a884f 450 if( ! fCaloPhotonsPHOS->GetEntriesFast() )
451 return;
452 else
453 LogSelection(6, fInternalRunNumber);
454
455
ada6b882 456 // Step 9: Consider pi0 (photon/cluster) pairs.
457 ConsiderPi0s();
c9fb2807 458 LogProgress(9);
ada6b882 459
460 // Step 10; Mixing
461 ConsiderPi0sMix();
c9fb2807 462 LogProgress(10);
2cde7444 463
393743be 464 // Step 11: MC
465 this->DoMC();
466
467 // Step 12: Update lists
ada6b882 468 UpdateLists();
c9fb2807 469 LogProgress(11);
ada6b882 470
2cde7444 471
ada6b882 472 // Post output data.
473 PostData(1, fOutputContainer);
474}
2cde7444 475
ada6b882 476//________________________________________________________________________
4c5399e6 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// }
ada6b882 485//________________________________________________________________________
85394e40 486void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed)
ada6b882 487{
107849ac 488 // Define centrality bins by their edges
85394e40 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");
ada6b882 494
495 fCentEdges = edges;
107849ac 496 fCentNMixed = nMixed;
497}
2cde7444 498
88014928 499//_____________________________________________________________________________
500void AliAnalysisTaskPi0Flow::SetEnablePHOSModule(int module, Bool_t enable)
501{
052d1b8e 502 if( module < 1 || kNMod < module )
503 AliFatal(Form("PHOS Module must be between 1 and %i", kNMod));
b628ce8b 504 else
505 fModuleEnabled[module-1] = enable;
88014928 506}
85394e40 507
508
ada6b882 509//________________________________________________________________________
510void AliAnalysisTaskPi0Flow::SetPHOSBadMap(Int_t mod, TH2I* badMapHist)
511 {
512 if(fPHOSBadMap[mod])
513 delete fPHOSBadMap[mod];
514
515 fPHOSBadMap[mod]=new TH2I(*badMapHist);
f5e08ac9 516 if(fDebug)
ada6b882 517 AliInfo(Form("Setting Bad Map Histogram %s",fPHOSBadMap[mod]->GetName()));
518 }
2cde7444 519
ada6b882 520//________________________________________________________________________
521Bool_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
2cde7444 524
ada6b882 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 ;
2cde7444 538 }
539 else{
ada6b882 540 AliError(Form("Can not find bad channels for detector %s ",det)) ;
2cde7444 541 }
2cde7444 542
ada6b882 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;
2cde7444 549
ada6b882 550 return kTRUE ;
551}
552//_____________________________________________________________________________
553void AliAnalysisTaskPi0Flow::FillPHOSCellQAHists()
554{
e2139d72 555 // Fill cell occupancy per module
556
ada6b882 557 AliVCaloCells * cells = fEvent->GetPHOSCells();
2cde7444 558
ada6b882 559 FillHistogram("hCenPHOSCells",fCentralityV0M,cells->GetNumberOfCells()) ;
560 FillHistogram("hCenTrack",fCentralityV0M,fEvent->GetNumberOfTracks()) ;
2cde7444 561
2cde7444 562
2cde7444 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);
ada6b882 566 Int_t relId[4] = {0,0,0,0};
2cde7444 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 }
da83c923 592 FillHistogram("hCellMultEvent",nCellModule[0]+nCellModule[1]+nCellModule[2]);
2cde7444 593 FillHistogram("hCellMultEventM1",nCellModule[0]);
594 FillHistogram("hCellMultEventM2",nCellModule[1]);
595 FillHistogram("hCellMultEventM3",nCellModule[2]);
596
ada6b882 597}
598//_____________________________________________________________________________
599void 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
91b112bd 610 AliVCaloCells* cells = dynamic_cast<AliVCaloCells*> (fEvent->GetPHOSCells());
ada6b882 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
b98b7f51 616
617
ada6b882 618 // check if cell/channel is good.
2cde7444 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] ;
88014928 627 if ( ! fModuleEnabled[mod-1] )
628 continue;
ada6b882 629 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
630 continue ; // reject if not.
631
f9dbbc02 632 Double_t distBC=clu->GetDistanceToBadChannel();
633 if(distBC<kMinBCDistance)
634 continue ;
635
ada6b882 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;
91b112bd 643 ecore = CoreEnergy(clu,cells);
ada6b882 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) );
b98b7f51 649 cluPHOS1.Recalibrate(fPHOSCalibData, static_cast<AliESDCaloCells*> (cells)); // modify the cell energies
ada6b882 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);
ada6b882 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
2cde7444 675 }
ada6b882 676 else if (fEventAOD) { // is ! ESD, AOD.
677 AliESDCaloCluster* aodCluster = (AliESDCaloCluster*) (clu);
678 aodCluster->GetMomentum(lorentzMomentum ,origo);
2cde7444 679 }
ada6b882 680 else {
681 AliError("(Calo)Cluster is neither ESD nor AOD");
682 continue;
2cde7444 683 }
2cde7444 684
ada6b882 685
b98b7f51 686 FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.);
ada6b882 687 if(lorentzMomentum.E()>1.5){
b98b7f51 688 FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.);
2cde7444 689 }
ada6b882 690
691 fCaloPhotonsPHOS->Add(new AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
692 AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
f8c12091 693 ph->SetCluster(clu);
ada6b882 694
2cde7444 695 ph->SetModule(mod) ;
ada6b882 696 lorentzMomentum*= ecore/lorentzMomentum.E() ;
697 ph->SetMomV2(&lorentzMomentum) ;
2cde7444 698 ph->SetNCells(clu->GetNCells());
699 ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
6ca8b4a1 700 //Evaluate CoreDispersion
701 Double_t m02=0.,m20=0. ;
702 EvalCoreLambdas(clu, cells, m02, m20) ;
3c8a190f 703 ph->SetDisp2Bit(TestCoreLambda(clu->E(),m20,m02)) ; //Correct order m20,m02
704// ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
2cde7444 705 if(ph->IsDispOK()){
b98b7f51 706 FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
2cde7444 707 }
ada6b882 708
709 // Track Matching
2cde7444 710 Double_t dx=clu->GetTrackDx() ;
711 Double_t dz=clu->GetTrackDz() ;
f369ddfc 712 Bool_t cpvBit=kTRUE ; //No track matched by default. True means: not from charged, according to veto.
261d581c 713 Bool_t cpvBit2=kTRUE ; //More Strict criterion
ada6b882 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 }
2cde7444 740 }
741 }
742 ph->SetCPVBit(cpvBit) ;
743 ph->SetCPV2Bit(cpvBit2) ;
744 if(cpvBit){
b98b7f51 745 FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.);
2cde7444 746 }
747 ph->SetEMCx(float(cellX)) ;
748 ph->SetEMCz(float(cellZ)) ;
ada6b882 749 // ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
750 ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolde
751
2cde7444 752 }
6ec40e8c 753 FillHistogram("hCenPHOS",fCentralityV0M, fCaloPhotonsPHOS->GetEntriesFast()) ;
ada6b882 754}
755//_____________________________________________________________________________
756void AliAnalysisTaskPi0Flow::FillSelectedClusterHistograms()
757{
6ec40e8c 758 for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
ada6b882 759 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
2cde7444 760
2cde7444 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() ;
da83c923 772
773 Double_t pt = ph1->Pt() ;
774 Double_t ptcore = ph1->GetMomV2()->Pt() ;
ada6b882 775
da83c923 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) ;
2cde7444 782 if(fHaveTPCRP)
da83c923 783 FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCentBin),ptcore,dphiT) ;
ada6b882 784
da83c923 785 FillHistogram(Form("hPhotAll_cen%d",fCentBin),pt) ;
786 FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ptcore) ;
2cde7444 787 if(ph1->IsntUnfolded()){
da83c923 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) ;
2cde7444 793 }
794 if(ph1->IsCPVOK()){
da83c923 795 FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),pt,dphiA) ;
796 FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),pt,dphiC) ;
2cde7444 797 if(fHaveTPCRP)
da83c923 798 FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),pt,dphiT) ;
2cde7444 799
da83c923 800 FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ptcore,dphiA) ;
801 FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ptcore,dphiC) ;
2cde7444 802 if(fHaveTPCRP)
da83c923 803 FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ptcore,dphiT) ;
2cde7444 804
da83c923 805 FillHistogram(Form("hPhotCPV_cen%d",fCentBin),pt) ;
806 FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ptcore) ;
2cde7444 807 }
808 if(ph1->IsCPV2OK()){
da83c923 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) ;
2cde7444 820 }
821 if(ph1->IsDispOK()){
da83c923 822 FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),pt,dphiA) ;
823 FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),pt,dphiC) ;
2cde7444 824 if(fHaveTPCRP)
da83c923 825 FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),pt,dphiT) ;
2cde7444 826
da83c923 827 FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ptcore,dphiA) ;
828 FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ptcore,dphiC) ;
2cde7444 829 if(fHaveTPCRP)
da83c923 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) ;
ada6b882 837
2cde7444 838 }
da83c923 839 FillHistogram(Form("hPhotDisp_cen%d",fCentBin),pt) ;
840 FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ptcore) ;
2cde7444 841 if(ph1->IsntUnfolded()){
da83c923 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) ;
2cde7444 857 }
da83c923 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) ;
2cde7444 871 if(ph1->IsCPVOK()){
da83c923 872 FillHistogram(Form("hPhotPhiV0ABoth2_cen%d",fCentBin),pt,dphiA) ;
873 FillHistogram(Form("hPhotPhiV0CBoth2_cen%d",fCentBin),pt,dphiC) ;
2cde7444 874 if(fHaveTPCRP)
da83c923 875 FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCentBin),pt,dphiT) ;
2cde7444 876
da83c923 877 FillHistogram(Form("hPhotPhiV0ABoth2core_cen%d",fCentBin),ptcore,dphiA) ;
878 FillHistogram(Form("hPhotPhiV0CBoth2core_cen%d",fCentBin),ptcore,dphiC) ;
2cde7444 879 if(fHaveTPCRP)
da83c923 880 FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCentBin),ptcore,dphiT) ;
2cde7444 881
da83c923 882 FillHistogram(Form("hPhotBoth2_cen%d",fCentBin),pt) ;
d97448f6 883 FillHistogram(Form("hPhotBoth2core_cen%d",fCentBin),ptcore) ;
2cde7444 884 }
ada6b882 885 }
2cde7444 886 }
ada6b882 887}
888//_____________________________________________________________________________
889void AliAnalysisTaskPi0Flow::ConsiderPi0s()
890{
891 char key[55];
6ec40e8c 892 for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
ada6b882 893 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
6ec40e8c 894 for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
ada6b882 895 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
2cde7444 896 TLorentzVector p12 = *ph1 + *ph2;
ada6b882 897 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
898 FillHistogram("hPHOSphi",fCentralityV0M,p12.Pt(),p12.Phi());
2cde7444 899 Double_t dphiA=p12.Phi()-fRPV0A ;
900 while(dphiA<0)dphiA+=TMath::Pi() ;
901 while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
ada6b882 902
2cde7444 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() ;
ada6b882 910
2cde7444 911 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
f0532ae0 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) ;
2cde7444 923 if(fHaveTPCRP)
f0532ae0 924 FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 925
f0532ae0 926 FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
927 FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
2cde7444 928 if(fHaveTPCRP)
f0532ae0 929 FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
ada6b882 930
931
f0532ae0 932 FillHistogram(Form("hPi0All_cen%d",fCentBin),m,pt) ;
933 FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
2cde7444 934 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
f0532ae0 935 FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),m,pt) ;
da83c923 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) ;
2cde7444 940 }
ada6b882 941
f0532ae0 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) ;
2cde7444 946 if(ph1->IsntUnfolded())
f0532ae0 947 FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt1) ;
2cde7444 948 if(ph2->IsntUnfolded())
f0532ae0 949 FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt2) ;
2cde7444 950 if(ph1->IsCPVOK()){
f0532ae0 951 FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt1) ;
952 FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 953 }
954 if(ph2->IsCPVOK()){
f0532ae0 955 FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt2) ;
956 FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 957 }
958 if(ph1->IsCPV2OK()){
f0532ae0 959 FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt1) ;
960 FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 961 }
962 if(ph2->IsCPV2OK()){
f0532ae0 963 FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt2) ;
964 FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
ada6b882 965 }
2cde7444 966 if(ph1->IsDispOK()){
f0532ae0 967 FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt1) ;
2cde7444 968 if(ph1->IsntUnfolded()){
f0532ae0 969 FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt1) ;
2cde7444 970 }
f0532ae0 971 FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 972 }
973 if(ph2->IsDispOK()){
f0532ae0 974 FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt2) ;
2cde7444 975 if(ph1->IsntUnfolded()){
f0532ae0 976 FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt2) ;
2cde7444 977 }
f0532ae0 978 FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 979 }
980 if(ph1->IsDisp2OK()){
f0532ae0 981 FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt1) ;
982 FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 983 }
984 if(ph2->IsDisp2OK()){
f0532ae0 985 FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt2) ;
986 FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 987 }
988 if(ph1->IsDispOK() && ph1->IsCPVOK()){
f0532ae0 989 FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt1) ;
990 FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 991 }
992 if(ph2->IsDispOK() && ph2->IsCPVOK()){
f0532ae0 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) ;
ada6b882 1003 }
1004
1005
1006 if(a<kAlphaCut){
f0532ae0 1007 FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1008 }
1009
1010 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
ada6b882 1011 snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
f0532ae0 1012 FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
1013 FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1014 if(fHaveTPCRP)
f0532ae0 1015 FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 1016
f0532ae0 1017 FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1018 FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
2cde7444 1019 if(fHaveTPCRP)
f0532ae0 1020 FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
2cde7444 1021
f0532ae0 1022 FillHistogram(Form("hPi0CPV_cen%d",fCentBin),m,pt) ;
1023 FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
ada6b882 1024
1025 if(a<kAlphaCut){
f0532ae0 1026 FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1027 }
ada6b882 1028 }
2cde7444 1029 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
da83c923 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
f0532ae0 1039 FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),m,pt) ;
1040 FillHistogram(Form("hPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
ada6b882 1041 if(a<kAlphaCut){
f0532ae0 1042 FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1043 }
ada6b882 1044 }
2cde7444 1045 if(ph1->IsDispOK() && ph2->IsDispOK()){
ada6b882 1046 snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
f0532ae0 1047 FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1048 FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1049 if(fHaveTPCRP)
f0532ae0 1050 FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
da83c923 1051
f0532ae0 1052 FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1053 FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
2cde7444 1054 if(fHaveTPCRP)
f0532ae0 1055 FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
ada6b882 1056
f0532ae0 1057 FillHistogram(Form("hPi0Disp_cen%d",fCentBin),m,pt) ;
1058 FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),mcore, ptcore) ;
1059
2cde7444 1060 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
f0532ae0 1061 FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),m,pt) ;
da83c923 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) ;
2cde7444 1067 }
ada6b882 1068
1069 if(a<kAlphaCut){
f0532ae0 1070 FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1071 }
2cde7444 1072 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
f0532ae0 1073 FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1074 FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1075 if(fHaveTPCRP)
f0532ae0 1076 FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 1077
f0532ae0 1078 FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1079 FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
2cde7444 1080 if(fHaveTPCRP)
f0532ae0 1081 FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
ada6b882 1082
f0532ae0 1083 FillHistogram(Form("hPi0Both_cen%d",fCentBin),m,pt) ;
1084 FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),mcore,ptcore) ;
ada6b882 1085
1086 if(a<kAlphaCut){
1087 snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
f0532ae0 1088 FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1089 }
2cde7444 1090 if(ph1->Module()==1 && ph2->Module()==1)
f0532ae0 1091 FillHistogram("hPi0M11",m,pt );
2cde7444 1092 else if(ph1->Module()==2 && ph2->Module()==2)
f0532ae0 1093 FillHistogram("hPi0M22",m,pt );
2cde7444 1094 else if(ph1->Module()==3 && ph2->Module()==3)
f0532ae0 1095 FillHistogram("hPi0M33",m,pt );
2cde7444 1096 else if(ph1->Module()==1 && ph2->Module()==2)
f0532ae0 1097 FillHistogram("hPi0M12",m,pt );
2cde7444 1098 else if(ph1->Module()==1 && ph2->Module()==3)
f0532ae0 1099 FillHistogram("hPi0M13",m,pt );
2cde7444 1100 else if(ph1->Module()==2 && ph2->Module()==3)
f0532ae0 1101 FillHistogram("hPi0M23",m,pt );
2cde7444 1102
1103 }
f0532ae0 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
2cde7444 1137 }
1138 } // end of loop i2
1139 } // end of loop i1
ada6b882 1140}
1141//_____________________________________________________________________________
1142void AliAnalysisTaskPi0Flow::ConsiderPi0sMix()
1143{
1144 char key[55];
1145
1146 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1147
6ec40e8c 1148 for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
ada6b882 1149 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
89e3079e 1150 for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
ada6b882 1151 TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
2cde7444 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() ;
ada6b882 1168
1169
2cde7444 1170 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
f0532ae0 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
ada6b882 1180
1181 snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
f0532ae0 1182 FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA) ;
1183 FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1184 if(fHaveTPCRP)
f0532ae0 1185 FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 1186
f0532ae0 1187 FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),mcore, ptcore, dphiA) ;
1188 FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),mcore, ptcore, dphiC) ;
2cde7444 1189 if(fHaveTPCRP)
f0532ae0 1190 FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),mcore, ptcore, dphiT) ;
2cde7444 1191
f0532ae0 1192 FillHistogram(Form("hMiPi0All_cen%d",fCentBin),m,pt) ;
1193 FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
2cde7444 1194 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
f0532ae0 1195 FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),m,pt) ;
da83c923 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) ;
2cde7444 1200 }
1201
f0532ae0 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) ;
2cde7444 1206 if(ph1->IsntUnfolded())
f0532ae0 1207 FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt1) ;
2cde7444 1208 if(ph2->IsntUnfolded())
f0532ae0 1209 FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt2) ;
2cde7444 1210 if(ph1->IsCPVOK()){
f0532ae0 1211 FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt1) ;
1212 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 1213 }
1214 if(ph2->IsCPVOK()){
f0532ae0 1215 FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt2) ;
1216 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 1217 }
1218 if(ph1->IsCPV2OK()){
f0532ae0 1219 FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt1) ;
1220 FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 1221 }
1222 if(ph2->IsCPV2OK()){
f0532ae0 1223 FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt2) ;
1224 FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
ada6b882 1225 }
2cde7444 1226 if(ph1->IsDispOK()){
f0532ae0 1227 FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt1) ;
2cde7444 1228 if(ph1->IsntUnfolded()){
f0532ae0 1229 FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt1) ;
2cde7444 1230 }
f0532ae0 1231 FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
2cde7444 1232 }
1233 if(ph2->IsDispOK()){
f0532ae0 1234 FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt2) ;
2cde7444 1235 if(ph1->IsntUnfolded()){
f0532ae0 1236 FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt2) ;
2cde7444 1237 }
f0532ae0 1238 FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 1239 }
1240 if(ph1->IsDisp2OK()){
f0532ae0 1241 FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt1) ;
1242 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
ada6b882 1243 }
2cde7444 1244 if(ph2->IsDisp2OK()){
f0532ae0 1245 FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt2) ;
1246 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore2) ;
2cde7444 1247 }
1248 if(ph1->IsDispOK() && ph1->IsCPVOK()){
ada6b882 1249 snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
f0532ae0 1250 FillHistogram(key,m,pt1) ;
ada6b882 1251 snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
f0532ae0 1252 FillHistogram(key,mcore,ptcore1) ;
2cde7444 1253 }
1254 if(ph2->IsDispOK() && ph2->IsCPVOK()){
ada6b882 1255 snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
f0532ae0 1256 FillHistogram(key,m,pt2) ;
ada6b882 1257 snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
f0532ae0 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) ;
ada6b882 1267 }
1268
1269
1270
1271 if(a<kAlphaCut){
f0532ae0 1272 FillHistogram(Form("hMiPi0All_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1273 }
1274 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
f0532ae0 1275 FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
1276 FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1277 if(fHaveTPCRP)
f0532ae0 1278 FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 1279
f0532ae0 1280 FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1281 FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
2cde7444 1282 if(fHaveTPCRP)
f0532ae0 1283 FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
ada6b882 1284
f0532ae0 1285 FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),m,pt) ;
1286 FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
2cde7444 1287
ada6b882 1288 if(a<kAlphaCut){
f0532ae0 1289 FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1290 }
1291 }
1292 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
f0532ae0 1293 FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),m,pt) ;
1294 FillHistogram(Form("hMiPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
2cde7444 1295
da83c923 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
ada6b882 1305 if(a<kAlphaCut){
f0532ae0 1306 FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1307 }
1308 }
1309 if(ph1->IsDispOK() && ph2->IsDispOK()){
f0532ae0 1310 FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1311 FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1312 if(fHaveTPCRP)
f0532ae0 1313 FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
ada6b882 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) ;
2cde7444 1317 if(fHaveTPCRP)
ada6b882 1318 FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1319
2cde7444 1320
f0532ae0 1321 FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),m,pt) ;
ada6b882 1322 FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
2cde7444 1323 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
f0532ae0 1324 FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),m,pt) ;
da83c923 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) ;
2cde7444 1329 }
ada6b882 1330
1331 if(a<kAlphaCut){
f0532ae0 1332 FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),m,pt) ;
ada6b882 1333 }
2cde7444 1334 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
f0532ae0 1335 FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1336 FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
2cde7444 1337 if(fHaveTPCRP)
f0532ae0 1338 FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
2cde7444 1339
ada6b882 1340 FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1341 FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
2cde7444 1342 if(fHaveTPCRP)
ada6b882 1343 FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
2cde7444 1344
f0532ae0 1345 FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),m,pt) ;
ada6b882 1346 FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1347
1348 if(a<kAlphaCut){
f0532ae0 1349 FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),m,pt) ;
2cde7444 1350 }
1351 }
1352 }
f0532ae0 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 }
2cde7444 1385 } // end of loop i2
1386 }
1387 } // end of loop i1
2cde7444 1388}
ada6b882 1389//_____________________________________________________________________________
1390void AliAnalysisTaskPi0Flow::UpdateLists()
2cde7444 1391{
ada6b882 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
107849ac 1394
ada6b882 1395 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
cdd3a269 1396 if( fDebug >= 2 )
107849ac 1397 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
6ec40e8c 1398 if(fCaloPhotonsPHOS->GetEntriesFast()>0){
ada6b882 1399 arrayList->AddFirst(fCaloPhotonsPHOS) ;
1400 fCaloPhotonsPHOS=0;
107849ac 1401 if(arrayList->GetEntries() > fCentNMixed[fCentBin]){ // Remove redundant events
ada6b882 1402 TObjArray * tmp = static_cast<TObjArray*>(arrayList->Last()) ;
1403 arrayList->RemoveLast() ;
1404 delete tmp ; // TODO: may conflict with delete done by list being owner.
2cde7444 1405 }
2cde7444 1406 }
ada6b882 1407 else
1408 fCaloPhotonsPHOS->Clear(); // TODO: redundant???
2cde7444 1409}
1410//_____________________________________________________________________________
1411void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x)const{
1412 //FillHistogram
9ed0694f 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)) ;
2cde7444 1418}
1419//_____________________________________________________________________________
1420void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y)const{
1421 //FillHistogram
9ed0694f 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)) ;
2cde7444 1427}
1428
1429//_____________________________________________________________________________
1430void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1431 //Fills 1D histograms with key
9ed0694f 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;
2cde7444 1438 }
9ed0694f 1439
1440 TH3 * th3 = dynamic_cast<TH3*> (obj);
1441 if(th3) {
1442 th3->Fill(x, y, z) ;
1443 return;
2cde7444 1444 }
9ed0694f 1445
1446 AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
2cde7444 1447}
1448
45c6dc2d 1449void 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
ada6b882 1463//_____________________________________________________________________________
1464AliVEvent* 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
2cde7444 1475//___________________________________________________________________________
393743be 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// }
2cde7444 1487
ada6b882 1488//___________________________________________________________________________
1489Int_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] ) {
cdd3a269 1499 if( fDebug >= 1 )
f5e08ac9 1500 AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
ada6b882 1501 return lastBinUpperIndex-1;
1502 }
1503 if( centralityV0M < fCentEdges[0] ) {
cdd3a269 1504 if( fDebug >= 1 )
f5e08ac9 1505 AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
ada6b882 1506 return 0;
1507 }
1508
1509 fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1510 return fCentBin;
1511}
1512
1513//___________________________________________________________________________
1514Int_t AliAnalysisTaskPi0Flow::GetRPBin()
1515{
1516 Double_t averageRP;
1517 if(fHaveTPCRP)
3c8a190f 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.;
ada6b882 1520 else
2b90a13e 1521 averageRP = (fRPV0A+fRPV0C) /2.;
ada6b882 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
cdd3a269 1530 if ( fDebug >= 2 )
f5e08ac9 1531 AliInfo(Form("Event Mixing Reaction Plane bin is: %d", fEMRPBin));
1532
ada6b882 1533 return fEMRPBin;
1534}
1535
1536
1537//_____________________________________________________________________________
c9fb2807 1538void AliAnalysisTaskPi0Flow::LogProgress(int step)
ada6b882 1539{
cdd3a269 1540 if(fDebug >= 2) {
ada6b882 1541 AliInfo(Form("step %d completed", step));
1542 }
1543 // the +0.5 is not realy neccisarry, but oh well... -henrik
c9fb2807 1544 //FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1545 //FillHistogram("hTotSelEvents", step+0.5);
1546}
1547
1548void 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
ada6b882 1554 FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1555 FillHistogram("hTotSelEvents", step+0.5);
1556}
1557
c9fb2807 1558
ada6b882 1559//___________________________________________________________________________
1560Int_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 }
90d2ba9f 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 }
2b90a13e 2063 if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) ) {
ada6b882 2064 AliWarning("Period not defined");
2065 }
2066 return 1;
2cde7444 2067}
2068//_____________________________________________________________________________
2069Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
ada6b882 2070
2cde7444 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) ;
ada6b882 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 ;
2cde7444 2079 return (R2<2.5*2.5) ;
ada6b882 2080
2cde7444 2081}
2082//_____________________________________________________________________________
2083Bool_t AliAnalysisTaskPi0Flow::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
ada6b882 2084
2cde7444 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) ;
ada6b882 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 ;
2cde7444 2093 return (R2<1.5*1.5) ;
ada6b882 2094
2cde7444 2095}
ada6b882 2096//____________________________________________________________________________
2097TList* 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
2cde7444 2116//____________________________________________________________________________
2117Double_t AliAnalysisTaskPi0Flow::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
2118 //Parameterization of LHC10h period
2119 //_true if neutral_
ada6b882 2120
2cde7444 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)+
ada6b882 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);
2cde7444 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) ;
53c2a02a 2126 Double_t mf = 0.; //
2b90a13e 2127 if(fEventAOD) mf = fEventAOD->GetMagneticField(); //Positive for ++ and negative for --
2128 else if(fEventESD) mf = fEventESD->GetMagneticField(); //Positive for ++ and negative for --
2129
2cde7444 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
ada6b882 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)) ;
2cde7444 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//____________________________________________________________________________
b509e3a5 2151void 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())) ;
2cde7444 2158 }
b509e3a5 2159 else{
2160 AliInfo(Form("Setting PHOS flattening with name %s \n",maps->GetName())) ;
9d29cbcb 2161 AliEPFlattener * h = (AliEPFlattener*)maps->At(0) ;
b509e3a5 2162 if(fTPCFlat) delete fTPCFlat ;
9d29cbcb 2163 fTPCFlat = new AliEPFlattener() ;
b509e3a5 2164 fTPCFlat = h ;
9d29cbcb 2165 h = (AliEPFlattener*)maps->At(1) ;
b509e3a5 2166 if(fV0AFlat) delete fV0AFlat ;
9d29cbcb 2167 fV0AFlat = new AliEPFlattener() ;
b509e3a5 2168 fV0AFlat = h ;
9d29cbcb 2169 h = (AliEPFlattener*)maps->At(2) ;
b509e3a5 2170 if(fV0CFlat) delete fV0CFlat ;
9d29cbcb 2171 fV0CFlat = new AliEPFlattener() ;
b509e3a5 2172 fV0CFlat = h ;
2173 }
2174
2cde7444 2175}
b509e3a5 2176 //____________________________________________________________________________
2177Double_t AliAnalysisTaskPi0Flow::ApplyFlattening(Double_t phi, Double_t c){
2178
2179 if(fTPCFlat)
2180 return fTPCFlat->MakeFlat(phi,c);
2181 return phi ;
2cde7444 2182
2cde7444 2183}
b509e3a5 2184//____________________________________________________________________________
2185Double_t AliAnalysisTaskPi0Flow::ApplyFlatteningV0A(Double_t phi, Double_t c){
2186
2187 if(fV0AFlat)
2188 return fV0AFlat->MakeFlat(phi,c);
2189 return phi ;
2cde7444 2190
2191}
2192//____________________________________________________________________________
b509e3a5 2193Double_t AliAnalysisTaskPi0Flow::ApplyFlatteningV0C(Double_t phi, Double_t c){
2194
2195 if(fV0CFlat)
2196 return fV0CFlat->MakeFlat(phi,c);
ada6b882 2197 return phi ;
2198
2199}
2cde7444 2200//____________________________________________________________________________
91b112bd 2201Double_t AliAnalysisTaskPi0Flow::CoreEnergy(AliVCluster * clu, AliVCaloCells * cells)
2202{
2cde7444 2203 //calculate energy of the cluster in the circle with radius distanceCut around the maximum
ada6b882 2204
2cde7444 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 ;
ada6b882 2209
f0532ae0 2210 const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2cde7444 2211// Calculates the center of gravity in the local PHOS-module coordinates
2212 Float_t wtot = 0;
f0532ae0 2213 const Int_t mulDigit=clu->GetNCells() ;
2214 Double_t xc[mulDigit] ;
2215 Double_t zc[mulDigit] ;
2216 Double_t ei[mulDigit] ;
2cde7444 2217 Double_t x = 0 ;
2218 Double_t z = 0 ;
2cde7444 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 ;
f0532ae0 2227 ei[iDigit]=elist[iDigit]*cells->GetCellAmplitude(clu->GetCellsAbsId()[iDigit]);
cdd3a269 2228 if( fDebug >= 3 )
f0532ae0 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() ) ) ;
2cde7444 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)
f0532ae0 2245 coreE += ei[iDigit] ;
2cde7444 2246 }
2247 //Apply non-linearity correction
2248 return fNonLinCorr->Eval(coreE) ;
2249}
2250//____________________________________________________________________________
2251Bool_t AliAnalysisTaskPi0Flow::AreNeibors(Int_t id1,Int_t id2){
ada6b882 2252 // return true if absId are "Neighbors" (adjacent, including diagornaly,)
2253 // false if not.
2cde7444 2254
ada6b882 2255 Int_t relid1[4] ;
2256 fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
2cde7444 2257
ada6b882 2258 Int_t relid2[4] ;
2259 fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
2cde7444 2260
ada6b882 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 }
2cde7444 2270
ada6b882 2271 // else false
2272 return false;
2cde7444 2273}
2274//____________________________________________________________________________
ada6b882 2275void AliAnalysisTaskPi0Flow::Reclusterize(AliVCluster * clu){
2cde7444 2276 //Re-clusterize to make continues cluster
ada6b882 2277
2cde7444 2278 const Int_t oldMulDigit=clu->GetNCells() ;
ada6b882 2279 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2cde7444 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 }
ada6b882 2310
2cde7444 2311 if(inClu==oldMulDigit) //no need to modify
ada6b882 2312 return ;
2cde7444 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++){
ada6b882 2319 tmpD[i]=dlist[i] ;
2cde7444 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 }
ada6b882 2328
2329
2330}
2331
2332//_____________________________________________________________________________
2333void 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) {
f5e08ac9 2339 if( fDebug )
ada6b882 2340 AliInfo(Form("no PHOS Geometric Misalignment Matrix for module %d", mod));
2341 continue;
2342 }
2343 else {
2344 fPHOSGeo->SetMisalMatrix(modMatrix, mod);
f5e08ac9 2345 if( fDebug )
cdd3a269 2346 AliInfo(Form("PHOS Geometric Misalignment Matrix set for module %d", mod));
ada6b882 2347 }
2348 }
2349 }
2cde7444 2350}
2351
2352//_____________________________________________________________________________
ada6b882 2353void AliAnalysisTaskPi0Flow::SetV0Calibration(){
2354 // assigns: fMultV0, fV0Cpol, fV0Apol, fMeanQ, and fWidthQ
2355
6a1a4e92 2356 if ( ! fManualV0EPCalc ) {
2357 if( fDebug >=2 )
2358 AliInfo("Not setting V0Calibration, only needed for manual V0 EP Calculation");
2359 return;
2360 }
2361
ada6b882 2362 int runNumber = this->fRunNumber;
2363
2cde7444 2364 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
2365 TFile *foadb = TFile::Open(oadbfilename.Data());
2366
2367 if(!foadb){
ada6b882 2368 AliError(Form("OADB file %s cannot be opened\n", oadbfilename.Data()));
2369 AliError("V0 Calibration not set !\n");
2cde7444 2370 return;
2371 }
2372
2373 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
2374 if(!cont){
ada6b882 2375 AliError("OADB object hMultV0BefCorr is not available in the file");
2376 AliError("V0 Calibration not set!\n");
2377 return;
2cde7444 2378 }
2379
ada6b882 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;
2cde7444 2388 }
2cde7444 2389
f5e08ac9 2390 if( fDebug ) AliInfo("Setting V0 calibration") ;
ada6b882 2391 fMultV0 = ((TH2F *) cont->GetObject(runNumber))->ProfileX();
2392
2393 TF1 *fpol0 = new TF1("fpol0","pol0");
e2139d72 2394 fMultV0->Fit(fpol0,"Q0","",0,31);
2cde7444 2395 fV0Cpol = fpol0->GetParameter(0);
e2139d72 2396 fMultV0->Fit(fpol0,"Q0","",32,64);
2cde7444 2397 fV0Apol = fpol0->GetParameter(0);
2398
2399 for(Int_t iside=0;iside<2;iside++){
2400 for(Int_t icoord=0;icoord<2;icoord++){
ada6b882 2401 for(Int_t i=0;i < kNCenBins;i++){
2cde7444 2402 char namecont[100];
2403 if(iside==0 && icoord==0)
53c2a02a 2404 snprintf(namecont,100,"hQxc2_%i",i);
2cde7444 2405 else if(iside==1 && icoord==0)
53c2a02a 2406 snprintf(namecont,100,"hQxa2_%i",i);
2cde7444 2407 else if(iside==0 && icoord==1)
53c2a02a 2408 snprintf(namecont,100,"hQyc2_%i",i);
2cde7444 2409 else if(iside==1 && icoord==1)
53c2a02a 2410 snprintf(namecont,100,"hQya2_%i",i);
2cde7444 2411
2412 cont = (AliOADBContainer*) foadb->Get(namecont);
2413 if(!cont){
ada6b882 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;
2cde7444 2417 }
2cde7444 2418
ada6b882 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;
2cde7444 2422 }
ada6b882 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;
2cde7444 2427 }
ada6b882 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// }
2cde7444 2457// fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2458// fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2459
2460 }
2461 }
2462 }
ada6b882 2463
2464 delete fpol0; fpol0=0;
2cde7444 2465}
ada6b882 2466
2467//_____________________________________________________________________________
2468void 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//_____________________________________________________________________________
2479void AliAnalysisTaskPi0Flow::SetGeometry()
2480{
2481 // Initialize the PHOS geometry
f5e08ac9 2482 if( kLHC10h == fPeriod && fEventESD ) {
ada6b882 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)) {
cdd3a269 2497 if( fDebug )
ada6b882 2498 AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
b98b7f51 2499 continue;
ada6b882 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//_____________________________________________________________________________
2511void AliAnalysisTaskPi0Flow::SetPHOSCalibData()
2512{
2513 if( fPHOSCalibData )
2514 delete fPHOSCalibData;
2515 fPHOSCalibData = 0;
2cde7444 2516
ada6b882 2517 // Calibration only needed for ESD
2518 if( fEventESD /*&& */ ) {
f5e08ac9 2519 if( kLHC10h == fPeriod && fEventESD ) {
ada6b882 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//_____________________________________________________________________________
2531void 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);
f5e08ac9 2547
cdd3a269 2548 if( fDebug >= 2 )
f5e08ac9 2549 AliInfo(Form("Vertex is set to (%.1f,%.1f,%.1f)", fVertex[0], fVertex[1], fVertex[2]));
ada6b882 2550
2551 fVtxBin=0 ;// No support for vtx binning implemented.
2552}
2553
2554//_____________________________________________________________________________
2555Bool_t AliAnalysisTaskPi0Flow::RejectEventVertex()
2556{
2557 if( ! fEvent->GetPrimaryVertex() )
2558 return true; // reject
c9fb2807 2559 LogSelection(1, fInternalRunNumber);
2560
10a684bf 2561 if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ )
ada6b882 2562 return true; // reject
c9fb2807 2563 LogSelection(2, fInternalRunNumber);
ada6b882 2564
2565 return false; // accept event.
2566}
2567
2568//_____________________________________________________________________________
2569void 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);
f5e08ac9 2581
cdd3a269 2582 if ( fDebug >= 2 )
f5e08ac9 2583 AliInfo(Form("Centrality (bin) is: %f (%d)", fCentralityV0M, fCentBin));
ada6b882 2584}
2585
2586//_____________________________________________________________________________
2587Bool_t AliAnalysisTaskPi0Flow::RejectCentrality()
2588{
2589 if( ! fEvent->GetCentrality() )
2590 return true; // reject
c9fb2807 2591 LogSelection(3, fInternalRunNumber);
2592
ada6b882 2593// if( fCentralityV0M <= 0. || fCentralityV0M>80. )
2594// return true; // reject
2595
2596 int lastBinUpperIndex = fCentEdges.GetSize() -1;
f5e08ac9 2597 if( fCentralityV0M > fCentEdges[lastBinUpperIndex] ) {
2598 if( fDebug )
2599 AliInfo("Rejecting due to centrality outside of binning.");
ada6b882 2600 return true; // reject
f5e08ac9 2601 }
c9fb2807 2602 LogSelection(4, fInternalRunNumber);
2603
f5e08ac9 2604 if( fCentralityV0M < fCentEdges[0] ) {
2605 if( fDebug )
2606 AliInfo("Rejecting due to centrality outside of binning.");
ada6b882 2607 return true; // reject
f5e08ac9 2608 }
c9fb2807 2609 LogSelection(5, fInternalRunNumber);
ada6b882 2610
2611 return false;
2612}
2613
2614
2615//_____________________________________________________________________________
2616void AliAnalysisTaskPi0Flow::EvalReactionPlane()
2617{
2618 // assigns: fHaveTPCRP and fRP
2619 // also does a few histogram fills
2620
2621 AliEventplane *eventPlane = fEvent->GetEventplane();
f5e08ac9 2622 if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2623
ada6b882 2624 Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
2625 FillHistogram("phiRP",reactionPlaneQ,fCentralityV0M) ;
2626
2627 if(reactionPlaneQ==999 || reactionPlaneQ < 0.){ //reaction plain was not defined
f5e08ac9 2628 if( fDebug ) AliInfo(Form("No Q Reaction Plane, value is %f", reactionPlaneQ));
2629 fHaveTPCRP = kFALSE;
ada6b882 2630 }
2631 else{
cdd3a269 2632 if( fDebug >= 2 ) AliInfo(Form("Q Reaction Plane is %f", reactionPlaneQ));
f5e08ac9 2633 fHaveTPCRP = kTRUE;
ada6b882 2634 }
f5e08ac9 2635
ada6b882 2636 if(fHaveTPCRP){
e945b860 2637 fRP = ApplyFlattening(reactionPlaneQ, fCentralityV0M) ;
ada6b882 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 }
f5e08ac9 2645 else
2646 fRP=0.;
ada6b882 2647}
2648
2649
2cde7444 2650//____________________________________________________________________________
ada6b882 2651void AliAnalysisTaskPi0Flow::EvalV0ReactionPlane(){
2652 // set: fRPV0A and fRPV0C
2cde7444 2653
6a1a4e92 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 }
2cde7444 2675
6a1a4e92 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));
6b1c1977 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) ;
2cde7444 2732
f5e08ac9 2733 // Flattening
e945b860 2734 fRPV0A=ApplyFlatteningV0A(fRPV0A,fCentralityV0M) ;
6b1c1977 2735 while (fRPV0A<0 ) fRPV0A+=TMath::Pi() ;
2736 while (fRPV0A>TMath::Pi()) fRPV0A-=TMath::Pi() ;
2737
e945b860 2738 fRPV0C=ApplyFlatteningV0C(fRPV0C,fCentralityV0M) ;
6b1c1977 2739 while (fRPV0C<0 ) fRPV0C+=TMath::Pi() ;
2740 while (fRPV0C>TMath::Pi()) fRPV0C-=TMath::Pi() ;
e945b860 2741
cdd3a269 2742 if( fDebug >= 2 )
6b1c1977 2743 AliInfo(Form("V0 Reaction Plane after flattening: A side: %f, C side: %f", fRPV0A, fRPV0C));
ada6b882 2744
2745 FillHistogram("phiRPV0Aflat",fRPV0A,fCentralityV0M) ;
ada6b882 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 }
2cde7444 2757}
6ca8b4a1 2758//____________________________________________________________________________
2759void 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
f0532ae0 2762 const Double_t rCut=4.5 ;
6ca8b4a1 2763
f0532ae0 2764 const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
6ca8b4a1 2765// Calculates the center of gravity in the local PHOS-module coordinates
2766 Float_t wtot = 0;
f0532ae0 2767 const Int_t mulDigit=clu->GetNCells() ;
2768 Double_t xc[mulDigit] ;
2769 Double_t zc[mulDigit] ;
3c8a190f 2770 Double_t wi[mulDigit] ;
6ca8b4a1 2771 Double_t x = 0 ;
2772 Double_t z = 0 ;
6ca8b4a1 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 ;
3c8a190f 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];
6ca8b4a1 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++) {
3c8a190f 2805 Double_t w=wi[iDigit];
2806 if (w>0.) {
6ca8b4a1 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//____________________________________________________________________________
2839Bool_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}