PWGJE/EMCALJetTasks/AliAnalysisTaskRhoAverage.cxx

   1 // $Id$
   2 //
   3 // Calculation of rho, method: sum of all particle pt / full acceptance area.
   4 //
   5 // Authors: S. Aiola
   6
   7 #include "AliAnalysisTaskRhoAverage.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TList.h>
  11 #include <TMath.h>
  12
  13 #include "AliAnalysisManager.h"
  14 #include "AliCentrality.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliLog.h"
  17 #include "AliRhoParameter.h"
  18 #include "AliVCluster.h"
  19 #include "AliVEventHandler.h"
  20 #include "AliVTrack.h"
  21
  22 ClassImp(AliAnalysisTaskRhoAverage)
  23
  24 //________________________________________________________________________
  25 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage() :
  26   AliAnalysisTaskRhoBase(),
  27   fTracksName(),
  28   fClustersName(),
  29   fJetsName(),
  30   fEtaMin(0),
  31   fEtaMax(0),
  32   fPhiMin(0),
  33   fPhiMax(0),
  34   fPtMin(0),
  35   fClusters(0),
  36   fJets(0),
  37   fTracks(0)
  38 {
  39   // Default constructor.
  40 }
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage(const char *name) :
  44   AliAnalysisTaskRhoBase(name),
  45   fTracksName("tracks"),
  46   fClustersName("caloClusters"),
  47   fJetsName("KtJets"),
  48   fEtaMin(-0.9),
  49   fEtaMax(0.9),
  50   fPhiMin(0),
  51   fPhiMax(2 * TMath::Pi()),
  52   fPtMin(0.15),
  53   fClusters(0),
  54   fJets(0),
  55   fTracks(0)
  56 {
  57   // Constructor.
  58 }
  59
  60 //________________________________________________________________________
  61 void AliAnalysisTaskRhoAverage::UserExec(Option_t *)
  62 {
  63   // Main loop, called for each event.
  64
  65   if (!fIsInit) {
  66     ExecOnce();
  67     fIsInit = 1;
  68   }
  69
  70   fRho->SetVal(-1);
  71
  72   Double_t rho = 0;
  73
  74   Int_t Ntracks = 0;
  75   if (fTracks)
  76     Ntracks = fTracks->GetEntriesFast();
  77
  78   Int_t Nclusters = 0;
  79   if (fClusters)
  80     Nclusters = fClusters->GetEntriesFast();
  81
  82   Int_t Njets = 0;
  83   if (fJets)
  84     Njets = fJets->GetEntriesFast();
  85
  86   Double_t maxJetPt = 0;
  87   Int_t maxJetId = -1;
  88   AliEmcalJet *maxJet = 0;
  89   for (Int_t ij = 0; ij < Njets; ij++) {
  90
  91     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
  92
  93     if (!jet) {
  94       AliError(Form("%s: Could not receive jet %d", GetName(), ij));
  95       continue;
  96     }
  97
  98     if (jet->Pt() > maxJetPt) {
  99       maxJetPt = jet->Pt();
 100       maxJetId = ij;
 101     }
 102   }
 103
 104   if (maxJetId >= 0)
 105     maxJet = static_cast<AliEmcalJet*>(fJets->At(maxJetId));
 106
 107   for (Int_t it = 0; it < Ntracks; ++it) {
 108
 109     AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
 110
 111     if (!track) {
 112       AliError(Form("%s: Could not receive track %d", GetName(), it));
 113       continue;
 114     }
 115
 116     if (track->Eta() < fEtaMin || track->Eta() > fEtaMax || track->Phi() < fPhiMin || track->Phi() > fPhiMax)
 117       continue;
 118
 119     if (track->Pt() < fPtMin)
 120       continue;
 121
 122     if (maxJet && IsJetTrack(maxJet, it))
 123       continue;
 124
 125     rho += track->Pt();
 126   }
 127
 128   Double_t vertex[] = {0, 0, 0};
 129   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 130
 131   for (Int_t ic = 0; ic < Nclusters; ++ic) {
 132
 133     AliVCluster *cluster = static_cast<AliVCluster*>(fClusters->At(ic));
 134
 135     if (!cluster) {
 136       AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
 137       continue;
 138     }
 139
 140     Float_t pos[3];
 141     cluster->GetPosition(pos);
 142     TVector3 clusVec(pos);
 143
 144     if (clusVec.Eta() < fEtaMin || clusVec.Eta() > fEtaMax ||
 145         clusVec.Phi() < fPhiMin || clusVec.Phi() > fPhiMax)
 146       continue;
 147
 148     TLorentzVector nPart;
 149     cluster->GetMomentum(nPart, const_cast<Double_t*>(vertex));
 150
 151     if (nPart.Et() < fPtMin)
 152       continue;
 153
 154     if (maxJet && IsJetCluster(maxJet, ic))
 155       continue;
 156
 157     rho += nPart.Et();
 158   }
 159
 160   Double_t area = (fEtaMax - fEtaMin) * (fPhiMax - fPhiMin);
 161
 162   if (maxJet)
 163     area -= maxJet->Area();
 164
 165   if (area>0) {
 166     rho /= area;
 167     fRho->SetVal(rho);
 168   } else {
 169     AliError(Form("%s: Area negative %f", GetName(), area));
 170   }
 171 }
 172
 173 //________________________________________________________________________
 174 void AliAnalysisTaskRhoAverage::ExecOnce()
 175 {
 176   // Initialize some settings that need to be determined in UserExec.
 177
 178   AliAnalysisTaskRhoBase::ExecOnce();
 179
 180   if (!fClustersName.IsNull()) {
 181     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fClustersName));
 182     if (!fClusters) {
 183       AliError(Form("%s: Pointer to jets %s == 0", GetName(), fClustersName.Data() ));
 184       return;
 185     }
 186   }
 187
 188   if (!fTracksName.IsNull()) {
 189     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 190     if (!fTracks) {
 191       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
 192       return;
 193     }
 194   }
 195
 196   if (!fJetsName.IsNull()) {
 197     fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
 198     if (!fJets) {
 199       AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
 200       return;
 201     }
 202   }
 203 }
 204
 205 //________________________________________________________________________
 206 Bool_t AliAnalysisTaskRhoAverage::IsJetTrack(AliEmcalJet* jet, Int_t itrack) const
 207 {
 208   // Return true if track is in jet.
 209
 210   for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
 211     Int_t ijettrack = jet->TrackAt(i);
 212     if (ijettrack == itrack)
 213       return kTRUE;
 214   }
 215   return kFALSE;
 216 }
 217
 218 //________________________________________________________________________
 219 Bool_t AliAnalysisTaskRhoAverage::IsJetCluster(AliEmcalJet* jet, Int_t iclus) const
 220 {
 221   // Return true if cluster is in jet.
 222
 223   for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
 224     Int_t ijetclus = jet->ClusterAt(i);
 225     if (ijetclus == iclus)
 226       return kTRUE;
 227   }
 228   return kFALSE;
 229 }