AliPhysics  b0c77bb (b0c77bb)
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  fUseTPCArea(0),
30  fExcludeAreaExcludedJets(0),
31  fHistOccCorrvsCent(0)
32 {
33  // Constructor.
34 }
35 
36 //________________________________________________________________________
38  AliAnalysisTaskRhoBase(name, histo),
39  fNExclLeadJets(0),
40  fExcludeOverlaps(0),
41  fRhoCMS(0),
42  fUseTPCArea(0),
43  fExcludeAreaExcludedJets(0),
44  fHistOccCorrvsCent(0)
45 {
46  // Constructor.
47 }
48 
49 //________________________________________________________________________
51 {
52  if (!fCreateHisto) return;
53 
55 
56  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
58 }
59 
60 //________________________________________________________________________
62 {
63  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
64  {
65  Int_t jet1Track = jet1->TrackAt(i);
66  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
67  {
68  Int_t jet2Track = jet2->TrackAt(j);
69  if (jet1Track == jet2Track)
70  return kTRUE;
71  }
72  }
73  return kFALSE;
74 }
75 
76 //________________________________________________________________________
78 {
79  if(jet->Pt()>5){
80  return kTRUE;
81  }else{
82  return kFALSE;
83  }
84 }
85 
86 //________________________________________________________________________
88 {
89  // Run the analysis.
90 
91  fOutRho->SetVal(0);
92  if (fOutRhoScaled)
93  fOutRhoScaled->SetVal(0);
94 
95  if (!fJets)
96  return kFALSE;
97  const Int_t Njets = fJets->GetEntries();
98 
99  AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
100  Int_t NjetsSig = 0;
101  if (sigjets) NjetsSig = sigjets->GetNJets();
102 
103  Int_t maxJetIds[] = {-1, -1};
104  Float_t maxJetPts[] = { 0, 0};
105 
106  if (fNExclLeadJets > 0) {
107  for (Int_t ij = 0; ij < Njets; ++ij) {
108  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
109  if (!jet) {
110  AliError(Form("%s: Could not receive jet %d", GetName(), ij));
111  continue;
112  }
113 
114  if (!AcceptJet(jet))
115  continue;
116 
117  //Search the two leading KT jets
118  if (jet->Pt() > maxJetPts[0]) {
119  maxJetPts[1] = maxJetPts[0];
120  maxJetIds[1] = maxJetIds[0];
121  maxJetPts[0] = jet->Pt();
122  maxJetIds[0] = ij;
123  } else if (jet->Pt() > maxJetPts[1]) {
124  maxJetPts[1] = jet->Pt();
125  maxJetIds[1] = ij;
126  }
127  }
128  if (fNExclLeadJets < 2) {
129  maxJetIds[1] = -1;
130  maxJetPts[1] = 0;
131  }
132  }
133 
134  static Double_t rhovec[999];
135  Int_t NjetAcc = 0;
136  Double_t TotaljetAreaPhys=0;
137  Double_t TotalAreaCovered=0;
138  Double_t TotalTPCArea=2*TMath::Pi()*0.9;
139 
140  // push all jets within selected acceptance into stack
141  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
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  // Take into account the area of real jets (no pure ghost jets)
150  // and also the area of later excluded jets.
151  // Some of these jets do not contribute to the rho calculation
152  // but to the factor with which this rho is scaled down
154  {
155  // Total area of physical jets that are used for the rho calculation
156  if(jet->GetNumberOfTracks()>0)
157  {
158  TotaljetAreaPhys+=jet->Area();
159  }
160  // Total area of all found jets ghost+phyical jets. This is a proxy
161  // for the available detector area in which the rho could have been calculated
162  TotalAreaCovered+=jet->Area();
163  }
164  // Exclude leading background jets (could be signal)
165  if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
166  continue;
167  // Exclude background jets that do not fullfill basic cuts defined in AliJetContainer
168  if (!AcceptJet(jet))
169  continue;
170 
171  // Search for overlap with signal jets
172  Bool_t isOverlapping = kFALSE;
173  if (sigjets) {
174  for(Int_t j=0;j<NjetsSig;j++)
175  {
176  AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
177  if(!signalJet)
178  continue;
179  if(!IsJetSignal(signalJet))
180  continue;
181 
182  if(IsJetOverlapping(signalJet, jet))
183  {
184  isOverlapping = kTRUE;
185  break;
186  }
187  }
188  }
189  // Exclude background jets that overlap with anti-kT signal jets
190  if(fExcludeOverlaps && isOverlapping)
191  continue;
192 
193  // Take into account only real jets (no pure ghost jets) that also
194  // contribute to the rho calculation
196  {
197  // Total area of physical jets that are used for the rho calculation
198  if(jet->GetNumberOfTracks()>0)
199  {
200  TotaljetAreaPhys+=jet->Area();
201  }
202  // Total area of all found jets ghost+phyical jets. This is a proxy
203  // for the available detector area in which the rho could have been calculated
204  TotalAreaCovered+=jet->Area();
205  }
206 
207  // Exclude pure ghost jets from the rho calculation.
208  // Use only jets that fulfill your background jet selection
209  // for the rho calculation.
210  // Eg. real signal jets should not bias the background rho
211  if(jet->GetNumberOfTracks()>0)
212  {
213  rhovec[NjetAcc] = jet->Pt() / jet->Area();
214  ++NjetAcc;
215  }
216  }
217 
218  Double_t OccCorr=1;
219  //Use the total TPC area in which rho is calculated as denominater
220  if(fUseTPCArea)
221  {
222  OccCorr = TotaljetAreaPhys/TotalTPCArea;
223  }
224  //Use the total area covered by all ghost+physical jets that pass the selection in the detector as denominator
225  else if (TotalAreaCovered>0)
226  {
227  OccCorr = TotaljetAreaPhys/TotalAreaCovered;
228  }
229 
230  if (fCreateHisto)
231  fHistOccCorrvsCent->Fill(fCent, OccCorr);
232 
233  if (NjetAcc > 0) {
234  //find median value
235  Double_t rho = TMath::Median(NjetAcc, rhovec);
236 
237  if(fRhoCMS){
238  rho = rho * OccCorr;
239  }
240 
241  fOutRho->SetVal(rho);
242 
243  if (fOutRhoScaled) {
244  Double_t rhoScaled = rho * GetScaleFactor(fCent);
245  fOutRhoScaled->SetVal(rhoScaled);
246  }
247  }
248 
249  return kTRUE;
250 }
256  const char *nTracks,
257  const char *nClusters,
258  const char *nRho,
259  Double_t jetradius,
260  UInt_t acceptance,
263  const Bool_t histo,
264  const char *nJetsSig,
265  const char *cutType,
266  Double_t jetptcut,
267  Double_t jetareacut,
268  Double_t emcareacut,
269  const char *suffix
270 )
271 {
272 
273  // Get the pointer to the existing analysis manager via the static access method.
274  //==============================================================================
275  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
276  if (!mgr)
277  {
278  ::Error("AddTaskRhoSparse", "No analysis manager to connect to.");
279  return NULL;
280  }
281 
282  // Check the analysis type using the event handlers connected to the analysis manager.
283  //==============================================================================
284  AliVEventHandler* handler = mgr->GetInputEventHandler();
285  if (!handler)
286  {
287  ::Error("AddTaskEmcalJetPerformance", "This task requires an input event handler");
288  return 0;
289  }
290 
291  enum EDataType_t {
292  kUnknown,
293  kESD,
294  kAOD
295  };
296 
297  EDataType_t dataType = kUnknown;
298 
299  if (handler->InheritsFrom("AliESDInputHandler")) {
300  dataType = kESD;
301  }
302  else if (handler->InheritsFrom("AliAODInputHandler")) {
303  dataType = kAOD;
304  }
305  //-------------------------------------------------------
306  // Init the task and do settings
307  //-------------------------------------------------------
308  TString trackName(nTracks);
309  TString clusName(nClusters);
310 
311  if (trackName == "usedefault") {
312  if (dataType == kESD) {
313  trackName = "Tracks";
314  }
315  else if (dataType == kAOD) {
316  trackName = "tracks";
317  }
318  else {
319  trackName = "";
320  }
321  }
322 
323  if (clusName == "usedefault") {
324  if (dataType == kESD) {
325  clusName = "CaloClusters";
326  }
327  else if (dataType == kAOD) {
328  clusName = "caloClusters";
329  }
330  else {
331  clusName = "";
332  }
333  }
334 
335  TString name("AliAnalysisTaskRhoSparse");
336  if (strcmp(suffix,"") != 0) {
337  name += "_";
338  name += suffix;
339  }
340 
341  AliAnalysisTaskRhoSparse* mgrTask = (AliAnalysisTaskRhoSparse*)(mgr->GetTask(name.Data()));
342  if (mgrTask) return mgrTask;
343 
344  AliAnalysisTaskRhoSparse *rhotask = new AliAnalysisTaskRhoSparse(name, histo);
345  rhotask->SetHistoBins(1000,-0.1,9.9);
346  rhotask->SetRhoCMS(1);
347  rhotask->SetSmallSystem(1);
348  rhotask->SetOutRhoName(nRho);
349 
350  AliTrackContainer* trackCont;
351  AliParticleContainer* partCont;
352  if (trackName == "mcparticles") {
353  AliMCParticleContainer* mcpartCont = rhotask->AddMCParticleContainer(trackName);
354  mcpartCont->SelectPhysicalPrimaries(kTRUE);
355  }
356  else if (trackName == "tracks" || trackName == "Tracks") {
357  trackCont = rhotask->AddTrackContainer(trackName);
358  }
359  else if (!trackName.IsNull()) {
360  rhotask->AddParticleContainer(trackName);
361  }
362  partCont = rhotask->GetParticleContainer(0);
363 
364  AliClusterContainer *clusterCont = rhotask->AddClusterContainer(clusName);
365  if (clusterCont) {
366  clusterCont->SetClusECut(0.);
367  clusterCont->SetClusPtCut(0.);
368  clusterCont->SetClusHadCorrEnergyCut(0.3);
369  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
370  }
371 
372  //
373  AliJetContainer *bkgJetCont = rhotask->AddJetContainer(jetType, AliJetContainer::kt_algorithm, rscheme, jetradius, acceptance, partCont, clusterCont);
374  if (bkgJetCont) {
375  //why?? bkgJetCont->SetJetAreaCut(jetareacut);
376  //why?? bkgJetCont->SetAreaEmcCut(emcareacut);
377  bkgJetCont->SetJetPtCut(0.);
378  }
379 
380  //This is only for excluding overlaps with signal jets in case signal jets are provided
381  AliJetContainer *sigJetCont;
382  if(nJetsSig!=0)
383  {
384  sigJetCont = rhotask->AddJetContainer(nJetsSig,cutType,jetradius);
385  if (sigJetCont)
386  {
387  sigJetCont->SetJetAreaCut(jetareacut);
388  sigJetCont->SetAreaEmcCut(emcareacut);
389  sigJetCont->SetJetPtCut(jetptcut);
390  sigJetCont->ConnectParticleContainer(trackCont);
391  sigJetCont->ConnectClusterContainer(clusterCont);
392  }
393  }
394  //-------------------------------------------------------
395  // Final settings, pass to manager and set the containers
396  //-------------------------------------------------------
397  mgr->AddTask(rhotask);
398 
399  // Create containers for input/output
400  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
401  if (histo) {
402  TString contname(name);
403  contname += "_histos";
404  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
405  TList::Class(),AliAnalysisManager::kOutputContainer,
406  Form("%s", AliAnalysisManager::GetCommonFileName()));
407  mgr->ConnectOutput(rhotask, 1, coutput1);
408  }
409 
410  return rhotask;
411 }
412 
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.
static AliAnalysisTaskRhoSparse * AddTaskRhoSparse(const char *nTracks="usedefault", const char *nClusters="usedefault", const char *nRho="Rho", Double_t jetradius=0.2, UInt_t acceptance=AliEmcalJet::kTPC, 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="")
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Bool_t fUseTPCArea
use the full TPC area for the denominator of the occupancy calculation
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
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)