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