PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx

   1 // $Id$
   2 //
   3 // Calculation of rho from a collection of jets.
   4 // If scale function is given the scaled rho will be exported
   5 // with the name as "fRhoName".apppend("_Scaled").
   6 //
   7 // Authors: R.Reed, S.Aiola
   8
   9 #include "AliAnalysisTaskRho.h"
  10
  11 #include <TClonesArray.h>
  12 #include <TF1.h>
  13 #include <TH1F.h>
  14 #include <TH2F.h>
  15 #include <TList.h>
  16 #include <TMath.h>
  17
  18 #include "AliAnalysisManager.h"
  19 #include "AliCentrality.h"
  20 #include "AliEmcalJet.h"
  21 #include "AliLog.h"
  22 #include "AliRhoParameter.h"
  23 #include "AliVCluster.h"
  24 #include "AliVEventHandler.h"
  25
  26 ClassImp(AliAnalysisTaskRho)
  27
  28 //________________________________________________________________________
  29 AliAnalysisTaskRho::AliAnalysisTaskRho() :
  30   AliAnalysisTaskRhoBase(),
  31   fTracksName("tracks"),
  32   fJetsName("KtJets"),
  33   fRhoScaledName(""),
  34   fPhiMin(0),
  35   fPhiMax(0),
  36   fEtaMin(0),
  37   fEtaMax(0),
  38   fAreaCut(0),
  39   fNExclLeadJets(0),
  40   fScaleFunction(0),
  41   fCreateHisto(kFALSE),
  42   fTracks(0),
  43   fJets(0),
  44   fOutputList(0),
  45   fHistCentrality(0),
  46   fHistJetPt(0),
  47   fHistJetArea(0),
  48   fHistRhovsCent(0),
  49   fHistDeltaRhovsCent(0),
  50   fHistDeltaRhoScalevsCent(0),
  51   fHistJetPtvsCent(0),
  52   fHistJetAreavsCent(0),
  53   fHistNjetvsCent(0),
  54   fHistRhovsNtrack(0),
  55   fHistDeltaRhovsNtrack(0),
  56   fHistDeltaRhoScalevsNtrack(0),
  57   fHistJetPtvsNtrack(0),
  58   fHistJetAreavsNtrack(0),
  59   fHistNjetvsNtrack(0),
  60   fRhoScaled(0)
  61 {
  62   // Constructor
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
  67   AliAnalysisTaskRhoBase(name),
  68   fTracksName("tracks"),
  69   fJetsName("KtJets"),
  70   fRhoScaledName(""),
  71   fPhiMin(0),
  72   fPhiMax(TMath::TwoPi()),
  73   fEtaMin(-0.5),
  74   fEtaMax(+0.5),
  75   fAreaCut(0.01),
  76   fNExclLeadJets(1),
  77   fScaleFunction(0),
  78   fCreateHisto(histo),
  79   fTracks(0),
  80   fJets(0),
  81   fOutputList(0),
  82   fHistCentrality(0),
  83   fHistJetPt(0),
  84   fHistJetArea(0),
  85   fHistRhovsCent(0),
  86   fHistDeltaRhovsCent(0),
  87   fHistDeltaRhoScalevsCent(0),
  88   fHistJetPtvsCent(0),
  89   fHistJetAreavsCent(0),
  90   fHistNjetvsCent(0),
  91   fHistRhovsNtrack(0),
  92   fHistDeltaRhovsNtrack(0),
  93   fHistDeltaRhoScalevsNtrack(0),
  94   fHistJetPtvsNtrack(0),
  95   fHistJetAreavsNtrack(0),
  96   fHistNjetvsNtrack(0),
  97   fRhoScaled(0)
  98 {
  99   // Constructor
 100
 101   if (fCreateHisto)
 102     DefineOutput(1, TList::Class());
 103 }
 104
 105 //________________________________________________________________________
 106 void AliAnalysisTaskRho::UserCreateOutputObjects()
 107 {
 108   // User create output objects, called at the beginning of the analysis.
 109
 110   AliAnalysisTaskRhoBase::UserCreateOutputObjects();
 111
 112   if (fScaleFunction) {
 113     fRhoScaledName = fRhoName;
 114     fRhoScaledName += "_Scaled";
 115     fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
 116   }
 117
 118   if (!fCreateHisto)
 119     return;
 120
 121   OpenFile(1);
 122   fOutputList = new TList();
 123   fOutputList->SetOwner();
 124
 125   fHistCentrality             = new TH1F("Centrality",            "Centrality",            101, -1,  100);
 126   fHistRhovsCent              = new TH2F("RhovsCent",             "RhovsCent",             101, -1,  100,   500,  0,   500);
 127   fHistDeltaRhovsCent         = new TH2F("DeltaRhovsCent",        "DetlaRhovsCent",        101, -1,  100,   500, -250, 250);
 128   fHistDeltaRhoScalevsCent    = new TH2F("DeltaRhoScalevsCent",   "DeltaRhoScalevsCent",   101, -1,  100,   500, -250, 250);
 129   fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   200,  0,   500);
 130   fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   100,  0,   1.0);
 131   fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   100,  0,   100);
 132
 133   fHistRhovsNtrack            = new TH2F("RhovsNtrack",           "RhovsNtrack",           500,  0,  2500,  500,  0,   500);
 134   fHistDeltaRhovsNtrack       = new TH2F("DeltaRhovsNtrack",      "DeltaRhovsNtrack",      500,  0,  2500,  500, -250, 250);
 135   fHistDeltaRhoScalevsNtrack  = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 500,  0,  2500,  500, -250, 250);
 136   fHistJetPtvsNtrack          = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         500,  0,  2500,  200,  0,   500);
 137   fHistJetAreavsNtrack        = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       500,  0,  2500,  100,  0,   1.0);
 138   fHistNjetvsNtrack           = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         500,  0,  2500,  100,  0,   100);
 139
 140   fHistJetPt                  = new TH1F("JetPt",                  "Jet Pt",               100,   0,    250);
 141   fHistJetArea                = new TH1F("JetArea",                "Jet Area",             100, 0.0,    1.0);
 142
 143   fOutputList->Add(fHistCentrality);
 144   fOutputList->Add(fHistRhovsCent);
 145   fOutputList->Add(fHistJetPt);
 146   fOutputList->Add(fHistJetArea);
 147   fOutputList->Add(fHistDeltaRhovsCent);
 148   fOutputList->Add(fHistDeltaRhoScalevsCent);
 149   fOutputList->Add(fHistJetPtvsCent);
 150   fOutputList->Add(fHistJetAreavsCent);
 151   fOutputList->Add(fHistNjetvsCent);
 152
 153   fOutputList->Add(fHistRhovsNtrack);
 154   fOutputList->Add(fHistDeltaRhovsNtrack);
 155   fOutputList->Add(fHistDeltaRhoScalevsNtrack);
 156   fOutputList->Add(fHistJetPtvsNtrack);
 157   fOutputList->Add(fHistJetAreavsNtrack);
 158   fOutputList->Add(fHistNjetvsNtrack);
 159
 160   PostData(1, fOutputList);
 161 }
 162
 163 //________________________________________________________________________
 164 void AliAnalysisTaskRho::UserExec(Option_t *)
 165 {
 166   // Main loop, called for each event.
 167
 168   if (!fIsInit) {
 169     ExecOnce();
 170     fIsInit = 1;
 171   }
 172
 173   if (!fJets)
 174     return;
 175
 176   fRho->SetVal(-1);
 177   if (fRhoScaled)
 178     fRhoScaled->SetVal(-1);
 179
 180   DetermineCent();
 181
 182   const Int_t Njets   = fJets->GetEntries();
 183
 184   Int_t maxJetIds[]   = {-1, -1};
 185   Float_t maxJetPts[] = { 0,  0};
 186
 187   if (fNExclLeadJets > 0) {
 188     for (Int_t ij = 0; ij < Njets; ++ij) {
 189       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 190       if (!jet) {
 191         AliError(Form("%s: Could not receive jet %d", GetName(), ij));
 192         continue;
 193       }
 194       if (jet->Pt() > maxJetPts[0]) {
 195         maxJetPts[1] = maxJetPts[0];
 196         maxJetIds[1] = maxJetIds[0];
 197         maxJetPts[0] = jet->Pt();
 198         maxJetIds[0] = ij;
 199       } else if (jet->Pt() > maxJetPts[1]) {
 200         maxJetPts[1] = jet->Pt();
 201         maxJetIds[1] = ij;
 202       }
 203     }
 204     if (fNExclLeadJets < 2) {
 205       maxJetIds[1] = -1;
 206       maxJetPts[1] = -1;
 207     }
 208   }
 209
 210   static Double_t rhovec[999];
 211   Int_t NjetAcc = 0;
 212
 213   // push all jets within selected acceptance into stack
 214   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 215
 216     // exlcuding lead jets
 217     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
 218       continue;
 219
 220     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 221     if (!jet)
 222       continue;
 223
 224     // applying some other cuts
 225     if (jet->Area() < fAreaCut)
 226       continue;
 227     if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
 228       continue;
 229     if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
 230       continue;
 231
 232     rhovec[NjetAcc] = jet->Pt() / jet->Area();
 233
 234     ++NjetAcc;
 235
 236     if (fCreateHisto) {
 237       // filling histograms
 238       const Int_t Ntracks = fTracks->GetEntriesFast();
 239       fHistJetPt->Fill(jet->Pt());
 240       fHistJetArea->Fill(jet->Area());
 241       fHistJetPtvsCent->Fill(fCent, jet->Pt());
 242       fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
 243       fHistJetAreavsCent->Fill(fCent, jet->Area());
 244       fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
 245     }
 246   }
 247
 248   if (NjetAcc > 0) {
 249     //find median value
 250     Double_t rho0      = TMath::Median(NjetAcc, rhovec);
 251     fRho->SetVal(rho0);
 252
 253     Double_t rhoScaled = rho0;
 254     if (fRhoScaled) {
 255       rhoScaled *= GetScaleFactor(fCent);
 256       fRhoScaled->SetVal(rhoScaled);
 257     }
 258
 259     if (fCreateHisto) {
 260       // filling other histograms
 261       Double_t rho        = GetRhoFactor(fCent);
 262       const Int_t Ntracks = fTracks->GetEntriesFast();
 263       fHistRhovsCent->Fill(fCent, rho0);
 264       fHistDeltaRhovsCent->Fill(fCent, rho0 - rho);
 265       fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rho);
 266       fHistRhovsNtrack->Fill(Ntracks, rho0);
 267       fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rho);
 268       fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rho);
 269     }
 270   }
 271
 272   if (fCreateHisto) {
 273     const Int_t Ntracks = fTracks->GetEntriesFast();
 274     fHistCentrality->Fill(fCent);
 275     fHistNjetvsCent->Fill(fCent, NjetAcc);
 276     fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
 277     PostData(1, fOutputList);
 278   }
 279 }
 280
 281 //________________________________________________________________________
 282 void AliAnalysisTaskRho::ExecOnce()
 283 {
 284   // Initialize some settings that need to be determined in UserExec.
 285
 286   AliAnalysisTaskRhoBase::ExecOnce();
 287
 288   if (fRhoScaled) {
 289     // add rho to event if not yet there
 290     if (!(InputEvent()->FindListObject(fRhoScaledName))) {
 291       InputEvent()->AddObject(fRhoScaled);
 292     } else {
 293       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
 294       return;
 295     }
 296   }
 297
 298   if (fCreateHisto) {
 299     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 300     if (!fTracks) {
 301       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
 302       return;
 303     }
 304   }
 305
 306   fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
 307   if (!fJets) {
 308     AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
 309     return;
 310   }
 311 }
 312
 313 //________________________________________________________________________
 314 Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
 315 {
 316   // Get scale factor.
 317
 318   Double_t scale = 1;
 319   if (fScaleFunction)
 320     scale = fScaleFunction->Eval(cent);
 321   return scale;
 322 }