PWGJE/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 <TH2F.h>
  12 #include <TLorentzVector.h>
  13 #include <TMath.h>
  14
  15 #include "AliEMCALGeometry.h"
  16 #include "AliLog.h"
  17 #include "AliVCluster.h"
  18 #include "AliVParticle.h"
  19 #include "AliParticleContainer.h"
  20
  21 ClassImp(AliAnalysisTaskScale)
  22
  23 //________________________________________________________________________
  24 AliAnalysisTaskScale::AliAnalysisTaskScale() :
  25   AliAnalysisTaskEmcalDev("AliAnalysisTaskScale", kTRUE),
  26   fScaleFunction(0),
  27   fEmcalArea(1),
  28   fTpcArea(1),
  29   fHistPtTPCvsCent(0),
  30   fHistPtEMCALvsCent(0),
  31   fHistEtvsCent(0),
  32   fHistScalevsCent(0),
  33   fHistDeltaScalevsCent(0),
  34   fHistScaleEmcalvsCent(0),
  35   fHistScale2EmcalvsCent(0),
  36   fHistChScalevsCent(0),
  37   fHistChScale2EmcalvsCent(0),
  38   fHistPtTPCvsNtrack(0),
  39   fHistPtEMCALvsNtrack(0),
  40   fHistEtvsNtrack(0),
  41   fHistScalevsNtrack(0),
  42   fHistDeltaScalevsNtrack(0),
  43   fHistScaleEmcalvsNtrack(0),
  44   fHistScale2EmcalvsNtrack(0),
  45   fHistChScalevsNtrack(0),
  46   fHistChScale2EmcalvsNtrack(0),
  47   fHistTrackPtvsCent(0),
  48   fHistClusterPtvsCent(0),
  49   fHistTrackEtaPhi(0),
  50   fHistClusterEtaPhi(0),
  51   fHistScalevsScale2Emcal(0),
  52   fHistScalevsScaleEmcal(0),
  53   fHistScaleEmcalvsScale2Emcal(0)
  54 {
  55   // Default constructor.
  56
  57   SetMakeGeneralHistograms(kTRUE);
  58 }
  59
  60 //________________________________________________________________________
  61 AliAnalysisTaskScale::AliAnalysisTaskScale(const char *name) :
  62   AliAnalysisTaskEmcalDev(name, kTRUE),
  63   fScaleFunction(0),
  64   fEmcalArea(1),
  65   fTpcArea(1),
  66   fHistPtTPCvsCent(0),
  67   fHistPtEMCALvsCent(0),
  68   fHistEtvsCent(0),
  69   fHistScalevsCent(0),
  70   fHistDeltaScalevsCent(0),
  71   fHistScaleEmcalvsCent(0),
  72   fHistScale2EmcalvsCent(0),
  73   fHistChScalevsCent(0),
  74   fHistChScale2EmcalvsCent(0),
  75   fHistPtTPCvsNtrack(0),
  76   fHistPtEMCALvsNtrack(0),
  77   fHistEtvsNtrack(0),
  78   fHistScalevsNtrack(0),
  79   fHistDeltaScalevsNtrack(0),
  80   fHistScaleEmcalvsNtrack(0),
  81   fHistScale2EmcalvsNtrack(0),
  82   fHistChScalevsNtrack(0),
  83   fHistChScale2EmcalvsNtrack(0),
  84   fHistTrackPtvsCent(0),
  85   fHistClusterPtvsCent(0),
  86   fHistTrackEtaPhi(0),
  87   fHistClusterEtaPhi(0),
  88   fHistScalevsScale2Emcal(0),
  89   fHistScalevsScaleEmcal(0),
  90   fHistScaleEmcalvsScale2Emcal(0)
  91 {
  92   // Constructor.
  93
  94   SetMakeGeneralHistograms(kTRUE);
  95 }
  96
  97 //________________________________________________________________________
  98 void AliAnalysisTaskScale::UserCreateOutputObjects()
  99 {
 100   // Create my user objects.
 101
 102   AliAnalysisTaskEmcalDev::UserCreateOutputObjects();
 103
 104   fHistPtTPCvsCent             = new TH2F("PtTPCvsCent","rho vs cent",            101, -1, 100,   500,   0, 1000);
 105   fHistPtEMCALvsCent           = new TH2F("PtEMCALvsCent","rho vs cent",          101, -1, 100,   500,   0, 1000);
 106   fHistEtvsCent                = new TH2F("EtvsCent","rho vs cent",               101, -1, 100,   500,   0, 1000);
 107   fHistScalevsCent             = new TH2F("ScalevsCent","rho vs cent",            101, -1, 100,   500,   0, 5);
 108   fHistDeltaScalevsCent        = new TH2F("DeltaScalevsCent","rho vs cent",       101, -1, 100,   500,  -2.5, 2.5);
 109   fHistScaleEmcalvsCent        = new TH2F("ScaleEmcalvsCent","",                  101, -1, 100,   500,   0, 5);
 110   fHistScale2EmcalvsCent       = new TH2F("Scale2EmcalvsCent","",                 101, -1, 100,   500,   0, 5);
 111   fHistChScalevsCent           = new TH2F("ChScalevsCent","",                     101, -1, 100,   500,   0, 5);
 112   fHistChScale2EmcalvsCent     = new TH2F("ChScale2EmcalvsCent","",               101, -1, 100,   500,   0, 5);
 113   fHistPtTPCvsNtrack           = new TH2F("PtTPCvsNtrack","rho vs cent",          500,  0, 2500,  500,   0, 1000);
 114   fHistPtEMCALvsNtrack         = new TH2F("PtEMCALvsNtrack","rho vs cent",        500,  0, 2500,  500,   0, 1000);
 115   fHistEtvsNtrack              = new TH2F("EtvsNtrack","rho vs cent",             500,  0, 2500,  500,   0, 1000);
 116   fHistScalevsNtrack           = new TH2F("ScalevsNtrack","rho vs cent",          500,  0, 2500,  500,   0, 5);
 117   fHistDeltaScalevsNtrack      = new TH2F("DeltaScalevsNtrack","rho vs cent",     500,  0, 2500,  500,  -2.5, 2.5);
 118   fHistScaleEmcalvsNtrack      = new TH2F("ScaleEmcalvsNtrack","",                500,  0, 2500,  500,   0, 5);
 119   fHistScale2EmcalvsNtrack     = new TH2F("Scale2EmcalvsNtrack","",               500,  0, 2500,  500,   0, 5);
 120   fHistChScalevsNtrack         = new TH2F("ChScalevsNtrack","",                   500,  0, 2500,  500,   0, 5);
 121   fHistChScale2EmcalvsNtrack   = new TH2F("ChScale2EmcalvsNtrack","",             500,  0, 2500,  500,   0, 5);
 122   fHistTrackPtvsCent           = new TH2F("TrackPtvsCent","Track pt vs cent",     101, -1, 100,   500,   0, 100);
 123   fHistClusterPtvsCent         = new TH2F("ClusterPtvsCent","Cluster pt vs cent", 101, -1, 100,   500,   0, 100);
 124   fHistTrackEtaPhi             = new TH2F("TrackEtaPhi","Track eta phi",          100, -1.0, 1.0, 101,   0, 2.02*TMath::Pi());
 125   fHistClusterEtaPhi           = new TH2F("ClusterEtaPhi","Cluster eta phi",      100, -1.0, 1.0, 101,   0, 2.02*TMath::Pi());
 126   fHistScalevsScale2Emcal      = new TH2F("ScalevsScale2Emcal","",                500,  0, 5,     500,   0, 5);
 127   fHistScalevsScaleEmcal       = new TH2F("ScalevsScaleEmcal","",                 500,  0, 5,     500,   0, 5);
 128   fHistScaleEmcalvsScale2Emcal = new TH2F("ScaleEmcalvsScale2Emcal","",           500,  0, 5,     500,   0, 5);
 129
 130   fOutput->Add(fHistPtTPCvsCent);
 131   fOutput->Add(fHistPtEMCALvsCent);
 132   fOutput->Add(fHistEtvsCent);
 133   fOutput->Add(fHistScalevsCent);
 134   fOutput->Add(fHistDeltaScalevsCent);
 135   fOutput->Add(fHistScaleEmcalvsCent);
 136   fOutput->Add(fHistScale2EmcalvsCent);
 137   fOutput->Add(fHistChScalevsCent);
 138   fOutput->Add(fHistChScale2EmcalvsCent);
 139   fOutput->Add(fHistPtTPCvsNtrack);
 140   fOutput->Add(fHistPtEMCALvsNtrack);
 141   fOutput->Add(fHistEtvsNtrack);
 142   fOutput->Add(fHistScalevsNtrack);
 143   fOutput->Add(fHistDeltaScalevsNtrack);
 144   fOutput->Add(fHistScaleEmcalvsNtrack);
 145   fOutput->Add(fHistScale2EmcalvsNtrack);
 146   fOutput->Add(fHistChScalevsNtrack);
 147   fOutput->Add(fHistChScale2EmcalvsNtrack);
 148   fOutput->Add(fHistTrackPtvsCent);
 149   fOutput->Add(fHistClusterPtvsCent);
 150   fOutput->Add(fHistTrackEtaPhi);
 151   fOutput->Add(fHistClusterEtaPhi);
 152   fOutput->Add(fHistScalevsScale2Emcal);
 153   fOutput->Add(fHistScalevsScaleEmcal);
 154   fOutput->Add(fHistScaleEmcalvsScale2Emcal);
 155
 156   PostData(1, fOutput);
 157 }
 158
 159 //________________________________________________________________________
 160 Double_t AliAnalysisTaskScale::GetScaleFactor(Double_t cent)
 161 {
 162   // Get scale function.
 163
 164   Double_t scale = -1;
 165   if (fScaleFunction)
 166     scale = fScaleFunction->Eval(cent);
 167   return scale;
 168 }
 169
 170 //________________________________________________________________________
 171 Bool_t AliAnalysisTaskScale::FillHistograms()
 172 {
 173   // Execute on each event.
 174
 175   const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 176   const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 177   const Double_t EmcalWidth = (EmcalMaxPhi-EmcalMinPhi)/2.0;
 178
 179   Double_t ptTPC   = 0;
 180   Double_t ptEMCAL = 0;
 181   Double_t ptEMCAL2 = 0;
 182
 183   const Int_t Ntracks = fTracks->GetEntries();
 184   for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 185     AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(iTracks));
 186
 187     if (!track)
 188       continue;
 189
 190     if (!AcceptTrack(track))
 191       continue;
 192
 193     if (TMath::Abs(track->Eta()) > 0.7)   // only accept tracks in the EMCal eta range
 194       continue;
 195
 196     fHistTrackPtvsCent->Fill(fCent,track->Pt());
 197     fHistTrackEtaPhi->Fill(track->Eta(),track->Phi());
 198     ptTPC += track->Pt();
 199     if ((track->Phi() > (EmcalMaxPhi+EmcalWidth)) || (track->Phi() < (EmcalMinPhi-EmcalWidth))) continue;
 200     ptEMCAL2 += track->Pt();
 201     if ((track->Phi() > EmcalMaxPhi) || (track->Phi() < EmcalMinPhi)) continue;
 202     ptEMCAL += track->Pt();
 203   }
 204
 205   if (ptTPC == 0)
 206     return kFALSE;
 207
 208   Double_t Et = 0;
 209   const Int_t Nclus = fCaloClusters->GetEntries();
 210   for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 211     AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(iClus));
 212     if (!c)
 213       continue;
 214
 215     if (!AcceptCluster(c))
 216       continue;
 217
 218     TLorentzVector nPart;
 219     c->GetMomentum(nPart, fVertex);
 220
 221     fHistClusterPtvsCent->Fill(fCent, nPart.Pt());
 222     fHistClusterEtaPhi->Fill(nPart.Eta(), nPart.Phi());
 223
 224     Et += nPart.Pt();
 225   }
 226
 227   Double_t scalecalc         = -1;
 228   if (ptEMCAL > 0 && Et > 0 && ptTPC > 0)
 229     scalecalc         =  ((Et + ptEMCAL) / fEmcalArea) * (fTpcArea / ptTPC);
 230   const Double_t scale             = GetScaleFactor(fCent);
 231   Double_t scalecalcemcal          = -1;
 232   if (ptEMCAL > 0)
 233     scalecalcemcal                 = (Et+ptEMCAL)/ptEMCAL;
 234   Double_t scalecalcemcal2         = -1;
 235   Double_t Chscalecalcemcal2       = -1;
 236   if (ptEMCAL2 > 0){
 237     scalecalcemcal2                = 2*(Et+ptEMCAL)/ptEMCAL2;
 238     Chscalecalcemcal2              = 2*ptEMCAL/ptEMCAL2;}
 239   const Double_t Chscalecalcemcal  = ((ptEMCAL) / fEmcalArea) * (fTpcArea / ptTPC);
 240
 241   fHistScaleEmcalvsCent->Fill(fCent,scalecalcemcal);
 242   fHistScale2EmcalvsCent->Fill(fCent,scalecalcemcal2);
 243   fHistChScalevsCent->Fill(fCent,Chscalecalcemcal);
 244   fHistChScale2EmcalvsCent->Fill(fCent,Chscalecalcemcal2);
 245   fHistScaleEmcalvsNtrack->Fill(Ntracks,scalecalcemcal);
 246   fHistScale2EmcalvsNtrack->Fill(Ntracks,scalecalcemcal2);
 247   fHistChScalevsNtrack->Fill(Ntracks,Chscalecalcemcal);
 248   fHistChScale2EmcalvsNtrack->Fill(Ntracks,Chscalecalcemcal2);
 249   fHistPtTPCvsCent->Fill(fCent, ptTPC);
 250   fHistPtEMCALvsCent->Fill(fCent, ptEMCAL);
 251   fHistEtvsCent->Fill(fCent, Et);
 252   fHistScalevsCent->Fill(fCent, scalecalc);
 253   fHistDeltaScalevsCent->Fill(fCent, scalecalc - scale);
 254   fHistPtTPCvsNtrack->Fill(Ntracks, ptTPC);
 255   fHistPtEMCALvsNtrack->Fill(Ntracks, ptEMCAL);
 256   fHistEtvsNtrack->Fill(Ntracks, Et);
 257   fHistScalevsNtrack->Fill(Ntracks, scalecalc);
 258   fHistDeltaScalevsNtrack->Fill(Ntracks, scalecalc - scale);
 259   fHistScalevsScale2Emcal->Fill(scalecalc,scalecalcemcal2);
 260   fHistScalevsScaleEmcal->Fill(scalecalc,scalecalcemcal);
 261   fHistScaleEmcalvsScale2Emcal->Fill(scalecalcemcal,scalecalcemcal2);
 262
 263   return kTRUE;
 264 }
 265
 266
 267 //________________________________________________________________________
 268 void AliAnalysisTaskScale::ExecOnce()
 269 {
 270   AliAnalysisTaskEmcalDev::ExecOnce();
 271
 272   const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 273   const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 274   const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 275   const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 276
 277   fEmcalArea  = (EmcalMaxPhi - EmcalMinPhi) * (EmcalMinEta - EmcalMaxEta);
 278
 279   AliParticleContainer *partCont = GetParticleContainer(0);
 280   if (!partCont) {
 281     AliError(Form("%s: No particle container found! Assuming tpc area = 1...",GetName()));
 282     fTpcArea = 1;
 283     return;
 284   }
 285
 286   Float_t TpcMaxPhi = partCont->GetParticlePhiMax();
 287   Float_t TpcMinPhi = partCont->GetParticlePhiMin();
 288
 289   if (TpcMaxPhi > TMath::Pi()*2) TpcMaxPhi = TMath::Pi()*2;
 290   if (TpcMinPhi < 0) TpcMinPhi = 0;
 291
 292   fTpcArea = (TpcMaxPhi - TpcMinPhi) * (EmcalMinEta - EmcalMaxEta);
 293 }