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