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