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