AliPhysics  2c6b7ad (2c6b7ad)
AliAnalysisTaskRhoAverage.cxx
Go to the documentation of this file.
1 //
2 // Calculation of rho, method: median all particle pt / multiplicity density.
3 //
4 // Authors: S. Aiola
5 
7 
8 #include <TClonesArray.h>
9 #include <TMath.h>
10 
11 #include "TLorentzVector.h"
12 #include "AliLog.h"
13 #include "AliRhoParameter.h"
14 #include "AliVCluster.h"
15 #include "AliVTrack.h"
16 #include "AliClusterContainer.h"
17 #include "AliParticleContainer.h"
18 
20 
21 //________________________________________________________________________
24  fRhoType(0),
25  fNExclLeadPart(0),
26  fUseMedian(kFALSE),
27  fTotalArea(1)
28 {
29  // Default constructor.
30 }
31 
32 //________________________________________________________________________
34  AliAnalysisTaskRhoBase(name, histo),
35  fRhoType(0),
36  fNExclLeadPart(0),
37  fUseMedian(kFALSE),
38  fTotalArea(1)
39 {
40  // Constructor.
41 }
42 
43 //________________________________________________________________________
45 {
46  // Run the analysis.
47 
48  const Int_t NMAX = 9999;
49  static Double_t rhovec[NMAX];
50  Int_t NpartAcc = 0;
51 
52  Int_t maxPartIds[] = {0, 0};
53  Float_t maxPartPts[] = {0, 0};
54 
55  // push all jets within selected acceptance into stack
56 
59 
60  if (fNExclLeadPart > 0) {
61 
62  if (tracks && (fRhoType == 0 || fRhoType == 1)) {
63 
64  AliVParticle *track = 0;
65  tracks->ResetCurrentID();
66  while ((track = tracks->GetNextAcceptParticle())) {
67 
68  if (track->Pt() > maxPartPts[0]) {
69  maxPartPts[1] = maxPartPts[0];
70  maxPartIds[1] = maxPartIds[0];
71  maxPartPts[0] = track->Pt();
72  maxPartIds[0] = tracks->GetCurrentID()+1;
73  }
74  else if (track->Pt() > maxPartPts[1]) {
75  maxPartPts[1] = track->Pt();
76  maxPartIds[1] = tracks->GetCurrentID()+1;
77  }
78  }
79  }
80 
81  if (clusters && (fRhoType == 0 || fRhoType == 2)) {
82 
83  AliVCluster *cluster = 0;
84  clusters->ResetCurrentID();
85  while ((cluster = clusters->GetNextAcceptCluster())) {
86  TLorentzVector nPart;
87  clusters->GetMomentum(nPart, clusters->GetCurrentID());
88 
89  if (nPart.Pt() > maxPartPts[0]) {
90  maxPartPts[1] = maxPartPts[0];
91  maxPartIds[1] = maxPartIds[0];
92  maxPartPts[0] = nPart.Pt();
93  maxPartIds[0] = -clusters->GetCurrentID()-1;
94  }
95  else if (nPart.Pt() > maxPartPts[1]) {
96  maxPartPts[1] = nPart.Pt();
97  maxPartIds[1] = -clusters->GetCurrentID()-1;
98  }
99  }
100  }
101 
102  if (fNExclLeadPart < 2) {
103  maxPartIds[1] = 0;
104  maxPartPts[1] = 0;
105  }
106  }
107 
108  if (tracks && (fRhoType == 0 || fRhoType == 1)) {
109  AliVParticle *track = 0;
110  tracks->ResetCurrentID();
111  while ((track = tracks->GetNextAcceptParticle()) && NpartAcc < NMAX) {
112 
113  // exlcuding lead particles
114  if (tracks->GetCurrentID() == maxPartIds[0]-1 || tracks->GetCurrentID() == maxPartIds[1]-1)
115  continue;
116 
117  rhovec[NpartAcc] = track->Pt();
118  ++NpartAcc;
119  }
120  }
121 
122  if (clusters && (fRhoType == 0 || fRhoType == 2)) {
123 
124  AliVCluster *cluster = 0;
125  clusters->ResetCurrentID();
126  while ((cluster = clusters->GetNextAcceptCluster()) && NpartAcc < NMAX) {
127  // exlcuding lead particles
128  if (clusters->GetCurrentID() == -maxPartIds[0]-1 || clusters->GetCurrentID() == -maxPartIds[1]-1)
129  continue;
130 
131  TLorentzVector nPart;
132  clusters->GetMomentum(nPart, clusters->GetCurrentID());
133 
134  rhovec[NpartAcc] = nPart.Pt();
135  ++NpartAcc;
136  }
137  }
138 
139  if (NpartAcc == NMAX) {
140  AliError(Form("%s: NpartAcc >= %d", GetName(), NMAX));
141  }
142 
143  Double_t rho = 0;
144 
145  if (NpartAcc > 0) {
146  if (fUseMedian)
147  rho = TMath::Median(NpartAcc, rhovec);
148  else
149  rho = TMath::Mean(NpartAcc, rhovec);
150 
151  rho *= NpartAcc / fTotalArea;
152  }
153 
154  fOutRho->SetVal(rho);
155 
156  if (fScaleFunction) {
157  Double_t rhoScaled = rho * GetScaleFactor(fCent);
158  fOutRhoScaled->SetVal(rhoScaled);
159  }
160 
161  return kTRUE;
162 }
163 
164 //________________________________________________________________________
166 {
168 
170  if (!partCont) {
171  AliError(Form("%s: No particle container found! Assuming area = 1...",GetName()));
172  fTotalArea = 1;
173  return;
174  }
175 
176  Float_t maxEta = partCont->GetParticleEtaMax();
177  Float_t minEta = partCont->GetParticleEtaMin();
178  Float_t maxPhi = partCont->GetParticlePhiMax();
179  Float_t minPhi = partCont->GetParticlePhiMin();
180 
181  if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
182  if (minPhi < 0) minPhi = 0;
183 
184  fTotalArea = (maxEta - minEta) * (maxPhi - minPhi);
185 
186  if (fTotalArea < 1e-6) {
187  AliError(Form("%s: Area = %f < 1e-6, assuming area = 1", GetName(), fTotalArea));
188  fTotalArea = 1;
189  }
190 }
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t GetParticleEtaMin() const
Container for particles within the EMCAL framework.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliRhoParameter * fOutRhoScaled
output rho object
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Double_t fCent
!event centrality
virtual Double_t GetScaleFactor(Double_t cent)
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Double_t GetParticlePhiMax() const
Double_t GetParticleEtaMax() const
Double_t GetParticlePhiMin() const
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
bool Bool_t
Definition: External.C:53
void ExecOnce()
Perform steps needed to initialize the analysis.
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()