AliPhysics  ff07904 (ff07904)
AliAnalysisTaskRhoSparse.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Calculation of rho from a collection of jets.
4 // If scale function is given the scaled rho will be exported
5 // with the name as "fOutRhoName".apppend("_Scaled").
6 //
7 // Authors: R.Reed, S.Aiola, M.Connors
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  fExcludeOverlaps(0),
27  fRhoCMS(0),
28  fHistOccCorrvsCent(0)
29 {
30  // Constructor.
31 }
32 
33 //________________________________________________________________________
35  AliAnalysisTaskRhoBase(name, histo),
36  fNExclLeadJets(0),
37  fExcludeOverlaps(0),
38  fRhoCMS(0),
39  fHistOccCorrvsCent(0)
40 {
41  // Constructor.
42 }
43 
44 //________________________________________________________________________
46 {
47  if (!fCreateHisto) return;
48 
50 
51  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
53 }
54 
55 //________________________________________________________________________
57 {
58  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
59  {
60  Int_t jet1Track = jet1->TrackAt(i);
61  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
62  {
63  Int_t jet2Track = jet2->TrackAt(j);
64  if (jet1Track == jet2Track)
65  return kTRUE;
66  }
67  }
68  return kFALSE;
69 }
70 
71 //________________________________________________________________________
73 {
74  if(jet->Pt()>5){
75  return kTRUE;
76  }else{
77  return kFALSE;
78  }
79 }
80 
81 
82 //________________________________________________________________________
84 {
85  // Run the analysis.
86 
87  fOutRho->SetVal(0);
88  if (fOutRhoScaled)
89  fOutRhoScaled->SetVal(0);
90 
91  if (!fJets)
92  return kFALSE;
93 
94  const Int_t Njets = fJets->GetEntries();
95 
96  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
97 
98  Int_t NjetsSig = 0;
99  if (sigjets) NjetsSig = sigjets->GetNJets();
100 
101  Int_t maxJetIds[] = {-1, -1};
102  Float_t maxJetPts[] = { 0, 0};
103 
104  if (fNExclLeadJets > 0) {
105  for (Int_t ij = 0; ij < Njets; ++ij) {
106  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
107  if (!jet) {
108  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
109  continue;
110  }
111 
112  if (!AcceptJet(jet))
113  continue;
114 
115  if (jet->Pt() > maxJetPts[0]) {
116  maxJetPts[1] = maxJetPts[0];
117  maxJetIds[1] = maxJetIds[0];
118  maxJetPts[0] = jet->Pt();
119  maxJetIds[0] = ij;
120  } else if (jet->Pt() > maxJetPts[1]) {
121  maxJetPts[1] = jet->Pt();
122  maxJetIds[1] = ij;
123  }
124  }
125  if (fNExclLeadJets < 2) {
126  maxJetIds[1] = -1;
127  maxJetPts[1] = 0;
128  }
129  }
130 
131  static Double_t rhovec[999];
132  Int_t NjetAcc = 0;
133  Double_t TotaljetArea=0;
134  Double_t TotaljetAreaPhys=0;
135  Double_t TotalTPCArea=2*TMath::Pi()*0.9;
136 
137  // push all jets within selected acceptance into stack
138  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
139 
140  // exlcuding lead jets
141  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
142  continue;
143 
144  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
145  if (!jet) {
146  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
147  continue;
148  }
149 
150  TotaljetArea+=jet->Area();
151 
152 
153  if (!AcceptJet(jet))
154  continue;
155 
156  // Search for overlap with signal jets
157  Bool_t isOverlapping = kFALSE;
158  if (sigjets) {
159  for(Int_t j=0;j<NjetsSig;j++)
160  {
161  AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
162  if(!signalJet)
163  continue;
164  if(!IsJetSignal(signalJet))
165  continue;
166 
167  if(IsJetOverlapping(signalJet, jet))
168  {
169  isOverlapping = kTRUE;
170  break;
171  }
172  }
173  }
174 
175  if(fExcludeOverlaps && isOverlapping)
176  continue;
177 
178  //This is to exclude pure ghost jets from the rho calculation
179  if(jet->GetNumberOfTracks()>0)
180  {
181  TotaljetAreaPhys+=jet->Area();
182  rhovec[NjetAcc] = jet->Pt() / jet->Area();
183  ++NjetAcc;
184  }
185  }
186 
187  Double_t OccCorr=TotaljetAreaPhys/TotalTPCArea;
188 
189  if (fCreateHisto)
190  fHistOccCorrvsCent->Fill(fCent, OccCorr);
191 
192  if (NjetAcc > 0) {
193  //find median value
194  Double_t rho = TMath::Median(NjetAcc, rhovec);
195 
196  if(fRhoCMS){
197  rho = rho * OccCorr;
198  }
199 
200  fOutRho->SetVal(rho);
201 
202  if (fOutRhoScaled) {
203  Double_t rhoScaled = rho * GetScaleFactor(fCent);
204  fOutRhoScaled->SetVal(rhoScaled);
205  }
206  }
207 
208  return kTRUE;
209 }
Double_t Area() const
Definition: AliEmcalJet.h:130
TH2F * fHistOccCorrvsCent
! occupancy correction vs. centrality
double Double_t
Definition: External.C:58
Definition: External.C:236
Bool_t fRhoCMS
flag to run CMS method
Bool_t fExcludeOverlaps
exclude background jets that overlap (share at least one track) with anti-KT signal jets ...
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Bool_t IsJetOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
AliRhoParameter * fOutRhoScaled
output rho object
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Int_t GetNJets() const
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation
Double_t fCent
!event centrality
TObjArray fJetCollArray
jet collection array
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:109
virtual Double_t GetScaleFactor(Double_t cent)
AliEmcalList * fOutput
!output list
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.
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
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
bool Bool_t
Definition: External.C:53
Bool_t IsJetSignal(AliEmcalJet *jet1)
Container for jet within the EMCAL jet framework.