AliPhysics  master (3d17d9d)
AliAnalysisTaskRhoDev.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
27 #include "AliAnalysisTaskRhoDev.h"
28 
29 #include <TClonesArray.h>
30 #include <TMath.h>
31 
32 #include <AliLog.h>
33 #include <AliVEventHandler.h>
34 #include <AliAnalysisManager.h>
35 
36 #include "AliEmcalJet.h"
37 #include "AliRhoParameter.h"
38 #include "AliJetContainer.h"
39 
40 ClassImp(AliAnalysisTaskRhoDev);
41 
44  fNExclLeadJets(0),
45  fRhoSparse(kFALSE),
46  fExclJetOverlap(),
47  fOccupancyFactor(0),
48  fHistOccCorrvsCent(nullptr)
49 {
50 }
51 
53  AliAnalysisTaskRhoBaseDev(name, histo),
54  fNExclLeadJets(0),
55  fRhoSparse(kFALSE),
59 {
60 }
61 
63 {
64  if (!fCreateHisto) return;
65 
67 
68  fHistOccCorrvsCent = new TH2F("fHistOccCorrvsCent", "fHistOccCorrvsCent;Centrality (%);#it{C}", 100, 0, 100, 2000, 0 , 2);
70 }
71 
72 std::pair<AliEmcalJet*, AliEmcalJet*> AliAnalysisTaskRhoDev::GetLeadingJets()
73 {
74  std::pair<AliEmcalJet*, AliEmcalJet*> maxJets = {nullptr, nullptr};
75  if (fNExclLeadJets <= 0) return maxJets;
76 
77  auto itJet = fSortedJets["Background"].begin();
78 
79  maxJets.first = *itJet;
80  if (fNExclLeadJets > 1) {
81  itJet++;
82  if (itJet != fSortedJets["Background"].end()) maxJets.second = *itJet;
83  }
84 
85  return maxJets;
86 }
87 
89 {
90  if (fJetCollArray.empty()) return;
91 
92  auto maxJets = GetLeadingJets();
93 
94  static Double_t rhovec[999];
95  Int_t NjetAcc = 0;
96  Double_t TotaljetArea = 0; // Total area of background jets (including ghost jets)
97  Double_t TotaljetAreaPhys = 0; // Total area of physical background jets (excluding ghost jets)
98  // Ghost jet is a jet made only of ghost particles
99 
100  AliJetContainer* bkgJetCont = fJetCollArray["Background"];
101  AliJetContainer* sigJetCont = nullptr;
102  if (!fExclJetOverlap.IsNull()) {
103  auto sigJetContIt = fJetCollArray.find(fExclJetOverlap.Data());
104  if (sigJetContIt != fJetCollArray.end()) sigJetCont = sigJetContIt->second;
105  }
106 
107  // push all jets within selected acceptance into stack
108  for (auto jet : bkgJetCont->accepted()) {
109 
110  TotaljetArea += jet->Area();
111 
112  if (jet->IsGhost()) continue;
113 
114  TotaljetAreaPhys += jet->Area();
115 
116  // excluding leading jets
117  if (jet == maxJets.first || jet == maxJets.second) continue;
118 
119  Bool_t overlapsWithSignal = kFALSE;
120  if (sigJetCont) {
121  for (auto sigJet : sigJetCont->accepted()) {
122  if (AreJetsOverlapping(jet, sigJet)) {
123  overlapsWithSignal = kTRUE;
124  break;
125  }
126  }
127  }
128 
129  if (overlapsWithSignal) continue;
130 
131  rhovec[NjetAcc] = jet->Pt() / jet->Area();
132  ++NjetAcc;
133  }
134 
135  // Occupancy correction for sparse event described in https://arxiv.org/abs/1207.2392
136  if (TotaljetArea > 0) {
137  fOccupancyFactor = TotaljetAreaPhys / TotaljetArea;
138  }
139  else {
140  fOccupancyFactor = 0;
141  }
142 
143  if (NjetAcc > 0) {
144  //find median value
145  Double_t rho = TMath::Median(NjetAcc, rhovec);
146 
147  if (fRhoSparse) rho = rho * fOccupancyFactor;
148 
149  fOutRho->SetVal(rho);
150  }
151 }
152 
154 {
156  if (!r) return kFALSE;
157 
159 
160  return kTRUE;
161 }
162 
164 {
165  if (fJetCollArray.count("Background") == 0) {
166  AliError("No signal jet collection found. Task will not run!");
167  return kFALSE;
168  }
169 
170  return kTRUE;
171 }
172 
174 {
175  // Get the pointer to the existing analysis manager via the static access method.
176  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
177  if (!mgr) {
178  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "No analysis manager to connect to.");
179  return nullptr;
180  }
181 
182  // Check the analysis type using the event handlers connected to the analysis manager.
183  AliVEventHandler* handler = mgr->GetInputEventHandler();
184  if (!handler) {
185  ::Error("AliAnalysisTaskRhoDev::AddTaskRhoDev", "This task requires an input event handler");
186  return nullptr;
187  }
188 
190 
191  if (handler->InheritsFrom("AliESDInputHandler")) {
192  dataType = kESD;
193  }
194  else if (handler->InheritsFrom("AliAODInputHandler")) {
195  dataType = kAOD;
196  }
197 
198  // Init the task and do settings
199  if (trackName == "usedefault") {
200  if (dataType == kESD) {
201  trackName = "Tracks";
202  }
203  else if (dataType == kAOD) {
204  trackName = "tracks";
205  }
206  else {
207  trackName = "";
208  }
209  }
210 
211  if (clusName == "usedefault") {
212  if (dataType == kESD) {
213  clusName = "CaloClusters";
214  }
215  else if (dataType == kAOD) {
216  clusName = "caloClusters";
217  }
218  else {
219  clusName = "";
220  }
221  }
222 
223  TString name(TString::Format("AliAnalysisTaskRhoDev_%s", nRho.Data()));
224  if (!suffix.IsNull()) {
225  name += "_";
226  name += suffix;
227  }
228 
229  AliAnalysisTaskRhoDev* mgrTask = dynamic_cast<AliAnalysisTaskRhoDev*>(mgr->GetTask(name.Data()));
230  if (mgrTask) {
231  ::Warning("AliAnalysisTaskRhoDev::AddTaskRhoDev", "Not adding the task again, since a task with the same name '%s' already exists", name.Data());
232  return mgrTask;
233  }
234 
235  AliAnalysisTaskRhoDev* rhotask = new AliAnalysisTaskRhoDev(name, histo);
236  rhotask->SetOutRhoName(nRho);
237 
238  AliParticleContainer* partCont = rhotask->AddParticleContainer(trackName.Data());
239  partCont->SetMinPt(trackPtCut);
240  AliClusterContainer *clusterCont = rhotask->AddClusterContainer(clusName.Data());
241  if (clusterCont) {
242  clusterCont->SetClusECut(0.);
243  clusterCont->SetClusPtCut(0.);
244  clusterCont->SetClusHadCorrEnergyCut(clusECut);
245  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
246  }
247 
248  AliJetContainer *jetCont = new AliJetContainer(jetType, AliJetContainer::kt_algorithm, rscheme, jetradius, partCont, clusterCont);
249  if (jetCont) {
250  jetCont->SetJetPtCut(0);
251  jetCont->SetJetAcceptanceType(acceptance);
252  jetCont->SetName("Background");
253  rhotask->AdoptJetContainer(jetCont);
254  }
255 
256  // Final settings, pass to manager and set the containers
257  mgr->AddTask(rhotask);
258 
259  // Create containers for input/output
260  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
261  if (histo) {
262  TString contname(name);
263  contname += "_histos";
264  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
265  TList::Class(), AliAnalysisManager::kOutputContainer,
266  Form("%s", AliAnalysisManager::GetCommonFileName()));
267  mgr->ConnectOutput(rhotask, 1, coutput1);
268  }
269 
270  return rhotask;
271 }
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
void SetName(const char *n)
Set the name of the class of the objets inside the underlying array.
Double_t fOccupancyFactor
!occupancy correction factor for sparse events
Base class for a task that calculates the UE.
AliClusterContainer * AddClusterContainer(EMCAL_STRINGVIEW branchName, EMCAL_STRINGVIEW contName="")
void SetMinPt(Double_t min)
TH2F * fHistOccCorrvsCent
!occupancy correction vs. centrality
Container for particles within the EMCAL framework.
AliParticleContainer * AddParticleContainer(EMCAL_STRINGVIEW branchName, EMCAL_STRINGVIEW contName="")
AliAnalysisTaskRhoDev()
Default constructor. Needed by ROOT I/O.
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="")
Create an instance of this class and add it to the analysis manager.
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" ...
Bool_t VerifyContainers()
Verify that the required particle, cluster and jet containers were provided.
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.
Bool_t FillHistograms()
Fill histograms.
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 UserCreateOutputObjects()
Creating user output.
void SetClusHadCorrEnergyCut(Double_t cut)
std::pair< AliEmcalJet *, AliEmcalJet * > GetLeadingJets()
UInt_t fNExclLeadJets
number of leading jets to be excluded from the median calculation