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