afa351a4329373d5485c7668f1acd5432fbe8a04
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskPi0V2.cxx
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 "TH3F.h"
36 #include "TCanvas.h"
37 #include "TList.h"
38
39 #include "AliAnalysisTaskSE.h"
40 #include "AliAnalysisManager.h"
41 #include "AliStack.h"
42 #include "AliESDtrackCuts.h"
43 #include "AliESDEvent.h"
44 #include "AliESDInputHandler.h"
45 #include "AliAODEvent.h"
46 #include "AliMCEvent.h"
47
48 #include "AliEventplane.h"
49 #include "AliEMCALGeometry.h"
50 #include "THnSparse.h"
51
52 using std::cout;
53 using std::endl;
54
55 ClassImp(AliAnalysisTaskPi0V2)
56
57 //________________________________________________________________________
58 AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2(const char *name) // All data members should be initialised here
59    :AliAnalysisTaskSE(name),
60     fOutput(0),
61     fESD(0),
62     fEvtSelect(1),
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.),
68     hEvtCount(0), hAllcentV0(0), hAllcentV0r(0), hAllcentV0A(0), hAllcentV0C(0), hAllcentTPC(0),
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
76 {
77     // Dummy constructor ALWAYS needed for I/O.
78     DefineInput(0, TChain::Class());
79     DefineOutput(1, TList::Class());                                            // for output list
80 }
81
82 //________________________________________________________________________
83 AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2() // All data members should be initialised here
84    :AliAnalysisTaskSE("default_name"),
85     fOutput(0),
86     fESD(0),
87     fEvtSelect(1),
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.),
93     hEvtCount(0), hAllcentV0(0), hAllcentV0r(0), hAllcentV0A(0), hAllcentV0C(0), hAllcentTPC(0),
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)
100 {
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
107 }
108
109 //________________________________________________________________________
110 AliAnalysisTaskPi0V2::~AliAnalysisTaskPi0V2()
111 {
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.
114         delete fOutput;
115 }
116 //_____________________________________________________________________
117 Double_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 //_____________________________________________________________________
143 Double_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
168   Double_t crossEnergy = 0.;
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 //_____________________________________________________________________
193 Bool_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
209   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
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 //______________________________________________________________________
245 Bool_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); 
259      if((1. - double(GetCrossEnergy(c,id))/maxE) > 0.97)
260     return kFALSE;
261
262
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
277 //  if(c->M20 >)
278
279   return kTRUE;
280
281 }
282 //_____________________________________________________________________
283 Bool_t AliAnalysisTaskPi0V2::IsGoodPion(const TLorentzVector &p1, const TLorentzVector &p2) const
284 {
285   // Is good pion?
286
287
288 //  Double_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
289 //  if (asym>0.7)
290 //    return kFALSE;
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 //_______________________________________________________________________
306 void 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
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();
333
334   //cout<<"cos2V0: "<<cos2phiV0<<"  cos2V0A: "<<cos2phiV0A<<"  cos2V0C: "<<cos2phiV0C<<endl;
335   //cout<<mass<<"  "<<pt<<"  "<<phi<<"  "<<endl;
336   //cout<<" dphiV0: "<<dphiV0<<"    dphiV0A: "<<dphiV0A<<"  dphiV0C: "<<dphiV0C<<"+++++++"<<endl;
337
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);
369
370
371 }
372 //_________________________________________________________________________________________________
373 void 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);
394
395 }
396 //________________________________________________________________________
397 void AliAnalysisTaskPi0V2::UserCreateOutputObjects()
398 {
399     // Create histograms
400     // Called once (on the worker node)
401         
402     fOutput = new TList();
403     fOutput->SetOwner();  // IMPORTANT!
404
405     hEvtCount = new TH1F("hEvtCount", " Event Plane", 10, 0.5, 10.5);
406     hEvtCount->GetXaxis()->SetBinLabel(1,"SemiMB");
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);
413     
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());
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
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.);    
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
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.);
456     fOutput->Add(hdifV0A_V0C);
457     fOutput->Add(hdifV0A_TPC);
458     fOutput->Add(hdifV0C_TPC);
459     fOutput->Add(hdifV0C_V0A);
460
461
462
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.);
466     fOutput->Add(hdifEMC_EP);
467     fOutput->Add(hdifful_EP);
468     fOutput->Add(hdifout_EP);
469
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
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
495     
496
497     PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
498 }
499
500 //________________________________________________________________________
501 void AliAnalysisTaskPi0V2::UserExec(Option_t *) 
502 {
503     // Main loop
504     // Called for each event
505         
506     // Create pointer to reconstructed event
507    AliVEvent *event = InputEvent();
508    if (!event) { Printf("ERROR: Could not retrieve event"); return; }
509
510   Bool_t isSelected =0;      
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 )
517         return; 
518
519     // create pointer to event
520     fESD = dynamic_cast<AliESDEvent*>(event);
521     if (!fESD) {
522         AliError("Cannot get the ESD event");
523         return;
524     }  
525
526     hEvtCount->Fill(1);
527     
528     const AliESDVertex* fvertex = fESD->GetPrimaryVertex();
529     if(TMath::Abs(fvertex->GetZ())>10.)
530       return;
531     Double_t vertex[3] = {fvertex->GetX(), fvertex->GetY(), fvertex->GetZ()};
532
533     hEvtCount->Fill(2);
534
535     if(fESD->GetCentrality()) {
536       fCentrality = 
537         fESD->GetCentrality()->GetCentralityPercentile("V0M");
538     } else{
539            return;
540     }
541
542     hEvtCount->Fill(3);
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);
572
573     }
574 //cout<<" fEPV0:"<<fEPV0<<" fEPV0A:"<<fEPV0A<<" fEPV0C:"<<fEPV0C<<" fEPV0Ar:"<<fEPV0Ar<<" fEPV0Cr:"<<fEPV0Cr<<" fEPV0r:"<<fEPV0AR4<<" fEPV0AR7:"<<fEPV0AR7<<" fEPV0CR0:"<<fEPV0CR0<<" fEPV0CR3:"<<fEPV0CR3<<"--------------------------------------------"<<endl;
575     
576 //cout<<" EPTPC: "<<fEPTPC<<" resoTPC: "<<fEPTPCreso<<"---------------"<<endl;
577
578     hEvtCount->Fill(4);
579
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;
584
585     hEvtCount->Fill(5);
586
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
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;
600
601    if(fEPTPC != -999.)
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);
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);
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)));
627         
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
633    
634
635     Int_t nCluster =  fESD->GetNumberOfCaloClusters(); 
636     for(Int_t i=0; i<nCluster; ++i){
637       AliESDCaloCluster *c1 = fESD->GetCaloCluster(i);
638       if(!c1->IsEMCAL()) continue;
639       if(!IsGoodCluster(c1)) continue;
640       for(Int_t j=i+1; j<nCluster; ++j){
641         AliESDCaloCluster *c2 = fESD->GetCaloCluster(j);
642         if(!c2->IsEMCAL()) continue;
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
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();
661
662         Double_t difTrack = TVector2::Phi_0_2pi(tPhi-fEPV0);  if(difTrack >TMath::Pi()) difTrack -= TMath::Pi();
663         if(esdtrack->IsEMCAL()){        
664           hdifEMC_EP->Fill(fCentrality, difTrack, tPt);
665         }else{
666           hdifout_EP->Fill(fCentrality, difTrack, tPt);
667         }
668         hdifful_EP->Fill(fCentrality,   difTrack, tPt);
669     }
670     hEvtCount->Fill(6);
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);
675 }
676
677
678 //________________________________________________________________________
679 void AliAnalysisTaskPi0V2::Terminate(Option_t *) 
680 {
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
695 }