AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
32 
38  fNExclLeadJets(0),
39  fRhoSparse(kFALSE),
40  fOccupancyFactor(1.),
41  fHistOccCorrvsCent(nullptr)
42 {
43 }
44 
52  AliAnalysisTaskRhoBaseDev(name, histo),
53  fNExclLeadJets(0),
54  fRhoSparse(kFALSE),
55  fOccupancyFactor(1.),
56  fHistOccCorrvsCent(nullptr)
57 {
58 }
59 
65 {
66  if (!fCreateHisto) return;
67 
69 
70  fHistOccCorrvsCent = new TH2F("fHistOccCorrvsCent", "fHistOccCorrvsCent", 100, 0, 100, 2000, 0 , 2);
72 }
73 
78 std::pair<AliEmcalJet*, AliEmcalJet*> AliAnalysisTaskRhoDev::GetLeadingJets()
79 {
80  std::pair<AliEmcalJet*, AliEmcalJet*> maxJets = {nullptr, nullptr};
81  if (fNExclLeadJets <= 0) return maxJets;
82 
83  auto itJet = fSortedJets["Signal"].begin();
84 
85  maxJets.first = *itJet;
86  if (fNExclLeadJets > 1) {
87  itJet++;
88  if (itJet != fSortedJets["Signal"].end()) maxJets.second = *itJet;
89  }
90 
91  return maxJets;
92 }
93 
101 {
102  for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i) {
103  Int_t jet1Track = jet1->TrackAt(i);
104  for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j) {
105  Int_t jet2Track = jet2->TrackAt(j);
106  if (jet1Track == jet2Track) return kTRUE;
107  }
108  }
109  return kFALSE;
110 }
111 
118 {
119  if (fJetCollArray.empty()) return;
120 
121  auto maxJets = GetLeadingJets();
122 
123  static Double_t rhovec[999];
124  Int_t NjetAcc = 0;
125  Double_t TotaljetArea = 0; // Total area of background jets (including ghost jets)
126  Double_t TotaljetAreaPhys = 0; // Total area of physical background jets (excluding ghost jets)
127  // Ghost jet is a jet made only of ghost particles
128 
129  AliJetContainer* bkgJetCont = fJetCollArray["Background"];
130  AliJetContainer* sigJetCont = nullptr;
131  auto sigJetContIt = fJetCollArray.find("Signal");
132  if (sigJetContIt != fJetCollArray.end()) sigJetCont = sigJetContIt->second;
133 
134  // push all jets within selected acceptance into stack
135  for (auto jet : bkgJetCont->accepted()) {
136 
137  TotaljetArea += jet->Area();
138 
139  if (jet->IsGhost()) continue;
140 
141  TotaljetAreaPhys += jet->Area();
142 
143  // excluding leading jets
144  if (jet == maxJets.first || jet == maxJets.second) continue;
145 
146  Bool_t overlapsWithSignal = kFALSE;
147  if (sigJetCont) {
148  for (auto sigJet : sigJetCont->accepted()) {
149  if (IsJetOverlapping(jet, sigJet)) {
150  overlapsWithSignal = kTRUE;
151  break;
152  }
153  }
154  }
155 
156  if (overlapsWithSignal) continue;
157 
158  rhovec[NjetAcc] = jet->Pt() / jet->Area();
159  ++NjetAcc;
160  }
161 
162 
163  if (NjetAcc > 0) {
164  //find median value
165  Double_t rho = TMath::Median(NjetAcc, rhovec);
166 
167  // Occupancy correction for sparse event described in https://arxiv.org/abs/1207.2392
168  if (TotaljetArea > 0) fOccupancyFactor = TotaljetAreaPhys / TotaljetArea;
169 
170  if (fRhoSparse) rho = rho * fOccupancyFactor;
171 
172  fOutRho->SetVal(rho);
173  }
174 }
175 
180 {
182  if (!r) return kFALSE;
183 
185 
186  return kTRUE;
187 }
188 
194 {
195  if (fJetCollArray.count("Background") == 0) {
196  AliError("No signal jet collection found. Task will not run!");
197  return kFALSE;
198  }
199 
200  return kTRUE;
201 }
202 
219 {
220  // Get the pointer to the existing analysis manager via the static access method.
221  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
222  if (!mgr) {
223  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "No analysis manager to connect to.");
224  return nullptr;
225  }
226 
227  // Check the analysis type using the event handlers connected to the analysis manager.
228  AliVEventHandler* handler = mgr->GetInputEventHandler();
229  if (!handler) {
230  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "This task requires an input event handler");
231  return nullptr;
232  }
233 
234  EDataType_t dataType = kUnknownDataType;
235 
236  if (handler->InheritsFrom("AliESDInputHandler")) {
237  dataType = kESD;
238  }
239  else if (handler->InheritsFrom("AliAODInputHandler")) {
240  dataType = kAOD;
241  }
242 
243  // Init the task and do settings
244  if (trackName == "usedefault") {
245  if (dataType == kESD) {
246  trackName = "Tracks";
247  }
248  else if (dataType == kAOD) {
249  trackName = "tracks";
250  }
251  else {
252  trackName = "";
253  }
254  }
255 
256  if (clusName == "usedefault") {
257  if (dataType == kESD) {
258  clusName = "CaloClusters";
259  }
260  else if (dataType == kAOD) {
261  clusName = "caloClusters";
262  }
263  else {
264  clusName = "";
265  }
266  }
267 
268  TString name(TString::Format("AliAnalysisTaskRhoDev_%s", nRho.Data()));
269  if (!suffix.IsNull()) {
270  name += "_";
271  name += suffix;
272  }
273 
274  AliAnalysisTaskRhoDev* mgrTask = dynamic_cast<AliAnalysisTaskRhoDev*>(mgr->GetTask(name.Data()));
275  if (mgrTask) {
276  ::Warning("AliAnalysisTaskRhoDev::AddTaskRhoDev", "Not adding the task again, since a task with the same name '%s' already exists", name.Data());
277  return mgrTask;
278  }
279 
280  AliAnalysisTaskRhoDev* rhotask = new AliAnalysisTaskRhoDev(name, histo);
281  rhotask->SetOutRhoName(nRho);
282 
283  AliParticleContainer* partCont = rhotask->AddParticleContainer(trackName.Data());
284  partCont->SetMinPt(trackPtCut);
285  AliClusterContainer *clusterCont = rhotask->AddClusterContainer(clusName.Data());
286  if (clusterCont) {
287  clusterCont->SetClusECut(0.);
288  clusterCont->SetClusPtCut(0.);
289  clusterCont->SetClusHadCorrEnergyCut(clusECut);
290  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
291  }
292 
293  AliJetContainer *jetCont = new AliJetContainer(jetType, AliJetContainer::kt_algorithm, rscheme, jetradius, partCont, clusterCont);
294  if (jetCont) {
295  jetCont->SetJetPtCut(0);
296  jetCont->SetJetAcceptanceType(acceptance);
297  jetCont->SetName("Background");
298  rhotask->AdoptJetContainer(jetCont);
299  }
300 
301  // Final settings, pass to manager and set the containers
302  mgr->AddTask(rhotask);
303 
304  // Create containers for input/output
305  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
306  if (histo) {
307  TString contname(name);
308  contname += "_histos";
309  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
310  TList::Class(), AliAnalysisManager::kOutputContainer,
311  Form("%s", AliAnalysisManager::GetCommonFileName()));
312  mgr->ConnectOutput(rhotask, 1, coutput1);
313  }
314 
315  return rhotask;
316 }
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
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
Bool_t IsJetOverlapping(AliEmcalJet *jet1, AliEmcalJet *jet2)
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
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
AliParticleContainer * AddParticleContainer(std::string branchName, std::string contName="")
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
std::map< std::string, std::list< AliEmcalJet * > > fSortedJets
!jets sorted by momentum
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
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.
void SetClusHadCorrEnergyCut(Double_t cut)
std::pair< AliEmcalJet *, AliEmcalJet * > GetLeadingJets()
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation