AliPhysics  a34469b (a34469b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskRhoMass.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Calculation of rho mass from a collection of jets.
4 // If scale function is given the scaled rho will be exported
5 // with the name as "fOutRhoMassName".Apppend("_Scaled").
6 //
7 // Authors: M. Verweij
8 
10 
11 #include <TClonesArray.h>
12 #include <TMath.h>
13 
14 #include "AliAnalysisManager.h"
15 #include "AliEmcalJet.h"
16 #include "AliLog.h"
17 #include "AliRhoParameter.h"
18 
20 
21 //________________________________________________________________________
24  fNExclLeadJets(0),
25  fJetRhoMassType(kMd),
26  fPionMassClusters(kFALSE),
27  fHistMdAreavsCent(0)
28 {
29  // Constructor.
30 }
31 
32 //________________________________________________________________________
34  AliAnalysisTaskRhoMassBase(name, histo),
35  fNExclLeadJets(0),
36  fJetRhoMassType(kMd),
37  fPionMassClusters(kFALSE),
38  fHistMdAreavsCent(0)
39 {
40  // Constructor.
41 }
42 
43 //________________________________________________________________________
45 {
46  // User create output objects, called at the beginning of the analysis.
47 
48  if (!fCreateHisto)
49  return;
50 
52 
53  fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt/2.);
54  fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
55  fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
57 }
58 
59 
60 //________________________________________________________________________
62 {
63  // Run the analysis.
64 
65  fOutRhoMass->SetVal(0);
67  fOutRhoMassScaled->SetVal(0);
68 
69  if (!fJets)
70  return kFALSE;
71 
72  const Int_t Njets = fJets->GetEntries();
73 
74  Int_t maxJetIds[] = {-1, -1};
75  Float_t maxJetPts[] = { 0, 0};
76 
77  if (fNExclLeadJets > 0) {
78  for (Int_t ij = 0; ij < Njets; ++ij) {
79  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
80  if (!jet) {
81  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
82  continue;
83  }
84 
85  if (!AcceptJet(jet))
86  continue;
87 
88  if (jet->Pt() > maxJetPts[0]) {
89  maxJetPts[1] = maxJetPts[0];
90  maxJetIds[1] = maxJetIds[0];
91  maxJetPts[0] = jet->Pt();
92  maxJetIds[0] = ij;
93  } else if (jet->Pt() > maxJetPts[1]) {
94  maxJetPts[1] = jet->Pt();
95  maxJetIds[1] = ij;
96  }
97  }
98  if (fNExclLeadJets < 2) {
99  maxJetIds[1] = -1;
100  maxJetPts[1] = 0;
101  }
102  }
103 
104  static Double_t rhomvec[999];
105  static Double_t Evec[999];
106  static Double_t Mvec[999];
107  Int_t NjetAcc = 0;
108 
109  // push all jets within selected acceptance into stack
110  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
111 
112  // exlcuding lead jets
113  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
114  continue;
115 
116  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
117  if (!jet) {
118  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
119  continue;
120  }
121 
122  if (!AcceptJet(jet))
123  continue;
124 
125  // Double_t sumM = GetSumMConstituents(jet);
126  // Double_t sumPt = GetSumPtConstituents(jet);
127  if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) {
128  //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area();
129  // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area();
130  rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
131  fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
132  Evec[NjetAcc] = jet->E();
133  Mvec[NjetAcc] = jet->M();
134  ++NjetAcc;
135  }
136  }
137 
138  if (NjetAcc > 0) {
139  //find median value
140  Double_t rhom = TMath::Median(NjetAcc, rhomvec);
141  fOutRhoMass->SetVal(rhom);
142 
143  Int_t Ntracks = fTracks->GetEntries();
144  Double_t meanM = TMath::Mean(NjetAcc, Mvec);
145  Double_t meanE = TMath::Mean(NjetAcc, Evec);
146  Double_t gamma = 0.;
147  if(meanM>0.) gamma = meanE/meanM;
148  fHistGammaVsNtrack->Fill(Ntracks,gamma);
149 
150  if (fOutRhoMassScaled) {
151  Double_t rhomScaled = rhom * GetScaleFactor(fCent);
152  fOutRhoMassScaled->SetVal(rhomScaled);
153  }
154  }
155 
156  return kTRUE;
157 }
158 
159 //________________________________________________________________________
161 
162  Double_t sum = 0.;
163 
164  AliVParticle *vp;
165  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
166  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
167  if(!vp) continue;
168  sum+=vp->M();
169  }
170  return sum;
171 }
172 
173 //________________________________________________________________________
175 
176  Double_t sum = 0.;
177 
178  AliVParticle *vp;
179  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
180  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
181  if(!vp) continue;
182  sum+=vp->Pt();
183  }
184  return sum;
185 }
186 
187 //________________________________________________________________________
189  //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
190  Double_t sum = 0.;
191  Double_t px = 0.;
192  Double_t py = 0.;
193  Double_t pz = 0.;
194  Double_t E = 0.;
195 
196  if (fTracks) {
197  AliVParticle *vp;
198  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
199  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
200  if(!vp) continue;
201  if(fJetRhoMassType==kMd) sum += TMath::Sqrt(vp->M()*vp->M() + vp->Pt()*vp->Pt()) - vp->Pt(); //sqrt(E^2-P^2+pt^2)=sqrt(E^2-pz^2)
202  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
203  else if(fJetRhoMassType==kMd4) {
204  px+=vp->Px();
205  py+=vp->Py();
206  pz+=vp->Pz();
207  E+=vp->E();
208  }
209  }
210  }
211 
212  if (fCaloClusters) {
213  AliVCluster *vp;
214  for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
215  vp = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
216  if(!vp) continue;
217  TLorentzVector nPart;
218  vp->GetMomentum(nPart, fVertex);
219  Double_t m = 0.;
220  if(fPionMassClusters) m = 0.13957;
221  if(fJetRhoMassType==kMd) sum += TMath::Sqrt(m*m + nPart.Pt()*nPart.Pt()) - nPart.Pt();
222  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(nPart.M()*nPart.M() + nPart.P()*nPart.P()) - nPart.P();
223  else if(fJetRhoMassType==kMd4) {
224  px+=nPart.Px();
225  py+=nPart.Py();
226  pz+=nPart.Pz();
227  E+=nPart.E();
228  }
229  }
230  }
231 
232  if(fJetRhoMassType==kMd4) {
233  Double_t pt = TMath::Sqrt(px*px + py*py);
234  Double_t m2 = E*E - pt*pt - pz*pz;
235  sum = TMath::Sqrt(m2 + pt*pt) - pt;
236  }
237  return sum;
238 }
Double_t Area() const
Definition: AliEmcalJet.h:117
double Double_t
Definition: External.C:58
Definition: External.C:236
Double_t GetSumMConstituents(AliEmcalJet *jet)
Double_t fMinBinPt
min pt in histograms
virtual Double_t GetScaleFactor(Double_t cent)
AliRhoParameter * fOutRhoMassScaled
output rho object
Double_t E() const
Definition: AliEmcalJet.h:106
Double_t GetSumPtConstituents(AliEmcalJet *jet)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
TClonesArray * fCaloClusters
!clusters
Short_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:124
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:125
float Float_t
Definition: External.C:68
Double_t GetMd(AliEmcalJet *jet)
Double_t fCent
!event centrality
TH2F * fHistGammaVsNtrack
rho mass scaled vs. no. of clusters
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:96
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TClonesArray * fTracks
!tracks
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
Double_t fVertex[3]
!event vertex
Bool_t fCreateHisto
whether or not create histograms
ClassImp(AliAnalysisTaskRhoMass) AliAnalysisTaskRhoMass
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
bool Bool_t
Definition: External.C:53
Double_t M() const
Definition: AliEmcalJet.h:107
Int_t fNbins
no. of pt bins