Getting rid of removed class AliPHOSEPFlattener
[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 "TProfile.h"
37 #include "TCanvas.h"
38 #include "TList.h"
39
40 #include "AliAnalysisTaskSE.h"
41 #include "AliAnalysisManager.h"
42 #include "AliStack.h"
43 #include "AliESDtrackCuts.h"
44 #include "AliESDEvent.h"
45 #include "AliESDInputHandler.h"
46 #include "AliAODEvent.h"
47 #include "AliMCEvent.h"
48 #include "AliVCluster.h"
49
50 #include "AliEventplane.h"
51 #include "AliEMCALGeometry.h"
52 #include "THnSparse.h"
53 #include "TClonesArray.h"
54 #include "TString.h"
55 #include "AliCaloPID.h"
56 #include "AliCalorimeterUtils.h"
57 #include "AliCaloTrackReader.h"
58 #include "AliEPFlattener.h"
59 #include "AliOADBContainer.h"
60
61 using std::cout;
62 using std::endl;
63
64 ClassImp(AliAnalysisTaskPi0V2)
65
66 //________________________________________________________________________
67 AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2(const char *name) // All data members should be initialised here
68    :AliAnalysisTaskSE(name),
69     fOutput(0),
70     fESD(0),fAOD(0),
71     fTracksName("PicoTrack"), fV1ClusName("CaloCluster"), fV2ClusName("CaloCluster"),
72     fTrigClass("CVLN_|CSEMI_|CCENT|CVHN"), type("AOD"),
73     fTracks(0), fV1Clus(0), fV2Clus(0),
74     fRunNumber(-999),fInterRunNumber(-999),
75     fVtxCut(15.),
76     fNcellCut(2.), fECut(1.), fEtaCut(0.65), fM02Cut(0.5),fDrCut(0.025), fPi0AsyCut(0), isV1Clus(1), isPhosCali(0),
77     fCentrality(99.),
78     fEPTPC(-999.),
79     fEPTPCreso(0.), 
80     fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
81     fEPV0AR4(-999.), fEPV0AR5(-999.), fEPV0AR6(-999.), fEPV0AR7(-999.), fEPV0CR0(-999.), fEPV0CR1(-999.), fEPV0CR2(-999.), fEPV0CR3(-999.),
82     hEvtCount(0), 
83     h2DcosV0A(0), h2DsinV0A(0), h2DcosV0C(0), h2DsinV0C(0), h2DcosTPC(0), h2DsinTPC(0), 
84     hEPTPC(0), hresoTPC(0),
85     hEPV0(0), hEPV0A(0), hEPV0C(0), hEPV0Ar(0), hEPV0Cr(0), hEPV0r(0), hEPV0AR4(0), hEPV0AR7(0), hEPV0CR0(0), hEPV0CR3(0),
86     hEPTPCCor(0), hEPV0ACor(0), hEPV0CCor(0),
87     hdifV0Ar_V0Cr(0), hdifV0A_V0CR0(0), hdifV0A_V0CR3(0), hdifV0ACR0_V0CR3(0), hdifV0C_V0AR4(0), hdifV0C_V0AR7(0), hdifV0AR4_V0AR7(0),
88     hdifV0A_V0C(0), hdifV0A_TPC(0), hdifV0C_TPC(0), hdifV0C_V0A(0), 
89     hM02vsPtA(0), hM02vsPtB(0), hClusDxDZA(0), hClusDxDZB(0),
90     hdifEMC_EPV0(0), hdifEMC_EPV0A(0), hdifEMC_EPV0C(0), hdifful_EPV0(0), hdifful_EPV0A(0), hdifful_EPV0C(0), 
91     hdifout_EPV0(0), hdifout_EPV0A(0), hdifout_EPV0C(0), 
92     fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"), fTPCFlat(0x0), fV0AFlat(0x0),  fV0CFlat(0x0),
93     fClusterPbV0(0), fClusterPbV0A(0), fClusterPbV0C(0), fClusterPbTPC(0),    
94     fHEPV0r(0x0), fHEPV0A(0x0), fHEPV0C(0x0), fHEPTPC(0x0)
95
96 {
97     // Dummy constructor ALWAYS needed for I/O.
98     DefineInput(0, TChain::Class());
99     DefineOutput(1, TList::Class());                                            // for output list
100 }
101
102 //________________________________________________________________________
103 AliAnalysisTaskPi0V2::AliAnalysisTaskPi0V2() // All data members should be initialised here
104    :AliAnalysisTaskSE("default_name"),
105     fOutput(0),
106     fESD(0),fAOD(0),
107     fTracksName("PicoTrack"), fV1ClusName("CaloCluster"), fV2ClusName("CaloCluster"),
108     fTrigClass("CVLN_|CSEMI_|CCENT|CVHN"), type("AOD"),
109     fTracks(0), fV1Clus(0), fV2Clus(0),
110     fRunNumber(-999),fInterRunNumber(-999),
111     fVtxCut(15.),
112     fNcellCut(2.), fECut(1.), fEtaCut(0.65), fM02Cut(0.5), fDrCut(0.025), fPi0AsyCut(0), isV1Clus(1),isPhosCali(0),
113     fCentrality(99.),
114     fEPTPC(-999.),
115     fEPTPCreso(0.),
116     fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
117     fEPV0AR4(-999.), fEPV0AR5(-999.), fEPV0AR6(-999.), fEPV0AR7(-999.), fEPV0CR0(-999.), fEPV0CR1(-999.), fEPV0CR2(-999.), fEPV0CR3(-999.),
118     hEvtCount(0), 
119     h2DcosV0A(0), h2DsinV0A(0), h2DcosV0C(0), h2DsinV0C(0), h2DcosTPC(0), h2DsinTPC(0),
120     hEPTPC(0), hresoTPC(0),
121     hEPV0(0), hEPV0A(0), hEPV0C(0), hEPV0Ar(0), hEPV0Cr(0), hEPV0r(0), hEPV0AR4(0), hEPV0AR7(0), hEPV0CR0(0), hEPV0CR3(0),
122     hEPTPCCor(0), hEPV0ACor(0), hEPV0CCor(0),
123     hdifV0Ar_V0Cr(0), hdifV0A_V0CR0(0), hdifV0A_V0CR3(0), hdifV0ACR0_V0CR3(0), hdifV0C_V0AR4(0), hdifV0C_V0AR7(0), hdifV0AR4_V0AR7(0),
124     hdifV0A_V0C(0), hdifV0A_TPC(0), hdifV0C_TPC(0), hdifV0C_V0A(0),
125     hM02vsPtA(0), hM02vsPtB(0), hClusDxDZA(0), hClusDxDZB(0),
126     hdifEMC_EPV0(0), hdifEMC_EPV0A(0), hdifEMC_EPV0C(0), hdifful_EPV0(0), hdifful_EPV0A(0), hdifful_EPV0C(0),
127     hdifout_EPV0(0), hdifout_EPV0A(0), hdifout_EPV0C(0), 
128     fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"), fTPCFlat(0x0), fV0AFlat(0x0),  fV0CFlat(0x0),
129     fClusterPbV0(0), fClusterPbV0A(0), fClusterPbV0C(0), fClusterPbTPC(0),    
130     fHEPV0r(0x0), fHEPV0A(0x0), fHEPV0C(0x0), fHEPTPC(0x0)
131 {
132     // Constructor
133     // Define input and output slots here (never in the dummy constructor)
134     // Input slot #0 works with a TChain - it is connected to the default input container
135     // Output slot #1 writes into a TH1 container
136     DefineInput(0, TChain::Class());
137     DefineOutput(1, TList::Class());                                            // for output list
138 }
139
140 //________________________________________________________________________
141 AliAnalysisTaskPi0V2::~AliAnalysisTaskPi0V2()
142 {
143     // Destructor. Clean-up the output list, but not the histograms that are put inside
144     // (the list is owner and will clean-up these histograms). Protect in PROOF case.
145      if(fTPCFlat)delete fTPCFlat;  fTPCFlat=0x0;
146      if(fV0AFlat)delete fV0AFlat;  fV0AFlat=0x0;
147      if(fV0CFlat)delete fV0CFlat;  fV0CFlat=0x0;
148      delete fOutput;
149 }
150 //_____________________________________________________________________
151 Double_t AliAnalysisTaskPi0V2::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
152 {
153   // Get maximum energy of attached cell.
154
155   id = -1;
156
157   AliVCaloCells *cells = 0;
158   if (type=="ESD"){
159     cells = fESD->GetEMCALCells();
160   } else if(type=="AOD"){
161     cells = fAOD->GetEMCALCells();
162   }
163   if (!cells)
164     return 0;
165
166   Double_t maxe = 0;
167   const Int_t ncells = cluster->GetNCells();
168   for (Int_t i=0; i<ncells; i++) {
169     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
170     if (e>maxe) {
171       maxe = e;
172       id   = cluster->GetCellAbsId(i);
173     }
174   }
175   return maxe;
176 }
177 //_____________________________________________________________________
178 Double_t AliAnalysisTaskPi0V2::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) const
179 {
180   // Calculate the energy of cross cells around the leading cell.
181
182   AliVCaloCells *cells;
183   if (type=="ESD"){
184     cells = fESD->GetEMCALCells();
185   } else if(type=="AOD"){
186     cells = fAOD->GetEMCALCells();
187   }
188   if (!cells)
189     return 0;
190
191   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
192   if (!geom)
193     return 0;
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   Int_t iphis   = -1;
202   Int_t ietas   = -1;
203
204   Double_t crossEnergy = 0.;
205
206   geom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
207   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
208
209   Int_t ncells = cluster->GetNCells();
210   for (Int_t i=0; i<ncells; i++) {
211     Int_t cellAbsId = cluster->GetCellAbsId(i);
212     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
213     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
214     Int_t aphidiff = TMath::Abs(iphi-iphis);
215     if (aphidiff>1)
216       continue;
217     Int_t aetadiff = TMath::Abs(ieta-ietas);
218     if (aetadiff>1)
219       continue;
220     if ( (aphidiff==1 && aetadiff==0) ||
221         (aphidiff==0 && aetadiff==1) ) {
222       crossEnergy += cells->GetCellAmplitude(cellAbsId);
223     }
224   }
225
226   return crossEnergy;
227 }
228 //_____________________________________________________________________
229 Bool_t AliAnalysisTaskPi0V2::IsWithinFiducialVolume(Short_t id) const
230 {
231   // Check if cell is within given fiducial volume.
232
233   Double_t fNFiducial = 1;
234
235   Int_t iSupMod = -1;
236   Int_t iTower  = -1;
237   Int_t iIphi   = -1;
238   Int_t iIeta   = -1;
239   Int_t iphi    = -1;
240   Int_t ieta    = -1;
241
242   Bool_t okrow = kFALSE;
243   Bool_t okcol = kFALSE;
244
245   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
246   if (!geom)
247     return kFALSE;
248
249   Int_t cellAbsId = id;
250   geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
251   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
252
253   // Check rows/phi
254   if (iSupMod < 10) {
255     if (iphi >= fNFiducial && iphi < 24-fNFiducial)
256       okrow = kTRUE;
257   } else {
258     if (iphi >= fNFiducial && iphi < 12-fNFiducial)
259       okrow = kTRUE;
260   }
261   // Check columns/eta
262   Bool_t noEMCALBorderAtEta0 = kTRUE;
263   if (!noEMCALBorderAtEta0) {
264     if (ieta > fNFiducial && ieta < 48-fNFiducial)
265       okcol = kTRUE;
266   } else {
267     if (iSupMod%2==0) {
268       if (ieta >= fNFiducial)
269         okcol = kTRUE;
270     } else {
271       if (ieta < 48-fNFiducial)
272         okcol = kTRUE;
273     }
274   }
275   if (okrow && okcol)
276      return kTRUE;
277
278   return kFALSE;
279 }
280 //______________________________________________________________________
281 Bool_t AliAnalysisTaskPi0V2::IsGoodCluster(const AliVCluster *c) const
282 {
283   if(!c)
284     return kFALSE;
285
286   if(c->GetNCells() < fNcellCut)
287    return kFALSE;
288
289   if(c->E() < fECut)
290    return kFALSE;
291   Short_t id = -1;
292   Double_t maxE = GetMaxCellEnergy(c, id); 
293      if((1. - double(GetCrossEnergy(c,id))/maxE) > 0.97)
294     return kFALSE;
295
296
297   Float_t pos1[3];
298   c->GetPosition(pos1);
299   TVector3 clsPos(pos1);
300   Double_t eta = clsPos.Eta();
301
302   if(TMath::Abs(eta) > fEtaCut)
303     return kFALSE;  
304
305   if (!IsWithinFiducialVolume(id))
306     return kFALSE;
307
308   if(c->GetM02() >fM02Cut)
309     return kFALSE;
310
311
312   return kTRUE;
313
314 }
315 //________________________________________________________________________________________________
316 Bool_t AliAnalysisTaskPi0V2::IsGoodClusterV1(const AliVCluster *c) const
317 {
318
319   if(!c)
320     return kFALSE;
321
322   if(c->GetNCells() < fNcellCut)
323    return kFALSE;
324
325   if(c->E() < fECut)
326    return kFALSE;
327
328   Short_t id = -1;
329   Double_t maxE = GetMaxCellEnergy(c, id);
330   if((1. - double(GetCrossEnergy(c,id))/maxE) > 0.97)
331     return kFALSE;
332
333
334   Float_t pos1[3];
335   c->GetPosition(pos1);
336   TVector3 clsPos(pos1);
337   Double_t eta = clsPos.Eta();
338
339   if(TMath::Abs(eta) > fEtaCut)
340     return kFALSE;
341
342   if (!IsWithinFiducialVolume(id))
343     return kFALSE;
344
345   if(c->GetM02() <fM02Cut)
346     return kFALSE;
347
348   Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
349   if(dr<fDrCut)
350     return kFALSE;
351
352   return kTRUE;
353
354 }
355 //_____________________________________________________________________
356 Bool_t AliAnalysisTaskPi0V2::IsGoodPion(const TLorentzVector &p1, const TLorentzVector &p2) const
357 {
358   // Is good pion?
359
360   if(fPi0AsyCut){
361     Double_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
362     if (asym>0.7)
363       return kFALSE;
364   }
365   TLorentzVector pion;
366   pion = p1 + p2;
367   Double_t eta = pion.Eta();
368   if(TMath::Abs(eta) > fEtaCut)
369     return kFALSE;
370
371   return kTRUE;
372 }
373 //_______________________________________________________________________
374 void AliAnalysisTaskPi0V2::FillPion(const TLorentzVector& p1, const TLorentzVector& p2, Double_t EPV0A, Double_t EPV0C, Double_t EPTPC)
375 {
376   // Fill histogram.
377
378   if (!IsGoodPion(p1,p2))
379     return;
380   TLorentzVector pion;
381   pion = p1 + p2;
382
383   Double_t mass = pion.M();
384   Double_t pt   = pion.Pt();
385   Double_t phi  = pion.Phi();
386
387   Double_t dphiV0A  = phi-EPV0A;
388   Double_t dphiV0C  = phi-EPV0C;
389   Double_t dphiTPC  = phi-EPTPC;
390
391   dphiV0A = TVector2::Phi_0_2pi(dphiV0A); if(dphiV0A >TMath::Pi())  dphiV0A -= TMath::Pi();
392   dphiV0C = TVector2::Phi_0_2pi(dphiV0C); if(dphiV0C >TMath::Pi())  dphiV0C -= TMath::Pi();
393   dphiTPC = TVector2::Phi_0_2pi(dphiTPC); if(dphiTPC >TMath::Pi())  dphiTPC -= TMath::Pi();
394
395   Double_t xV0A[4]; // Match ndims in fH V0A EP
396   xV0A[0]       = mass;
397   xV0A[1]       = pt;
398   xV0A[2]       = fCentrality;
399   xV0A[3]       = dphiV0A;
400   fHEPV0A->Fill(xV0A);
401
402   Double_t xV0C[4]; // Match ndims in fH V0C EP
403   xV0C[0]       = mass;
404   xV0C[1]       = pt;
405   xV0C[2]       = fCentrality;
406   xV0C[3]       = dphiV0C;
407   fHEPV0C->Fill(xV0C);
408
409   Double_t xTPC[4]; // Match ndims in fH TPC EP
410   xTPC[0]       = mass;
411   xTPC[1]       = pt;
412   xTPC[2]       = fCentrality;
413   xTPC[3]       = dphiTPC;
414   fHEPTPC->Fill(xTPC);
415
416 }
417 //________________________________________________________________________________________________________________________________
418 void AliAnalysisTaskPi0V2::FillCluster(const TLorentzVector& p1, Double_t EPV0A, Double_t EPV0C, Double_t EPTPC, AliVCluster *c)
419 {
420   //cluster(photon) v2 method
421 //  Double_t Pt   = p1.Pt();
422   Double_t Et   = p1.Et();
423   Double_t Phi  = p1.Phi();
424   Double_t M02  = c->GetM02();
425   Double_t DxClus = c->GetTrackDx();
426   Double_t DzClus = c->GetTrackDz();
427   Double_t dr = TMath::Sqrt(DxClus*DxClus + DzClus*DzClus);
428
429   Double_t difClusV0A = TVector2::Phi_0_2pi(Phi-EPV0A);  if(difClusV0A >TMath::Pi()) difClusV0A -= TMath::Pi();
430   Double_t difClusV0C = TVector2::Phi_0_2pi(Phi-EPV0C);  if(difClusV0C >TMath::Pi()) difClusV0C -= TMath::Pi();
431   Double_t difClusTPC = TVector2::Phi_0_2pi(Phi-EPTPC);  if(difClusTPC >TMath::Pi()) difClusTPC -= TMath::Pi();
432
433   Double_t DataV0A[5];
434   DataV0A[0] = Et;
435   DataV0A[1] = M02;
436   DataV0A[2] = fCentrality;
437   DataV0A[3] = difClusV0A;
438   DataV0A[4] = dr;
439   fClusterPbV0A->Fill(DataV0A);
440
441   Double_t DataV0C[5];
442   DataV0C[0] = Et;
443   DataV0C[1] = M02;
444   DataV0C[2] = fCentrality;
445   DataV0C[3] = difClusV0C;
446   DataV0C[4] = dr;
447   fClusterPbV0C->Fill(DataV0C);
448
449   Double_t DataTPC[5];
450   DataTPC[0] = Et;
451   DataTPC[1] = M02;
452   DataTPC[2] = fCentrality;
453   DataTPC[3] = difClusTPC;
454   DataTPC[4] = dr;
455   fClusterPbTPC->Fill(DataTPC);
456
457 }
458 //_________________________________________________________________________________________________
459 void AliAnalysisTaskPi0V2::GetMom(TLorentzVector& p, const AliVCluster *c, Double_t *vertex)
460 {
461   // Calculate momentum.
462   Float_t posMom[3];
463   c->GetPosition(posMom);
464   TVector3 clsPos2(posMom);
465
466   Double_t e   = c->E();
467   Double_t r   = clsPos2.Perp();
468   Double_t eta = clsPos2.Eta();
469   Double_t phi = clsPos2.Phi();
470
471   TVector3 pos;
472   pos.SetPtEtaPhi(r,eta,phi);
473
474   if (vertex) { //calculate direction relative to vertex
475     pos -= vertex;
476   }
477
478   Double_t rad = pos.Mag();
479   p.SetPxPyPzE(e*pos.x()/rad, e*pos.y()/rad, e*pos.z()/rad, e);
480
481 }
482 //________________________________________________________________________
483 void AliAnalysisTaskPi0V2::UserCreateOutputObjects()
484 {
485     // Create histograms
486     // Called once (on the worker node)
487         
488     fOutput = new TList();
489     fOutput->SetOwner();  // IMPORTANT!
490
491     hEvtCount = new TH1F("hEvtCount", " Event Plane", 9, 0.5, 9.5);
492     hEvtCount->GetXaxis()->SetBinLabel(1,"All");
493     hEvtCount->GetXaxis()->SetBinLabel(2,"Evt");
494     hEvtCount->GetXaxis()->SetBinLabel(3,"Trg Class");
495     hEvtCount->GetXaxis()->SetBinLabel(4,"Vtx");
496     hEvtCount->GetXaxis()->SetBinLabel(5,"Cent");
497     hEvtCount->GetXaxis()->SetBinLabel(6,"EPtask");
498     hEvtCount->GetXaxis()->SetBinLabel(7,"ClusterTask");
499     hEvtCount->GetXaxis()->SetBinLabel(8,"Pass");
500     fOutput->Add(hEvtCount);
501     
502     hEPTPC   = new TH2F("hEPTPC",   "EPTPC     vs cent", 100, 0., 100., 100, 0., TMath::Pi());
503     hresoTPC = new TH2F("hresoTPC", "TPc reso  vs cent", 100, 0., 100., 100, 0., 1.);
504     hEPV0    = new TH2F("hEPV0",    "EPV0      vs cent", 100, 0., 100., 100, 0., TMath::Pi());
505     hEPV0A   = new TH2F("hEPV0A",   "EPV0A     vs cent", 100, 0., 100., 100, 0., TMath::Pi());
506     hEPV0C   = new TH2F("hEPV0C",   "EPV0C     vs cent", 100, 0., 100., 100, 0., TMath::Pi());
507     hEPV0Ar  = new TH2F("hEPV0Ar",  "EPV0Ar    vs cent", 100, 0., 100., 100, 0., TMath::Pi());
508     hEPV0Cr  = new TH2F("hEPV0Cr",  "EPV0Cr    vs cent", 100, 0., 100., 100, 0., TMath::Pi());
509     hEPV0r   = new TH2F("hEPV0r",   "EPV0r     vs cent", 100, 0., 100., 100, 0., TMath::Pi());
510     hEPV0AR4 = new TH2F("hEPV0AR4", "EPV0AR4   vs cent", 100, 0., 100., 100, 0., TMath::Pi());
511     hEPV0AR7 = new TH2F("hEPV0AR7", "EPV0AR7   vs cent", 100, 0., 100., 100, 0., TMath::Pi());
512     hEPV0CR0 = new TH2F("hEPV0CR0", "EPV0CR0   vs cent", 100, 0., 100., 100, 0., TMath::Pi());
513     hEPV0CR3 = new TH2F("hEPV0CR3", "EPV0CR3   vs cent", 100, 0., 100., 100, 0., TMath::Pi());
514     fOutput->Add(hEPTPC);
515     fOutput->Add(hresoTPC);
516     fOutput->Add(hEPV0);
517     fOutput->Add(hEPV0A);
518     fOutput->Add(hEPV0C);
519     fOutput->Add(hEPV0Ar);
520     fOutput->Add(hEPV0Cr);
521     fOutput->Add(hEPV0r);
522     fOutput->Add(hEPV0AR4);
523     fOutput->Add(hEPV0AR7);
524     fOutput->Add(hEPV0CR0);
525     fOutput->Add(hEPV0CR3);
526
527     hEPTPCCor  = new TH2F("hEPTPCCor",   "EPTPC  vs cent after PHOS Correct", 100, 0., 100., 100, 0., TMath::Pi());
528     hEPV0ACor  = new TH2F("hEPV0ACor",   "EPV0A  vs cent after PHOS Correct", 100, 0., 100., 100, 0., TMath::Pi());
529     hEPV0CCor  = new TH2F("hEPV0CCor",   "EPV0C  vs cent after PHOS Correct", 100, 0., 100., 100, 0., TMath::Pi());
530     fOutput->Add(hEPTPCCor);
531     fOutput->Add(hEPV0ACor);
532     fOutput->Add(hEPV0CCor);
533
534     hdifV0Ar_V0Cr    = new TH2F("hdifV0Ar_V0Cr",    "EP Ar-Cr ", 100, 0., 100., 100, -1., 1.);    
535     hdifV0A_V0CR0    = new TH2F("hdifV0A_V0CR0",    "EP A-R0 ",  100, 0., 100., 100, -1., 1.);    
536     hdifV0A_V0CR3    = new TH2F("hdifV0A_V0CR3",    "EP A-R3 ",  100, 0., 100., 100, -1., 1.);    
537     hdifV0ACR0_V0CR3 = new TH2F("hdifV0ACR0_V0CR3", "EP R0-R3 ", 100, 0., 100., 100, -1., 1.);    
538     hdifV0C_V0AR4    = new TH2F("hdifV0C_V0AR4",    "EP C-R4 ",  100, 0., 100., 100, -1., 1.);    
539     hdifV0C_V0AR7    = new TH2F("hdifV0C_V0AR7",    "EP C-R7 ",  100, 0., 100., 100, -1., 1.);    
540     hdifV0AR4_V0AR7  = new TH2F("hdifV0AR4_V0AR7",  "EP R4-R7 ", 100, 0., 100., 100, -1., 1.);    
541     fOutput->Add(hdifV0Ar_V0Cr);
542     fOutput->Add(hdifV0A_V0CR0);
543     fOutput->Add(hdifV0A_V0CR3);
544     fOutput->Add(hdifV0ACR0_V0CR3);
545     fOutput->Add(hdifV0C_V0AR4);
546     fOutput->Add(hdifV0C_V0AR7);
547     fOutput->Add(hdifV0AR4_V0AR7);
548
549     hdifV0A_V0C = new TH2F("hdifV0A_V0C", "EP A-C  ", 100, 0., 100., 100, -1., 1.);
550     hdifV0A_TPC = new TH2F("hdifV0A_TPC", "EP A-TPC", 100, 0., 100., 100, -1., 1.);
551     hdifV0C_TPC = new TH2F("hdifV0C_TPC", "EP C-TPC", 100, 0., 100., 100, -1., 1.);
552     hdifV0C_V0A = new TH2F("hdifV0C_V0A", "EP C-A  ", 100, 0., 100., 100, -1., 1.);
553     fOutput->Add(hdifV0A_V0C);
554     fOutput->Add(hdifV0A_TPC);
555     fOutput->Add(hdifV0C_TPC);
556     fOutput->Add(hdifV0C_V0A);
557
558
559
560     hdifEMC_EPV0  = new TH3F("hdifEMC_EPV0",  "dif phi in EMC with EP",  100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
561     hdifEMC_EPV0A = new TH3F("hdifEMC_EPV0A", "dif phi in EMC with EP",  100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
562     hdifEMC_EPV0C = new TH3F("hdifEMC_EPV0C", "dif phi in EMC with EP",  100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
563     fOutput->Add(hdifEMC_EPV0);
564     fOutput->Add(hdifEMC_EPV0A);
565     fOutput->Add(hdifEMC_EPV0C);
566
567     hdifful_EPV0 = new TH3F("hdifful_EPV0",    "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
568     hdifful_EPV0A = new TH3F("hdifful_EPV0A",  "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
569     hdifful_EPV0C = new TH3F("hdifful_EPV0C",  "dif phi in full with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
570     fOutput->Add(hdifful_EPV0);
571     fOutput->Add(hdifful_EPV0A);
572     fOutput->Add(hdifful_EPV0C);
573
574     hdifout_EPV0  = new TH3F("hdifout_EPV0",  "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
575     hdifout_EPV0A = new TH3F("hdifout_EPV0A", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
576     hdifout_EPV0C = new TH3F("hdifout_EPV0C", "dif phi NOT in EMC with EP", 100, 0., 100., 100, 0., TMath::Pi(), 15, 0., 15.);
577     fOutput->Add(hdifout_EPV0);
578     fOutput->Add(hdifout_EPV0A);
579     fOutput->Add(hdifout_EPV0C);
580
581   if(isV1Clus){
582                       //  Et   M02  spdcent DeltaPhi    Dr  
583     Int_t    bins[5] = {  500, 350,  100,     100,      100}; // binning
584     Double_t min[5]  = {  0.0, 0.0,    0,     0.0,      0 }; // min x
585     Double_t max[5]  = { 50.0, 3.5,  100,  TMath::Pi(), 0.1}; // max x
586
587     fClusterPbV0A = new THnSparseF("fClusterPbV0A","",5,bins,min,max);
588     fClusterPbV0A->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0A->GetAxis(1)->SetTitle("M02"); fClusterPbV0A->GetAxis(2)->SetTitle("V0M Centrality");
589     fClusterPbV0A->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0A->GetAxis(4)->SetTitle("Dr"); 
590     fOutput->Add(fClusterPbV0A);
591
592     fClusterPbV0C = new THnSparseF("fClusterPbV0C","",5,bins,min,max);
593     fClusterPbV0C->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0C->GetAxis(1)->SetTitle("M02"); fClusterPbV0C->GetAxis(2)->SetTitle("V0M Centrality");
594     fClusterPbV0C->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0C->GetAxis(4)->SetTitle("Dr");
595     fOutput->Add(fClusterPbV0C);
596
597     fClusterPbTPC = new THnSparseF("fClusterPbTPC","",5,bins,min,max);
598     fClusterPbTPC->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbTPC->GetAxis(1)->SetTitle("M02"); fClusterPbTPC->GetAxis(2)->SetTitle("V0M Centrality");
599     fClusterPbTPC->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbTPC->GetAxis(4)->SetTitle("Dr");
600     fOutput->Add(fClusterPbTPC);
601    }
602   
603     h2DcosV0A = new TProfile("h2DcosV0A", "cos(Phi) V0r vs Run NUmber", 200, 0., 200.);
604     h2DsinV0A = new TProfile("h2DsinV0A", "sin(Phi) V0r vs Run NUmber", 200, 0., 200.);
605     h2DcosV0C = new TProfile("h2DcosV0C", "cos(Phi) V0r vs Run NUmber", 200, 0., 200.);
606     h2DsinV0C = new TProfile("h2DsinV0C", "sin(Phi) V0r vs Run NUmber", 200, 0., 200.);
607     h2DcosTPC = new TProfile("h2DcosTPC", "cos(Phi) V0r vs Run NUmber", 200, 0., 200.);
608     h2DsinTPC = new TProfile("h2DsinTPC", "sin(Phi) V0r vs Run NUmber", 200, 0., 200.);
609     fOutput->Add(h2DcosV0A);
610     fOutput->Add(h2DsinV0A);
611     fOutput->Add(h2DcosV0C);
612     fOutput->Add(h2DsinV0C);
613     fOutput->Add(h2DcosTPC);
614     fOutput->Add(h2DsinTPC);
615
616     if(isV1Clus){
617       hM02vsPtA = new TH2F("hM02vsPtA", "M02 vs Et before cut", 5000, 0, 50, 400, 0, 4.);
618       hM02vsPtB = new TH2F("hM02vsPtB", "M02 vs Et before cut", 5000, 0, 50, 400, 0, 4.);
619       fOutput->Add(hM02vsPtA);
620       fOutput->Add(hM02vsPtB);
621      }
622       hClusDxDZA = new TH2F("hClusDxDZA", "clus Dx vs Dz", 1000, -1., 1., 1000, -1., 1);  
623       hClusDxDZB = new TH2F("hClusDxDZB", "clus Dx vs Dz", 1000, -1., 1., 1000, -1., 1);
624       fOutput->Add(hClusDxDZA);
625       fOutput->Add(hClusDxDZB);
626     
627   if(!isV1Clus){
628     const Int_t ndims = 4;
629     Int_t nMgg=500, nPt=40, nCent=20, nDeltaPhi=315;
630     Int_t binsv1[ndims] = {nMgg, nPt, nCent, nDeltaPhi};
631     Double_t xmin[ndims] = { 0,   0.,  0,   0.};
632     Double_t xmax[ndims] = { 0.5, 20., 100, 3.15};
633     fHEPV0A = new THnSparseF("fHEPV0A",   "Flow histogram EPV0A", ndims, binsv1, xmin, xmax);
634     fHEPV0C = new THnSparseF("fHEPV0C",   "Flow histogram EPV0C", ndims, binsv1, xmin, xmax);
635     fHEPTPC = new THnSparseF("fHEPTPC",   "Flow histogram EPTPC", ndims, binsv1, xmin, xmax);
636     fHEPV0A->GetAxis(0)->SetTitle("m_{#gamma#gamma} "); fHEPV0A->GetAxis(1)->SetTitle("p_{T}[GeV]"); fHEPV0A->GetAxis(2)->SetTitle("centrality");fHEPV0A->GetAxis(3)->SetTitle("#delta #phi");
637     fHEPV0C->GetAxis(0)->SetTitle("m_{#gamma#gamma} "); fHEPV0C->GetAxis(1)->SetTitle("p_{T}[GeV]"); fHEPV0C->GetAxis(2)->SetTitle("centrality");fHEPV0C->GetAxis(3)->SetTitle("#delta #phi");
638     fHEPTPC->GetAxis(0)->SetTitle("m_{#gamma#gamma} "); fHEPTPC->GetAxis(1)->SetTitle("p_{T}[GeV]"); fHEPTPC->GetAxis(2)->SetTitle("centrality");fHEPTPC->GetAxis(3)->SetTitle("#delta #phi");
639     fOutput->Add(fHEPV0A);
640     fOutput->Add(fHEPV0C);
641     fOutput->Add(fHEPTPC);
642   }
643     
644
645     PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
646 }
647
648 //________________________________________________________________________
649 void AliAnalysisTaskPi0V2::UserExec(Option_t *) 
650 {
651   // Main loop
652   // Called for each event
653
654   hEvtCount->Fill(1);
655   // Create pointer to reconstructed event
656   AliVEvent *event = InputEvent();
657    if (!event) { Printf("ERROR: Could not retrieve event"); return; }
658    // create pointer to event
659   type = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetDataType();
660   if(!type){
661      AliError("Cannot get the event");
662      return;
663    }
664
665   if(type=="ESD"){
666     fESD = dynamic_cast<AliESDEvent*>(event);
667     if (!fESD) {
668             AliError("Cannot get the ESD event");
669             return;
670     }
671   } else if (type=="AOD"){
672     fAOD = dynamic_cast<AliAODEvent*>(event);
673     if (!fAOD) {
674             AliError("Cannot get the AOD event");
675             return;
676     }
677   }
678
679   if(type=="ESD") fRunNumber = fESD->GetRunNumber();
680   else if(type=="AOD") fRunNumber = fAOD->GetRunNumber();
681   fInterRunNumber = ConvertToInternalRunNumber(fRunNumber);
682
683   hEvtCount->Fill(2);
684   if(!fTrigClass.IsNull()){
685     TString fired;
686     if(type=="ESD"){
687     fired = fESD->GetFiredTriggerClasses();
688     } else if(type=="AOD"){
689     fired = fAOD->GetFiredTriggerClasses();
690     }
691     if (!fired.Contains("-B-"))
692       return;
693     TObjArray *arr = fTrigClass.Tokenize("|");
694     if (!arr)
695       return;
696     Bool_t match = 0;
697     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
698       TObject *obj = arr->At(i);
699       if (!obj)
700         continue;
701       if (fired.Contains(obj->GetName())) {
702         match = 1;
703         break;
704       }
705     }
706     delete arr;
707     if (
708         !match && //select by Trigger classes in KCentral and KSemiCentral
709         !(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kAnyINT ))  // always accept MB
710         )
711       return; //Not match skip this event
712   }
713     hEvtCount->Fill(3);
714
715   if(isPhosCali){
716     SetFlatteningData();
717   }
718     const AliVVertex* fvertex;
719     if(type=="ESD"){
720     fvertex = fESD->GetPrimaryVertex();
721     }else if(type=="AOD"){
722     fvertex = fAOD->GetPrimaryVertex();
723     }
724     if(TMath::Abs(fvertex->GetZ())>fVtxCut)
725       return;
726     Double_t vertex[3] = {fvertex->GetX(), fvertex->GetY(), fvertex->GetZ()};
727
728     hEvtCount->Fill(4);
729
730     if(type=="ESD"){
731       if(fESD->GetCentrality()) {
732         fCentrality = 
733                 fESD->GetCentrality()->GetCentralityPercentile("CL1"); //spd vertex
734     } 
735    } else if(type=="AOD"){
736       if(fAOD->GetCentrality()) {
737         fCentrality = 
738                 fAOD->GetCentrality()->GetCentralityPercentile("CL1"); //spd vertex
739     }  
740    }
741     hEvtCount->Fill(5);
742     if(type=="ESD"){
743       AliEventplane *ep = fESD->GetEventplane();
744     if (ep) {
745       if (ep->GetQVector())
746         fEPTPC    = ep->GetQVector()->Phi()/2. ;
747       else
748         fEPTPC = -999.;
749       if (ep->GetQsub1()&&ep->GetQsub2())
750         fEPTPCreso  = TMath::Cos(2.*(ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.));
751       else
752         fEPTPCreso = -1;
753
754       fEPV0    = ep->GetEventplane("V0",  fESD);
755       fEPV0A   = ep->GetEventplane("V0A", fESD);
756       fEPV0C   = ep->GetEventplane("V0C", fESD);
757       Double_t qx=0, qy=0;
758       Double_t qxr=0, qyr=0;
759       fEPV0Ar  = ep->CalculateVZEROEventPlane(fESD, 4, 5, 2, qxr, qyr);
760       fEPV0Cr  = ep->CalculateVZEROEventPlane(fESD, 2, 3, 2, qx,  qy);
761       qxr += qx;
762       qyr += qy;
763       fEPV0r   = TMath::ATan2(qyr,qxr)/2.;
764       fEPV0AR4 = ep->CalculateVZEROEventPlane(fESD, 4, 2, qx, qy);
765       fEPV0AR5 = ep->CalculateVZEROEventPlane(fESD, 5, 2, qx, qy);
766       fEPV0AR6 = ep->CalculateVZEROEventPlane(fESD, 6, 2, qx, qy);
767       fEPV0AR7 = ep->CalculateVZEROEventPlane(fESD, 7, 2, qx, qy);
768       fEPV0CR0 = ep->CalculateVZEROEventPlane(fESD, 0, 2, qx, qy);
769       fEPV0CR1 = ep->CalculateVZEROEventPlane(fESD, 1, 2, qx, qy);
770       fEPV0CR2 = ep->CalculateVZEROEventPlane(fESD, 2, 2, qx, qy);
771       fEPV0CR3 = ep->CalculateVZEROEventPlane(fESD, 3, 2, qx, qy);
772     }
773     } else if(type=="AOD"){
774       AliEventplane *ep = fAOD->GetEventplane();
775     if (ep) {
776       if (fAOD->GetHeader()){
777         fEPTPC  = fAOD->GetHeader()->GetEventplane();
778       }
779       else
780         fEPTPC = -999.;
781       if (ep->GetQsub1()&&ep->GetQsub2())
782         fEPTPCreso  = TMath::Cos(2.*(ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.));
783       else
784         fEPTPCreso = -1;
785
786       fEPV0    = ep->GetEventplane("V0",  fAOD);
787       fEPV0A   = ep->GetEventplane("V0A", fAOD);
788       fEPV0C   = ep->GetEventplane("V0C", fAOD);
789       Double_t qx=0, qy=0;
790       Double_t qxr=0, qyr=0;
791       fEPV0Ar  = ep->CalculateVZEROEventPlane(fAOD, 4, 5, 2, qxr, qyr);
792       fEPV0Cr  = ep->CalculateVZEROEventPlane(fAOD, 2, 3, 2, qx,  qy);
793       qxr += qx;
794       qyr += qy;
795       fEPV0r   = TMath::ATan2(qyr,qxr)/2.;
796       fEPV0AR4 = ep->CalculateVZEROEventPlane(fAOD, 4, 2, qx, qy);
797       fEPV0AR5 = ep->CalculateVZEROEventPlane(fAOD, 5, 2, qx, qy);
798       fEPV0AR6 = ep->CalculateVZEROEventPlane(fAOD, 6, 2, qx, qy);
799       fEPV0AR7 = ep->CalculateVZEROEventPlane(fAOD, 7, 2, qx, qy);
800       fEPV0CR0 = ep->CalculateVZEROEventPlane(fAOD, 0, 2, qx, qy);
801       fEPV0CR1 = ep->CalculateVZEROEventPlane(fAOD, 1, 2, qx, qy);
802       fEPV0CR2 = ep->CalculateVZEROEventPlane(fAOD, 2, 2, qx, qy);
803       fEPV0CR3 = ep->CalculateVZEROEventPlane(fAOD, 3, 2, qx, qy);
804     }
805   }
806     FillEPQA(); //Fill the EP QA
807
808     hEvtCount->Fill(6);
809
810     if(isPhosCali){
811      // PHOS Flattening
812      fEPV0A = ApplyFlatteningV0A(fEPV0A, fCentrality); //V0A after Phos flatten
813      fEPV0C = ApplyFlatteningV0C(fEPV0C, fCentrality); //V0C after Phos flatten
814      fEPTPC = ApplyFlattening(fEPTPC, fCentrality);    //TPC after Phos flatten
815     }
816
817     fEPV0   = TVector2::Phi_0_2pi(fEPV0);    if(fEPV0>TMath::Pi())   fEPV0  = fEPV0  - TMath::Pi();
818     fEPV0r  = TVector2::Phi_0_2pi(fEPV0r);   if(fEPV0r>TMath::Pi())  fEPV0r = fEPV0r - TMath::Pi();
819     fEPV0A  = TVector2::Phi_0_2pi(fEPV0A);   if(fEPV0A>TMath::Pi())  fEPV0A = fEPV0A - TMath::Pi();
820     fEPV0C  = TVector2::Phi_0_2pi(fEPV0C);   if(fEPV0C>TMath::Pi())  fEPV0C = fEPV0C - TMath::Pi();
821     fEPV0Ar = TVector2::Phi_0_2pi(fEPV0Ar);  if(fEPV0Ar>TMath::Pi()) fEPV0Ar = fEPV0Ar - TMath::Pi();
822     fEPV0Cr = TVector2::Phi_0_2pi(fEPV0Cr);  if(fEPV0Cr>TMath::Pi()) fEPV0Cr = fEPV0Cr - TMath::Pi();
823     fEPV0AR4   = TVector2::Phi_0_2pi(fEPV0AR4);    if(fEPV0AR4>TMath::Pi())   fEPV0AR4  = fEPV0AR4 - TMath::Pi();
824     fEPV0AR7   = TVector2::Phi_0_2pi(fEPV0AR7);    if(fEPV0AR7>TMath::Pi())   fEPV0AR7  = fEPV0AR7 - TMath::Pi();
825     fEPV0CR0   = TVector2::Phi_0_2pi(fEPV0CR0);    if(fEPV0CR0>TMath::Pi())   fEPV0CR0  = fEPV0CR0 - TMath::Pi();
826     fEPV0CR3   = TVector2::Phi_0_2pi(fEPV0CR3);    if(fEPV0CR3>TMath::Pi())   fEPV0CR3  = fEPV0CR3 - TMath::Pi();
827
828    if(fEPTPC != -999. &&  fEPTPC != -1)
829    hEPTPC->Fill(fCentrality,  fEPTPC); 
830    if(fEPTPCreso!=-1) hresoTPC->Fill(fCentrality, fEPTPCreso);
831    hEPV0->Fill(fCentrality,   fEPV0);
832    hEPV0A->Fill(fCentrality,  fEPV0A);
833    hEPV0C->Fill(fCentrality,  fEPV0C);
834    hEPV0Ar->Fill(fCentrality, fEPV0Ar);
835    hEPV0Cr->Fill(fCentrality, fEPV0Cr);
836    hEPV0r->Fill(fCentrality,  fEPV0r);
837    hEPV0AR4->Fill(fCentrality, fEPV0AR4);
838    hEPV0AR7->Fill(fCentrality, fEPV0AR7);
839    hEPV0CR0->Fill(fCentrality, fEPV0CR0);
840    hEPV0CR3->Fill(fCentrality, fEPV0CR3);
841
842    if(!isPhosCali){
843      SetFlatteningData();
844      hEPTPCCor->Fill(fCentrality, ApplyFlattening(fEPTPC, fCentrality));
845      hEPV0ACor->Fill(fCentrality, ApplyFlatteningV0A(fEPV0A, fCentrality));
846      hEPV0CCor->Fill(fCentrality, ApplyFlatteningV0C(fEPV0C, fCentrality));
847    } else {
848      hEPTPCCor->Fill(fCentrality, fEPTPC);
849      hEPV0ACor->Fill(fCentrality, fEPV0A);
850      hEPV0CCor->Fill(fCentrality, fEPV0C);
851   } 
852
853    hdifV0Ar_V0Cr->Fill(fCentrality, TMath::Cos(2.*(fEPV0Ar - fEPV0Cr)));
854    hdifV0A_V0CR0->Fill(fCentrality, TMath::Cos(2.*(fEPV0A - fEPV0CR0)));
855    hdifV0A_V0CR3->Fill(fCentrality, TMath::Cos(2.*(fEPV0A - fEPV0CR3)));
856    hdifV0ACR0_V0CR3->Fill(fCentrality, TMath::Cos(2*(fEPV0CR0 - fEPV0CR3)));
857    hdifV0C_V0AR4->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0AR4)));
858    hdifV0C_V0AR7->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0AR7)));
859    hdifV0AR4_V0AR7->Fill(fCentrality, TMath::Cos(2*(fEPV0AR4 - fEPV0AR7)));
860         
861    hdifV0A_V0C->Fill(fCentrality, TMath::Cos(2*(fEPV0A - fEPV0C)));
862    if(fEPTPC!=-1 && fEPTPC!=-999.){
863      hdifV0A_TPC->Fill(fCentrality, TMath::Cos(2*(fEPV0A - fEPTPC)));
864      hdifV0C_TPC->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPTPC)));
865    }
866    hdifV0C_V0A->Fill(fCentrality, TMath::Cos(2*(fEPV0C - fEPV0A)));
867     // Cluster loop for reconstructed event
868
869 //================ for v2 clusterize analysis==============================================
870   if(!isV1Clus){
871     if (!fV2ClusName.IsNull() && !fV2Clus) {
872       fV2Clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fV2ClusName));
873       if (!fV2Clus) {
874         AliError(Form("%s: Could not retrieve v2 cluster name %s!", GetName(), fV2ClusName.Data()));
875         return;
876       }
877     }
878     Int_t nCluster =  fV2Clus->GetEntries(); 
879     for(Int_t i=0; i<nCluster; ++i){
880       AliVCluster *c1 = static_cast<AliVCluster*>(fV2Clus->At(i));      
881       if(!c1) continue;
882       hClusDxDZA->Fill(c1->GetTrackDz(), c1->GetTrackDx());
883       if(!c1->IsEMCAL()) continue;
884       if(!IsGoodCluster(c1)) continue;
885       hClusDxDZB->Fill(c1->GetTrackDz(), c1->GetTrackDx());
886       TLorentzVector p1;
887       GetMom(p1, c1, vertex);
888       for(Int_t j=i+1; j<nCluster; ++j){
889           AliVCluster *c2 = static_cast<AliVCluster*>(fV2Clus->At(j));      
890         if(!c2) continue;
891         if(!c2->IsEMCAL()) continue;
892         if(!IsGoodCluster(c2)) continue;
893         TLorentzVector p2;
894         GetMom(p2, c2, vertex);
895         FillPion(p1, p2, fEPV0A, fEPV0C, fEPTPC);
896       }
897     }
898   }
899   //================ for v1 clusterize analysis==============================================
900   if(isV1Clus){
901     if (!fV2ClusName.IsNull() && !fV1Clus) {
902       fV1Clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fV1ClusName));
903       if (!fV1Clus) {
904         AliError(Form("%s: Could not retrieve v1 cluster name %s!", GetName(), fV1ClusName.Data()));
905         return;
906       }
907     }
908     Int_t nClusterV1 = fV1Clus->GetEntries();
909     for(Int_t i=0; i<nClusterV1; ++i){
910       AliVCluster *c3 = dynamic_cast<AliVCluster*>(fV1Clus->At(i));      
911       if(!c3) continue;
912       if(!c3->IsEMCAL()) continue;
913       Double_t M02c3 = c3->GetM02();
914       Double_t Dxc3  = c3->GetTrackDx();
915       Double_t Dzc3  = c3->GetTrackDz(); 
916
917       hClusDxDZA->Fill(Dzc3, Dxc3);
918       Float_t clsPosEt[3] = {0,0,0};
919       c3->GetPosition(clsPosEt);
920       TVector3 clsVec(clsPosEt);
921       Double_t Et = c3->E()*TMath::Sin(clsVec.Theta());
922       hM02vsPtA->Fill(Et, M02c3);
923       if(!IsGoodClusterV1(c3)) continue;
924       hM02vsPtB->Fill(Et, M02c3);
925       hClusDxDZB->Fill(Dzc3, Dxc3);
926       TLorentzVector p3;
927       GetMom(p3, c3, vertex);
928       FillCluster(p3, fEPV0A, fEPV0C, fEPTPC, c3);
929     }
930   }
931
932     hEvtCount->Fill(7);
933
934    if (!fTracksName.IsNull() && !fTracks) {
935      fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
936      if (!fTracks) {
937        AliError(Form("%s: Could not retrieve tracks %s!", GetName(), fTracksName.Data())); 
938        return;
939      }
940   }
941
942    Int_t ntracks = fTracks->GetEntries();
943    for(Int_t i=0; i<ntracks; ++i){
944      AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(i));
945      if(!track) continue;
946      Double_t tPhi = track->Phi();
947      Double_t tPt  = track->Pt();
948      Double_t Eta  = track->Eta();
949
950      Double_t difTrackV0  = TVector2::Phi_0_2pi(tPhi-fEPV0);   if(difTrackV0  >TMath::Pi()) difTrackV0  -= TMath::Pi();
951      Double_t difTrackV0A = TVector2::Phi_0_2pi(tPhi-fEPV0A);  if(difTrackV0A >TMath::Pi()) difTrackV0A -= TMath::Pi();
952      Double_t difTrackV0C = TVector2::Phi_0_2pi(tPhi-fEPV0C);  if(difTrackV0C >TMath::Pi()) difTrackV0C -= TMath::Pi();
953      Double_t difTrackTPC = TVector2::Phi_0_2pi(tPhi-fEPTPC);  if(difTrackTPC >TMath::Pi()) difTrackTPC -= TMath::Pi();
954      if(tPhi*TMath::RadToDeg()>80. && tPhi*TMath::RadToDeg()<180. && Eta <0.7 && Eta >(-0.7)){  
955        hdifEMC_EPV0->Fill(fCentrality, difTrackV0, tPt);
956        hdifEMC_EPV0A->Fill(fCentrality, difTrackV0A, tPt);
957        hdifEMC_EPV0C->Fill(fCentrality, difTrackV0C, tPt);
958      }else{
959        hdifout_EPV0->Fill(fCentrality, difTrackV0, tPt);
960        hdifout_EPV0A->Fill(fCentrality, difTrackV0A, tPt);
961        hdifout_EPV0C->Fill(fCentrality, difTrackV0C, tPt);
962      }
963      hdifful_EPV0->Fill(fCentrality,    difTrackV0, tPt);
964      hdifful_EPV0A->Fill(fCentrality,   difTrackV0A, tPt);
965      hdifful_EPV0C->Fill(fCentrality,   difTrackV0C, tPt);
966    } 
967     hEvtCount->Fill(8);
968
969     // NEW HISTO should be filled before this point, as PostData puts the
970     // information for this iteration of the UserExec in the container
971     PostData(1, fOutput);
972 }
973 //____________________________________________________________________
974 Int_t AliAnalysisTaskPi0V2::ConvertToInternalRunNumber(Int_t n)
975 {
976     switch(n){
977     case  170593 : return 179 ;
978     case  170572 : return 178 ;
979     case  170556 : return 177 ;
980     case  170552 : return 176 ;
981     case  170546 : return 175 ;
982     case  170390 : return 174 ;
983     case  170389 : return 173 ;
984     case  170388 : return 172 ;
985     case  170387 : return 171 ;
986     case  170315 : return 170 ;
987     case  170313 : return 169 ;
988     case  170312 : return 168 ;
989     case  170311 : return 167 ;
990     case  170309 : return 166 ;
991     case  170308 : return 165 ;
992     case  170306 : return 164 ;
993     case  170270 : return 163 ;
994     case  170269 : return 162 ;
995     case  170268 : return 161 ;
996     case  170267 : return 160 ;
997     case  170264 : return 159 ;
998     case  170230 : return 158 ;
999     case  170228 : return 157 ;
1000     case  170208 : return 156 ;
1001     case  170207 : return 155 ;
1002     case  170205 : return 154 ;
1003     case  170204 : return 153 ;
1004     case  170203 : return 152 ;
1005     case  170195 : return 151 ;
1006     case  170193 : return 150 ;
1007     case  170163 : return 149 ;
1008     case  170162 : return 148 ;
1009     case  170159 : return 147 ;
1010     case  170155 : return 146 ;
1011     case  170152 : return 145 ;
1012     case  170091 : return 144 ;
1013     case  170089 : return 143 ;
1014     case  170088 : return 142 ;
1015     case  170085 : return 141 ;
1016     case  170084 : return 140 ;
1017     case  170083 : return 139 ;
1018     case  170081 : return 138 ;
1019     case  170040 : return 137 ;
1020     case  170038 : return 136 ;
1021     case  170036 : return 135 ;
1022     case  170027 : return 134 ;
1023     case  169981 : return 133 ;
1024     case  169975 : return 132 ;
1025     case  169969 : return 131 ;
1026     case  169965 : return 130 ;
1027     case  169961 : return 129 ;
1028     case  169956 : return 128 ;
1029     case  169926 : return 127 ;
1030     case  169924 : return 126 ;
1031     case  169923 : return 125 ;
1032     case  169922 : return 124 ;
1033     case  169919 : return 123 ;
1034     case  169918 : return 122 ;
1035     case  169914 : return 121 ;
1036     case  169859 : return 120 ;
1037     case  169858 : return 119 ;
1038     case  169855 : return 118 ;
1039     case  169846 : return 117 ;
1040     case  169838 : return 116 ;
1041     case  169837 : return 115 ;
1042     case  169835 : return 114 ;
1043     case  169683 : return 113 ;
1044     case  169628 : return 112 ;
1045     case  169591 : return 111 ;
1046     case  169590 : return 110 ;
1047     case  169588 : return 109 ;
1048     case  169587 : return 108 ;
1049     case  169586 : return 107 ;
1050     case  169584 : return 106 ;
1051     case  169557 : return 105 ;
1052     case  169555 : return 104 ;
1053     case  169554 : return 103 ;
1054     case  169553 : return 102 ;
1055     case  169550 : return 101 ;
1056     case  169515 : return 100 ;
1057     case  169512 : return 99 ;
1058     case  169506 : return 98 ;
1059     case  169504 : return 97 ;
1060     case  169498 : return 96 ;
1061     case  169475 : return 95 ;
1062     case  169420 : return 94 ;
1063     case  169419 : return 93 ;
1064     case  169418 : return 92 ;
1065     case  169417 : return 91 ;
1066     case  169415 : return 90 ;
1067     case  169411 : return 89 ;
1068     case  169238 : return 88 ;
1069     case  169236 : return 87 ;
1070     case  169167 : return 86 ;
1071     case  169160 : return 85 ;
1072     case  169156 : return 84 ;
1073     case  169148 : return 83 ;
1074     case  169145 : return 82 ;
1075     case  169144 : return 81 ;
1076     case  169143 : return 80 ;
1077     case  169138 : return 79 ;
1078     case  169099 : return 78 ;
1079     case  169094 : return 77 ;
1080     case  169091 : return 76 ;
1081     case  169045 : return 75 ;
1082     case  169044 : return 74 ;
1083     case  169040 : return 73 ;
1084     case  169035 : return 72 ;
1085     case  168992 : return 71 ;
1086     case  168988 : return 70 ;
1087     case  168984 : return 69 ;
1088     case  168826 : return 68 ;
1089     case  168777 : return 67 ;
1090     case  168514 : return 66 ;
1091     case  168512 : return 65 ;
1092     case  168511 : return 64 ;
1093     case  168467 : return 63 ;
1094     case  168464 : return 62 ;
1095     case  168461 : return 61 ;
1096     case  168460 : return 60 ;
1097     case  168458 : return 59 ;
1098     case  168362 : return 58 ;
1099     case  168361 : return 57 ;
1100     case  168356 : return 56 ;
1101     case  168342 : return 55 ;
1102     case  168341 : return 54 ;
1103     case  168325 : return 53 ;
1104     case  168322 : return 52 ;
1105     case  168318 : return 51 ;
1106     case  168311 : return 50 ;
1107     case  168310 : return 49 ;
1108     case  168213 : return 48 ;
1109     case  168212 : return 47 ;
1110     case  168208 : return 46 ;
1111     case  168207 : return 45 ;
1112     case  168206 : return 44 ;
1113     case  168205 : return 43 ;
1114     case  168204 : return 42 ;
1115     case  168203 : return 41 ;
1116     case  168181 : return 40 ;
1117     case  168177 : return 39 ;
1118     case  168175 : return 38 ;
1119     case  168173 : return 37 ;
1120     case  168172 : return 36 ;
1121     case  168171 : return 35 ;
1122     case  168115 : return 34 ;
1123     case  168108 : return 33 ;
1124     case  168107 : return 32 ;
1125     case  168105 : return 31 ;
1126     case  168104 : return 30 ;
1127     case  168103 : return 29 ;
1128     case  168076 : return 28 ;
1129     case  168069 : return 27 ;
1130     case  168068 : return 26 ;
1131     case  168066 : return 25 ;
1132     case  167988 : return 24 ;
1133     case  167987 : return 23 ;
1134     case  167986 : return 22 ;
1135     case  167985 : return 21 ;
1136     case  167921 : return 20 ;
1137     case  167920 : return 19 ;
1138     case  167915 : return 18 ;
1139     case  167909 : return 17 ;
1140     case  167903 : return 16 ;
1141     case  167902 : return 15 ;
1142     case  167818 : return 14 ;
1143     case  167814 : return 13 ;
1144     case  167813 : return 12 ;
1145     case  167808 : return 11 ;
1146     case  167807 : return 10 ;
1147     case  167806 : return 9 ;
1148     case  167713 : return 8 ;
1149     case  167712 : return 7 ;
1150     case  167711 : return 6 ;
1151     case  167706 : return 5 ;
1152     case  167693 : return 4 ;
1153     case  166532 : return 3 ;
1154     case  166530 : return 2 ;
1155     case  166529 : return 1 ;
1156
1157     default : return 199;
1158     }
1159 }
1160 //_______________________________________________________________________
1161 void AliAnalysisTaskPi0V2::FillEPQA()
1162 {
1163   
1164   h2DcosV0A->Fill(fInterRunNumber, TMath::Cos(fEPV0A));
1165   h2DsinV0A->Fill(fInterRunNumber, TMath::Sin(fEPV0A));
1166   h2DcosV0C->Fill(fInterRunNumber, TMath::Cos(fEPV0C));
1167   h2DsinV0C->Fill(fInterRunNumber, TMath::Sin(fEPV0C));
1168   h2DcosTPC->Fill(fInterRunNumber, TMath::Cos(fEPTPC));
1169   h2DsinTPC->Fill(fInterRunNumber, TMath::Sin(fEPTPC));
1170
1171
1172 }
1173 //_________________________________________________________________________________
1174 void AliAnalysisTaskPi0V2::SetFlatteningData(){
1175   //Read objects with flattening parameters
1176   AliOADBContainer flatContainer("phosFlat");
1177   flatContainer.InitFromFile(fEPcalibFileName.Data(),"phosFlat");
1178   TObjArray *maps = (TObjArray*)flatContainer.GetObject(fRunNumber,"phosFlat");
1179   if(!maps){
1180       AliError(Form("Can not read Flattening for run %d. \n From file >%s<\n",fRunNumber,fEPcalibFileName.Data())) ;    
1181   }
1182   else{
1183     AliInfo(Form("Setting PHOS flattening with name %s \n",maps->GetName())) ;
1184     AliEPFlattener * h = (AliEPFlattener*)maps->At(0) ;  
1185     if(fTPCFlat) delete fTPCFlat ;
1186     fTPCFlat = new AliEPFlattener();
1187     fTPCFlat = h ;
1188     h = (AliEPFlattener*)maps->At(1);  
1189     if(fV0AFlat) delete fV0AFlat ;
1190     fV0AFlat = new AliEPFlattener();
1191     fV0AFlat = h ;
1192     h = (AliEPFlattener*)maps->At(2);  
1193     if(fV0CFlat) delete fV0CFlat ;
1194     fV0CFlat = new AliEPFlattener();
1195     fV0CFlat = h ;
1196   }    
1197  
1198 }
1199  //____________________________________________________________________________
1200 Double_t  AliAnalysisTaskPi0V2::ApplyFlattening(Double_t phi, Double_t c){
1201  
1202   if(fTPCFlat)
1203     return fTPCFlat->MakeFlat(phi,c);
1204   return phi ;
1205
1206 }
1207 //____________________________________________________________________________
1208 Double_t  AliAnalysisTaskPi0V2::ApplyFlatteningV0A(Double_t phi, Double_t c){
1209  
1210   if(fV0AFlat)
1211     return fV0AFlat->MakeFlat(phi,c);
1212   return phi ;
1213
1214 }
1215 //____________________________________________________________________________
1216 Double_t  AliAnalysisTaskPi0V2::ApplyFlatteningV0C(Double_t phi, Double_t c){
1217  
1218   if(fV0CFlat)
1219     return fV0CFlat->MakeFlat(phi,c);
1220   return phi ;
1221
1222 }
1223 //________________________________________________________________________
1224 void AliAnalysisTaskPi0V2::Terminate(Option_t *) 
1225 {
1226     // Draw result to screen, or perform fitting, normalizations
1227     // Called once at the end of the query
1228 //    fOutput = dynamic_cast<TList*> (GetOutputData(1));
1229  //   if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
1230         
1231     // Get the physics selection histograms with the selection statistics
1232     //AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1233     //AliESDInputHandler *inputH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
1234     //TH2F *histStat = (TH2F*)inputH->GetStatistics();
1235    
1236    
1237     // NEW HISTO should be retrieved from the TList container in the above way,
1238     // so it is available to draw on a canvas such as below
1239
1240 }