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