PWGJE/EMCALJetTasks/AliAnalysisTaskRhoAverage.cxx

   1 // $Id$
   2 //
   3 // Calculation of rho, method: median all particle pt / multiplicity density.
   4 //
   5 // Authors: S. Aiola
   6
   7 #include "AliAnalysisTaskRhoAverage.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TMath.h>
  11
  12 #include "AliLog.h"
  13 #include "AliRhoParameter.h"
  14 #include "AliVCluster.h"
  15 #include "AliVTrack.h"
  16 #include "AliParticleContainer.h"
  17
  18 ClassImp(AliAnalysisTaskRhoAverage)
  19
  20 //________________________________________________________________________
  21 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage() :
  22   AliAnalysisTaskRhoBase("AliAnalysisTaskRhoAverage"),
  23   fRhoType(0),
  24   fNExclLeadPart(0),
  25   fUseMedian(kFALSE),
  26   fTotalArea(1)
  27 {
  28   // Default constructor.
  29 }
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage(const char *name, Bool_t histo) :
  33   AliAnalysisTaskRhoBase(name, histo),
  34   fRhoType(0),
  35   fNExclLeadPart(0),
  36   fUseMedian(kFALSE),
  37   fTotalArea(1)
  38 {
  39   // Constructor.
  40 }
  41
  42 //________________________________________________________________________
  43 Bool_t AliAnalysisTaskRhoAverage::Run()
  44 {
  45   // Run the analysis.
  46
  47   const Int_t NMAX = 9999;
  48   static Double_t rhovec[NMAX];
  49   Int_t NpartAcc = 0;
  50
  51   Int_t   maxPartIds[] = {0, 0};
  52   Float_t maxPartPts[] = {0, 0};
  53
  54   // push all jets within selected acceptance into stack
  55
  56   if (fNExclLeadPart > 0) {
  57
  58     if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
  59
  60       const Int_t Ntracks = fTracks->GetEntriesFast();
  61
  62       for (Int_t it = 0; it < Ntracks; ++it) {
  63
  64         AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
  65
  66         if (!track) {
  67           AliError(Form("%s: Could not receive track %d", GetName(), it));
  68           continue;
  69         }
  70
  71         if (!AcceptTrack(track))
  72           continue;
  73
  74         if (track->Pt() > maxPartPts[0]) {
  75           maxPartPts[1] = maxPartPts[0];
  76           maxPartIds[1] = maxPartIds[0];
  77           maxPartPts[0] = track->Pt();
  78           maxPartIds[0] = it+1;
  79         }
  80         else if (track->Pt() > maxPartPts[1]) {
  81           maxPartPts[1] = track->Pt();
  82           maxPartIds[1] = it+1;
  83         }
  84       }
  85     }
  86
  87     if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
  88
  89       const Int_t Nclusters = fCaloClusters->GetEntriesFast();
  90
  91       for (Int_t ic = 0; ic < Nclusters; ++ic) {
  92
  93         AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
  94
  95         if (!cluster) {
  96           AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
  97           continue;
  98         }
  99
 100         if (!AcceptCluster(cluster))
 101           continue;
 102
 103         TLorentzVector nPart;
 104         cluster->GetMomentum(nPart, fVertex);
 105
 106         if (nPart.Pt() > maxPartPts[0]) {
 107           maxPartPts[1] = maxPartPts[0];
 108           maxPartIds[1] = maxPartIds[0];
 109           maxPartPts[0] = nPart.Pt();
 110           maxPartIds[0] = -ic-1;
 111         }
 112         else if (nPart.Pt() > maxPartPts[1]) {
 113           maxPartPts[1] = nPart.Pt();
 114           maxPartIds[1] = -ic-1;
 115         }
 116       }
 117     }
 118
 119     if (fNExclLeadPart < 2) {
 120       maxPartIds[1] = 0;
 121       maxPartPts[1] = 0;
 122     }
 123   }
 124
 125   if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
 126
 127     const Int_t Ntracks = fTracks->GetEntriesFast();
 128
 129     for (Int_t it = 0; it < Ntracks && NpartAcc < NMAX; ++it) {
 130
 131       // exlcuding lead particles
 132       if (it == maxPartIds[0]-1 || it == maxPartIds[1]-1)
 133         continue;
 134
 135       AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
 136
 137       if (!track) {
 138         AliError(Form("%s: Could not receive track %d", GetName(), it));
 139         continue;
 140       }
 141
 142       if (!AcceptTrack(track))
 143         continue;
 144
 145       rhovec[NpartAcc] = track->Pt();
 146       ++NpartAcc;
 147     }
 148   }
 149
 150   if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
 151
 152     const Int_t Nclusters = fCaloClusters->GetEntriesFast();
 153
 154     for (Int_t ic = 0; ic < Nclusters && NpartAcc < NMAX; ++ic) {
 155
 156       // exlcuding lead particles
 157       if (ic == -maxPartIds[0]-1 || ic == -maxPartIds[1]-1)
 158         continue;
 159
 160       AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
 161
 162       if (!cluster) {
 163         AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
 164         continue;
 165       }
 166
 167       if (!AcceptCluster(cluster))
 168         continue;
 169
 170       TLorentzVector nPart;
 171       cluster->GetMomentum(nPart, fVertex);
 172
 173       rhovec[NpartAcc] = nPart.Pt();
 174       ++NpartAcc;
 175     }
 176   }
 177
 178   if (NpartAcc == NMAX) {
 179     AliError(Form("%s: NpartAcc >= %d", GetName(), NMAX));
 180   }
 181
 182   Double_t rho = 0;
 183
 184   if (NpartAcc > 0) {
 185     if (fUseMedian)
 186       rho = TMath::Median(NpartAcc, rhovec);
 187     else
 188       rho = TMath::Mean(NpartAcc, rhovec);
 189
 190     rho *= NpartAcc / fTotalArea;
 191   }
 192
 193   fOutRho->SetVal(rho);
 194
 195   if (fScaleFunction) {
 196     Double_t rhoScaled = rho * GetScaleFactor(fCent);
 197     fOutRhoScaled->SetVal(rhoScaled);
 198   }
 199
 200   return kTRUE;
 201 }
 202
 203 //________________________________________________________________________
 204 void AliAnalysisTaskRhoAverage::ExecOnce()
 205 {
 206   AliAnalysisTaskRhoBase::ExecOnce();
 207
 208   AliParticleContainer *partCont = GetParticleContainer(0);
 209   if (!partCont) {
 210     AliError(Form("%s: No particle container found! Assuming area = 1...",GetName()));
 211     fTotalArea = 1;
 212     return;
 213   }
 214
 215   Float_t maxEta = partCont->GetParticleEtaMax();
 216   Float_t minEta = partCont->GetParticleEtaMin();
 217   Float_t maxPhi = partCont->GetParticlePhiMax();
 218   Float_t minPhi = partCont->GetParticlePhiMin();
 219
 220   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 221   if (minPhi < 0) minPhi = 0;
 222
 223   fTotalArea = (maxEta - minEta) * (maxPhi - minPhi);
 224
 225   if (fTotalArea < 1e-6) {
 226     AliError(Form("%s: Area = %f < 1e-6, assuming area = 1", GetName(), fTotalArea));
 227     fTotalArea = 1;
 228   }
 229 }