PWGGA/EMCALJetTasks/AliAnalysisTaskScale.cxx

   1 // $Id$
   2 //
   3 // Scale task.
   4 //
   5 // Author: R.Reed, M.Connors
   6
   7 #include "AliAnalysisTaskScale.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TF1.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TList.h>
  14 #include <TLorentzVector.h>
  15
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliLog.h"
  20 #include "AliVCluster.h"
  21 #include "AliVEvent.h"
  22 #include "AliVTrack.h"
  23
  24 ClassImp(AliAnalysisTaskScale)
  25
  26 //________________________________________________________________________
  27 AliAnalysisTaskScale::AliAnalysisTaskScale() :
  28   AliAnalysisTaskSE(),
  29   fTracksName(),
  30   fClustersName(),
  31   fMinTrackPt(0.15),
  32   fMinClusterPt(0.15),
  33   fScaleFunction(0),
  34   fGeom(0),
  35   fOutputList(0),
  36   fHistCentrality(0),
  37   fHistPtTPCvsCent(0),
  38   fHistPtEMCALvsCent(0),
  39   fHistEtvsCent(0),
  40   fHistScalevsCent(0),
  41   fHistDeltaScalevsCent(0),
  42   fHistPtTPCvsNtrack(0),
  43   fHistPtEMCALvsNtrack(0),
  44   fHistEtvsNtrack(0),
  45   fHistScalevsNtrack(0),
  46   fHistDeltaScalevsNtrack(0),
  47   fHistTrackPtvsCent(0),
  48   fHistClusterPtvsCent(0),
  49   fHistTrackEtaPhi(0),
  50   fHistClusterEtaPhi(0)
  51 {
  52   // Default constructor.
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskScale::AliAnalysisTaskScale(const char *name) :
  57   AliAnalysisTaskSE(name),
  58   fTracksName("Tracks"),
  59   fClustersName("CaloClusters"),
  60   fMinTrackPt(0.15),
  61   fMinClusterPt(0.15),
  62   fScaleFunction(0),
  63   fGeom(0),
  64   fOutputList(0),
  65   fHistCentrality(0),
  66   fHistPtTPCvsCent(0),
  67   fHistPtEMCALvsCent(0),
  68   fHistEtvsCent(0),
  69   fHistScalevsCent(0),
  70   fHistDeltaScalevsCent(0),
  71   fHistPtTPCvsNtrack(0),
  72   fHistPtEMCALvsNtrack(0),
  73   fHistEtvsNtrack(0),
  74   fHistScalevsNtrack(0),
  75   fHistDeltaScalevsNtrack(0),
  76   fHistTrackPtvsCent(0),
  77   fHistClusterPtvsCent(0),
  78   fHistTrackEtaPhi(0),
  79   fHistClusterEtaPhi(0)
  80 {
  81   // Constructor.
  82
  83   DefineOutput(1, TList::Class());
  84 }
  85
  86 //________________________________________________________________________
  87 void AliAnalysisTaskScale::UserCreateOutputObjects()
  88 {
  89   // Create my user objects.
  90
  91   OpenFile(1);
  92   fOutputList = new TList();
  93   fOutputList->SetOwner();
  94
  95   fHistCentrality         = new TH1F("Centrality","Centrality",              101, -1, 100);
  96   fHistPtTPCvsCent        = new TH2F("PtTPCvsCent","rho vs cent",            101, -1, 100, 500,   0, 1000);
  97   fHistPtEMCALvsCent      = new TH2F("PtEMCALvsCent","rho vs cent",          101, -1, 100, 500,   0, 1000);
  98   fHistEtvsCent           = new TH2F("EtvsCent","rho vs cent",               101, -1, 100, 500,   0, 1000);
  99   fHistScalevsCent        = new TH2F("ScalevsCent","rho vs cent",            101, -1, 100, 400,   0, 4);
 100   fHistDeltaScalevsCent   = new TH2F("DeltaScalevsCent","rho vs cent",       101, -1, 100, 400,  -2, 2);
 101   fHistPtTPCvsNtrack      = new TH2F("PtTPCvsNtrack","rho vs cent",          500,  0, 2500, 500,  0, 1000);
 102   fHistPtEMCALvsNtrack    = new TH2F("PtEMCALvsNtrack","rho vs cent",        500,  0, 2500, 500,  0, 1000);
 103   fHistEtvsNtrack         = new TH2F("EtvsNtrack","rho vs cent",             500,  0, 2500, 500,  0, 1000);
 104   fHistScalevsNtrack      = new TH2F("ScalevsNtrack","rho vs cent",          500,  0, 2500, 400,  0, 4);
 105   fHistDeltaScalevsNtrack = new TH2F("DeltaScalevsNtrack","rho vs cent",     500,  0, 2500, 400, -2, 2);
 106   fHistTrackPtvsCent      = new TH2F("TrackPtvsCent","Track pt vs cent",     101, -1, 100,  500,  0, 100);
 107   fHistClusterPtvsCent    = new TH2F("ClusterPtvsCent","Cluster pt vs cent", 101, -1, 100,  500,  0, 100);
 108   fHistTrackEtaPhi        = new TH2F("TrackEtaPhi","Track eta phi",          100, -1.0, 1.0, 64,  0, 6.4);
 109   fHistClusterEtaPhi      = new TH2F("ClusterEtaPhi","Cluster eta phi",      100, -1.0, 1.0, 64, -3.2, 3.2);
 110
 111   fOutputList->Add(fHistCentrality);
 112   fOutputList->Add(fHistPtTPCvsCent);
 113   fOutputList->Add(fHistPtEMCALvsCent);
 114   fOutputList->Add(fHistEtvsCent);
 115   fOutputList->Add(fHistScalevsCent);
 116   fOutputList->Add(fHistDeltaScalevsCent);
 117   fOutputList->Add(fHistPtTPCvsNtrack);
 118   fOutputList->Add(fHistPtEMCALvsNtrack);
 119   fOutputList->Add(fHistEtvsNtrack);
 120   fOutputList->Add(fHistScalevsNtrack);
 121   fOutputList->Add(fHistDeltaScalevsNtrack);
 122   fOutputList->Add(fHistTrackPtvsCent);
 123   fOutputList->Add(fHistClusterPtvsCent);
 124   fOutputList->Add(fHistTrackEtaPhi);
 125   fOutputList->Add(fHistClusterEtaPhi);
 126
 127   PostData(1, fOutputList);
 128 }
 129
 130 //________________________________________________________________________
 131 Double_t AliAnalysisTaskScale::GetScaleFactor(Double_t cent)
 132 {
 133   // Get scale function.
 134
 135   Double_t scale = -1;
 136   if (fScaleFunction)
 137     scale = fScaleFunction->Eval(cent);
 138   return scale;
 139 }
 140
 141 //________________________________________________________________________
 142 void AliAnalysisTaskScale::UserExec(Option_t *)
 143 {
 144   // Execute on each event.
 145
 146   if (!fGeom)
 147     fGeom = AliEMCALGeometry::GetInstance();
 148
 149   if (!fGeom) {
 150     AliFatal("Can not create geometry");
 151     return;
 152   }
 153
 154   // get input collections
 155   TClonesArray *tracks = 0;
 156   TClonesArray *clusters = 0;
 157   TList *l = InputEvent()->GetList();
 158   tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
 159   if (!tracks) {
 160     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 161     return;
 162   }
 163   clusters = dynamic_cast<TClonesArray*>(l->FindObject(fClustersName));
 164   if (!clusters){
 165     AliError(Form("Pointer to clusters %s == 0", fClustersName.Data() ));
 166     return;
 167   }
 168
 169   // get centrality
 170   Double_t cent = -1;
 171   AliCentrality *centrality = InputEvent()->GetCentrality();
 172   if (centrality)
 173     cent = centrality->GetCentralityPercentile("V0M");
 174   if (cent < 0) {
 175     AliError(Form("Centrality negative: %f", cent));
 176     return;
 177   }
 178   fHistCentrality->Fill(cent);
 179
 180   // get vertex
 181   Double_t vertex[3] = {0, 0, 0};
 182   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 183
 184   const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 185   const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 186   const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 187   const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 188
 189   const Double_t TpcMinPhi   = 0;
 190   const Double_t TpcMaxPhi   = 2 * TMath::Pi();
 191
 192   const Double_t TpcArea     = (TpcMaxPhi - TpcMinPhi) * (EmcalMinEta - EmcalMaxEta);
 193   const Double_t EmcalArea   = (EmcalMaxPhi - EmcalMinPhi) * (EmcalMinEta - EmcalMaxEta);
 194
 195   Double_t ptTPC   = 0;
 196   Double_t ptEMCAL = 0;
 197
 198   const Int_t Ntracks = tracks->GetEntries();
 199   for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 200     AliVTrack *track = static_cast<AliVTrack*>(tracks->At(iTracks));
 201
 202     if (!track)
 203       continue;
 204
 205     if (TMath::Abs(track->Eta()) > 0.7)   // only accept tracks in the EMCal eta range
 206       continue;
 207
 208     fHistTrackPtvsCent->Fill(cent,track->Pt());
 209     fHistTrackEtaPhi->Fill(track->Eta(),track->Phi());
 210
 211     if (track->Pt()< fMinTrackPt)
 212       continue;
 213
 214     ptTPC += track->Pt();
 215     if ((track->Phi() > EmcalMaxPhi) || (track->Phi() < EmcalMinPhi))
 216       continue;
 217
 218     ptEMCAL += track->Pt();
 219   }
 220
 221   Double_t Et = 0;
 222   const Int_t Nclus = clusters->GetEntries();
 223   for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 224     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(iClus));
 225     if (!c)
 226       continue;
 227     if (!c->IsEMCAL())
 228       continue;
 229     TLorentzVector nPart;
 230     c->GetMomentum(nPart, vertex);
 231
 232     fHistClusterPtvsCent->Fill(cent,nPart.Pt());
 233     fHistClusterEtaPhi->Fill(nPart.Eta(),nPart.Phi());
 234
 235     if (nPart.Pt()< fMinClusterPt)
 236       continue;
 237
 238     Et += nPart.Pt();
 239   }
 240
 241   if (ptTPC == 0)
 242     return;
 243
 244   const Double_t scalecalc = ((Et + ptEMCAL) / EmcalArea) * (TpcArea / ptTPC);
 245
 246   const Double_t scale     = GetScaleFactor(cent);
 247
 248   fHistPtTPCvsCent->Fill(cent, ptTPC);
 249   fHistPtEMCALvsCent->Fill(cent, ptEMCAL);
 250   fHistEtvsCent->Fill(cent, Et);
 251   fHistScalevsCent->Fill(cent, scalecalc);
 252   fHistDeltaScalevsCent->Fill(cent, scalecalc - scale);
 253   fHistPtTPCvsNtrack->Fill(Ntracks, ptTPC);
 254   fHistPtEMCALvsNtrack->Fill(Ntracks, ptEMCAL);
 255   fHistEtvsNtrack->Fill(Ntracks, Et);
 256   fHistScalevsNtrack->Fill(Ntracks, scalecalc);
 257   fHistDeltaScalevsNtrack->Fill(Ntracks, scalecalc - scale);
 258
 259   PostData(1, fOutputList);
 260 }
 261
 262 //________________________________________________________________________
 263 void AliAnalysisTaskScale::Terminate(Option_t *)
 264 {
 265   // Called once at the end of the analysis.
 266 }