AliPhysics  master (3d17d9d)
AliAnalysisTaskRhoMassSparse.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2014, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
28 
29 #include <TClonesArray.h>
30 #include <TMath.h>
31 
32 #include "AliAnalysisManager.h"
33 #include "AliEmcalJet.h"
34 #include "AliLog.h"
35 #include "AliRhoParameter.h"
36 #include "AliJetContainer.h"
37 
39 
42  fNExclLeadJets(0),
43  fJetRhoMassType(kMd),
44  fPionMassClusters(kFALSE),
45  fHistMdAreavsCent(0)
46 {
47 }
48 
50  AliAnalysisTaskRhoMassBase(name, histo),
51  fNExclLeadJets(0),
52  fJetRhoMassType(kMd),
53  fPionMassClusters(kFALSE),
54  fHistMdAreavsCent(0)
55 {
56 }
57 
59 {
60  if (!fCreateHisto)
61  return;
62 
64 
65  fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt/2.);
66  fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
67  fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
69 
70  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
72 
73 }
74 
76 {
77  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
78  {
79  Int_t jet1Track = jet1->TrackAt(i);
80  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
81  {
82  Int_t jet2Track = jet2->TrackAt(j);
83  if (jet1Track == jet2Track)
84  return kTRUE;
85  }
86  }
87  return kFALSE;
88 }
89 
91 {
92  if(jet->Pt()>5){
93  return kTRUE;
94  }else{
95  return kFALSE;
96  }
97 }
98 
99 
101 {
102  fOutRhoMass->SetVal(0);
103  if (fOutRhoMassScaled)
104  fOutRhoMassScaled->SetVal(0);
105 
106  if (!fJets)
107  return kFALSE;
108 
109  const Int_t Njets = fJets->GetEntries();
110 
111  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
112 
113  Int_t NjetsSig = 0;
114  if (sigjets) NjetsSig = sigjets->GetNJets();
115 
116  Int_t maxJetIds[] = {-1, -1};
117  Float_t maxJetPts[] = { 0, 0};
118 
119  if (fNExclLeadJets > 0) {
120  for (Int_t ij = 0; ij < Njets; ++ij) {
121  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
122  if (!jet) {
123  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
124  continue;
125  }
126 
127  if (!AcceptJet(jet))
128  continue;
129 
130  if (jet->Pt() > maxJetPts[0]) {
131  maxJetPts[1] = maxJetPts[0];
132  maxJetIds[1] = maxJetIds[0];
133  maxJetPts[0] = jet->Pt();
134  maxJetIds[0] = ij;
135  } else if (jet->Pt() > maxJetPts[1]) {
136  maxJetPts[1] = jet->Pt();
137  maxJetIds[1] = ij;
138  }
139  }
140  if (fNExclLeadJets < 2) {
141  maxJetIds[1] = -1;
142  maxJetPts[1] = 0;
143  }
144  }
145 
146  static Double_t rhomvec[999];
147  static Double_t Evec[999];
148  static Double_t Mvec[999];
149  Int_t NjetAcc = 0;
150  Double_t TotaljetArea=0;
151  Double_t TotaljetAreaPhys=0;
152 
153  // push all jets within selected acceptance into stack
154  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
155 
156  // exlcuding lead jets
157  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
158  continue;
159 
160  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
161  if (!jet) {
162  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
163  continue;
164  }
165 
166  TotaljetArea+=jet->Area();
167 
168  if(jet->Pt()>0.1){
169  TotaljetAreaPhys+=jet->Area();
170  }
171 
172  if (!AcceptJet(jet))
173  continue;
174 
175  // Search for overlap with signal jets
176  Bool_t isOverlapping = kFALSE;
177  if (sigjets) {
178  for(Int_t j=0;j<NjetsSig;j++)
179  {
180  AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
181  if(!signalJet)
182  continue;
183  if(!IsJetSignal(signalJet))
184  continue;
185 
186  if(IsJetOverlapping(signalJet, jet))
187  {
188  isOverlapping = kTRUE;
189  break;
190  }
191  }
192  }
193 
194  if(isOverlapping)
195  continue;
196 
197  // Double_t sumM = GetSumMConstituents(jet);
198  // Double_t sumPt = GetSumPtConstituents(jet);
199 
200  if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) {
201  //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area();
202  // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area();
203  rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
204  fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
205  Evec[NjetAcc] = jet->E();
206  Mvec[NjetAcc] = jet->M();
207  ++NjetAcc;
208  }
209  }
210 
211  Double_t OccCorr=0.0;
212  if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
213 
214  if (fCreateHisto)
215  fHistOccCorrvsCent->Fill(fCent, OccCorr);
216 
217 
218  if (NjetAcc > 0) {
219  //find median value
220  Double_t rhom = TMath::Median(NjetAcc, rhomvec);
221  if(fRhoCMS){
222  rhom = rhom * OccCorr;
223  }
224 
225  fOutRhoMass->SetVal(rhom);
226 
227  Int_t Ntracks = fTracks->GetEntries();
228  Double_t meanM = TMath::Mean(NjetAcc, Mvec);
229  Double_t meanE = TMath::Mean(NjetAcc, Evec);
230  Double_t gamma = 0.;
231  if(meanM>0.) gamma = meanE/meanM;
232  fHistGammaVsNtrack->Fill(Ntracks,gamma);
233 
234  if (fOutRhoMassScaled) {
235  Double_t rhomScaled = rhom * GetScaleFactor(fCent);
236  fOutRhoMassScaled->SetVal(rhomScaled);
237  }
238  }
239 
240  return kTRUE;
241 }
242 
244 
245  Double_t sum = 0.;
246 
247  AliVParticle *vp;
248  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
249  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
250  if(!vp) continue;
251  sum+=vp->M();
252  }
253  return sum;
254 }
255 
257 
258  Double_t sum = 0.;
259 
260  AliVParticle *vp;
261  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
262  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
263  if(!vp) continue;
264  sum+=vp->Pt();
265  }
266  return sum;
267 }
268 
270  Double_t sum = 0.;
271  Double_t px = 0.;
272  Double_t py = 0.;
273  Double_t pz = 0.;
274  Double_t E = 0.;
275 
276  if (fTracks) {
277  AliVParticle *vp;
278  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
279  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
280  if(!vp) continue;
281  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)
282  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
283  else if(fJetRhoMassType==kMd4) {
284  px+=vp->Px();
285  py+=vp->Py();
286  pz+=vp->Pz();
287  E+=vp->E();
288  }
289  }
290  }
291 
292  if (fCaloClusters) {
293  AliVCluster *vp;
294  for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
295  vp = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
296  if(!vp) continue;
297  TLorentzVector nPart;
298  vp->GetMomentum(nPart, fVertex);
299  Double_t m = 0.;
300  if(fPionMassClusters) m = 0.13957;
301  if(fJetRhoMassType==kMd) sum += TMath::Sqrt(m*m + nPart.Pt()*nPart.Pt()) - nPart.Pt();
302  else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(nPart.M()*nPart.M() + nPart.P()*nPart.P()) - nPart.P();
303  else if(fJetRhoMassType==kMd4) {
304  px+=nPart.Px();
305  py+=nPart.Py();
306  pz+=nPart.Pz();
307  E+=nPart.E();
308  }
309  }
310  }
311 
312  if(fJetRhoMassType==kMd4) {
313  Double_t pt = TMath::Sqrt(px*px + py*py);
314  Double_t m2 = E*E - pt*pt - pz*pz;
315  sum = TMath::Sqrt(m2 + pt*pt) - pt;
316  }
317  return sum;
318 }
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Double_t GetSumMConstituents(AliEmcalJet *jet)
Definition: External.C:236
Calculation of rho mass from a collection of jets.
Double_t fMinBinPt
min pt in histograms
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:137
virtual Double_t GetScaleFactor(Double_t cent)
Get scale factor.
AliRhoParameter * fOutRhoMassScaled
! output scaled rho object
Double_t E() const
Definition: AliEmcalJet.h:119
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation
Bool_t fRhoCMS
flag to run CMS method
TH2F * fHistOccCorrvsCent
! occupancy correction vs. centrality
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
Double_t GetSumPtConstituents(AliEmcalJet *jet)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Bool_t fPionMassClusters
assume pion mass for clusters
JetRhoMassType fJetRhoMassType
method for rho_m calculation
TClonesArray * fCaloClusters
!clusters
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
float Float_t
Definition: External.C:68
Int_t GetNJets() const
Double_t fCent
!event centrality
void UserCreateOutputObjects()
User create output objects, called at the beginning of the analysis.
rho_m using addition of 4-vectors
Double_t GetMd(AliEmcalJet *jet)
Get md as defined in http://arxiv.org/pdf/1211.2811.pdf.
TH2F * fHistGammaVsNtrack
! Gamma(<E>/<M>) vs Ntrack
TH2F * fHistMdAreavsCent
! Md/Area vs cent for all kt clusters
TObjArray fJetCollArray
jet collection array
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:109
Bool_t IsJetOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TClonesArray * fTracks
!tracks
Base class for rho mass calculation.
Double_t fVertex[3]
!event vertex
AliRhoParameter * fOutRhoMass
! output rho object
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:51
void UserCreateOutputObjects()
User create output objects, called at the beginning of the analysis.
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
bool Bool_t
Definition: External.C:53
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins