]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskPi0V2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added name argument to AddTaskPHOSPi0Flow.C
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskPi0V2.cxx
CommitLineData
04b116e8 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
3c40321c 16/* $Id: AliAnalysisTaskPi0V2.cxx 55404 2012-03-29 10:10:19Z fca $ */
17
04b116e8 18/* AliAnalysisTaskPi0V2.cxx
19 *
20 * Template task producing a P_t spectrum and pseudorapidity distribution.
21 * Includes explanations of physics and primary track selections
22 *
23 * Instructions for adding histograms can be found below, starting with NEW HISTO
24 *
25 * Based on tutorial example from offline pages
26 * Edited by Arvinder Palaha
27 */
3c40321c 28#include "AliAnalysisTaskPi0V2.h"
29
04b116e8 30#include "Riostream.h"
31#include "TChain.h"
32#include "TTree.h"
33#include "TH1F.h"
34#include "TH2F.h"
35#include "TH3F.h"
36#include "TCanvas.h"
37#include "TList.h"
3c40321c 38
70d53162 39#include "AliAnalysisTaskSE.h"
04b116e8 40#include "AliAnalysisManager.h"
41#include "AliStack.h"
42#include "AliESDtrackCuts.h"
3c40321c 43#include "AliESDEvent.h"
44#include "AliESDInputHandler.h"
04b116e8 45#include "AliAODEvent.h"
70d53162 46#include "AliMCEvent.h"
04b116e8 47
48#include "AliEventplane.h"
49#include "AliEMCALGeometry.h"
50#include "THnSparse.h"
3c40321c 51
9d6804f4 52using std::cout;
53using std::endl;
54
3c40321c 55ClassImp(AliAnalysisTaskPi0V2)
56
57//________________________________________________________________________
965c985f 58AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2(const char *name) // All data members should be initialised here
59 :AliAnalysisTaskSE(name),
04b116e8 60 fOutput(0),
04b116e8 61 fESD(0),
670ffa5c 62 fTrackCuts(0),
8071d5b2 63 fEvtSelect(1),
670ffa5c 64 fVtxCut(10.),
65 fNcellCut(2), fECut(1), fEtaCut(0.65), fM02Cut(0.5), fPi0AsyCut(0),
04b116e8 66 fCentrality(99.),
67 fEPTPC(-999.),
68 fEPTPCreso(0.),
69 fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
70 fEPV0AR4(-999.), fEPV0AR5(-999.), fEPV0AR6(-999.), fEPV0AR7(-999.), fEPV0CR0(-999.), fEPV0CR1(-999.), fEPV0CR2(-999.), fEPV0CR3(-999.),
ef7e23cf 71 hEvtCount(0), hAllcentV0(0), hAllcentV0r(0), hAllcentV0A(0), hAllcentV0C(0), hAllcentTPC(0),
04b116e8 72 hEPTPC(0), hresoTPC(0),
73 hEPV0(0), hEPV0A(0), hEPV0C(0), hEPV0Ar(0), hEPV0Cr(0), hEPV0r(0), hEPV0AR4(0), hEPV0AR7(0), hEPV0CR0(0), hEPV0CR3(0),
74 hdifV0A_V0CR0(0), hdifV0A_V0CR3(0), hdifV0ACR0_V0CR3(0), hdifV0C_V0AR4(0), hdifV0C_V0AR7(0), hdifV0AR4_V0AR7(0),
75 hdifV0A_V0C(0), hdifV0A_TPC(0), hdifV0C_TPC(0), hdifV0C_V0A(0),
93df010a 76 hdifEMC_EPV0(0), hdifEMC_EPV0A(0), hdifEMC_EPV0C(0), hdifful_EPV0(0), hdifful_EPV0A(0), hdifful_EPV0C(0),
77 hdifout_EPV0(0), hdifout_EPV0A(0), hdifout_EPV0C(0), hdifEMC_EPTPC(0), hdifful_EPTPC(0), hdifout_EPTPC(0),
04b116e8 78 fHEPV0r(0), fHEPV0A(0), fHEPV0C(0), fHEPTPC(0)
79
3c40321c 80{
04b116e8 81 // Dummy constructor ALWAYS needed for I/O.
670ffa5c 82 fTrackCuts = new AliESDtrackCuts();
04b116e8 83 DefineInput(0, TChain::Class());
84 DefineOutput(1, TList::Class()); // for output list
3c40321c 85}
86
87//________________________________________________________________________
965c985f 88AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2() // All data members should be initialised here
8f40bd27 89 :AliAnalysisTaskSE("default_name"),
04b116e8 90 fOutput(0),
04b116e8 91 fESD(0),
670ffa5c 92 fTrackCuts(0),
8071d5b2 93 fEvtSelect(1),
670ffa5c 94 fVtxCut(10.),
95 fNcellCut(2), fECut(1), fEtaCut(0.65), fM02Cut(0.5), fPi0AsyCut(0),
04b116e8 96 fCentrality(99.),
97 fEPTPC(-999.),
98 fEPTPCreso(0.),
99 fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
100 fEPV0AR4(-999.), fEPV0AR5(-999.), fEPV0AR6(-999.), fEPV0AR7(-999.), fEPV0CR0(-999.), fEPV0CR1(-999.), fEPV0CR2(-999.), fEPV0CR3(-999.),
ef7e23cf 101 hEvtCount(0), hAllcentV0(0), hAllcentV0r(0), hAllcentV0A(0), hAllcentV0C(0), hAllcentTPC(0),
04b116e8 102 hEPTPC(0), hresoTPC(0),
103 hEPV0(0), hEPV0A(0), hEPV0C(0), hEPV0Ar(0), hEPV0Cr(0), hEPV0r(0), hEPV0AR4(0), hEPV0AR7(0), hEPV0CR0(0), hEPV0CR3(0),
104 hdifV0A_V0CR0(0), hdifV0A_V0CR3(0), hdifV0ACR0_V0CR3(0), hdifV0C_V0AR4(0), hdifV0C_V0AR7(0), hdifV0AR4_V0AR7(0),
105 hdifV0A_V0C(0), hdifV0A_TPC(0), hdifV0C_TPC(0), hdifV0C_V0A(0),
93df010a 106 hdifEMC_EPV0(0), hdifEMC_EPV0A(0), hdifEMC_EPV0C(0), hdifful_EPV0(0), hdifful_EPV0A(0), hdifful_EPV0C(0),
107 hdifout_EPV0(0), hdifout_EPV0A(0), hdifout_EPV0C(0), hdifEMC_EPTPC(0), hdifful_EPTPC(0), hdifout_EPTPC(0),
04b116e8 108 fHEPV0r(0), fHEPV0A(0), fHEPV0C(0), fHEPTPC(0)
3c40321c 109{
04b116e8 110 // Constructor
111 // Define input and output slots here (never in the dummy constructor)
112 // Input slot #0 works with a TChain - it is connected to the default input container
113 // Output slot #1 writes into a TH1 container
670ffa5c 114 fTrackCuts = new AliESDtrackCuts();
04b116e8 115 DefineInput(0, TChain::Class());
116 DefineOutput(1, TList::Class()); // for output list
3c40321c 117}
118
119//________________________________________________________________________
120AliAnalysisTaskPi0V2::~AliAnalysisTaskPi0V2()
121{
04b116e8 122 // Destructor. Clean-up the output list, but not the histograms that are put inside
123 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
670ffa5c 124 delete fTrackCuts;
125 delete fOutput;
3c40321c 126}
127//_____________________________________________________________________
128Double_t AliAnalysisTaskPi0V2::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
129{
130 // Get maximum energy of attached cell.
131
132 id = -1;
133
134 AliVCaloCells *cells = 0;
135 if (fESD)
136 cells = fESD->GetEMCALCells();
3c40321c 137 if (!cells)
138 return 0;
139
140 Double_t maxe = 0;
141 const Int_t ncells = cluster->GetNCells();
142 for (Int_t i=0; i<ncells; i++) {
143 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
144 if (e>maxe) {
145 maxe = e;
146 id = cluster->GetCellAbsId(i);
147 }
148 }
149 return maxe;
150}
151//_____________________________________________________________________
152Double_t AliAnalysisTaskPi0V2::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) const
153{
154 // Calculate the energy of cross cells around the leading cell.
155
156 AliVCaloCells *cells = 0;
157 if (fESD)
158 cells = fESD->GetEMCALCells();
3c40321c 159 if (!cells)
160 return 0;
161
162 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
163 if (!geom)
164 return 0;
165
166 Int_t iSupMod = -1;
167 Int_t iTower = -1;
168 Int_t iIphi = -1;
169 Int_t iIeta = -1;
170 Int_t iphi = -1;
171 Int_t ieta = -1;
172 Int_t iphis = -1;
173 Int_t ietas = -1;
174
965c985f 175 Double_t crossEnergy = 0.;
3c40321c 176
177 geom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
178 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
179
180 Int_t ncells = cluster->GetNCells();
181 for (Int_t i=0; i<ncells; i++) {
182 Int_t cellAbsId = cluster->GetCellAbsId(i);
183 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
184 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
185 Int_t aphidiff = TMath::Abs(iphi-iphis);
186 if (aphidiff>1)
187 continue;
188 Int_t aetadiff = TMath::Abs(ieta-ietas);
189 if (aetadiff>1)
190 continue;
191 if ( (aphidiff==1 && aetadiff==0) ||
192 (aphidiff==0 && aetadiff==1) ) {
193 crossEnergy += cells->GetCellAmplitude(cellAbsId);
194 }
195 }
196
197 return crossEnergy;
198}
199//_____________________________________________________________________
200Bool_t AliAnalysisTaskPi0V2::IsWithinFiducialVolume(Short_t id) const
201{
202 // Check if cell is within given fiducial volume.
203
204 Double_t fNFiducial = 1;
205
206 Int_t iSupMod = -1;
207 Int_t iTower = -1;
208 Int_t iIphi = -1;
209 Int_t iIeta = -1;
210 Int_t iphi = -1;
211 Int_t ieta = -1;
212
213 Bool_t okrow = kFALSE;
214 Bool_t okcol = kFALSE;
215
04b116e8 216 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
3c40321c 217 if (!geom)
218 return kFALSE;
219
220 Int_t cellAbsId = id;
221 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
222 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
223
224 // Check rows/phi
225 if (iSupMod < 10) {
226 if (iphi >= fNFiducial && iphi < 24-fNFiducial)
227 okrow = kTRUE;
228 } else {
229 if (iphi >= fNFiducial && iphi < 12-fNFiducial)
230 okrow = kTRUE;
231 }
232 // Check columns/eta
233 Bool_t noEMCALBorderAtEta0 = kTRUE;
234 if (!noEMCALBorderAtEta0) {
235 if (ieta > fNFiducial && ieta < 48-fNFiducial)
236 okcol = kTRUE;
237 } else {
238 if (iSupMod%2==0) {
239 if (ieta >= fNFiducial)
240 okcol = kTRUE;
241 } else {
242 if (ieta < 48-fNFiducial)
243 okcol = kTRUE;
244 }
245 }
246 if (okrow && okcol)
247 return kTRUE;
248
249 return kFALSE;
250}
251//______________________________________________________________________
252Bool_t AliAnalysisTaskPi0V2::IsGoodCluster(const AliESDCaloCluster *c) const
253{
254
255 if(!c)
256 return kFALSE;
257
670ffa5c 258 if(c->GetNCells() < fNcellCut)
3c40321c 259 return kFALSE;
260
670ffa5c 261 if(c->E() < fECut)
3c40321c 262 return kFALSE;
263
264 Short_t id = -1;
265 Double_t maxE = GetMaxCellEnergy(c, id);
965c985f 266 if((1. - double(GetCrossEnergy(c,id))/maxE) > 0.97)
3c40321c 267 return kFALSE;
268
965c985f 269
3c40321c 270 Float_t pos1[3];
271 c->GetPosition(pos1);
272 TVector3 clsPos(pos1);
273 Double_t eta = clsPos.Eta();
274
670ffa5c 275 if(TMath::Abs(eta) > fEtaCut)
3c40321c 276 return kFALSE;
277
278 if (!IsWithinFiducialVolume(id))
279 return kFALSE;
280
670ffa5c 281 if(c->GetM02() >fM02Cut)
3c40321c 282 return kFALSE;
283
04b116e8 284// if(c->M20 >)
285
3c40321c 286 return kTRUE;
70d53162 287
04b116e8 288}
3c40321c 289//_____________________________________________________________________
290Bool_t AliAnalysisTaskPi0V2::IsGoodPion(const TLorentzVector &p1, const TLorentzVector &p2) const
291{
292 // Is good pion?
293
670ffa5c 294 if(fPi0AsyCut){
295 Double_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
296 if (asym>0.7)
297 return kFALSE;
298 }
3c40321c 299 TLorentzVector pion;
300 pion = p1 + p2;
301 Double_t eta = pion.Eta();
670ffa5c 302 if(TMath::Abs(eta) > fEtaCut)
3c40321c 303 return kFALSE;
304
305 return kTRUE;
306}
307//_______________________________________________________________________
308void AliAnalysisTaskPi0V2::FillPion(const TLorentzVector& p1, const TLorentzVector& p2, Double_t EPV0r, Double_t EPV0A, Double_t EPV0C, Double_t EPTPC)
309{
310 // Fill histogram.
311
312 if (!IsGoodPion(p1,p2))
313 return;
314 TLorentzVector pion;
315 pion = p1 + p2;
316
317 Double_t mass = pion.M();
318 Double_t pt = pion.Pt();
319 Double_t phi = pion.Phi();
320
321 Double_t dphiV0 = phi-EPV0r;
322 Double_t dphiV0A = phi-EPV0A;
323 Double_t dphiV0C = phi-EPV0C;
324 Double_t dphiTPC = phi-EPTPC;
325
326 Double_t cos2phiV0 = TMath::Cos(2.*(dphiV0));
327 Double_t cos2phiV0A = TMath::Cos(2.*(dphiV0A));
328 Double_t cos2phiV0C = TMath::Cos(2.*(dphiV0C));
329 Double_t cos2phiTPC = TMath::Cos(2.*(dphiTPC));
330
c6ee6f73 331 dphiV0 = TVector2::Phi_0_2pi(dphiV0); if(dphiV0 >TMath::Pi()) dphiV0 -= TMath::Pi();
332 dphiV0A = TVector2::Phi_0_2pi(dphiV0A); if(dphiV0A >TMath::Pi()) dphiV0A -= TMath::Pi();
333 dphiV0C = TVector2::Phi_0_2pi(dphiV0C); if(dphiV0C >TMath::Pi()) dphiV0C -= TMath::Pi();
334 dphiTPC = TVector2::Phi_0_2pi(dphiTPC); if(dphiTPC >TMath::Pi()) dphiTPC -= TMath::Pi();
3c40321c 335
336 Double_t xV0[5]; // Match ndims in fH V0 EP
337 xV0[0] = mass;
338 xV0[1] = pt;
339 xV0[2] = fCentrality;
340 xV0[3] = dphiV0;
341 xV0[4] = cos2phiV0;
342 fHEPV0r->Fill(xV0);
343
344 Double_t xV0A[5]; // Match ndims in fH V0A EP
345 xV0A[0] = mass;
346 xV0A[1] = pt;
347 xV0A[2] = fCentrality;
348 xV0A[3] = dphiV0A;
349 xV0A[4] = cos2phiV0A;
350 fHEPV0A->Fill(xV0A);
351
352 Double_t xV0C[5]; // Match ndims in fH V0C EP
353 xV0C[0] = mass;
354 xV0C[1] = pt;
355 xV0C[2] = fCentrality;
356 xV0C[3] = dphiV0C;
357 xV0C[4] = cos2phiV0C;
358 fHEPV0C->Fill(xV0C);
359
360 Double_t xTPC[5]; // Match ndims in fH TPC EP
361 xTPC[0] = mass;
362 xTPC[1] = pt;
363 xTPC[2] = fCentrality;
364 xTPC[3] = dphiTPC;
365 xTPC[4] = cos2phiTPC;
366 fHEPTPC->Fill(xTPC);
70d53162 367
04b116e8 368
369}
3c40321c 370//_________________________________________________________________________________________________
371void AliAnalysisTaskPi0V2::GetMom(TLorentzVector& p, const AliESDCaloCluster *c, Double_t *vertex)
372{
373 // Calculate momentum.
374 Float_t posMom[3];
375 c->GetPosition(posMom);
376 TVector3 clsPos2(posMom);
377
378 Double_t e = c->E();
379 Double_t r = clsPos2.Perp();
380 Double_t eta = clsPos2.Eta();
381 Double_t phi = clsPos2.Phi();
382
383 TVector3 pos;
384 pos.SetPtEtaPhi(r,eta,phi);
385
386 if (vertex) { //calculate direction relative to vertex
387 pos -= vertex;
388 }
389
390 Double_t rad = pos.Mag();
391 p.SetPxPyPzE(e*pos.x()/rad, e*pos.y()/rad, e*pos.z()/rad, e);
70d53162 392
04b116e8 393}
3c40321c 394//________________________________________________________________________
395void AliAnalysisTaskPi0V2::UserCreateOutputObjects()
396{
04b116e8 397 // Create histograms
398 // Called once (on the worker node)
3c40321c 399
04b116e8 400 fOutput = new TList();
401 fOutput->SetOwner(); // IMPORTANT!
ef7e23cf 402
403 hEvtCount = new TH1F("hEvtCount", " Event Plane", 10, 0.5, 10.5);
0c842c93 404 hEvtCount->GetXaxis()->SetBinLabel(1,"SemiMB");
ef7e23cf 405 hEvtCount->GetXaxis()->SetBinLabel(2,"vert");
406 hEvtCount->GetXaxis()->SetBinLabel(3,"cent");
407 hEvtCount->GetXaxis()->SetBinLabel(4,"EPtask");
408 hEvtCount->GetXaxis()->SetBinLabel(5,"EPvalue");
409 hEvtCount->GetXaxis()->SetBinLabel(6,"Pass");
410 fOutput->Add(hEvtCount);
04b116e8 411
ebaf288d 412 hEPTPC = new TH2F("hEPTPC", "EPTPC vs cent", 100, 0., 100., 100, 0., TMath::Pi());
413 hresoTPC = new TH2F("hresoTPC", "TPc reso vs cent", 100, 0., 100., 100, 0., 1.);
414 hEPV0 = new TH2F("hEPV0", "EPV0 vs cent", 100, 0., 100., 100, 0., TMath::Pi());
415 hEPV0A = new TH2F("hEPV0A", "EPV0A vs cent", 100, 0., 100., 100, 0., TMath::Pi());
416 hEPV0C = new TH2F("hEPV0C", "EPV0C vs cent", 100, 0., 100., 100, 0., TMath::Pi());
417 hEPV0Ar = new TH2F("hEPV0Ar", "EPV0Ar vs cent", 100, 0., 100., 100, 0., TMath::Pi());
418 hEPV0Cr = new TH2F("hEPV0Cr", "EPV0Cr vs cent", 100, 0., 100., 100, 0., TMath::Pi());
419 hEPV0r = new TH2F("hEPV0r", "EPV0r vs cent", 100, 0., 100., 100, 0., TMath::Pi());
420 hEPV0AR4 = new TH2F("hEPV0AR4", "EPV0AR4 vs cent", 100, 0., 100., 100, 0., TMath::Pi());
421 hEPV0AR7 = new TH2F("hEPV0AR7", "EPV0AR7 vs cent", 100, 0., 100., 100, 0., TMath::Pi());
422 hEPV0CR0 = new TH2F("hEPV0CR0", "EPV0CR0 vs cent", 100, 0., 100., 100, 0., TMath::Pi());
423 hEPV0CR3 = new TH2F("hEPV0CR3", "EPV0CR3 vs cent", 100, 0., 100., 100, 0., TMath::Pi());
04b116e8 424 fOutput->Add(hEPTPC);
425 fOutput->Add(hresoTPC);
426 fOutput->Add(hEPV0);
427 fOutput->Add(hEPV0A);
428 fOutput->Add(hEPV0C);
429 fOutput->Add(hEPV0Ar);
430 fOutput->Add(hEPV0Cr);
431 fOutput->Add(hEPV0r);
432 fOutput->Add(hEPV0AR4);
433 fOutput->Add(hEPV0AR7);
434 fOutput->Add(hEPV0CR0);
435 fOutput->Add(hEPV0CR3);
436
ebaf288d 437 hdifV0A_V0CR0 = new TH2F("hdifV0A_V0CR0", "EP A-R0 ", 100, 0., 100., 100, -1., 1.);
438 hdifV0A_V0CR3 = new TH2F("hdifV0A_V0CR3", "EP A-R3 ", 100, 0., 100., 100, -1., 1.);
439 hdifV0ACR0_V0CR3 = new TH2F("hdifV0ACR0_V0CR3", "EP R0-R3 ", 100, 0., 100., 100, -1., 1.);
440 hdifV0C_V0AR4 = new TH2F("hdifV0C_V0AR4", "EP C-R4 ", 100, 0., 100., 100, -1., 1.);
441 hdifV0C_V0AR7 = new TH2F("hdifV0C_V0AR7", "EP C-R7 ", 100, 0., 100., 100, -1., 1.);
442 hdifV0AR4_V0AR7 = new TH2F("hdifV0AR4_V0AR7", "EP R4-R7 ", 100, 0., 100., 100, -1., 1.);
04b116e8 443 fOutput->Add(hdifV0A_V0CR0);
444 fOutput->Add(hdifV0A_V0CR3);
445 fOutput->Add(hdifV0ACR0_V0CR3);
446 fOutput->Add(hdifV0C_V0AR4);
447 fOutput->Add(hdifV0C_V0AR7);
448 fOutput->Add(hdifV0AR4_V0AR7);
449
ebaf288d 450 hdifV0A_V0C = new TH2F("hdifV0A_V0C", "EP A-C ", 100, 0., 100., 100, -1., 1.);
451 hdifV0A_TPC = new TH2F("hdifV0A_TPC", "EP A-TPC", 100, 0., 100., 100, -1., 1.);
452 hdifV0C_TPC = new TH2F("hdifV0C_TPC", "EP C-TPC", 100, 0., 100., 100, -1., 1.);
453 hdifV0C_V0A = new TH2F("hdifV0C_V0A", "EP C-A ", 100, 0., 100., 100, -1., 1.);
04b116e8 454 fOutput->Add(hdifV0A_V0C);
455 fOutput->Add(hdifV0A_TPC);
456 fOutput->Add(hdifV0C_TPC);
457 fOutput->Add(hdifV0C_V0A);
458
459
460
93df010a 461 hdifEMC_EPV0 = new TH3F("hdifEMC_EPV0", "dif phi in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
462 hdifEMC_EPV0A = new TH3F("hdifEMC_EPV0A", "dif phi in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
463 hdifEMC_EPV0C = new TH3F("hdifEMC_EPV0C", "dif phi in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
464 fOutput->Add(hdifEMC_EPV0);
465 fOutput->Add(hdifEMC_EPV0A);
466 fOutput->Add(hdifEMC_EPV0C);
467
468 hdifful_EPV0 = new TH3F("hdifful_EPV0", "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
469 hdifful_EPV0A = new TH3F("hdifful_EPV0A", "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
470 hdifful_EPV0C = new TH3F("hdifful_EPV0C", "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
471 fOutput->Add(hdifful_EPV0);
472 fOutput->Add(hdifful_EPV0A);
473 fOutput->Add(hdifful_EPV0C);
474
475 hdifout_EPV0 = new TH3F("hdifout_EPV0", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
476 hdifout_EPV0A = new TH3F("hdifout_EPV0A", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
477 hdifout_EPV0C = new TH3F("hdifout_EPV0C", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
478 fOutput->Add(hdifout_EPV0);
479 fOutput->Add(hdifout_EPV0A);
480 fOutput->Add(hdifout_EPV0C);
481
482 hdifEMC_EPTPC = new TH3F("hdifEMC_EPTPC", "dif phi in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
483 hdifful_EPTPC = new TH3F("hdifful_EPTPC", "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
484 hdifout_EPTPC = new TH3F("hdifout_EPTPC", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
485 fOutput->Add(hdifEMC_EPTPC);
486 fOutput->Add(hdifful_EPTPC);
487 fOutput->Add(hdifout_EPTPC);
04b116e8 488
ebaf288d 489 hAllcentV0 = new TH1F("hAllcentV0", "All cent EP V0", 100, 0., TMath::Pi());
490 hAllcentV0r = new TH1F("hAllcentV0r", "All cent EP V0r", 100, 0., TMath::Pi());
491 hAllcentV0A = new TH1F("hAllcentV0A", "All cent EP V0A", 100, 0., TMath::Pi());
492 hAllcentV0C = new TH1F("hAllcentV0C", "All cent EP V0C", 100, 0., TMath::Pi());
493 hAllcentTPC = new TH1F("hAllcentTPC", "All cent EP TPC", 100, 0., TMath::Pi());
494 fOutput->Add(hAllcentV0);
495 fOutput->Add(hAllcentV0r);
496 fOutput->Add(hAllcentV0A);
497 fOutput->Add(hAllcentV0C);
498 fOutput->Add(hAllcentTPC);
499
04b116e8 500 const Int_t ndims = 5;
501 Int_t nMgg=500, nPt=40, nCent=20, nDeltaPhi=315, ncos2phi=500;
502 Int_t bins[ndims] = {nMgg, nPt, nCent, nDeltaPhi, ncos2phi};
503 Double_t xmin[ndims] = { 0, 0., 0, 0., -1.};
504 Double_t xmax[ndims] = { 0.5, 20., 100, 3.15, 1.};
505 fHEPV0r = new THnSparseF("fHEPV0r", "Flow histogram EPV0", ndims, bins, xmin, xmax);
506 fHEPV0A = new THnSparseF("fHEPV0A", "Flow histogram EPV0A", ndims, bins, xmin, xmax);
507 fHEPV0C = new THnSparseF("fHEPV0C", "Flow histogram EPV0C", ndims, bins, xmin, xmax);
508 fHEPTPC = new THnSparseF("fHEPTPC", "Flow histogram EPTPC", ndims, bins, xmin, xmax);
509 fOutput->Add(fHEPV0r);
510 fOutput->Add(fHEPV0A);
511 fOutput->Add(fHEPV0C);
512 fOutput->Add(fHEPTPC);
513
3c40321c 514
04b116e8 515
516 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
3c40321c 517}
518
519//________________________________________________________________________
520void AliAnalysisTaskPi0V2::UserExec(Option_t *)
521{
04b116e8 522 // Main loop
523 // Called for each event
524
04b116e8 525 // Create pointer to reconstructed event
526 AliVEvent *event = InputEvent();
527 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
528
0c842c93 529 Bool_t isSelected =0;
670ffa5c 530 if(fEvtSelect == 1){ //MB+SemiCentral
8071d5b2 531 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral));
670ffa5c 532 } else if (fEvtSelect == 2){ //MB+Central+SemiCentral
8071d5b2 533 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral | AliVEvent::kCentral));
670ffa5c 534 } else if(fEvtSelect == 3){ //MB
535 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB ));
8071d5b2 536 }
537 if(!isSelected )
0c842c93 538 return;
539
04b116e8 540 // create pointer to event
541 fESD = dynamic_cast<AliESDEvent*>(event);
542 if (!fESD) {
543 AliError("Cannot get the ESD event");
544 return;
545 }
04b116e8 546
ef7e23cf 547 hEvtCount->Fill(1);
548
549 const AliESDVertex* fvertex = fESD->GetPrimaryVertex();
670ffa5c 550 if(TMath::Abs(fvertex->GetZ())>fVtxCut)
c6ee6f73 551 return;
04b116e8 552 Double_t vertex[3] = {fvertex->GetX(), fvertex->GetY(), fvertex->GetZ()};
70d53162 553
ef7e23cf 554 hEvtCount->Fill(2);
555
04b116e8 556 if(fESD->GetCentrality()) {
557 fCentrality =
670ffa5c 558 fESD->GetCentrality()->GetCentralityPercentile("CL1"); //spd vertex
ef7e23cf 559 } else{
560 return;
04b116e8 561 }
3c40321c 562
ef7e23cf 563 hEvtCount->Fill(3);
04b116e8 564 AliEventplane *ep = fESD->GetEventplane();
565 if (ep) {
566 if (ep->GetQVector())
567 fEPTPC = ep->GetQVector()->Phi()/2. ;
568 else
569 fEPTPC = -999.;
570 if (ep->GetQsub1()&&ep->GetQsub2())
571 fEPTPCreso = TMath::Cos(2.*(ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.));
572 else
573 fEPTPCreso = -1;
574
575 fEPV0 = ep->GetEventplane("V0", fESD);
576 fEPV0A = ep->GetEventplane("V0A", fESD);
577 fEPV0C = ep->GetEventplane("V0C", fESD);
578 Double_t qx=0, qy=0;
579 Double_t qxr=0, qyr=0;
580 fEPV0Ar = ep->CalculateVZEROEventPlane(fESD, 4, 5, 2, qxr, qyr);
581 fEPV0Cr = ep->CalculateVZEROEventPlane(fESD, 2, 3, 2, qx, qy);
582 qxr += qx;
583 qyr += qy;
584 fEPV0r = TMath::ATan2(qyr,qxr)/2.;
585 fEPV0AR4 = ep->CalculateVZEROEventPlane(fESD, 4, 2, qx, qy);
586 fEPV0AR5 = ep->CalculateVZEROEventPlane(fESD, 5, 2, qx, qy);
587 fEPV0AR6 = ep->CalculateVZEROEventPlane(fESD, 6, 2, qx, qy);
588 fEPV0AR7 = ep->CalculateVZEROEventPlane(fESD, 7, 2, qx, qy);
589 fEPV0CR0 = ep->CalculateVZEROEventPlane(fESD, 0, 2, qx, qy);
590 fEPV0CR1 = ep->CalculateVZEROEventPlane(fESD, 1, 2, qx, qy);
591 fEPV0CR2 = ep->CalculateVZEROEventPlane(fESD, 2, 2, qx, qy);
592 fEPV0CR3 = ep->CalculateVZEROEventPlane(fESD, 3, 2, qx, qy);
3c40321c 593
04b116e8 594 }
965c985f 595
ef7e23cf 596 hEvtCount->Fill(4);
597
670ffa5c 598 if( fEPV0A<-2. || fEPV0C<-2. || fEPV0AR4<-2.
599 || fEPV0AR7<-2. || fEPV0CR0<-2. || fEPV0CR3<-2.
600 || fEPTPC<-2. || fEPV0r<-2. || fEPV0Ar<-2.
601 || fEPV0Cr<-2.) return;
ef7e23cf 602
603 hEvtCount->Fill(5);
604
ebaf288d 605 fEPV0 = TVector2::Phi_0_2pi(fEPV0); if(fEPV0>TMath::Pi()) fEPV0 = fEPV0 - TMath::Pi();
606 fEPV0r = TVector2::Phi_0_2pi(fEPV0r); if(fEPV0r>TMath::Pi()) fEPV0r = fEPV0r - TMath::Pi();
607 fEPV0A = TVector2::Phi_0_2pi(fEPV0A); if(fEPV0A>TMath::Pi()) fEPV0A = fEPV0A - TMath::Pi();
608 fEPV0C = TVector2::Phi_0_2pi(fEPV0C); if(fEPV0C>TMath::Pi()) fEPV0C = fEPV0C - TMath::Pi();
609 fEPV0Ar = TVector2::Phi_0_2pi(fEPV0Ar); if(fEPV0Ar>TMath::Pi()) fEPV0Ar = fEPV0Ar - TMath::Pi();
610 fEPV0Cr = TVector2::Phi_0_2pi(fEPV0Cr); if(fEPV0Cr>TMath::Pi()) fEPV0Cr = fEPV0Cr - TMath::Pi();
611 fEPV0AR4 = TVector2::Phi_0_2pi(fEPV0AR4); if(fEPV0AR4>TMath::Pi()) fEPV0AR4 = fEPV0AR4 - TMath::Pi();
612 fEPV0AR7 = TVector2::Phi_0_2pi(fEPV0AR7); if(fEPV0AR7>TMath::Pi()) fEPV0AR7 = fEPV0AR7 - TMath::Pi();
613 fEPV0CR0 = TVector2::Phi_0_2pi(fEPV0CR0); if(fEPV0CR0>TMath::Pi()) fEPV0CR0 = fEPV0CR0 - TMath::Pi();
614 fEPV0CR3 = TVector2::Phi_0_2pi(fEPV0CR3); if(fEPV0CR3>TMath::Pi()) fEPV0CR3 = fEPV0CR3 - TMath::Pi();
615
8071d5b2 616 if(fEPTPC != -999.)
04b116e8 617 hEPTPC->Fill(fCentrality, fEPTPC);
618 if(fEPTPCreso!=-1) hresoTPC->Fill(fCentrality, fEPTPCreso);
619 hEPV0->Fill(fCentrality, fEPV0);
620 hEPV0A->Fill(fCentrality, fEPV0A);
621 hEPV0C->Fill(fCentrality, fEPV0C);
622 hEPV0Ar->Fill(fCentrality, fEPV0Ar);
623 hEPV0Cr->Fill(fCentrality, fEPV0Cr);
624 hEPV0r->Fill(fCentrality, fEPV0r);
ebaf288d 625 hEPV0AR4->Fill(fCentrality, fEPV0AR4);
626 hEPV0AR7->Fill(fCentrality, fEPV0AR7);
627 hEPV0CR0->Fill(fCentrality, fEPV0CR0);
628 hEPV0CR3->Fill(fCentrality, fEPV0CR3);
629
630 hAllcentV0->Fill(fEPV0);
631 hAllcentV0r->Fill(fEPV0r);
632 hAllcentV0A->Fill(fEPV0A);
633 hAllcentV0C->Fill(fEPV0C);
634 hAllcentTPC->Fill(fEPTPC);
04b116e8 635
636 hdifV0A_V0CR0->Fill(fCentrality, TMath::Cos(2.*(fEPV0A - fEPV0CR0)));
637 hdifV0A_V0CR3->Fill(fCentrality, TMath::Cos(2.*(fEPV0A - fEPV0CR3)));
638 hdifV0ACR0_V0CR3->Fill(fCentrality, TMath::Cos(2*(fEPV0CR0 - fEPV0CR3)));
639 hdifV0C_V0AR4->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0AR4)));
640 hdifV0C_V0AR7->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0AR7)));
641 hdifV0AR4_V0AR7->Fill(fCentrality, TMath::Cos(2*(fEPV0AR4 - fEPV0AR7)));
3c40321c 642
04b116e8 643 hdifV0A_V0C->Fill(fCentrality, TMath::Cos(2*(fEPV0A - fEPV0C)));
644 hdifV0A_TPC->Fill(fCentrality, TMath::Cos(2*(fEPV0A - fEPTPC)));
645 hdifV0C_TPC->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPTPC)));
646 hdifV0C_V0A->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0A)));
647 // Cluster loop for reconstructed event
04b116e8 648
670ffa5c 649 Int_t nCluster = fESD->GetNumberOfCaloClusters();
650 for(Int_t i=0; i<nCluster; ++i){
651 AliESDCaloCluster *c1 = fESD->GetCaloCluster(i);
652 if(!c1->IsEMCAL()) continue;
653 if(!IsGoodCluster(c1)) continue;
654 for(Int_t j=i+1; j<nCluster; ++j){
655 AliESDCaloCluster *c2 = fESD->GetCaloCluster(j);
656 if(!c2->IsEMCAL()) continue;
657 if(!IsGoodCluster(c2)) continue;
658 TLorentzVector p1;
659 GetMom(p1, c1, vertex);
660 TLorentzVector p2;
661 GetMom(p2, c2, vertex);
662 FillPion(p1, p2, fEPV0r, fEPV0A, fEPV0C, fEPTPC);
663 }
664 }
665
666 //for track analysis.
667 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
668 fTrackCuts->SetRequireTPCRefit(kTRUE);
669 fTrackCuts->SetRequireITSRefit(kTRUE);
670 fTrackCuts->SetEtaRange(-0.7,0.7);
671 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
672 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
673 fTrackCuts->SetMinNClustersTPC(100);
674
675 Int_t nTrack = fESD->GetNumberOfTracks();
676 for(Int_t i=0; i<nTrack; ++i){
677 AliESDtrack* esdtrack = fESD->GetTrack(i); // pointer to reconstructed to track
678 if(!fTrackCuts->AcceptTrack(esdtrack))
679 continue;
680 if(!esdtrack) {
681 AliError(Form("ERROR: Could not retrieve esdtrack %d",i));
04b116e8 682 continue;
683 }
684 Double_t tPhi = esdtrack->Phi();
685 Double_t tPt = esdtrack->Pt();
ebaf288d 686
93df010a 687 Double_t difTrackV0 = TVector2::Phi_0_2pi(tPhi-fEPV0); if(difTrackV0 >TMath::Pi()) difTrackV0 -= TMath::Pi();
688 Double_t difTrackV0A = TVector2::Phi_0_2pi(tPhi-fEPV0A); if(difTrackV0A >TMath::Pi()) difTrackV0A -= TMath::Pi();
689 Double_t difTrackV0C = TVector2::Phi_0_2pi(tPhi-fEPV0C); if(difTrackV0C >TMath::Pi()) difTrackV0C -= TMath::Pi();
690 Double_t difTrackTPC = TVector2::Phi_0_2pi(tPhi-fEPTPC); if(difTrackTPC >TMath::Pi()) difTrackTPC -= TMath::Pi();
04b116e8 691 if(esdtrack->IsEMCAL()){
93df010a 692 hdifEMC_EPV0->Fill(fCentrality, difTrackV0, tPt);
693 hdifEMC_EPV0A->Fill(fCentrality, difTrackV0A, tPt);
694 hdifEMC_EPV0C->Fill(fCentrality, difTrackV0C, tPt);
695 hdifEMC_EPTPC->Fill(fCentrality, difTrackTPC, tPt);
04b116e8 696 }else{
93df010a 697 hdifout_EPV0->Fill(fCentrality, difTrackV0, tPt);
698 hdifout_EPV0A->Fill(fCentrality, difTrackV0A, tPt);
699 hdifout_EPV0C->Fill(fCentrality, difTrackV0C, tPt);
700 hdifout_EPTPC->Fill(fCentrality, difTrackTPC, tPt);
04b116e8 701 }
93df010a 702 hdifful_EPV0->Fill(fCentrality, difTrackV0, tPt);
703 hdifful_EPV0A->Fill(fCentrality, difTrackV0A, tPt);
704 hdifful_EPV0C->Fill(fCentrality, difTrackV0C, tPt);
705 hdifful_EPTPC->Fill(fCentrality, difTrackTPC, tPt);
04b116e8 706 }
ef7e23cf 707 hEvtCount->Fill(6);
04b116e8 708
709 // NEW HISTO should be filled before this point, as PostData puts the
710 // information for this iteration of the UserExec in the container
711 PostData(1, fOutput);
3c40321c 712}
713
04b116e8 714
3c40321c 715//________________________________________________________________________
716void AliAnalysisTaskPi0V2::Terminate(Option_t *)
717{
04b116e8 718 // Draw result to screen, or perform fitting, normalizations
719 // Called once at the end of the query
720// fOutput = dynamic_cast<TList*> (GetOutputData(1));
721 // if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
722
723 // Get the physics selection histograms with the selection statistics
724 //AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
725 //AliESDInputHandler *inputH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
726 //TH2F *histStat = (TH2F*)inputH->GetStatistics();
727
728
729 // NEW HISTO should be retrieved from the TList container in the above way,
730 // so it is available to draw on a canvas such as below
731
3c40321c 732}