AliPhysics  c7b8e89 (c7b8e89)
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 <AliVEventHandler.h>
16 #include "AliEmcalJet.h"
17 #include "AliLog.h"
18 #include "AliRhoParameter.h"
19 #include "AliJetContainer.h"
20 
22 
23 //________________________________________________________________________
26  fNExclLeadJets(0),
27  fExcludeOverlaps(0),
28  fRhoCMS(0),
29  fHistOccCorrvsCent(0)
30 {
31  // Constructor.
32 }
33 
34 //________________________________________________________________________
36  AliAnalysisTaskRhoBase(name, histo),
37  fNExclLeadJets(0),
38  fExcludeOverlaps(0),
39  fRhoCMS(0),
40  fHistOccCorrvsCent(0)
41 {
42  // Constructor.
43 }
44 
45 //________________________________________________________________________
47 {
48  if (!fCreateHisto) return;
49 
51 
52  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
54 }
55 
56 //________________________________________________________________________
58 {
59  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
60  {
61  Int_t jet1Track = jet1->TrackAt(i);
62  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
63  {
64  Int_t jet2Track = jet2->TrackAt(j);
65  if (jet1Track == jet2Track)
66  return kTRUE;
67  }
68  }
69  return kFALSE;
70 }
71 
72 //________________________________________________________________________
74 {
75  if(jet->Pt()>5){
76  return kTRUE;
77  }else{
78  return kFALSE;
79  }
80 }
81 
82 
83 //________________________________________________________________________
85 {
86  // Run the analysis.
87 
88  fOutRho->SetVal(0);
89  if (fOutRhoScaled)
90  fOutRhoScaled->SetVal(0);
91 
92  if (!fJets)
93  return kFALSE;
94  const Int_t Njets = fJets->GetEntries();
95 
96  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
97  Int_t NjetsSig = 0;
98  if (sigjets) NjetsSig = sigjets->GetNJets();
99 
100  Int_t maxJetIds[] = {-1, -1};
101  Float_t maxJetPts[] = { 0, 0};
102 
103  if (fNExclLeadJets > 0) {
104  for (Int_t ij = 0; ij < Njets; ++ij) {
105  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
106  if (!jet) {
107  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
108  continue;
109  }
110 
111  if (!AcceptJet(jet))
112  continue;
113 
114  //Search the two leading KT jets
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 TotaljetAreaPhys=0;
134  Double_t TotalTPCArea=2*TMath::Pi()*0.9;
135 
136  // push all jets within selected acceptance into stack
137  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
138 
139  // exlcuding leading background jets (could be signal)
140  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
141  continue;
142 
143  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
144  if (!jet) {
145  AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
146  continue;
147  }
148 
149  if (!AcceptJet(jet))
150  continue;
151 
152  // Search for overlap with signal jets
153  Bool_t isOverlapping = kFALSE;
154  if (sigjets) {
155  for(Int_t j=0;j<NjetsSig;j++)
156  {
157  AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
158  if(!signalJet)
159  continue;
160  if(!IsJetSignal(signalJet))
161  continue;
162 
163  if(IsJetOverlapping(signalJet, jet))
164  {
165  isOverlapping = kTRUE;
166  break;
167  }
168  }
169  }
170 
171  if(fExcludeOverlaps && isOverlapping)
172  continue;
173 
174  //This is to exclude pure ghost jets from the rho calculation
175  if(jet->GetNumberOfTracks()>0)
176  {
177  TotaljetAreaPhys+=jet->Area();
178  rhovec[NjetAcc] = jet->Pt() / jet->Area();
179  ++NjetAcc;
180  }
181  }
182 
183  Double_t OccCorr=TotaljetAreaPhys/TotalTPCArea;
184 
185  if (fCreateHisto)
186  fHistOccCorrvsCent->Fill(fCent, OccCorr);
187 
188  if (NjetAcc > 0) {
189  //find median value
190  Double_t rho = TMath::Median(NjetAcc, rhovec);
191 
192  if(fRhoCMS){
193  rho = rho * OccCorr;
194  }
195 
196  fOutRho->SetVal(rho);
197 
198  if (fOutRhoScaled) {
199  Double_t rhoScaled = rho * GetScaleFactor(fCent);
200  fOutRhoScaled->SetVal(rhoScaled);
201  }
202  }
203 
204  return kTRUE;
205 }
211  const char *nTracks,
212  const char *nClusters,
213  const char *nRho,
214  Double_t jetradius,
215  UInt_t acceptance,
218  const Bool_t histo,
219  const char *nJetsSig,
220  const char *cutType,
221  Double_t jetptcut,
222  Double_t jetareacut,
223  Double_t emcareacut,
224  const char *suffix
225  )
226 {
227 
228  // Get the pointer to the existing analysis manager via the static access method.
229  //==============================================================================
230  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
231  if (!mgr)
232  {
233  ::Error("AddTaskRhoSparse", "No analysis manager to connect to.");
234  return NULL;
235  }
236 
237  // Check the analysis type using the event handlers connected to the analysis manager.
238  //==============================================================================
239  AliVEventHandler* handler = mgr->GetInputEventHandler();
240  if (!handler)
241  {
242  ::Error("AddTaskEmcalJetPerformance", "This task requires an input event handler");
243  return 0;
244  }
245 
246  enum EDataType_t {
247  kUnknown,
248  kESD,
249  kAOD
250  };
251 
252  EDataType_t dataType = kUnknown;
253 
254  if (handler->InheritsFrom("AliESDInputHandler")) {
255  dataType = kESD;
256  }
257  else if (handler->InheritsFrom("AliAODInputHandler")) {
258  dataType = kAOD;
259  }
260  //-------------------------------------------------------
261  // Init the task and do settings
262  //-------------------------------------------------------
263  TString trackName(nTracks);
264  TString clusName(nClusters);
265 
266  if (trackName == "usedefault") {
267  if (dataType == kESD) {
268  trackName = "Tracks";
269  }
270  else if (dataType == kAOD) {
271  trackName = "tracks";
272  }
273  else {
274  trackName = "";
275  }
276  }
277 
278  if (clusName == "usedefault") {
279  if (dataType == kESD) {
280  clusName = "CaloClusters";
281  }
282  else if (dataType == kAOD) {
283  clusName = "caloClusters";
284  }
285  else {
286  clusName = "";
287  }
288  }
289 
290  TString name("AliAnalysisTaskRhoSparse");
291  if (strcmp(suffix,"") != 0) {
292  name += "_";
293  name += suffix;
294  }
295 
296  AliAnalysisTaskRhoSparse* mgrTask = (AliAnalysisTaskRhoSparse*)(mgr->GetTask(name.Data()));
297  if (mgrTask) return mgrTask;
298 
299  AliAnalysisTaskRhoSparse *rhotask = new AliAnalysisTaskRhoSparse(name, histo);
300  rhotask->SetHistoBins(1000,-0.1,9.9);
301  rhotask->SetRhoCMS(1);
302  rhotask->SetSmallSystem(1);
303  rhotask->SetOutRhoName(nRho);
304 
305  AliTrackContainer* trackCont;
306  AliParticleContainer* partCont;
307  if (trackName == "mcparticles") {
308  AliMCParticleContainer* mcpartCont = rhotask->AddMCParticleContainer(trackName);
309  mcpartCont->SelectPhysicalPrimaries(kTRUE);
310  }
311  else if (trackName == "tracks" || trackName == "Tracks") {
312  trackCont = rhotask->AddTrackContainer(trackName);
313  }
314  else if (!trackName.IsNull()) {
315  rhotask->AddParticleContainer(trackName);
316  }
317  partCont = rhotask->GetParticleContainer(0);
318 
319  AliClusterContainer *clusterCont = rhotask->AddClusterContainer(clusName);
320  if (clusterCont) {
321  clusterCont->SetClusECut(0.);
322  clusterCont->SetClusPtCut(0.);
323  clusterCont->SetClusHadCorrEnergyCut(0.3);
324  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
325  }
326 
327  //
328  AliJetContainer *bkgJetCont = rhotask->AddJetContainer(jetType, AliJetContainer::kt_algorithm, rscheme, jetradius, acceptance, partCont, clusterCont);
329  if (bkgJetCont) {
330  //why?? bkgJetCont->SetJetAreaCut(jetareacut);
331  //why?? bkgJetCont->SetAreaEmcCut(emcareacut);
332  bkgJetCont->SetJetPtCut(0.);
333  }
334 
335  //This is only for excluding overlaps with signal jets in case signal jets are provided
336  AliJetContainer *sigJetCont;
337  if(nJetsSig!=0)
338  {
339  sigJetCont = rhotask->AddJetContainer(nJetsSig,cutType,jetradius);
340  if (sigJetCont)
341  {
342  sigJetCont->SetJetAreaCut(jetareacut);
343  sigJetCont->SetAreaEmcCut(emcareacut);
344  sigJetCont->SetJetPtCut(jetptcut);
345  sigJetCont->ConnectParticleContainer(trackCont);
346  sigJetCont->ConnectClusterContainer(clusterCont);
347  }
348  }
349  //-------------------------------------------------------
350  // Final settings, pass to manager and set the containers
351  //-------------------------------------------------------
352  mgr->AddTask(rhotask);
353 
354  // Create containers for input/output
355  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
356  if (histo) {
357  TString contname(name);
358  contname += "_histos";
359  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
360  TList::Class(),AliAnalysisManager::kOutputContainer,
361  Form("%s", AliAnalysisManager::GetCommonFileName()));
362  mgr->ConnectOutput(rhotask, 1, coutput1);
363  }
364 
365  return rhotask;
366 }
367 
void SetAreaEmcCut(Double_t a=0.99)
Double_t Area() const
Definition: AliEmcalJet.h:130
Calculation of background rho for sparse events (pp and pPb).
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
void SetOutRhoName(const char *name)
Container with name, TClonesArray and cuts for particles.
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Bool_t fExcludeOverlaps
exclude background jets that overlap (share at least one track) with anti-KT signal jets ...
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Bool_t IsJetOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
AliRhoParameter * fOutRhoScaled
output rho object
void SetSmallSystem(Bool_t setter=kTRUE)
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Int_t GetNJets() const
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation
Double_t fCent
!event centrality
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
TObjArray fJetCollArray
jet collection array
TClonesArray * fJets
! jets
static AliAnalysisTaskRhoSparse * AddTaskRhoSparse(const char *nTracks="usedefault", const char *nClusters="usedefault", const char *nRho="Rho", Double_t jetradius=0.2, UInt_t acceptance=AliEmcalJet::kTPCfid, AliJetContainer::EJetType_t jetType=AliJetContainer::kChargedJet, AliJetContainer::ERecoScheme_t rscheme=AliJetContainer::pt_scheme, const Bool_t histo=kFALSE, const char *nJetsSig="", const char *cutType="TPC", Double_t jetptcut=0.0, Double_t jetareacut=0.01, Double_t emcareacut=0, const char *suffix="")
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
virtual Double_t GetScaleFactor(Double_t cent)
AliEmcalList * fOutput
!output list
void SetHistoBins(Int_t nbins, Double_t min, Double_t max)
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.
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
void SelectPhysicalPrimaries(Bool_t s)
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 SetClusPtCut(Double_t cut)
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
bool Bool_t
Definition: External.C:53
EDataType_t
Switch for the data type.
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
void ConnectClusterContainer(AliClusterContainer *c)
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
Bool_t IsJetSignal(AliEmcalJet *jet1)
Container for jet within the EMCAL jet framework.
void SetJetAreaCut(Float_t cut)
void SetClusHadCorrEnergyCut(Double_t cut)