AliPhysics  2c6b7ad (2c6b7ad)
AliAnalysisTaskRhoDev.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2017, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include "AliAnalysisTaskRhoDev.h"
17 
18 #include <TClonesArray.h>
19 #include <TMath.h>
20 
21 #include <AliLog.h>
22 #include <AliVEventHandler.h>
23 #include <AliAnalysisManager.h>
24 
25 #include "AliEmcalJet.h"
26 #include "AliRhoParameter.h"
27 #include "AliJetContainer.h"
28 
30 ClassImp(AliAnalysisTaskRhoDev);
32 
38  fNExclLeadJets(0),
39  fRhoSparse(kFALSE),
40  fExclJetOverlap(),
41  fOccupancyFactor(0),
42  fHistOccCorrvsCent(nullptr)
43 {
44 }
45 
53  AliAnalysisTaskRhoBaseDev(name, histo),
54  fNExclLeadJets(0),
55  fRhoSparse(kFALSE),
59 {
60 }
61 
67 {
68  if (!fCreateHisto) return;
69 
71 
72  fHistOccCorrvsCent = new TH2F("fHistOccCorrvsCent", "fHistOccCorrvsCent;Centrality (%);#it{C}", 100, 0, 100, 2000, 0 , 2);
74 }
75 
80 std::pair<AliEmcalJet*, AliEmcalJet*> AliAnalysisTaskRhoDev::GetLeadingJets()
81 {
82  std::pair<AliEmcalJet*, AliEmcalJet*> maxJets = {nullptr, nullptr};
83  if (fNExclLeadJets <= 0) return maxJets;
84 
85  auto itJet = fSortedJets["Background"].begin();
86 
87  maxJets.first = *itJet;
88  if (fNExclLeadJets > 1) {
89  itJet++;
90  if (itJet != fSortedJets["Background"].end()) maxJets.second = *itJet;
91  }
92 
93  return maxJets;
94 }
95 
102 {
103  if (fJetCollArray.empty()) return;
104 
105  auto maxJets = GetLeadingJets();
106 
107  static Double_t rhovec[999];
108  Int_t NjetAcc = 0;
109  Double_t TotaljetArea = 0; // Total area of background jets (including ghost jets)
110  Double_t TotaljetAreaPhys = 0; // Total area of physical background jets (excluding ghost jets)
111  // Ghost jet is a jet made only of ghost particles
112 
113  AliJetContainer* bkgJetCont = fJetCollArray["Background"];
114  AliJetContainer* sigJetCont = nullptr;
115  if (!fExclJetOverlap.IsNull()) {
116  auto sigJetContIt = fJetCollArray.find(fExclJetOverlap.Data());
117  if (sigJetContIt != fJetCollArray.end()) sigJetCont = sigJetContIt->second;
118  }
119 
120  // push all jets within selected acceptance into stack
121  for (auto jet : bkgJetCont->accepted()) {
122 
123  TotaljetArea += jet->Area();
124 
125  if (jet->IsGhost()) continue;
126 
127  TotaljetAreaPhys += jet->Area();
128 
129  // excluding leading jets
130  if (jet == maxJets.first || jet == maxJets.second) continue;
131 
132  Bool_t overlapsWithSignal = kFALSE;
133  if (sigJetCont) {
134  for (auto sigJet : sigJetCont->accepted()) {
135  if (AreJetsOverlapping(jet, sigJet)) {
136  overlapsWithSignal = kTRUE;
137  break;
138  }
139  }
140  }
141 
142  if (overlapsWithSignal) continue;
143 
144  rhovec[NjetAcc] = jet->Pt() / jet->Area();
145  ++NjetAcc;
146  }
147 
148  // Occupancy correction for sparse event described in https://arxiv.org/abs/1207.2392
149  if (TotaljetArea > 0) {
150  fOccupancyFactor = TotaljetAreaPhys / TotaljetArea;
151  }
152  else {
153  fOccupancyFactor = 0;
154  }
155 
156  if (NjetAcc > 0) {
157  //find median value
158  Double_t rho = TMath::Median(NjetAcc, rhovec);
159 
160  if (fRhoSparse) rho = rho * fOccupancyFactor;
161 
162  fOutRho->SetVal(rho);
163  }
164 }
165 
170 {
172  if (!r) return kFALSE;
173 
175 
176  return kTRUE;
177 }
178 
184 {
185  if (fJetCollArray.count("Background") == 0) {
186  AliError("No signal jet collection found. Task will not run!");
187  return kFALSE;
188  }
189 
190  return kTRUE;
191 }
192 
209 {
210  // Get the pointer to the existing analysis manager via the static access method.
211  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
212  if (!mgr) {
213  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "No analysis manager to connect to.");
214  return nullptr;
215  }
216 
217  // Check the analysis type using the event handlers connected to the analysis manager.
218  AliVEventHandler* handler = mgr->GetInputEventHandler();
219  if (!handler) {
220  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "This task requires an input event handler");
221  return nullptr;
222  }
223 
225 
226  if (handler->InheritsFrom("AliESDInputHandler")) {
227  dataType = kESD;
228  }
229  else if (handler->InheritsFrom("AliAODInputHandler")) {
230  dataType = kAOD;
231  }
232 
233  // Init the task and do settings
234  if (trackName == "usedefault") {
235  if (dataType == kESD) {
236  trackName = "Tracks";
237  }
238  else if (dataType == kAOD) {
239  trackName = "tracks";
240  }
241  else {
242  trackName = "";
243  }
244  }
245 
246  if (clusName == "usedefault") {
247  if (dataType == kESD) {
248  clusName = "CaloClusters";
249  }
250  else if (dataType == kAOD) {
251  clusName = "caloClusters";
252  }
253  else {
254  clusName = "";
255  }
256  }
257 
258  TString name(TString::Format("AliAnalysisTaskRhoDev_%s", nRho.Data()));
259  if (!suffix.IsNull()) {
260  name += "_";
261  name += suffix;
262  }
263 
264  AliAnalysisTaskRhoDev* mgrTask = dynamic_cast<AliAnalysisTaskRhoDev*>(mgr->GetTask(name.Data()));
265  if (mgrTask) {
266  ::Warning("AliAnalysisTaskRhoDev::AddTaskRhoDev", "Not adding the task again, since a task with the same name '%s' already exists", name.Data());
267  return mgrTask;
268  }
269 
270  AliAnalysisTaskRhoDev* rhotask = new AliAnalysisTaskRhoDev(name, histo);
271  rhotask->SetOutRhoName(nRho);
272 
273  AliParticleContainer* partCont = rhotask->AddParticleContainer(trackName.Data());
274  partCont->SetMinPt(trackPtCut);
275  AliClusterContainer *clusterCont = rhotask->AddClusterContainer(clusName.Data());
276  if (clusterCont) {
277  clusterCont->SetClusECut(0.);
278  clusterCont->SetClusPtCut(0.);
279  clusterCont->SetClusHadCorrEnergyCut(clusECut);
280  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
281  }
282 
283  AliJetContainer *jetCont = new AliJetContainer(jetType, AliJetContainer::kt_algorithm, rscheme, jetradius, partCont, clusterCont);
284  if (jetCont) {
285  jetCont->SetJetPtCut(0);
286  jetCont->SetJetAcceptanceType(acceptance);
287  jetCont->SetName("Background");
288  rhotask->AdoptJetContainer(jetCont);
289  }
290 
291  // Final settings, pass to manager and set the containers
292  mgr->AddTask(rhotask);
293 
294  // Create containers for input/output
295  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
296  if (histo) {
297  TString contname(name);
298  contname += "_histos";
299  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
300  TList::Class(), AliAnalysisManager::kOutputContainer,
301  Form("%s", AliAnalysisManager::GetCommonFileName()));
302  mgr->ConnectOutput(rhotask, 1, coutput1);
303  }
304 
305  return rhotask;
306 }
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
Bool_t AreJetsOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
double Double_t
Definition: External.C:58
EDataType_t
Switch for the data type.
Definition: External.C:236
Declaration of class AliAnalysisTaskRhoDev.
Double_t fOccupancyFactor
!occupancy correction factor for sparse events
Base class for a task that calculates the UE.
TH2F * fHistOccCorrvsCent
!occupancy correction vs. centrality
Container for particles within the EMCAL framework.
Bool_t fCreateHisto
whether or not create histograms
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
void SetOutRhoName(const char *name)
unsigned int UInt_t
Definition: External.C:33
static AliAnalysisTaskRhoDev * AddTaskRhoDev(TString nTracks="usedefault", Double_t trackPtCut=0.15, TString nClusters="usedefault", Double_t clusECut=0.30, TString 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, Bool_t histo=kTRUE, TString suffix="")
Bool_t fRhoSparse
flag to run CMS method as described in https://arxiv.org/abs/1207.2392
TString fExclJetOverlap
name of the jet collection that should be used to reject jets that are considered "signal" ...
AliParticleContainer * AddParticleContainer(std::string branchName, std::string contName="")
Class for a task that calculates the UE.
void SetClusPtCut(Double_t cut)
void AdoptJetContainer(AliJetContainer *cont)
std::map< std::string, AliJetContainer * > fJetCollArray
jet collection array
void SetJetAcceptanceType(UInt_t type)
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
Container structure for EMCAL clusters.
AliRhoParameter * fOutRho
!output rho object
Double_t fCent
!event centrality
Container for jet within the EMCAL jet framework.
std::map< std::string, std::list< AliEmcalJet * > > fSortedJets
!jets sorted by momentum
void SetClusHadCorrEnergyCut(Double_t cut)
std::pair< AliEmcalJet *, AliEmcalJet * > GetLeadingJets()
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation