PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx

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