AliPhysics  59e0e03 (59e0e03)
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  fRhoCMS(0),
27  fHistOccCorrvsCent(0)
28 {
29  // Constructor.
30 }
31 
32 //________________________________________________________________________
34  AliAnalysisTaskRhoBase(name, histo),
35  fNExclLeadJets(0),
36  fRhoCMS(0),
37  fHistOccCorrvsCent(0)
38 {
39  // Constructor.
40 }
41 
42 //________________________________________________________________________
44 {
45  if (!fCreateHisto) return;
46 
48 
49  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
51 }
52 
53 //________________________________________________________________________
55 {
56  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
57  {
58  Int_t jet1Track = jet1->TrackAt(i);
59  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
60  {
61  Int_t jet2Track = jet2->TrackAt(j);
62  if (jet1Track == jet2Track)
63  return kTRUE;
64  }
65  }
66  return kFALSE;
67 }
68 
69 //________________________________________________________________________
71 {
72  if(jet->Pt()>5){
73  return kTRUE;
74  }else{
75  return kFALSE;
76  }
77 }
78 
79 
80 //________________________________________________________________________
82 {
83  // Run the analysis.
84 
85  fOutRho->SetVal(0);
86  if (fOutRhoScaled)
87  fOutRhoScaled->SetVal(0);
88 
89  if (!fJets)
90  return kFALSE;
91 
92  const Int_t Njets = fJets->GetEntries();
93 
94  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
95 
96  Int_t NjetsSig = 0;
97  if (sigjets) NjetsSig = sigjets->GetNJets();
98 
99  Int_t maxJetIds[] = {-1, -1};
100  Float_t maxJetPts[] = { 0, 0};
101 
102  if (fNExclLeadJets > 0) {
103  for (Int_t ij = 0; ij < Njets; ++ij) {
104  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
105  if (!jet) {
106  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
107  continue;
108  }
109 
110  if (!AcceptJet(jet))
111  continue;
112 
113  if (jet->Pt() > maxJetPts[0]) {
114  maxJetPts[1] = maxJetPts[0];
115  maxJetIds[1] = maxJetIds[0];
116  maxJetPts[0] = jet->Pt();
117  maxJetIds[0] = ij;
118  } else if (jet->Pt() > maxJetPts[1]) {
119  maxJetPts[1] = jet->Pt();
120  maxJetIds[1] = ij;
121  }
122  }
123  if (fNExclLeadJets < 2) {
124  maxJetIds[1] = -1;
125  maxJetPts[1] = 0;
126  }
127  }
128 
129  static Double_t rhovec[999];
130  Int_t NjetAcc = 0;
131  Double_t TotaljetArea=0;
132  Double_t TotaljetAreaPhys=0;
133 
134  // push all jets within selected acceptance into stack
135  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
136 
137  // exlcuding lead jets
138  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
139  continue;
140 
141  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
142  if (!jet) {
143  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
144  continue;
145  }
146 
147  TotaljetArea+=jet->Area();
148 
149  if(jet->Pt()>0.1){
150  TotaljetAreaPhys+=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(isOverlapping)
176  continue;
177 
178  if(jet->Pt()>0.1){
179  rhovec[NjetAcc] = jet->Pt() / jet->Area();
180  ++NjetAcc;
181  }
182  }
183 
184  Double_t OccCorr=0.0;
185  if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
186 
187  if (fCreateHisto)
188  fHistOccCorrvsCent->Fill(fCent, OccCorr);
189 
190  if (NjetAcc > 0) {
191  //find median value
192  Double_t rho = TMath::Median(NjetAcc, rhovec);
193 
194  if(fRhoCMS){
195  rho = rho * OccCorr;
196  }
197 
198  fOutRho->SetVal(rho);
199 
200  if (fOutRhoScaled) {
201  Double_t rhoScaled = rho * GetScaleFactor(fCent);
202  fOutRhoScaled->SetVal(rhoScaled);
203  }
204  }
205 
206  return kTRUE;
207 }
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Definition: External.C:236
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
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.