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