PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0V2ShSh.cxx

   1 // $Id: AliAnalysisTaskEMCALPi0V2ShSh.cxx$
   2
   3 #include "AliAnalysisTaskEMCALPi0V2ShSh.h"
   4
   5 //Root include files
   6 //#include <Riostream.h>
   7 #include <TParticle.h>
   8 #include <TRefArray.h>
   9 #include <TCanvas.h>
  10 #include <TChain.h>
  11 #include <TFile.h>
  12 #include <TH1F.h>
  13 #include <TH1D.h>
  14 #include <TH2F.h>
  15 #include <THnSparse.h>
  16 #include <TList.h>
  17 #include <TMath.h>
  18 #include <TVirtualFFT.h>
  19
  20 //AliRoot include files
  21 #include "AliAnalysisTaskSE.h"
  22 #include "AliRunLoader.h"
  23 #include "AliAnalysisManager.h"
  24 #include "AliAnalysisTask.h"
  25 #include "AliStack.h"
  26 #include "AliESDEvent.h"
  27 #include "AliESDVertex.h"
  28 #include "AliESDCaloCluster.h"
  29 #include "AliESDEvent.h"
  30 #include "AliESDHeader.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDtrack.h"
  33 #include "AliKFParticle.h"
  34 #include "AliVCluster.h"
  35 #include "AliCentrality.h"
  36 #include "AliEventplane.h"
  37
  38
  39
  40 ClassImp(AliAnalysisTaskEMCALPi0V2ShSh)
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh() :
  44   AliAnalysisTaskSE(),
  45   fCaloClusters(0), fEventPlane(0), fCentralityV0M(99.), fESD(0), fOutputList(0),
  46   fEPTPC(-999.), fEPTPCResolution(0.),
  47   fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
  48   fEPV0A4r(-999.), fEPV0A5r(-999.), fEPV0A6r(-999.), fEPV0A7r(-999.), fEPV0C0r(-999.), fEPV0C1r(-999.), fEPV0C2r(-999.), fEPV0C3r(-999.),
  49   fHistAllcentV0(0), fHistAllcentV0r(0), fHistAllcentV0A(0), fHistAllcentV0C(0), fHistAllcentTPC(0),
  50   fHistEPTPC(0), fHistEPTPCResolution(0),
  51   fHistEPV0(0), fHistEPV0A(0), fHistEPV0C(0), fHistEPV0Ar(0), fHistEPV0Cr(0), fHistEPV0r(0), fHistEPV0A4r(0), fHistEPV0A7r(0), fHistEPV0C0r(0), fHistEPV0C3r(0),
  52   fHistdifV0A_V0C0r(0), fHistdifV0A_V0C3r(0), fHistdifV0C0r_V0C3r(0), fHistdifV0C_V0A4r(0), fHistdifV0C_V0A7r(0), fHistdifV0A4r_V0A7r(0), fHistdifV0Ar_V0Cr(0),
  53   fHistClusterEta(0), fHistClusterPhi(0), fHistClusterE(0), fHistClusterEt(0), fHistClusterN(0), fHistClusterM02(0),
  54   fHistClusterEN(0), fHistClusterEM02(0), fHistClusterPhiEta(0), fHistClusterEtN(0), fHistClusterEtM02(0), fHistClusterdphiV0(0),
  55   fHistTrackPt(0), fHistTrackEta(0), fHistTrackPhi(0), fHistTrackPhiEta(0),
  56   fClusterPbV0(0), fClusterPbV0A(0), fClusterPbV0C(0), fClusterPbTPC(0)
  57 {
  58   // Default constructor.
  59 }
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskEMCALPi0V2ShSh::AliAnalysisTaskEMCALPi0V2ShSh(const char *name) :
  63   AliAnalysisTaskSE(name),
  64   fCaloClusters(0), fEventPlane(0), fCentralityV0M(99.), fESD(0), fOutputList(0),
  65   fEPTPC(-999.), fEPTPCResolution(0.),
  66   fEPV0(-999.), fEPV0A(-999.), fEPV0C(-999.), fEPV0Ar(-999.), fEPV0Cr(-999.), fEPV0r(-999.),
  67   fEPV0A4r(-999.), fEPV0A5r(-999.), fEPV0A6r(-999.), fEPV0A7r(-999.), fEPV0C0r(-999.), fEPV0C1r(-999.), fEPV0C2r(-999.), fEPV0C3r(-999.),
  68   fHistAllcentV0(0), fHistAllcentV0r(0), fHistAllcentV0A(0), fHistAllcentV0C(0), fHistAllcentTPC(0),
  69   fHistEPTPC(0), fHistEPTPCResolution(0),
  70   fHistEPV0(0), fHistEPV0A(0), fHistEPV0C(0), fHistEPV0Ar(0), fHistEPV0Cr(0), fHistEPV0r(0), fHistEPV0A4r(0), fHistEPV0A7r(0), fHistEPV0C0r(0), fHistEPV0C3r(0),
  71   fHistdifV0A_V0C0r(0), fHistdifV0A_V0C3r(0), fHistdifV0C0r_V0C3r(0), fHistdifV0C_V0A4r(0), fHistdifV0C_V0A7r(0), fHistdifV0A4r_V0A7r(0), fHistdifV0Ar_V0Cr(0),
  72   fHistClusterEta(0), fHistClusterPhi(0), fHistClusterE(0), fHistClusterEt(0), fHistClusterN(0), fHistClusterM02(0),
  73   fHistClusterEN(0), fHistClusterEM02(0), fHistClusterPhiEta(0), fHistClusterEtN(0), fHistClusterEtM02(0), fHistClusterdphiV0(0),
  74   fHistTrackPt(0), fHistTrackEta(0), fHistTrackPhi(0), fHistTrackPhiEta(0),
  75   fClusterPbV0(0), fClusterPbV0A(0), fClusterPbV0C(0), fClusterPbTPC(0)
  76 {
  77   // Constructor
  78   // Define input and output slots here
  79   // Input slot #0 works with a TChain
  80   DefineInput(0, TChain::Class());
  81   // Output slot #0 id reserved by the base class for AOD
  82   // Output slot #1 writes into a TH1 container
  83   DefineOutput(1, TList::Class());
  84 }
  85
  86 //________________________________________________________________________
  87 void AliAnalysisTaskEMCALPi0V2ShSh::UserCreateOutputObjects()
  88 {
  89   // Create histograms, called once.
  90
  91   fCaloClusters = new TRefArray();
  92
  93   fOutputList = new TList();
  94   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
  95
  96   fHistTrackPt = new TH1F("fHistTrackPt","Track Transverse Momentum Distribution of Pb+Pb",100,0.0,30.0);
  97   fHistTrackPt->GetYaxis()->SetTitle("Entries"); fHistTrackPt->GetXaxis()->SetTitle("P_{t} [GeV/c]");
  98   fOutputList->Add(fHistTrackPt);
  99
 100   fHistTrackEta = new TH1F("fHistTrackEta","Track Pseudorapidity Distribution of Pb+Pb",100,-1.0,1.0);
 101   fHistTrackEta->GetYaxis()->SetTitle("Entries"); fHistTrackEta->GetXaxis()->SetTitle("#eta");
 102   fOutputList->Add(fHistTrackEta);
 103
 104   fHistTrackPhi = new TH1F("fHistTrackPhi","Track #phi Distribution of Pb+Pb",100,0.0,6.29);
 105   fHistTrackPhi->GetYaxis()->SetTitle("Entries"); fHistTrackPhi->GetXaxis()->SetTitle("#phi [rad]");
 106   fOutputList->Add(fHistTrackPhi);
 107
 108   fHistTrackPhiEta = new TH2F("fHistTrackPhiEta","Track Pseudorapidity vs #phi of Pb+Pb",100,-1.0,1.0,100,0.0,6.29);
 109   fHistTrackPhiEta->GetXaxis()->SetTitle("#eta"); fHistTrackPhiEta->GetYaxis()->SetTitle("#phi [rad]");
 110   fOutputList->Add(fHistTrackPhiEta);
 111
 112   fHistClusterEta = new TH1F("fHistClusterEta","Cluster Pseudorapidity Distribution of Pb+Pb",100,-1.0,1.0);
 113   fHistClusterEta->GetYaxis()->SetTitle("Entries"); fHistClusterEta->GetXaxis()->SetTitle("#eta");
 114   fOutputList->Add(fHistClusterEta);
 115
 116   fHistClusterPhi = new TH1F("fHistClusterPhi","Cluster #phi Distribution of Pb+Pb",100,0.0,6.29);
 117   fHistClusterPhi->GetYaxis()->SetTitle("Entries"); fHistClusterPhi->GetXaxis()->SetTitle("#phi [rad]");
 118   fOutputList->Add(fHistClusterPhi);
 119
 120   fHistClusterPhiEta = new TH2F("fHistClusterPhiEta","Cluster Pseudorapidity vs #phi of Pb+Pb",100,-1.0,1.0,100,0.0,6.29);
 121   fHistClusterPhiEta->GetXaxis()->SetTitle("#eta"); fHistClusterPhiEta->GetYaxis()->SetTitle("#phi [rad]");
 122   fOutputList->Add(fHistClusterPhiEta);
 123
 124   fHistClusterM02 = new TH1F("fHistClusterM02","Cluster M02 Distribution of Pb+Pb",100,0.0,3.0);
 125   fHistClusterM02->GetYaxis()->SetTitle("Entries"); fHistClusterM02->GetXaxis()->SetTitle("M02");
 126   fOutputList->Add(fHistClusterM02);
 127
 128   fHistClusterE = new TH1F("fHistClusterE","Cluster Energy Distribution of Pb+Pb",100,0.0,20.0);
 129   fHistClusterE->GetYaxis()->SetTitle("Entries"); fHistClusterE->GetXaxis()->SetTitle("Energy [GeV]");
 130   fOutputList->Add(fHistClusterE);
 131
 132   fHistClusterEt = new TH1F("fHistClusterEt","Cluster Transverse Energy Distribution of Pb+Pb",100,0.0,20.0);
 133   fHistClusterEt->GetYaxis()->SetTitle("Entries"); fHistClusterEt->GetXaxis()->SetTitle("Transverse Energy [GeV]");
 134   fOutputList->Add(fHistClusterEt);
 135
 136   fHistClusterEM02 = new TH2F("fHistClusterEM02","Cluster Energy vs M02 of Pb+Pb",100,0.0,20.0,100,0.0,3.0);
 137   fHistClusterEM02->GetYaxis()->SetTitle("M02"); fHistClusterEM02->GetXaxis()->SetTitle("Energy [GeV]");
 138   fOutputList->Add(fHistClusterEM02);
 139
 140   fHistClusterEtM02 = new TH2F("fHistClusterEtM02","Cluster Transverse Energy vs M02 of Pb+Pb",100,0.0,20.0,100,0.0,3.0);
 141   fHistClusterEtM02->GetYaxis()->SetTitle("M02"); fHistClusterEtM02->GetXaxis()->SetTitle("Transverse Energy [GeV]");
 142   fOutputList->Add(fHistClusterEtM02);
 143
 144   fHistClusterN = new TH1F("fHistClusterN","Cluster N Distribution of Pb+Pb",30,0.0,30.0);
 145   fHistClusterN->GetYaxis()->SetTitle("Entries"); fHistClusterN->GetXaxis()->SetTitle("N");
 146   fOutputList->Add(fHistClusterN);
 147
 148   fHistClusterEN = new TH2F("fHistClusterEN","N vs Cluster Energy of Pb+Pb",100,0.0,20.0,30,0.0,30.0);
 149   fHistClusterEN->GetYaxis()->SetTitle("N"); fHistClusterEN->GetXaxis()->SetTitle("Energy [GeV]");
 150   fOutputList->Add(fHistClusterEN);
 151
 152   fHistClusterEtN = new TH2F("fHistClusterEtN","N vs Cluster Transverse Energy of Pb+Pb",100,0.0,20.0,30,0.0,30.0);
 153   fHistClusterEtN->GetYaxis()->SetTitle("N"); fHistClusterEtN->GetXaxis()->SetTitle("Transverse Energy [GeV]");
 154   fOutputList->Add(fHistClusterEtN);
 155
 156   fHistClusterdphiV0 = new TH1D("fHistClusterdphiV0","Cluster dphiV0 Distribution of Pb+Pb",100,0.0,TMath::Pi());
 157   fHistClusterdphiV0->GetYaxis()->SetTitle("Entries"); fHistClusterdphiV0->GetXaxis()->SetTitle("dphiV0 [rad]");
 158   fOutputList->Add(fHistClusterdphiV0);
 159
 160   fHistEPV0 = new TH2F("fHistEPV0","V0 Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 161   fHistEPV0->GetYaxis()->SetTitle("V0 Event Plane Angle [rad]"); fHistEPV0->GetXaxis()->SetTitle("V0M Centrality");
 162   fOutputList->Add(fHistEPV0);
 163
 164   fHistEPV0A = new TH2F("fHistEPV0A","V0A Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 165   fHistEPV0A->GetYaxis()->SetTitle("V0A Event Plane Angle [rad]"); fHistEPV0A->GetXaxis()->SetTitle("V0M Centrality");
 166   fOutputList->Add(fHistEPV0A);
 167
 168   fHistEPV0C = new TH2F("fHistEPV0","V0C Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 169   fHistEPV0C->GetYaxis()->SetTitle("V0C Event Plane Angle [rad]"); fHistEPV0C->GetXaxis()->SetTitle("V0M Centrality");
 170   fOutputList->Add(fHistEPV0C);
 171
 172   fHistEPV0r = new TH2F("fHistEPV0r","V0r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 173   fHistEPV0r->GetYaxis()->SetTitle("V0r Event Plane Angle [rad]"); fHistEPV0r->GetXaxis()->SetTitle("V0M Centrality");
 174   fOutputList->Add(fHistEPV0r);
 175
 176   fHistEPV0Ar = new TH2F("fHistEPV0","V0Ar Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 177   fHistEPV0Ar->GetYaxis()->SetTitle("V0Ar Event Plane Angle [rad]"); fHistEPV0Ar->GetXaxis()->SetTitle("V0M Centrality");
 178   fOutputList->Add(fHistEPV0Ar);
 179
 180   fHistEPV0Cr = new TH2F("fHistEPV0Cr","V0Cr Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 181   fHistEPV0Cr->GetYaxis()->SetTitle("V0Cr Event Plane Angle [rad]"); fHistEPV0Cr->GetXaxis()->SetTitle("V0M Centrality");
 182   fOutputList->Add(fHistEPV0Cr);
 183
 184   fHistEPV0A4r = new TH2F("fHistEPV0Cr","V0A4r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 185   fHistEPV0A4r->GetYaxis()->SetTitle("V0A4r Event Plane Angle [rad]"); fHistEPV0A4r->GetXaxis()->SetTitle("V0M Centrality");
 186   fOutputList->Add(fHistEPV0A4r);
 187
 188   fHistEPV0A7r = new TH2F("fHistEPV0Cr","V0A7r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 189   fHistEPV0A7r->GetYaxis()->SetTitle("V0A7r Event Plane Angle [rad]"); fHistEPV0A7r->GetXaxis()->SetTitle("V0M Centrality");
 190   fOutputList->Add(fHistEPV0A7r);
 191
 192   fHistEPV0C0r = new TH2F("fHistEPV0Cr","V0C0r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 193   fHistEPV0C0r->GetYaxis()->SetTitle("V0C0r Event Plane Angle [rad]"); fHistEPV0C0r->GetXaxis()->SetTitle("V0M Centrality");
 194   fOutputList->Add(fHistEPV0C0r);
 195
 196   fHistEPV0C3r = new TH2F("fHistEPV0Cr","V0C3r Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 197   fHistEPV0C3r->GetYaxis()->SetTitle("V0C3r Event Plane Angle [rad]"); fHistEPV0C3r->GetXaxis()->SetTitle("V0M Centrality");
 198   fOutputList->Add(fHistEPV0C3r);
 199
 200   fHistdifV0A_V0C0r = new TH2F("fHistdifV0A_V0C0r","(V0A - V0C0r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 201   fHistdifV0A_V0C0r->GetYaxis()->SetTitle("Cos[2*(V0A - V0C0r)]"); fHistdifV0A_V0C0r->GetXaxis()->SetTitle("V0M Centrality");
 202   fOutputList->Add(fHistdifV0A_V0C0r);
 203
 204   fHistdifV0A_V0C3r = new TH2F("fHistdifV0A_V0C3r","(V0A - V0C3r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 205   fHistdifV0A_V0C3r->GetYaxis()->SetTitle("Cos[2*(V0A - V0C3r)]"); fHistdifV0A_V0C3r->GetXaxis()->SetTitle("V0M Centrality");
 206   fOutputList->Add(fHistdifV0A_V0C3r);
 207
 208   fHistdifV0C0r_V0C3r = new TH2F("fHistdifV0C0r_V0C3r","(V0C0r - V0C3r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 209   fHistdifV0C0r_V0C3r->GetYaxis()->SetTitle("Cos[2*(V0C0r - V0C3r)]"); fHistdifV0C0r_V0C3r->GetXaxis()->SetTitle("V0M Centrality");
 210   fOutputList->Add(fHistdifV0C0r_V0C3r);
 211
 212   fHistdifV0C_V0A4r = new TH2F("fHistdifV0C_V0A4r","(V0C - V0A4r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 213   fHistdifV0C_V0A4r->GetYaxis()->SetTitle("Cos[2*(V0C - V0C4r)]"); fHistdifV0C_V0A4r->GetXaxis()->SetTitle("V0M Centrality");
 214   fOutputList->Add(fHistdifV0C_V0A4r);
 215
 216   fHistdifV0C_V0A7r = new TH2F("fHistdifV0C_V0A3r","(V0C - V0A7r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 217   fHistdifV0C_V0A7r->GetYaxis()->SetTitle("Cos[2*(V0C - V0A7r)]"); fHistdifV0C_V0A7r->GetXaxis()->SetTitle("V0M Centrality");
 218   fOutputList->Add(fHistdifV0C_V0A7r);
 219
 220   fHistdifV0A4r_V0A7r = new TH2F("fHistdifV0A4r_V0A7r","(V0A4r - V0A7r) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 221   fHistdifV0A4r_V0A7r->GetYaxis()->SetTitle("Cos[2*(V0A4r - V0A7r)]"); fHistdifV0A4r_V0A7r->GetXaxis()->SetTitle("V0M Centrality");
 222   fOutputList->Add(fHistdifV0A4r_V0A7r);
 223
 224   fHistdifV0Ar_V0Cr = new TH2F("fHistdifV0Ar_V0Cr","(V0Ar - V0Cr) vs Centrality of Pb+Pb",100,0,100,100,-1.0,1.0);
 225   fHistdifV0Ar_V0Cr->GetYaxis()->SetTitle("Cos[2*(V0Ar - V0Cr)]"); fHistdifV0Ar_V0Cr->GetXaxis()->SetTitle("V0M Centrality");
 226   fOutputList->Add(fHistdifV0Ar_V0Cr);
 227
 228   fHistAllcentV0 = new TH1F("fHistAllcentV0","V0 Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
 229   fHistAllcentV0->GetXaxis()->SetTitle("V0 Event Plane Angle [rad]"); fHistAllcentV0->GetYaxis()->SetTitle("");
 230   fOutputList->Add(fHistAllcentV0);
 231
 232   fHistAllcentV0r = new TH1F("fHistAllcentV0r","V0r Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
 233   fHistAllcentV0r->GetXaxis()->SetTitle("V0r Event Plane Angle [rad]"); fHistAllcentV0r->GetYaxis()->SetTitle("");
 234   fOutputList->Add(fHistAllcentV0r);
 235
 236   fHistAllcentV0A = new TH1F("fHistAllcentV0A","V0A Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
 237   fHistAllcentV0A->GetXaxis()->SetTitle("V0A Event Plane Angle [rad]"); fHistAllcentV0A->GetYaxis()->SetTitle("");
 238   fOutputList->Add(fHistAllcentV0A);
 239
 240   fHistAllcentV0C = new TH1F("fHistAllcentV0C","V0C Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
 241   fHistAllcentV0C->GetXaxis()->SetTitle("V0C Event Plane Angle [rad]"); fHistAllcentV0C->GetYaxis()->SetTitle("");
 242   fOutputList->Add(fHistAllcentV0C);
 243
 244   fHistAllcentTPC = new TH1F("fHistAllcentTPC","TPC Event Plane Angle of Pb+Pb",100,0.0,TMath::Pi());
 245   fHistAllcentTPC->GetXaxis()->SetTitle("TPC Event Plane Angle [rad]"); fHistAllcentTPC->GetYaxis()->SetTitle("");
 246   fOutputList->Add(fHistAllcentTPC);
 247
 248   fHistEPTPC = new TH2F("fHistEPTPC","TPC Event Plane Angle vs Centrality of Pb+Pb",100,0,100,100,0.0,TMath::Pi());
 249   fHistEPTPC->GetYaxis()->SetTitle("Event Plane Angle [rad]"); fHistEPTPC->GetXaxis()->SetTitle("V0M Centrality");
 250   fOutputList->Add(fHistEPTPC);
 251
 252   fHistEPTPCResolution = new TH2F("fHistEPTPCResolution","TPC Resolution vs Centrality of Pb+Pb",100,0,100,100,0.0,1.0);
 253   fHistEPTPCResolution->GetYaxis()->SetTitle("TPC Resolution"); fHistEPTPCResolution->GetXaxis()->SetTitle("V0M Centrality");
 254   fOutputList->Add(fHistEPTPCResolution);
 255
 256
 257                     //  Et   M02  V0Mcent DeltaPhi    Cos[2*DeltaPhi]
 258   Int_t    bins[5] = {  500, 250,  100,     100,          100  }; // binning
 259   Double_t min[5]  = {  0.0, 0.0,    0,     0.0,         -1.0}; // min x
 260   Double_t max[5]  = { 50.0, 2.5,  100,  TMath::Pi(),     1.0}; // max x
 261
 262   fClusterPbV0 = new THnSparseF("fClusterPbV0","",5,bins,min,max);
 263   fClusterPbV0->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0->GetAxis(1)->SetTitle("M02"); fClusterPbV0->GetAxis(2)->SetTitle("V0M Centrality");
 264   fClusterPbV0->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
 265   fOutputList->Add(fClusterPbV0);
 266
 267   fClusterPbV0A = new THnSparseF("fClusterPbV0A","",5,bins,min,max);
 268   fClusterPbV0A->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0A->GetAxis(1)->SetTitle("M02"); fClusterPbV0A->GetAxis(2)->SetTitle("V0M Centrality");
 269   fClusterPbV0A->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0A->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
 270   fOutputList->Add(fClusterPbV0A);
 271
 272   fClusterPbV0C = new THnSparseF("fClusterPbV0C","",5,bins,min,max);
 273   fClusterPbV0C->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbV0C->GetAxis(1)->SetTitle("M02"); fClusterPbV0C->GetAxis(2)->SetTitle("V0M Centrality");
 274   fClusterPbV0C->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbV0C->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
 275   fOutputList->Add(fClusterPbV0C);
 276
 277   fClusterPbTPC = new THnSparseF("fClusterPbTPC","",5,bins,min,max);
 278   fClusterPbTPC->GetAxis(0)->SetTitle("Transverse Energy [GeV]"); fClusterPbTPC->GetAxis(1)->SetTitle("M02"); fClusterPbTPC->GetAxis(2)->SetTitle("V0M Centrality");
 279   fClusterPbTPC->GetAxis(3)->SetTitle("Delta(#phi) [rad]"); fClusterPbTPC->GetAxis(4)->SetTitle("Cos[2*Delta(#phi)]");
 280   fOutputList->Add(fClusterPbTPC);
 281
 282
 283   PostData(1, fOutputList);
 284 }
 285
 286 //________________________________________________________________________
 287 void AliAnalysisTaskEMCALPi0V2ShSh::UserExec(Option_t *)
 288 {
 289   // Main loop, called for each event.
 290
 291   // Create pointer to reconstructed event
 292   AliVEvent *event = InputEvent();
 293   if (!event) { Printf("ERROR: Could not retrieve event\n"); return; }
 294
 295   fESD = dynamic_cast<AliESDEvent*>(event);
 296   if (!fESD) {
 297     printf("ERROR: Could not retrieve the ESD event\n");
 298     return;
 299   }
 300
 301   //Bool_t isSelected =0;
 302   //isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral));
 303   //if(!isSelected) { return; }
 304
 305   if(fESD->GetCentrality()){
 306      fCentralityV0M = fESD->GetCentrality()->GetCentralityPercentile("V0M");
 307   } else {
 308     printf("ERROR: Could not retrieve the Centrality\n");
 309     return;
 310   }
 311
 312   fEventPlane = fESD->GetEventplane();
 313   if (fEventPlane) {
 314      VZEROEventPlane();
 315   } else {
 316      printf("ERROR: Could not retrieve the Centrality\n");
 317      return;
 318   }
 319
 320   if(fESD->GetEMCALClusters(fCaloClusters)) { //select EMCAL clusters only
 321      FillClusterHists();
 322   } else{
 323      printf("ERROR: Could not get EMCAL Cluster\n");
 324      return;
 325   }
 326
 327   FillTrackHists();
 328
 329   fCaloClusters->Clear();
 330   PostData(1, fOutputList);
 331 }
 332
 333 //________________________________________________________________________
 334 void AliAnalysisTaskEMCALPi0V2ShSh::VZEROEventPlane()
 335 {// Calculate the V0 Event Plane
 336
 337       if (fEventPlane->GetQVector()) {
 338         fEPTPC = TVector2::Phi_0_2pi(fEventPlane->GetQVector()->Phi())/2.0; //if(fEPTPC>TMath::Pi()) {fEPTPC-=TMath::Pi();}
 339       } else { fEPTPC = -999.; }
 340
 341       if (fEventPlane->GetQsub1()&&fEventPlane->GetQsub2()) {
 342         fEPTPCResolution = TMath::Cos(2.0*(fEventPlane->GetQsub1()->Phi()/2.0-fEventPlane->GetQsub2()->Phi()/2.0)); }
 343       else { fEPTPCResolution = -1; }
 344
 345       fEPV0  = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0",  fESD)); if(fEPV0>TMath::Pi()) {fEPV0-=TMath::Pi();}
 346       fEPV0A = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0A", fESD)); if(fEPV0A>TMath::Pi()) {fEPV0A-=TMath::Pi();}
 347       fEPV0C = TVector2::Phi_0_2pi(fEventPlane->GetEventplane("V0C", fESD)); if(fEPV0C>TMath::Pi()) {fEPV0C-=TMath::Pi();}
 348
 349       Double_t qx=0, qy=0, qxr=0, qyr=0;
 350       fEPV0Ar = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 4, 5, 2, qxr, qyr)); if(fEPV0Ar>TMath::Pi()) {fEPV0Ar-=TMath::Pi();}
 351       fEPV0Cr = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 2, 3, 2, qx,  qy)); if(fEPV0Cr>TMath::Pi()) {fEPV0Cr-=TMath::Pi();}
 352       qxr += qx; qyr += qy;
 353       fEPV0r   = TVector2::Phi_0_2pi(TMath::ATan2(qyr,qxr))/2.0;
 354
 355       fEPV0A4r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 4, 2, qx, qy)); if(fEPV0A4r>TMath::Pi()) {fEPV0A4r-=TMath::Pi();}
 356       fEPV0A5r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 5, 2, qx, qy)); if(fEPV0A5r>TMath::Pi()) {fEPV0A5r-=TMath::Pi();}
 357       fEPV0A6r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 6, 2, qx, qy)); if(fEPV0A6r>TMath::Pi()) {fEPV0A6r-=TMath::Pi();}
 358       fEPV0A7r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 7, 2, qx, qy)); if(fEPV0A7r>TMath::Pi()) {fEPV0A7r-=TMath::Pi();}
 359       fEPV0C0r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 0, 2, qx, qy)); if(fEPV0C0r>TMath::Pi()) {fEPV0C0r-=TMath::Pi();}
 360       fEPV0C1r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 1, 2, qx, qy)); if(fEPV0C1r>TMath::Pi()) {fEPV0C1r-=TMath::Pi();}
 361       fEPV0C2r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 2, 2, qx, qy)); if(fEPV0C2r>TMath::Pi()) {fEPV0C2r-=TMath::Pi();}
 362       fEPV0C3r = TVector2::Phi_0_2pi(fEventPlane->CalculateVZEROEventPlane(fESD, 3, 2, qx, qy)); if(fEPV0C3r>TMath::Pi()) {fEPV0C3r-=TMath::Pi();}
 363
 364       fHistEPTPC->Fill(fCentralityV0M,  fEPTPC);
 365       if(fEPTPCResolution!=-1) { fHistEPTPCResolution->Fill(fCentralityV0M, fEPTPCResolution); }
 366       fHistEPV0->Fill(fCentralityV0M, fEPV0);
 367       fHistEPV0A->Fill(fCentralityV0M, fEPV0A);
 368       fHistEPV0C->Fill(fCentralityV0M, fEPV0C);
 369       fHistEPV0Ar->Fill(fCentralityV0M, fEPV0Ar);
 370       fHistEPV0Cr->Fill(fCentralityV0M, fEPV0Cr);
 371       fHistEPV0r->Fill(fCentralityV0M, fEPV0r);
 372       fHistEPV0A4r->Fill(fCentralityV0M, fEPV0A4r);
 373       fHistEPV0A7r->Fill(fCentralityV0M, fEPV0A7r);
 374       fHistEPV0C0r->Fill(fCentralityV0M, fEPV0C0r);
 375       fHistEPV0C3r->Fill(fCentralityV0M, fEPV0C3r);
 376
 377       fHistAllcentV0->Fill(fEPV0);
 378       fHistAllcentV0r->Fill(fEPV0r);
 379       fHistAllcentV0A->Fill(fEPV0A);
 380       fHistAllcentV0C->Fill(fEPV0C);
 381       fHistAllcentTPC->Fill(fEPTPC);
 382
 383       fHistdifV0A_V0C0r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A - fEPV0C0r)));
 384       fHistdifV0A_V0C3r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A - fEPV0C3r)));
 385       fHistdifV0C0r_V0C3r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C0r - fEPV0C3r)));
 386       fHistdifV0C_V0A4r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C - fEPV0A4r)));
 387       fHistdifV0C_V0A7r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0C - fEPV0A7r)));
 388       fHistdifV0A4r_V0A7r->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0A4r - fEPV0A7r)));
 389       fHistdifV0Ar_V0Cr->Fill(fCentralityV0M, TMath::Cos(2.0*(fEPV0Ar - fEPV0Cr)));
 390
 391 }
 392
 393 //________________________________________________________________________
 394 void AliAnalysisTaskEMCALPi0V2ShSh::FillClusterHists()
 395 {// Fill cluster histograms.
 396
 397   Float_t pos[3] ;
 398
 399   if(!fCaloClusters)
 400     return;
 401
 402   const Int_t nclus = fCaloClusters->GetEntries();
 403
 404   if(nclus==0)
 405     return;
 406
 407   for(Int_t iclus=0;iclus<nclus;iclus++){
 408     AliESDCaloCluster *clus = (AliESDCaloCluster *) fCaloClusters->At(iclus); // retrieve cluster from esd
 409     if(!clus)
 410       continue;
 411     if(!clus->IsEMCAL()){
 412       printf("ERROR: Current Cluster is not an EMCAL Cluster\n");
 413       continue;
 414     }
 415     clus->GetPosition(pos);
 416     TVector3 vpos(pos[0],pos[1],pos[2]);
 417     Double_t Transverse_Energy = ((clus->E())/ (TMath::CosH(vpos.Eta())));
 418     Double_t dphiV0 = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0r); if(dphiV0>TMath::Pi()) {dphiV0-=TMath::Pi();}
 419     Double_t dphiV0A = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0A); if(dphiV0A>TMath::Pi()) {dphiV0A-=TMath::Pi();}
 420     Double_t dphiV0C = TVector2::Phi_0_2pi(vpos.Phi()-fEPV0C); if(dphiV0C>TMath::Pi()) {dphiV0C-=TMath::Pi();}
 421     Double_t dphiTPC = TVector2::Phi_0_2pi(vpos.Phi()-fEPTPC); if(dphiTPC>TMath::Pi()) {dphiTPC-=TMath::Pi();}
 422
 423     Double_t DataV0[5];
 424     DataV0[0] = Transverse_Energy;
 425     DataV0[1] = clus->GetM02();
 426     DataV0[2] = fCentralityV0M;
 427     DataV0[3] = dphiV0;
 428     DataV0[4] = TMath::Cos(2.0*(dphiV0));
 429     fClusterPbV0->Fill(DataV0);
 430
 431     Double_t DataV0A[5];
 432     DataV0A[0] = Transverse_Energy;
 433     DataV0A[1] = clus->GetM02();
 434     DataV0A[2] = fCentralityV0M;
 435     DataV0A[3] = dphiV0A;
 436     DataV0A[4] = TMath::Cos(2.0*(dphiV0A));
 437     fClusterPbV0A->Fill(DataV0A);
 438
 439     Double_t DataV0C[5];
 440     DataV0C[0] = Transverse_Energy;
 441     DataV0C[1] = clus->GetM02();
 442     DataV0C[2] = fCentralityV0M;
 443     DataV0C[3] = dphiV0C;
 444     DataV0C[4] = TMath::Cos(2.0*(dphiV0C));
 445     fClusterPbV0C->Fill(DataV0C);
 446
 447     Double_t DataTPC[5];
 448     DataTPC[0] = Transverse_Energy;
 449     DataTPC[1] = clus->GetM02();
 450     DataTPC[2] = fCentralityV0M;
 451     DataTPC[3] = dphiTPC;
 452     DataTPC[4] = TMath::Cos(2.0*(dphiTPC));
 453     fClusterPbTPC->Fill(DataTPC);
 454
 455     fHistClusterE->Fill(clus->E());
 456     fHistClusterEt->Fill(Transverse_Energy);
 457     fHistClusterM02->Fill(clus->GetM02());
 458     fHistClusterN->Fill(clus->GetNCells());
 459     fHistClusterPhi->Fill(vpos.Phi());
 460     fHistClusterEta->Fill(vpos.Eta());
 461     fHistClusterPhiEta->Fill(vpos.Eta(),vpos.Phi());
 462     fHistClusterEN->Fill(clus->E(),clus->GetNCells());
 463     fHistClusterEM02->Fill(clus->E(),clus->GetM02());
 464     fHistClusterEtN->Fill(Transverse_Energy,clus->GetNCells());
 465     fHistClusterEtM02->Fill(Transverse_Energy,clus->GetM02());
 466     fHistClusterdphiV0->Fill(dphiV0);
 467
 468   }
 469 }
 470
 471 //________________________________________________________________________
 472 void AliAnalysisTaskEMCALPi0V2ShSh::FillTrackHists()
 473 {// Fill track histograms.
 474
 475   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 476     AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 477     if (!track){
 478       printf("ERROR: Could not retreive esdtrack\n");
 479       continue;
 480     }
 481     fHistTrackPt->Fill(track->Pt());
 482     fHistTrackEta->Fill(track->Eta());
 483     fHistTrackPhi->Fill(track->Phi());
 484     fHistTrackPhiEta->Fill(track->Eta(),track->Phi());
 485
 486   }
 487 }
 488
 489 //________________________________________________________________________
 490 void AliAnalysisTaskEMCALPi0V2ShSh::Terminate(Option_t *)
 491 {
 492   // Draw result to screen, or perform fitting, normalizations
 493   // Called once at the end of the query
 494
 495  /* fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 496   if(!fOutputList) { Printf("ERROR: could not retrieve TList fOutputList"); return; }
 497
 498   fHistTrackPt = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackPt"));
 499   if (!fHistTrackPt) { Printf("ERROR: could not retrieve fHistTrackPt"); return;}
 500   fHistTrackEta = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackEta"));
 501   if (!fHistTrackEta) { Printf("ERROR: could not retrieve fHistTrackEta"); return;}
 502   fHistTrackPhi = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistTrackPhi"));
 503   if (!fHistTrackPhi) { Printf("ERROR: could not retrieve fHistTrackPhi"); return;}
 504   fHistClusterEta = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterEta"));
 505   if (!fHistClusterEta) { Printf("ERROR: could not retrieve fHistClusterEta"); return;}
 506   fHistClusterPhi = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterPhi"));
 507   if (!fHistClusterPhi) { Printf("ERROR: could not retrieve fHistClusterPhi"); return;}
 508   fHistClusterE = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterE"));
 509   if (!fHistClusterE) { Printf("ERROR: could not retrieve fHistClusterE"); return;}
 510   fHistClusterN = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterN"));
 511   if (!fHistClusterN) { Printf("ERROR: could not retrieve fHistClusterN"); return;}
 512   fHistClusterM02 = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistClusterM02"));
 513   if (!fHistClusterM02) { Printf("ERROR: could not retrieve fHistClusterM02"); return;}
 514
 515     // Get the physics selection histograms with the selection statistics
 516     //AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 517     //AliESDInputHandler *inputH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
 518     //TH2F *histStat = (TH2F*)inputH->GetStatistics();
 519
 520
 521   TCanvas *canvas1 = new TCanvas("canvas1","Track P_{T} & #eta & #phi",10,10,1020,510);
 522   canvas1->Divide(3,1);
 523   canvas1->cd(1)->SetLogy();
 524   fHistTrackPt->DrawCopy("");
 525   canvas1->cd(2);
 526   fHistTrackEta->DrawCopy("");
 527   canvas1->cd(3);
 528   fHistTrackPhi->DrawCopy("");
 529
 530   TCanvas *canvas2 = new TCanvas("canvas2","Cluster #eta & #phi",10,10,1020,510);
 531   canvas2->Divide(2,1);
 532   canvas2->cd(1);
 533   fHistClusterEta->DrawCopy("");
 534   canvas2->cd(2);
 535   fHistClusterPhi->DrawCopy("");
 536
 537   TCanvas *canvas3 = new TCanvas("canvas3","Cluster E & N & M02",10,10,1020,510);
 538   canvas3->Divide(3,1);
 539   canvas3->cd(1)->SetLogy();
 540   fHistClusterE->DrawCopy("");
 541   canvas3->cd(2)->SetLogy();
 542   fHistClusterN->DrawCopy("");
 543   canvas3->cd(3)->SetLogy();
 544   fHistClusterM02->DrawCopy("");
 545
 546   TCanvas *canvas4 = new TCanvas("canvas4","Track #phi vs #eta &  Cluster #phi vs #eta",10,10,1020,510);
 547   canvas4->Divide(2,1);
 548   canvas4->cd(1)->SetLogz();
 549   fHistTrackPhiEta->DrawCopy("colorz");
 550   canvas4->cd(2)->SetLogz();
 551   fHistClusterPhiEta->DrawCopy("colorz");
 552
 553   TCanvas *canvas5 = new TCanvas("canvas5","Cluster E vs N &  E vs M02",10,10,1020,510);
 554   canvas5->Divide(2,1);
 555   canvas5->cd(1)->SetLogz();
 556   fHistClusterEN->DrawCopy("colorz");
 557   canvas5->cd(2)->SetLogz();
 558   fHistClusterEM02->DrawCopy("colorz");
 559
 560   TCanvas *canvas6 = new TCanvas("canvas6","Cluster Et & Et vs N & Et vs M02",10,10,1020,510);
 561   canvas6->cd(1)->SetLogy();
 562   fHistClusterEt->DrawCopy("");
 563   canvas6->cd(2)->SetLogz();
 564   fHistClusterEtN->DrawCopy("colorz");
 565   canvas6->cd(3)->SetLogz();
 566   fHistClusterEtM02->DrawCopy("colorz");
 567
 568   canvas1->SaveAs("lhc11h_2_Track_PT_Eta_Phi.jpg");
 569   canvas2->SaveAs("lhc11h_2_Cluster_Eta_Phi.jpg");
 570   canvas3->SaveAs("lhc11h_2_Cluster_E_N_M02.jpg");
 571   canvas4->SaveAs("lhc11h_2_Track_PhivsEta_Cluster_PhivsEta.jpg");
 572   canvas5->SaveAs("lhc11h_2_Cluster_EvsN_EvsM02.jpg");
 573   canvas6->SaveAs("lhc11h_2_Cluster_Et_EtvsN_EtvsM02.jpg"); */
 574 }