AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskRhoMassSparse.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 #include "AliJetContainer.h"
19 
21 
22 //________________________________________________________________________
25  fNExclLeadJets(0),
26  fJetRhoMassType(kMd),
27  fPionMassClusters(kFALSE),
28  fHistMdAreavsCent(0)
29 {
30  // Constructor.
31 }
32 
33 //________________________________________________________________________
35  AliAnalysisTaskRhoMassBase(name, histo),
36  fNExclLeadJets(0),
37  fJetRhoMassType(kMd),
38  fPionMassClusters(kFALSE),
39  fHistMdAreavsCent(0)
40 {
41  // Constructor.
42 }
43 
44 //________________________________________________________________________
46 {
47  // User create output objects, called at the beginning of the analysis.
48 
49  if (!fCreateHisto)
50  return;
51 
53 
54  fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt/2.);
55  fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
56  fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
58 
59  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
61 
62 }
63 
64 //________________________________________________________________________
66 {
67  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
68  {
69  Int_t jet1Track = jet1->TrackAt(i);
70  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
71  {
72  Int_t jet2Track = jet2->TrackAt(j);
73  if (jet1Track == jet2Track)
74  return kTRUE;
75  }
76  }
77  return kFALSE;
78 }
79 
80 //________________________________________________________________________
82 {
83  if(jet->Pt()>5){
84  return kTRUE;
85  }else{
86  return kFALSE;
87  }
88 }
89 
90 
91 //________________________________________________________________________
93 {
94  // Run the analysis.
95 
96  fOutRhoMass->SetVal(0);
98  fOutRhoMassScaled->SetVal(0);
99 
100  if (!fJets)
101  return kFALSE;
102 
103  const Int_t Njets = fJets->GetEntries();
104 
105  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
106 
107  Int_t NjetsSig = 0;
108  if (sigjets) NjetsSig = sigjets->GetNJets();
109 
110  Int_t maxJetIds[] = {-1, -1};
111  Float_t maxJetPts[] = { 0, 0};
112 
113  if (fNExclLeadJets > 0) {
114  for (Int_t ij = 0; ij < Njets; ++ij) {
115  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
116  if (!jet) {
117  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
118  continue;
119  }
120 
121  if (!AcceptJet(jet))
122  continue;
123 
124  if (jet->Pt() > maxJetPts[0]) {
125  maxJetPts[1] = maxJetPts[0];
126  maxJetIds[1] = maxJetIds[0];
127  maxJetPts[0] = jet->Pt();
128  maxJetIds[0] = ij;
129  } else if (jet->Pt() > maxJetPts[1]) {
130  maxJetPts[1] = jet->Pt();
131  maxJetIds[1] = ij;
132  }
133  }
134  if (fNExclLeadJets < 2) {
135  maxJetIds[1] = -1;
136  maxJetPts[1] = 0;
137  }
138  }
139 
140  static Double_t rhomvec[999];
141  static Double_t Evec[999];
142  static Double_t Mvec[999];
143  Int_t NjetAcc = 0;
144  Double_t TotaljetArea=0;
145  Double_t TotaljetAreaPhys=0;
146 
147  // push all jets within selected acceptance into stack
148  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
149 
150  // exlcuding lead jets
151  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
152  continue;
153 
154  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
155  if (!jet) {
156  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
157  continue;
158  }
159 
160  TotaljetArea+=jet->Area();
161 
162  if(jet->Pt()>0.1){
163  TotaljetAreaPhys+=jet->Area();
164  }
165 
166  if (!AcceptJet(jet))
167  continue;
168 
169  // Search for overlap with signal jets
170  Bool_t isOverlapping = kFALSE;
171  if (sigjets) {
172  for(Int_t j=0;j<NjetsSig;j++)
173  {
174  AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
175  if(!signalJet)
176  continue;
177  if(!IsJetSignal(signalJet))
178  continue;
179 
180  if(IsJetOverlapping(signalJet, jet))
181  {
182  isOverlapping = kTRUE;
183  break;
184  }
185  }
186  }
187 
188  if(isOverlapping)
189  continue;
190 
191  // Double_t sumM = GetSumMConstituents(jet);
192  // Double_t sumPt = GetSumPtConstituents(jet);
193 
194  if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) {
195  //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area();
196  // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area();
197  rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
198  fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
199  Evec[NjetAcc] = jet->E();
200  Mvec[NjetAcc] = jet->M();
201  ++NjetAcc;
202  }
203  }
204 
205  Double_t OccCorr=0.0;
206  if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
207 
208  if (fCreateHisto)
209  fHistOccCorrvsCent->Fill(fCent, OccCorr);
210 
211 
212  if (NjetAcc > 0) {
213  //find median value
214  Double_t rhom = TMath::Median(NjetAcc, rhomvec);
215  if(fRhoCMS){
216  rhom = rhom * OccCorr;
217  }
218 
219  fOutRhoMass->SetVal(rhom);
220 
221  Int_t Ntracks = fTracks->GetEntries();
222  Double_t meanM = TMath::Mean(NjetAcc, Mvec);
223  Double_t meanE = TMath::Mean(NjetAcc, Evec);
224  Double_t gamma = 0.;
225  if(meanM>0.) gamma = meanE/meanM;
226  fHistGammaVsNtrack->Fill(Ntracks,gamma);
227 
228  if (fOutRhoMassScaled) {
229  Double_t rhomScaled = rhom * GetScaleFactor(fCent);
230  fOutRhoMassScaled->SetVal(rhomScaled);
231  }
232  }
233 
234  return kTRUE;
235 }
236 
237 //________________________________________________________________________
239 
240  Double_t sum = 0.;
241 
242  AliVParticle *vp;
243  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
244  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
245  if(!vp) continue;
246  sum+=vp->M();
247  }
248  return sum;
249 }
250 
251 //________________________________________________________________________
253 
254  Double_t sum = 0.;
255 
256  AliVParticle *vp;
257  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
258  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
259  if(!vp) continue;
260  sum+=vp->Pt();
261  }
262  return sum;
263 }
264 
265 //________________________________________________________________________
267  //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
268  Double_t sum = 0.;
269  Double_t px = 0.;
270  Double_t py = 0.;
271  Double_t pz = 0.;
272  Double_t E = 0.;
273 
274  if (fTracks) {
275  AliVParticle *vp;
276  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
277  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
278  if(!vp) continue;
279  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)
280  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
281  else if(fJetRhoMassType==kMd4) {
282  px+=vp->Px();
283  py+=vp->Py();
284  pz+=vp->Pz();
285  E+=vp->E();
286  }
287  }
288  }
289 
290  if (fCaloClusters) {
291  AliVCluster *vp;
292  for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
293  vp = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
294  if(!vp) continue;
295  TLorentzVector nPart;
296  vp->GetMomentum(nPart, fVertex);
297  Double_t m = 0.;
298  if(fPionMassClusters) m = 0.13957;
299  if(fJetRhoMassType==kMd) sum += TMath::Sqrt(m*m + nPart.Pt()*nPart.Pt()) - nPart.Pt();
300  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(nPart.M()*nPart.M() + nPart.P()*nPart.P()) - nPart.P();
301  else if(fJetRhoMassType==kMd4) {
302  px+=nPart.Px();
303  py+=nPart.Py();
304  pz+=nPart.Pz();
305  E+=nPart.E();
306  }
307  }
308  }
309 
310  if(fJetRhoMassType==kMd4) {
311  Double_t pt = TMath::Sqrt(px*px + py*py);
312  Double_t m2 = E*E - pt*pt - pz*pz;
313  sum = TMath::Sqrt(m2 + pt*pt) - pt;
314  }
315  return sum;
316 }
Double_t Area() const
Definition: AliEmcalJet.h:97
ClassImp(AliAnalysisTaskRhoMassSparse) AliAnalysisTaskRhoMassSparse
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:86
TList * fOutput
!output list
TH2F * fHistOccCorrvsCent
Md/Area vs cent for all kt clusters.
Double_t GetSumPtConstituents(AliEmcalJet *jet)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:106
TClonesArray * fCaloClusters
!clusters
Short_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:104
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:105
Int_t GetNJets() const
Double_t fCent
!event centrality
TH2F * fHistGammaVsNtrack
rho mass scaled vs. no. of clusters
TObjArray fJetCollArray
jet collection array
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:76
Bool_t IsJetOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
Double_t fMaxBinPt
max pt in histograms
TClonesArray * fTracks
!tracks
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:127
Double_t fVertex[3]
!event vertex
Bool_t fCreateHisto
whether or not create histograms
AliEmcalJet * GetAcceptJet(Int_t i) const
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
Double_t M() const
Definition: AliEmcalJet.h:87
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins