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