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