PWGJE/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 <TMath.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliLog.h"
  17 #include "AliRhoParameter.h"
  18
  19 ClassImp(AliAnalysisTaskRho)
  20
  21 //________________________________________________________________________
  22 AliAnalysisTaskRho::AliAnalysisTaskRho() :
  23   AliAnalysisTaskRhoBase("AliAnalysisTaskRho"),
  24   fNExclLeadJets(0)
  25 {
  26   // Constructor.
  27 }
  28
  29 //________________________________________________________________________
  30 AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
  31   AliAnalysisTaskRhoBase(name, histo),
  32   fNExclLeadJets(0)
  33 {
  34   // Constructor.
  35 }
  36
  37
  38 //________________________________________________________________________
  39 Bool_t AliAnalysisTaskRho::Run()
  40 {
  41   // Run the analysis.
  42
  43   fRho->SetVal(0);
  44   if (fRhoScaled)
  45     fRhoScaled->SetVal(0);
  46
  47   if (!fJets)
  48     return kFALSE;
  49
  50   const Int_t Njets   = fJets->GetEntries();
  51
  52   Int_t maxJetIds[]   = {-1, -1};
  53   Float_t maxJetPts[] = { 0,  0};
  54
  55   if (fNExclLeadJets > 0) {
  56     for (Int_t ij = 0; ij < Njets; ++ij) {
  57       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
  58       if (!jet) {
  59         AliError(Form("%s: Could not receive jet %d", GetName(), ij));
  60         continue;
  61       }
  62
  63       if (!AcceptJet(jet))
  64         continue;
  65
  66       if (jet->Pt() > maxJetPts[0]) {
  67         maxJetPts[1] = maxJetPts[0];
  68         maxJetIds[1] = maxJetIds[0];
  69         maxJetPts[0] = jet->Pt();
  70         maxJetIds[0] = ij;
  71       } else if (jet->Pt() > maxJetPts[1]) {
  72         maxJetPts[1] = jet->Pt();
  73         maxJetIds[1] = ij;
  74       }
  75     }
  76     if (fNExclLeadJets < 2) {
  77       maxJetIds[1] = -1;
  78       maxJetPts[1] = 0;
  79     }
  80   }
  81
  82   static Double_t rhovec[999];
  83   Int_t NjetAcc = 0;
  84
  85   // push all jets within selected acceptance into stack
  86   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
  87
  88     // exlcuding lead jets
  89     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
  90       continue;
  91
  92     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
  93     if (!jet) {
  94       AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
  95       continue;
  96     }
  97
  98     if (!AcceptJet(jet))
  99       continue;
 100
 101     rhovec[NjetAcc] = jet->Pt() / jet->Area();
 102     ++NjetAcc;
 103   }
 104
 105
 106   if (NjetAcc > 0) {
 107     //find median value
 108     Double_t rho = TMath::Median(NjetAcc, rhovec);
 109     fRho->SetVal(rho);
 110
 111     if (fRhoScaled) {
 112       Double_t rhoScaled = rho * GetScaleFactor(fCent);
 113       fRhoScaled->SetVal(rhoScaled);
 114     }
 115   }
 116
 117   return kTRUE;
 118 }