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