3 // Calculation of rho, method: median all particle pt / multiplicity density.
7 #include "AliAnalysisTaskRhoAverage.h"
9 #include <TClonesArray.h>
13 #include "AliRhoParameter.h"
14 #include "AliVCluster.h"
15 #include "AliVTrack.h"
16 #include "AliParticleContainer.h"
18 ClassImp(AliAnalysisTaskRhoAverage)
20 //________________________________________________________________________
21 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage() :
22 AliAnalysisTaskRhoBase("AliAnalysisTaskRhoAverage"),
28 // Default constructor.
31 //________________________________________________________________________
32 AliAnalysisTaskRhoAverage::AliAnalysisTaskRhoAverage(const char *name, Bool_t histo) :
33 AliAnalysisTaskRhoBase(name, histo),
42 //________________________________________________________________________
43 Bool_t AliAnalysisTaskRhoAverage::Run()
47 const Int_t NMAX = 9999;
48 static Double_t rhovec[NMAX];
51 Int_t maxPartIds[] = {0, 0};
52 Float_t maxPartPts[] = {0, 0};
54 // push all jets within selected acceptance into stack
56 if (fNExclLeadPart > 0) {
58 if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
60 const Int_t Ntracks = fTracks->GetEntriesFast();
62 for (Int_t it = 0; it < Ntracks; ++it) {
64 AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
67 AliError(Form("%s: Could not receive track %d", GetName(), it));
71 if (!AcceptTrack(track))
74 if (track->Pt() > maxPartPts[0]) {
75 maxPartPts[1] = maxPartPts[0];
76 maxPartIds[1] = maxPartIds[0];
77 maxPartPts[0] = track->Pt();
80 else if (track->Pt() > maxPartPts[1]) {
81 maxPartPts[1] = track->Pt();
87 if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
89 const Int_t Nclusters = fCaloClusters->GetEntriesFast();
91 for (Int_t ic = 0; ic < Nclusters; ++ic) {
93 AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
96 AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
100 if (!AcceptCluster(cluster))
103 TLorentzVector nPart;
104 cluster->GetMomentum(nPart, fVertex);
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;
112 else if (nPart.Pt() > maxPartPts[1]) {
113 maxPartPts[1] = nPart.Pt();
114 maxPartIds[1] = -ic-1;
119 if (fNExclLeadPart < 2) {
125 if (fTracks && (fRhoType == 0 || fRhoType == 1)) {
127 const Int_t Ntracks = fTracks->GetEntriesFast();
129 for (Int_t it = 0; it < Ntracks && NpartAcc < NMAX; ++it) {
131 // exlcuding lead particles
132 if (it == maxPartIds[0]-1 || it == maxPartIds[1]-1)
135 AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
138 AliError(Form("%s: Could not receive track %d", GetName(), it));
142 if (!AcceptTrack(track))
145 rhovec[NpartAcc] = track->Pt();
150 if (fCaloClusters && (fRhoType == 0 || fRhoType == 2)) {
152 const Int_t Nclusters = fCaloClusters->GetEntriesFast();
154 for (Int_t ic = 0; ic < Nclusters && NpartAcc < NMAX; ++ic) {
156 // exlcuding lead particles
157 if (ic == -maxPartIds[0]-1 || ic == -maxPartIds[1]-1)
160 AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
163 AliError(Form("%s: Could not receive cluster %d", GetName(), ic));
167 if (!AcceptCluster(cluster))
170 TLorentzVector nPart;
171 cluster->GetMomentum(nPart, fVertex);
173 rhovec[NpartAcc] = nPart.Pt();
178 if (NpartAcc == NMAX) {
179 AliError(Form("%s: NpartAcc >= %d", GetName(), NMAX));
186 rho = TMath::Median(NpartAcc, rhovec);
188 rho = TMath::Mean(NpartAcc, rhovec);
190 rho *= NpartAcc / fTotalArea;
193 fOutRho->SetVal(rho);
195 if (fScaleFunction) {
196 Double_t rhoScaled = rho * GetScaleFactor(fCent);
197 fOutRhoScaled->SetVal(rhoScaled);
203 //________________________________________________________________________
204 void AliAnalysisTaskRhoAverage::ExecOnce()
206 AliAnalysisTaskRhoBase::ExecOnce();
208 AliParticleContainer *partCont = GetParticleContainer(0);
210 AliError(Form("%s: No particle container found! Assuming area = 1...",GetName()));
215 Float_t maxEta = partCont->GetParticleEtaMax();
216 Float_t minEta = partCont->GetParticleEtaMin();
217 Float_t maxPhi = partCont->GetParticlePhiMax();
218 Float_t minPhi = partCont->GetParticlePhiMin();
220 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
221 if (minPhi < 0) minPhi = 0;
223 fTotalArea = (maxEta - minEta) * (maxPhi - minPhi);
225 if (fTotalArea < 1e-6) {
226 AliError(Form("%s: Area = %f < 1e-6, assuming area = 1", GetName(), fTotalArea));