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