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