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