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