PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx

   1 // $Id: $
   2 //
   3 // Calculation of rho
   4 //
   5 // Authors: R.Reed, S.Aiola
   6
   7 #include <TList.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TClonesArray.h>
  11 #include <TF1.h>
  12 #include <TMath.h>
  13
  14 #include "AliLog.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliVEventHandler.h"
  17 #include "AliCentrality.h"
  18 #include "AliEmcalJet.h"
  19 #include "AliVCluster.h"
  20
  21 #include "AliAnalysisTaskRho.h"
  22
  23 ClassImp(AliAnalysisTaskRho)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskRho::AliAnalysisTaskRho() :
  27   AliAnalysisTaskRhoBase(),
  28   fTracksName("tracks"),
  29   fJetsName("KtJets"),
  30   fClustersName("caloClusters"),
  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 AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name) :
  62   AliAnalysisTaskRhoBase(name),
  63   fTracksName("tracks"),
  64   fJetsName("KtJets"),
  65   fClustersName("caloClusters"),
  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   DefineOutput(1, TList::Class());
  96 }
  97
  98 //________________________________________________________________________
  99 void AliAnalysisTaskRho::UserCreateOutputObjects()
 100 {
 101   AliAnalysisTaskRhoBase::UserCreateOutputObjects();
 102
 103   fRhoScaledName = fRhoName;
 104   fRhoScaledName += "_Scaled";
 105   fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);
 106
 107   OpenFile(1);
 108   fOutputList = new TList();
 109   fOutputList->SetOwner();
 110
 111   if (!fCreateHisto) {
 112     PostData(1, fOutputList);
 113     return;
 114   }
 115
 116   fHistCentrality             = new TH1F("Centrality",            "Centrality",            101, -1,  100);
 117   fHistRhovsCent              = new TH2F("RhovsCent",             "RhovsCent",             101, -1,  100,   500,  0,   500);
 118   fHistDeltaRhovsCent         = new TH2F("DeltaRhovsCent",        "DetlaRhovsCent",        101, -1,  100,   500, -250, 250);
 119   fHistDeltaRhoScalevsCent    = new TH2F("DeltaRhoScalevsCent",   "DeltaRhoScalevsCent",   101, -1,  100,   500, -250, 250);
 120   fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   200,  0,   500);
 121   fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   100,  0,   1.0);
 122   fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   100,  0,   100);
 123
 124   fHistRhovsNtrack            = new TH2F("RhovsNtrack",           "RhovsNtrack",           500,  0,  2500,  500,  0,   500);
 125   fHistDeltaRhovsNtrack       = new TH2F("DeltaRhovsNtrack",      "DeltaRhovsNtrack",      500,  0,  2500,  500, -250, 250);
 126   fHistDeltaRhoScalevsNtrack  = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 500,  0,  2500,  500, -250, 250);
 127   fHistJetPtvsNtrack          = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         500,  0,  2500,  200,  0,   500);
 128   fHistJetAreavsNtrack        = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       500,  0,  2500,  100,  0,   1.0);
 129   fHistNjetvsNtrack           = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         500,  0,  2500,  100,  0,   100);
 130
 131   fHistJetPt                  = new TH1F("JetPt",                  "Jet Pt",               100,   0,    250);
 132   fHistJetArea                = new TH1F("JetArea",                "Jet Area",             100, 0.0,    1.0);
 133
 134   fOutputList->Add(fHistCentrality);
 135   fOutputList->Add(fHistRhovsCent);
 136   fOutputList->Add(fHistJetPt);
 137   fOutputList->Add(fHistJetArea);
 138   fOutputList->Add(fHistDeltaRhovsCent);
 139   fOutputList->Add(fHistDeltaRhoScalevsCent);
 140   fOutputList->Add(fHistJetPtvsCent);
 141   fOutputList->Add(fHistJetAreavsCent);
 142   fOutputList->Add(fHistNjetvsCent);
 143
 144   fOutputList->Add(fHistRhovsNtrack);
 145   fOutputList->Add(fHistDeltaRhovsNtrack);
 146   fOutputList->Add(fHistDeltaRhoScalevsNtrack);
 147   fOutputList->Add(fHistJetPtvsNtrack);
 148   fOutputList->Add(fHistJetAreavsNtrack);
 149   fOutputList->Add(fHistNjetvsNtrack);
 150
 151   PostData(1, fOutputList);
 152
 153 }
 154
 155 //________________________________________________________________________
 156 Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
 157 {
 158   Double_t scale = 1;
 159   if (fScaleFunction)
 160     scale = fScaleFunction->Eval(cent);
 161   return scale;
 162 }
 163
 164
 165 //________________________________________________________________________
 166 void AliAnalysisTaskRho::UserExec(Option_t *)
 167 {
 168   // Main loop, called for each event.
 169
 170   AliAnalysisTaskRhoBase::UserExec("");
 171
 172   fRho->SetVal(-1);
 173
 174    // add rho to event if not yet there
 175   if (!(InputEvent()->FindListObject(fRhoScaledName))) {
 176     new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, -1);
 177     InputEvent()->AddObject(fRhoScaled);
 178   }
 179   else {
 180     fRhoScaled->SetVal(-1);
 181   }
 182
 183   // optimization in case autobranch loading is off
 184   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 185   if (fTracksName == "Tracks")
 186     am->LoadBranch("Tracks");
 187
 188   TClonesArray *jets = 0;
 189   TClonesArray *tracks = 0;
 190
 191   tracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 192   if (!tracks) {
 193   AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 194    return;
 195   }
 196
 197   jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
 198   if (!jets) {
 199     AliError(Form("Pointer to jets %s == 0", fJetsName.Data() ));
 200     return;
 201   }
 202
 203   if (fCreateHisto)
 204     fHistCentrality->Fill(fCent);
 205
 206   const Int_t Ntracks = tracks->GetEntries();
 207   const Int_t Njets = jets->GetEntries();
 208
 209   Int_t maxJetIds[] = {-1, -1};
 210   Int_t maxJetPts[] = {0, 0};
 211   if (fNExclLeadJets > 0) {
 212     for (Int_t ij = 0; ij < Njets; ij++) {
 213
 214       AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ij));
 215
 216       if (!jet) {
 217         AliError(Form("Could not receive jet %d", ij));
 218         continue;
 219       }
 220
 221       if (jet->Pt() > maxJetPts[0]) {
 222         maxJetPts[1] = maxJetPts[0];
 223         maxJetIds[1] = maxJetPts[0];
 224         maxJetPts[0] = jet->Pt();
 225         maxJetIds[0] = ij;
 226       }
 227
 228       if (jet->Pt() > maxJetPts[1]) {
 229         maxJetPts[1] = jet->Pt();
 230         maxJetIds[1] = ij;
 231       }
 232     }
 233     if (fNExclLeadJets < 2) {
 234       maxJetIds[1] = -1;
 235       maxJetPts[1] = -1;
 236     }
 237   }
 238
 239   static Double_t rhovec[999];
 240   Int_t NjetAcc = 0;
 241
 242   // push all jets within selected acceptance into stack
 243   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 244
 245     // exlcuding lead jets
 246     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
 247       continue;
 248
 249     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
 250
 251     if (!jet)
 252       continue;
 253
 254     // applying some other cuts
 255     if (jet->Area() < fAreaCut)
 256       continue;
 257     if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
 258       continue;
 259     if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
 260       continue;
 261     if (jet->Area() == 0)
 262       continue;
 263
 264     rhovec[NjetAcc] = jet->Pt() / jet->Area();
 265
 266     NjetAcc++;
 267
 268     if (fCreateHisto) {
 269       // filling histograms
 270       fHistJetPt->Fill(jet->Pt());
 271       fHistJetArea->Fill(jet->Area());
 272       fHistJetPtvsCent->Fill(fCent, jet->Pt());
 273       fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
 274       fHistJetAreavsCent->Fill(fCent, jet->Area());
 275       fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
 276     }
 277   }
 278
 279   if (fCreateHisto) {
 280     fHistNjetvsCent->Fill(fCent, NjetAcc);
 281     fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
 282   }
 283
 284   Double_t scale = GetScaleFactor(fCent);
 285   Double_t rhochem = GetRhoFactor(fCent);
 286
 287   Double_t rho0 = -1;
 288
 289   if (NjetAcc > 0){
 290     //find median value
 291     rho0 = TMath::Median(NjetAcc, rhovec);
 292
 293     Double_t rhoScaled = rho0 * scale;
 294
 295     fRho->SetVal(rho0);
 296     fRhoScaled->SetVal(rhoScaled);
 297
 298     if (fCreateHisto) {
 299       // filling other histograms
 300       fHistRhovsCent->Fill(fCent, rho0);
 301       fHistDeltaRhovsCent->Fill(fCent, rho0 - rhochem);
 302       fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rhochem);
 303       fHistRhovsNtrack->Fill(Ntracks, rho0);
 304       fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rhochem);
 305       fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rhochem);
 306     }
 307   }
 308
 309   if (fCreateHisto)
 310     PostData(1, fOutputList);
 311 }
 312
 313
 314 //________________________________________________________________________
 315 void AliAnalysisTaskRho::Terminate(Option_t *)
 316 {
 317
 318 }