AliPhysics  ee192bd (ee192bd)
AliAnalysisTaskEmcalJetSubstructureTree.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 <algorithm>
28 #include <array>
29 #include <bitset>
30 #include <iostream>
31 #include <string>
32 #include <set>
33 #include <sstream>
34 #include <vector>
35 
36 #include <fastjet/ClusterSequence.hh>
37 #include <fastjet/contrib/Nsubjettiness.hh>
38 #include <fastjet/contrib/SoftDrop.hh>
39 
40 #include <THistManager.h>
41 #include <TLinearBinning.h>
42 #include <TLorentzVector.h>
43 #include <TMath.h>
44 #include <TObjString.h>
45 #include <TString.h>
46 #include <TVector3.h>
47 
48 #include "AliAODEvent.h"
49 #include "AliAODInputHandler.h"
50 #include "AliAnalysisManager.h"
51 #include "AliAnalysisDataSlot.h"
52 #include "AliAnalysisDataContainer.h"
54 #include "AliCDBEntry.h"
55 #include "AliCDBManager.h"
56 #include "AliClusterContainer.h"
57 #include "AliJetContainer.h"
60 #include "AliEmcalJet.h"
61 #include "AliEmcalList.h"
64 #include "AliLog.h"
65 #include "AliParticleContainer.h"
66 #include "AliRhoParameter.h"
67 #include "AliTrackContainer.h"
68 #include "AliTriggerCluster.h"
69 #include "AliTriggerConfiguration.h"
70 #include "AliVCluster.h"
71 #include "AliVParticle.h"
72 
73 #ifdef EXPERIMENTAL_JETCONSTITUENTS
76 #endif
77 
78 
82 
83 namespace EmcalTriggerJets {
84 
87  fJetSubstructureTree(nullptr),
88  fGlobalTreeParams(nullptr),
89  fSoftDropMeasured(nullptr),
90  fSoftDropTrue(nullptr),
91  fNSubMeasured(nullptr),
92  fNSubTrue(nullptr),
93  fKineRec(nullptr),
94  fKineSim(nullptr),
95  fJetStructureMeasured(nullptr),
96  fJetStructureTrue(nullptr),
97  fQAHistos(nullptr),
98  fLumiMonitor(nullptr),
99  fSDZCut(0.1),
100  fSDBetaCut(0),
101  fReclusterizer(kCAAlgo),
102  fHasRecEvent(false),
103  fHasTrueEvent(false),
104  fTriggerSelectionBits(AliVEvent::kAny),
105  fTriggerSelectionString(""),
106  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
107  fUseTriggerSelectionForData(false),
108  fUseDownscaleWeight(false),
109  fUseChargedConstituents(true),
110  fUseNeutralConstituents(true),
111  fFillPart(true),
112  fFillRho(true),
113  fFillSoftDrop(true),
114  fFillNSub(true),
115  fFillStructGlob(true)
116 {
117 }
118 
120  AliAnalysisTaskEmcalJet(name, kTRUE),
127  fKineRec(nullptr),
128  fKineSim(nullptr),
133  fSDZCut(0.1),
134  fSDBetaCut(0),
136  fHasRecEvent(false),
137  fHasTrueEvent(false),
138  fTriggerSelectionBits(AliVEvent::kAny),
140  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
142  fUseDownscaleWeight(false),
145  fFillPart(true),
146  fFillRho(true),
147  fFillSoftDrop(true),
148  fFillNSub(true),
149  fFillStructGlob(true)
150 {
151  DefineOutput(2, TTree::Class());
152  SetUseAliAnaUtils(true);
153 }
154 
158  if(fSoftDropTrue) delete fSoftDropTrue;
159  if(fNSubMeasured) delete fNSubMeasured;
160  if(fNSubTrue) delete fNSubTrue;
161  if(fKineRec) delete fKineRec;
162  if(fKineSim) delete fKineSim;
165 }
166 
169 
170  // Make QA for constituent clusters
171  TLinearBinning jetptbinning(9, 20, 200),
172  clusterenergybinning(200, 0., 200),
173  timebinning(1000, -500., 500.),
174  m02binning(100, 0., 1.),
175  ncellbinning(101, -0.5, 100.5);
176  fQAHistos = new THistManager("QAhistos");
177  fQAHistos->CreateTH1("hEventCounter", "Event counter", 1, 0.5, 1.5);
178  fQAHistos->CreateTH2("hClusterConstE", "EMCAL cluster energy vs jet pt; p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
179  fQAHistos->CreateTH2("hClusterConstTime", "EMCAL cluster time vs. jet pt; p_{t, jet} (GeV/c); t_{cl} (ns)", jetptbinning, timebinning);
180  fQAHistos->CreateTH2("hClusterConstM02", "EMCAL cluster M02 vs. jet pt; p{t, jet} (GeV/c); M02", jetptbinning, m02binning);
181  fQAHistos->CreateTH2("hClusterConstNcell", "EMCAL cluster ncell vs. jet pt; p{t, jet} (GeV/c); Number of cells", jetptbinning, ncellbinning);
182 
183  // Test of constituent QA
184 #ifdef EXPERIMENTAL_JETCONSTITUENTS
185  fQAHistos->CreateTH2("hChargedConstituentPt", "charged constituent pt vs jet pt (via constituent map); p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
186  fQAHistos->CreateTH2("hChargedIndexPt", "charged constituent pt vs jet pt (via index map); p_{t, jet} (GeV/c); p_{t, ch} (GeV/c)", jetptbinning, clusterenergybinning);
187 
188  fQAHistos->CreateTH2("hClusterConstituentEDefault", "cluster constituent default energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
189  fQAHistos->CreateTH2("hClusterConstituentENLC", "cluster constituent non-linearity-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
190  fQAHistos->CreateTH2("hClusterConstituentEHC", "cluster constituent hadronic-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
191  fQAHistos->CreateTH2("hClusterIndexENLC", "cluster constituent non-linearity-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
192  fQAHistos->CreateTH2("hClusterIndexEHC", "cluster constituent hadronic-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
193  fQAHistos->CreateTH2("hLeadingChargedConstituentPt", "Pt of the leading charged constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
194  fQAHistos->CreateTH2("hLeadingClusterConstituentPt", "Pt of the leading cluster constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
195 #endif
196  for(auto h : *(fQAHistos->GetListOfHistograms())) fOutput->Add(h);
197 
198  OpenFile(2);
199  std::string treename = this->GetOutputSlot(2)->GetContainer()->GetName();
200  fJetSubstructureTree = new TTree(treename.data(), "Tree with jet substructure information");
201 
203  fGlobalTreeParams->LinkJetTreeBranches(fJetSubstructureTree, fFillRho);
206  if(fFillPart) {
209  }
210  if(fFillSoftDrop) {
213  if(fFillPart) {
216  }
217  }
218  if(fFillNSub) {
221  if(fFillPart) {
224  }
225  }
226 
227  if(fFillStructGlob){
230  if(fFillPart){
233  }
234  }
235 
236  PostData(1, fOutput);
237  PostData(2, fJetSubstructureTree);
238 }
239 
243  }
244 }
245 
249  AliParticleContainer *particles = GetParticleContainer("mcparticles");
250 
251  AliJetContainer *mcjets = GetJetContainer("mcjets");
252  AliJetContainer *datajets = GetJetContainer("datajets");
253 
254  FillLuminosity(); // Makes only sense in data
255 
256  // for(auto e : *(fInputEvent->GetList())) std::cout << e->GetName() << std::endl;
257 
258  std::stringstream rhoTagData, rhoTagMC;
259  if(datajets) rhoTagData << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(datajets->GetJetRadius() * 10.);
260  if(mcjets) rhoTagMC << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(mcjets->GetJetRadius() * 10.);
261 
262  if(fFillRho){
263  std::string rhoSparseData = "RhoSparse_Full_" + rhoTagData.str(), rhoSparseMC = "RhoSparse_Full_" + rhoTagMC.str(),
264  rhoMassData = "RhoMassSparse_Full_" + rhoTagData.str(), rhoMassMC = "RhoMassSparse_Full_" + rhoTagMC.str();
265  AliRhoParameter *rhoPtRec = GetRhoFromEvent(rhoSparseData.data()),
266  *rhoMassRec = GetRhoFromEvent(rhoMassData.data()),
267  *rhoPtSim = GetRhoFromEvent(rhoSparseMC.data()),
268  *rhoMassSim = GetRhoFromEvent(rhoMassMC.data());
269  AliDebugStream(2) << "Found rho parameter for reconstructed pt: " << (rhoPtRec ? "yes" : "no") << ", value: " << (rhoPtRec ? rhoPtRec->GetVal() : 0.) << std::endl;
270  AliDebugStream(2) << "Found rho parameter for sim pt: " << (rhoPtSim ? "yes" : "no") << ", value: " << (rhoPtSim ? rhoPtSim->GetVal() : 0.) << std::endl;
271  AliDebugStream(2) << "Found rho parameter for reconstructed Mass: " << (rhoMassRec ? "yes" : "no") << ", value: " << (rhoMassRec ? rhoMassRec->GetVal() : 0.) << std::endl;
272  AliDebugStream(2) << "Found rho parameter for sim Mass: " << (rhoMassSim ? "yes" : "no") << ", value: " << (rhoMassSim ? rhoMassSim->GetVal() : 0.) << std::endl;
273  Double_t rhopars[4] = {
274  rhoPtRec ? rhoPtRec->GetVal() : 0.,
275  rhoPtSim ? rhoPtSim->GetVal() : 0.,
276  rhoMassRec ? rhoMassRec->GetVal() : 0.,
277  rhoMassSim ? rhoMassSim->GetVal() : 0.
278  };
279  memcpy(this->fGlobalTreeParams->fRhoParamters, rhopars, sizeof(Double_t) * 4);
280  }
281 
282  AliDebugStream(1) << "Inspecting jet radius " << (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius()) << std::endl;
283  this->fGlobalTreeParams->fJetRadius = (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius());
284  fGlobalTreeParams->fTriggerClusterIndex = -1; // Reset trigger cluster index
285 
286  if(datajets && !mcjets){
287  // decode trigger string in order to determine the trigger clusters
288  std::vector<std::string> clusternames;
289  auto triggerinfos = PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data());
290  for(auto t : triggerinfos) {
291  if(std::find(clusternames.begin(), clusternames.end(), t.Triggercluster()) == clusternames.end()) clusternames.emplace_back(t.Triggercluster());
292  }
293  bool isCENT = (std::find(clusternames.begin(), clusternames.end(), "CENT") != clusternames.end()),
294  isCENTNOTRD = (std::find(clusternames.begin(), clusternames.end(), "CENTNOTRD") != clusternames.end()),
295  isCALO = (std::find(clusternames.begin(), clusternames.end(), "CALO") != clusternames.end()),
296  isCALOFAST = (std::find(clusternames.begin(), clusternames.end(), "CALOFAST") != clusternames.end());
297  if(isCENT) fGlobalTreeParams->fTriggerClusterIndex = 0;
298  else if(isCENTNOTRD) fGlobalTreeParams->fTriggerClusterIndex = 1;
299  else if(isCALO) fGlobalTreeParams->fTriggerClusterIndex = 2;
300  else if(isCALOFAST) fGlobalTreeParams->fTriggerClusterIndex = 3;
301  }
302 
303  double weight = 1.;
305  AliDebugStream(2) << "Trigger selection string: " << fTriggerSelectionString << std::endl;
306  TString selectionString = (fTriggerSelectionBits & AliVEvent::kINT7) ? "INT7" : fTriggerSelectionString;
307  auto triggerstring = MatchTrigger(selectionString.Data());
308  AliDebugStream(2) << "Getting downscale correction factor for trigger string " << triggerstring << std::endl;
310  }
311  AliDebugStream(1) << "Using downscale weight " << weight << std::endl;
312  this->fGlobalTreeParams->fEventWeight = weight;
313 
314 
315  // Count events (for spectrum analysis)
316  fQAHistos->FillTH1("hEventCounter", 1);
317 
318  AliSoftdropDefinition softdropSettings;
319  softdropSettings.fBeta = fSDBetaCut;
320  softdropSettings.fZ = fSDZCut;
321  switch(fReclusterizer) {
322  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
323  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
324  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
325  };
326 
327  AliNSubjettinessDefinition nsubjettinessSettings;
328  nsubjettinessSettings.fBeta = 1.;
329  nsubjettinessSettings.fRadius = 0.4;
330 
331  if(datajets) {
332  AliDebugStream(1) << "In data jets branch: found " << datajets->GetNJets() << " jets, " << datajets->GetNAcceptedJets() << " were accepted\n";
333  AliDebugStream(1) << "Having MC information: " << (mcjets ? TString::Format("yes, with %d jets", mcjets->GetNJets()) : "no") << std::endl;
334  if(mcjets) {
335  AliDebugStream(1) << "In MC jets branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
336  }
337  for(auto jet : datajets->accepted()) {
338  double pt = jet->Pt(), pz = jet->Pz(), E = jet->E(), M = TMath::Sqrt(E*E - pt*pt - pz*pz);
339  AliDebugStream(2) << "Next jet: pt:" << jet->Pt() << ", E: " << E << ", pz: " << pz << ", M(self): " << M << "M(fj)" << jet->M() << std::endl;
340  AliEmcalJet *associatedJet = jet->ClosestJet();
341 
342  if(mcjets) {
343  if(!associatedJet) {
344  AliDebugStream(2) << "Not found associated jet" << std::endl;
345  continue;
346  }
347  if(!(SelectJet(*jet, tracks) && SelectJet(*associatedJet, particles))) continue;
348  try {
349  DoConstituentQA(jet, tracks, clusters);
350  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, datajets->GetJetRadius() * 2., tracks, clusters, {softdropSettings, nsubjettinessSettings}),
351  structureMC = fFillPart ? MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings}) : AliJetSubstructureData();
353  if(fKineSim) *fKineSim = MakeJetKineParameters(*associatedJet);
354  if(fSoftDropMeasured) *fSoftDropMeasured = structureData.fSoftDrop;
355  if(fSoftDropTrue) *fSoftDropTrue = structureMC.fSoftDrop;
356  if(fNSubMeasured) *fNSubMeasured = structureData.fNsubjettiness;
357  if(fNSubTrue) *fNSubTrue = structureMC.fNsubjettiness;
358  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
359  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*associatedJet, particles, nullptr), MakePtD(*associatedJet, particles, nullptr)};
360  fJetSubstructureTree->Fill();
361  } catch(ReclusterizerException &e) {
362  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
363  } catch(SubstructureException &e) {
364  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
365  }
366  } else {
367  if(!SelectJet(*jet, tracks)) continue;
368  try {
369  DoConstituentQA(jet, tracks, clusters);
370  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
373  if(fNSubMeasured) *fNSubMeasured = structure.fNsubjettiness;
374  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
375  fJetSubstructureTree->Fill();
376  } catch(ReclusterizerException &e) {
377  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
378  } catch(SubstructureException &e) {
379  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
380  }
381  }
382  }
383  } else {
384  if(mcjets) {
385  // for MCgen train
386  AliDebugStream(1) << "In MC pure jet branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
387  for(auto j : mcjets->accepted()){
388  AliEmcalJet *mcjet = static_cast<AliEmcalJet *>(j);
389  try {
390  AliJetSubstructureData structure = MakeJetSubstructure(*mcjet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
391  if(this->fKineSim) *fKineSim = MakeJetKineParameters(*mcjet);
392  if(fSoftDropTrue) *fSoftDropTrue = structure.fSoftDrop;
393  if(fNSubTrue) *fNSubTrue = structure.fNsubjettiness;
394  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*mcjet, particles, nullptr), MakePtD(*mcjet, particles, nullptr)};
395  fJetSubstructureTree->Fill();
396  } catch (ReclusterizerException &e) {
397  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
398  } catch (SubstructureException &e) {
399  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
400  }
401  }
402  }
403  }
404 
405  return true;
406 }
407 
409  AliDebugStream(1) << "Trigger selection called\n";
410  if(!(fHasRecEvent || fHasTrueEvent)){
411  AliErrorStream() << "Impossible combination: Neither rec nor true event available. Rejecting ..." << std::endl;
412  return false;
413  }
414  // Run trigger selection (not on pure MCgen train - pure MCgen train has no rec event, ESD event is fake there)
415  if(fHasRecEvent){
416  if(!fHasTrueEvent){
417  // Pure data - do EMCAL trigger selection from selection string
418  AliDebugStream(1) << "Applying trigger selection for trigger bits " << std::bitset<sizeof(decltype(fTriggerSelectionBits)) * 8>(fTriggerSelectionBits) << "and trigger selection string " << fTriggerSelectionString << std::endl;
419  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
420  AliDebugStream(1) << "Passed trigger bit selection" << std::endl;
421  if(fTriggerSelectionString.Length()) {
422  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
423  AliDebugStream(1) << "Passed trigger string section" << std::endl;
425  auto trgselresult = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
426  AliDebugStream(1) << "Found trigger decision object: " << (trgselresult ? "yes" : "no") << std::endl;
427  if(!trgselresult){
428  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
429  return false;
430  }
431  if(!trgselresult->IsEventSelected(fTriggerSelectionString)) return false;
432  AliDebugStream(1) << "Data event selected" << std::endl;
433  }
434  }
435  } else {
436  // Simulation - do EMCAL trigger selection from trigger selection object
437  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false; // Require INT7 trigger - EMCAL triggers will be a subset
439  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
440  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
441  if(!mctrigger){
442  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
443  return false;
444  }
445  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
446  }
447  }
448  }
449  return true; // trigger selected or pure MCgen information
450 }
451 
453  AliCDBManager * cdb = AliCDBManager::Instance();
454  if(!fMCEvent && cdb){
455  // Get List of trigger clusters
456  AliCDBEntry *en = cdb->Get("GRP/CTP/Config");
457  AliTriggerConfiguration *trg = static_cast<AliTriggerConfiguration *>(en->GetObject());
458  std::vector<std::string> clusternames;
459  for(auto c : trg->GetClusters()) {
460  AliTriggerCluster *clust = static_cast<AliTriggerCluster *>(c);
461  std::string clustname = clust->GetName();
462  auto iscent = clustname.find("CENT") != std::string::npos, iscalo = clustname.find("CALO") != std::string::npos;
463  if(!(iscalo || iscent)) continue;
464  clusternames.emplace_back(clustname);
465  }
466 
467  // Set the x-axis of the luminosity monitor histogram
468  fLumiMonitor = new TH1F("hLumiMonitor", "Luminosity monitor", clusternames.size(), 0, clusternames.size());
469  int currentbin(1);
470  for(auto c : clusternames) {
471  fLumiMonitor->GetXaxis()->SetBinLabel(currentbin++, c.data());
472  }
473  fOutput->Add(fLumiMonitor);
474  }
475 }
476 
479  auto downscalefactors = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance();
480  if(fInputEvent->GetFiredTriggerClasses().Contains("INT7")) {
481  for(auto trigger : PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data())){
482  auto int7trigger = trigger.IsTriggerClass("INT7");
483  auto bunchcrossing = trigger.BunchCrossing() == "B";
484  auto nopf = trigger.PastFutureProtection() == "NOPF";
485  AliDebugStream(4) << "Full name: " << trigger.ExpandClassName() << ", INT7 trigger: " << (int7trigger ? "Yes" : "No") << ", bunch crossing: " << (bunchcrossing ? "Yes" : "No") << ", no past-future protection: " << (nopf ? "Yes" : "No") << ", Cluster: " << trigger.Triggercluster() << std::endl;
486  if(int7trigger && bunchcrossing && nopf) {
487  double downscale = downscalefactors->GetDownscaleFactorForTriggerClass(trigger.ExpandClassName());
488  AliDebugStream(5) << "Using downscale " << downscale << std::endl;
489  fLumiMonitor->Fill(trigger.Triggercluster().data(), 1./downscale);
490  }
491  }
492  }
493  }
494 }
495 
497  AliJetKineParameters result;
498  result.fPt = TMath::Abs(jet.Pt());
499  result.fE = jet.E();
500  result.fEta = jet.Eta();
501  result.fPhi = jet.Phi();
502  result.fArea = jet.Area();
503  result.fMass = jet.M();
504  result.fNEF = jet.NEF();
505  result.fNCharged = jet.GetNumberOfTracks();
506  result.fNNeutral = jet.GetNumberOfClusters();
507  return result;
508 }
509 
511  const int kClusterOffset = 30000; // In order to handle tracks and clusters in the same index space the cluster index needs and offset, large enough so that there is no overlap with track indices
512  std::vector<fastjet::PseudoJet> constituents;
513  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
514  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
515  if(tracks && (fUseChargedConstituents || isMC)){ // Neutral particles part of particle container in case of MC
516  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
517  auto track = jet.TrackAt(itrk, tracks->GetArray());
518  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
519  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
520  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
521  constituentTrack.set_user_index(jet.TrackAt(itrk));
522  constituents.push_back(constituentTrack);
523  }
524  }
525 
526  if(clusters && fUseNeutralConstituents){
527  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
528  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
529  TLorentzVector clustervec;
530  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
531  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
532  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
533  constituents.push_back(constituentCluster);
534  }
535  }
536 
537  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
538  if(!constituents.size()){
539  AliErrorStream() << "Jet has 0 constituents." << std::endl;
540  throw ReclusterizerException();
541  }
542  // Redo jet finding on constituents with a
543  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
544  std::vector<fastjet::PseudoJet> outputjets;
545  try {
546  fastjet::ClusterSequence jetfinder(constituents, jetdef);
547  outputjets = jetfinder.inclusive_jets(0);
549  return result;
550  } catch (fastjet::Error &e) {
551  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
552  throw ReclusterizerException();
553  } catch (SoftDropException &e) {
554  AliErrorStream() << "Softdrop exception caught: " << e.what() << std::endl;
555  throw ReclusterizerException();
556  }
557 }
558 
560  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
561  softdropAlgorithm.set_verbose_structure(kTRUE);
562  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
563  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
564  AliDebugStream(4) << "Jet has " << jet.constituents().size() << " constituents" << std::endl;
565  auto groomed = softdropAlgorithm(jet);
566  try {
567  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
568 
569  AliSoftDropParameters result({softdropstruct.symmetry(),
570  groomed.m(),
571  softdropstruct.delta_R(),
572  groomed.perp(),
573  softdropstruct.delta_R(),
574  softdropstruct.mu(),
575  softdropstruct.dropped_count()});
576  return result;
577  } catch(std::bad_cast &e) {
578  throw SoftDropException();
579  }
580 }
581 
584  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet),
585  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet)
586  });
587  return result;
588 }
589 
591  if(!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
592  throw SubstructureException();
593  TVector3 jetvec(jet.Px(), jet.Py(), jet.Pz());
594  Double_t den(0.), num(0.);
595  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
596  if(tracks && (fUseChargedConstituents || isMC)){
597  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
598  auto track = jet.TrackAt(itrk, tracks->GetArray());
599  if(!track){
600  AliErrorStream() << "Associated constituent particle / track not found\n";
601  continue;
602  }
603  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
604  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
605  TVector3 trackvec(track->Px(), track->Py(), track->Pz());
606 
607  num += track->Pt() * trackvec.DrEtaPhi(jetvec);
608  den += +track->Pt();
609  }
610  }
611  if(clusters && fUseNeutralConstituents) {
612  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
613  auto clust = jet.ClusterAt(icl, clusters->GetArray());
614  if(!clust) {
615  AliErrorStream() << "Associated constituent cluster not found\n";
616  continue;
617  }
618  TLorentzVector clusterp;
619  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
620 
621  num += clusterp.Pt() * clusterp.Vect().DrEtaPhi(jetvec);
622  den += clusterp.Pt();
623  }
624  }
625  return num/den;
626 }
627 
629  if (!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
630  throw SubstructureException();
631  Double_t den(0.), num(0.);
632  bool isMC = dynamic_cast<const AliMCParticleContainer *>(particles);
633  if(particles && (fUseChargedConstituents || isMC)){
634  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
635  auto trk = jet.TrackAt(itrk, particles->GetArray());
636  if(!trk){
637  AliErrorStream() << "Associated constituent particle / track not found\n";
638  continue;
639  }
640  if(!trk->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
641  if(trk->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
642  num += trk->Pt() * trk->Pt();
643  den += trk->Pt();
644  }
645  }
646  if(clusters && fUseNeutralConstituents){
647  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
648  auto clust = jet.ClusterAt(icl, clusters->GetArray());
649  if(!clust) {
650  AliErrorStream() << "Associated constituent cluster not found\n";
651  continue;
652  }
653  TLorentzVector clusterp;
654  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
655  num += clusterp.Pt() * clusterp.Pt();
656  den += clusterp.Pt();
657  }
658  }
659  return TMath::Sqrt(num)/den;
660 }
661 
663  for(int icl = 0; icl < jet->GetNumberOfClusters(); icl++){
664  auto clust = jet->ClusterAt(icl, clusters->GetArray());
665  AliDebugStream(3) << "cluster time " << clust->GetTOF() << std::endl;
666  fQAHistos->FillTH2("hClusterConstE", jet->Pt(),clust->GetUserDefEnergy(clusters->GetDefaultClusterEnergy()));
667  fQAHistos->FillTH2("hClusterConstTime", jet->Pt(), clust->GetTOF());
668  fQAHistos->FillTH2("hClusterConstM02", jet->Pt(), clust->GetM02());
669  fQAHistos->FillTH2("hClusterConstNcell", jet->Pt(), clust->GetNCells());
670 
671 #ifdef EXPERIMENTAL_JETCONSTITUENTS
672  fQAHistos->FillTH2("hClusterIndexENLC", jet->Pt(), clust->GetNonLinCorrEnergy());
673  fQAHistos->FillTH2("hClusterIndexEHC", jet->Pt(), clust->GetHadCorrEnergy());
674 #endif
675  }
676 
677 #ifdef EXPERIMENTAL_JETCONSTITUENTS
678  // Loop over charged particles - fill test histogram
679  for(int itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
680  auto part = jet->TrackAt(itrk, cont->GetArray());
681  fQAHistos->FillTH2("hChargedIndexPt", jet->Pt(), part->Pt());
682  }
683 
684  // Look over charged constituents
685  AliDebugStream(2) << "Jet: Number of particle constituents: " << jet->GetParticleConstituents().size() << std::endl;
686  for(auto part : jet->GetParticleConstituents()) {
687  //auto part = static_cast<PWG::JETFW::AliEmcalParticleJetConstituent *>(pconst);
688  AliDebugStream(3) << "Found particle constituent with pt " << part.Pt() << ", from VParticle " << part.GetParticle()->Pt() << std::endl;
689  fQAHistos->FillTH2("hChargedConstituentPt", jet->Pt(), part.Pt());
690  }
691 
692  // Loop over neutral constituents
693  AliDebugStream(2) << "Jet: Number of cluster constituents: " << jet->GetClusterConstituents().size() << std::endl;
694  for(auto clust : jet->GetClusterConstituents()){
695  //auto clust = static_cast<PWG::JETFW::AliEmcalClusterJetConstituent *>(cconst);
696  AliDebugStream(3) << "Found cluster constituent with energy " << clust.E() << " using energy definition " << static_cast<int>(clust.GetDefaultEnergyType()) << std::endl;
697  fQAHistos->FillTH2("hClusterConstituentEDefault", jet->Pt(), clust.E());
698  fQAHistos->FillTH2("hClusterConstituentENLC", jet->Pt(), clust.GetCluster()->GetNonLinCorrEnergy());
699  fQAHistos->FillTH2("hClusterConstituentEHC", jet->Pt(), clust.GetCluster()->GetHadCorrEnergy());
700  }
701 
702  // Fill global observables: Leading charged and cluster constituents
703  auto leadingcharged = jet->GetLeadingParticleConstituent();
704  auto leadingcluster = jet->GetLeadingClusterConstituent();
705  if(leadingcluster){
706  fQAHistos->FillTH1("hLeadingClusterConstituentPt", jet->Pt(), leadingcluster->GetCluster()->GetHadCorrEnergy());
707  }
708  if(leadingcharged) {
709  fQAHistos->FillTH1("hLeadingChargedConstituentPt", jet->Pt(), leadingcharged->GetParticle()->Pt());
710  }
711 #endif
712 }
713 
715  int ncharged = 0, nneutral = jet.GetNumberOfClusters();
716  if(particles) {
717  for(decltype(jet.GetNumberOfTracks()) ipart = 0; ipart < jet.GetNumberOfTracks(); ipart++){
718  auto part = jet.TrackAt(ipart, particles->GetArray());
719  if(!part) continue;
720  if(part->Charge()) ncharged++;
721  else nneutral++;
722  }
723  }
724  // check if the jet has at least one consituent for jet substructure
725  int nallowed = 0;
726  nallowed += fUseChargedConstituents ? ncharged : 0;
727  nallowed += fUseNeutralConstituents ? nneutral : 0;
728  return nallowed > 0;
729 }
730 
731 std::string AliAnalysisTaskEmcalJetSubstructureTree::MatchTrigger(const std::string &triggertoken) const {
732  std::vector<std::string> tokens;
733  std::string result;
734  std::stringstream decoder(fInputEvent->GetFiredTriggerClasses().Data());
735  while(std::getline(decoder, result, ' ')) tokens.emplace_back(result);
736  result.clear();
737  for(auto t : tokens) {
738  if(t.find(triggertoken) != std::string::npos) {
739  // take first occurrence - downscale factor should normally be the same
740  result = t;
741  break;
742  }
743  }
744  return result;
745 }
746 
747 bool AliAnalysisTaskEmcalJetSubstructureTree::IsSelectEmcalTriggers(const std::string &triggerstring) const {
748  const std::array<std::string, 8> kEMCALTriggers = {
749  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
750  };
751  bool isEMCAL = false;
752  for(auto emcaltrg : kEMCALTriggers) {
753  if(triggerstring.find(emcaltrg) != std::string::npos) {
754  isEMCAL = true;
755  break;
756  }
757  }
758  return isEMCAL;
759 }
760 
762  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
763 
764  Bool_t isAOD(kFALSE);
765  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
766  if(inputhandler) {
767  if(inputhandler->IsA() == AliAODInputHandler::Class()){
768  std::cout << "Analysing AOD events\n";
769  isAOD = kTRUE;
770  } else {
771  std::cout << "Analysing ESD events\n";
772  }
773  }
774 
775  std::stringstream taskname;
776  taskname << "JetSubstructureTreemaker_R" << std::setw(2) << std::setfill('0') << int(jetradius*10) << trigger;
778  mgr->AddTask(treemaker);
779  treemaker->SetMakeGeneralHistograms(kTRUE);
780  if(isMC) treemaker->SetHasTrueEvent(true);
781  if(isData) treemaker->SetHasRecEvent(true);
782 
783  // Adding containers
784  if(isMC) {
785  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
786  particles->SetMinPt(0.);
787 
788  AliJetContainer *mcjets = treemaker->AddJetContainer(
789  jettype,
791  recombinationScheme,
792  jetradius,
794  particles, nullptr);
795  mcjets->SetName("mcjets");
796  mcjets->SetJetPtCut(20.);
797  }
798 
799  if(isData) {
800  AliTrackContainer *tracks(nullptr);
801  if((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet)){
803  std::cout << "Track container name: " << tracks->GetName() << std::endl;
804  tracks->SetMinPt(0.15);
805  }
806  AliClusterContainer *clusters(nullptr);
807  if((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet)){
808  std::cout << "Using full or neutral jets ..." << std::endl;
810  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
811  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
812  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
813  } else {
814  std::cout << "Using charged jets ... " << std::endl;
815  }
816 
817  AliJetContainer *datajets = treemaker->AddJetContainer(
818  jettype,
820  recombinationScheme,
821  jetradius,
823  tracks, clusters);
824  datajets->SetName("datajets");
825  datajets->SetJetPtCut(20.);
826 
827  treemaker->SetUseAliAnaUtils(true, true);
828  treemaker->SetVzRange(-10., 10);
829 
830  // configure trigger selection
831  std::string triggerstring(trigger);
832  if(triggerstring.find("INT7") != std::string::npos) {
833  treemaker->SetTriggerBits(AliVEvent::kINT7);
834  } else if(triggerstring.find("EJ1") != std::string::npos) {
835  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
836  treemaker->SetTriggerString("EJ1");
837  } else if(triggerstring.find("EJ2") != std::string::npos) {
838  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
839  treemaker->SetTriggerString("EJ2");
840  } else if(triggerstring.find("EG1") != std::string::npos) {
841  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
842  treemaker->SetTriggerString("EG1");
843  } else if(triggerstring.find("EG2") != std::string::npos) {
844  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
845  treemaker->SetTriggerString("EG2");
846  }
847  }
848 
849  std::string jettypestring;
850  switch(jettype) {
851  case AliJetContainer::kFullJet: jettypestring = "FullJets"; break;
852  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; break;
853  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; break;
854  default: jettypestring = "Undef";
855  };
856 
857  // Connecting containers
858  std::stringstream outputfile, histname, treename;
859  outputfile << mgr->GetCommonFileName() << ":JetSubstructure_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
860  histname << "JetSubstructureHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
861  treename << "JetSubstructureTree_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
862  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
863  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
864  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer(treename.str().data(), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()));
865 
866  return treemaker;
867 }
868 
869 void AliSoftDropParameters::LinkJetTreeBranches(TTree *jettree, const char *tag) {
870  LinkBranch(jettree, &fZg, Form("Zg%s", tag), "D");
871  LinkBranch(jettree, &fRg, Form("Rg%s", tag), "D");
872  LinkBranch(jettree, &fMg, Form("Mg%s", tag), "D");
873  LinkBranch(jettree, &fPtg, Form("Ptg%s", tag), "D");
874  LinkBranch(jettree, &fMug, Form("Mug%s", tag), "D");
875  LinkBranch(jettree, &fDeltaR, Form("DeltaRg%s", tag), "D");
876  LinkBranch(jettree, &fNDropped, Form("NDropped%s", tag), "I");
877 };
878 
880  LinkBranch(jettree, &fOneSubjettiness, Form("OneSubjettiness%s", tag), "D");
881  LinkBranch(jettree, &fTwoSubjettiness, Form("TwoSubjettiness%s", tag), "D");
882 }
883 
885  LinkBranch(jettree, &fAngularity, Form("Angularity%s", tag), "D");
886  LinkBranch(jettree, &fPtD, Form("PtD%s", tag), "D");
887 }
888 
889 void AliJetKineParameters::LinkJetTreeBranches(TTree *jettree, const char *tag){
890  LinkBranch(jettree, &fPt, Form("PtJet%s", tag), "D");
891  LinkBranch(jettree, &fE, Form("EJet%s", tag), "D");
892  LinkBranch(jettree, &fEta, Form("Eta%s", tag), "D");
893  LinkBranch(jettree, &fPhi, Form("Phi%s", tag), "D");
894  LinkBranch(jettree, &fArea, Form("Area%s", tag), "D");
895  LinkBranch(jettree, &fMass, Form("Mass%s", tag), "D");
896  LinkBranch(jettree, &fNEF, Form("NEF%s", tag), "D");
897  LinkBranch(jettree, &fNCharged, Form("NCharged%s", tag), "I");
898  LinkBranch(jettree, &fNNeutral, Form("NNeutral%s", tag), "I");
899 }
900 
902  LinkBranch(jettree, &fJetRadius, "Radius", "D");
903  LinkBranch(jettree, &fEventWeight, "EventWeight", "D");
904  LinkBranch(jettree, &fTriggerClusterIndex, "TriggerClusterIndex", "I");
905  if(fillRho) {
906  std::string varnames[] = {"RhoPtRec", "RhoPtSim", "RhoMassRec", "RhoMassSim"};
907  for(int i = 0; i < 4; i++){
908  LinkBranch(jettree, fRhoParamters + i, varnames[i].data(), "D");
909  }
910  }
911 }
912 
913 void LinkBranch(TTree *jettree, void *data, const char *branchname, const char *type) {
914  jettree->Branch(branchname, data, Form("%s/%s", branchname, type));
915 }
916 } /* namespace EmcalTriggerJets */
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Bool_t fHasRecEvent
Has reconstructed event (for trigger selection)
Class creating a linear binning, used in the histogram manager.
Double_t MakePtD(const AliEmcalJet &jet, const AliParticleContainer *const particles, const AliClusterContainer *const clusters) const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
AliNSubjettinessParameters * fNSubTrue
! Data field for true n-subjettiness parameters in jet tree
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
AliJetSubstructureData MakeJetSubstructure(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters, const AliJetSubstructureSettings &settings) const
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:137
static AliEmcalDownscaleFactorsOCDB * Instance()
Bool_t fHasTrueEvent
Has Monte-Carlo truth (for trigger selection)
AliSoftDropParameters * fSoftDropMeasured
! Data field for measured soft drop parameters in jet tree
TCanvas * c
Definition: TestFitELoss.C:172
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
const std::vector< PWG::JETFW::AliEmcalParticleJetConstituent > & GetParticleConstituents() const
Get container with particle (track / MC particle) constituents.
Definition: AliEmcalJet.h:184
Int_t fTriggerClusterIndex
Index of the trigger cluster (0 - CENT, 1 - CENTNOTRD)
Structure for results from the soft drop algorithm.
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
virtual void UserExecOnce()
Task initializations handled in user tasks.
Double_t E() const
Definition: AliEmcalJet.h:119
const PWG::JETFW::AliEmcalParticleJetConstituent * GetLeadingParticleConstituent() const
Get the leading particle constituent.
virtual bool Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
void SetVzRange(Double_t min, Double_t max)
AliNSubjettinessParameters * fNSubMeasured
! Data field for measured n-subjettiness parameters in jet tree
static std::vector< PWG::EMCAL::Triggerinfo > DecodeTriggerString(const std::string &triggerstring)
Decoding trigger string.
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
fastjet::JetAlgorithm fRecluserAlgo
Reclusterization algorithm.
Container for particles within the EMCAL framework.
Int_t GetDefaultClusterEnergy() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliJetKineParameters MakeJetKineParameters(const AliEmcalJet &jet) const
Double_t Px() const
Definition: AliEmcalJet.h:106
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:68
AliRhoParameter * GetRhoFromEvent(const char *name)
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
void SetJetPtCut(Float_t cut)
TPC acceptance.
Definition: AliEmcalJet.h:67
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
static AliAnalysisTaskEmcalJetSubstructureTree * AddEmcalJetSubstructureTreeMaker(Bool_t isMC, Bool_t isData, Double_t jetradius, AliJetContainer::EJetType_t jettype, AliJetContainer::ERecoScheme_t recombinationScheme, Bool_t useDCAL, const char *name)
void LinkBranch(TTree *jettree, void *data, const char *branchname, const char *type)
Helper function linking struct members to branches in the jet substructure tree.
AliJetStructureParameters * fJetStructureTrue
! True jet substructure paramteres
Int_t GetNJets() const
AliNSubjettinessParameters MakeNsubjettinessParameters(const fastjet::PseudoJet &jet, const AliNSubjettinessDefinition &cut) const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Definition for the algorithm obtaining the softdrop parameters.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
virtual void RunChanged(Int_t newrun)
Process tasks relevant when a file with a different run number is processed.
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t MakeAngularity(const AliEmcalJet &jet, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const
AliJetStructureParameters * fJetStructureMeasured
! Measured jet substructure parameters
Double_t Pt() const
Definition: AliEmcalJet.h:109
Bool_t isMC
AliJetKineParameters * fKineRec
! Detector level jet kinematics
const PWG::JETFW::AliEmcalClusterJetConstituent * GetLeadingClusterConstituent() const
Get the leading cluster constituent.
void DoConstituentQA(const AliEmcalJet *jet, const AliParticleContainer *tracks, const AliClusterContainer *clusters)
Float_t GetJetRadius() const
const std::vector< PWG::JETFW::AliEmcalClusterJetConstituent > & GetClusterConstituents() const
Get container with cluster constituents.
Definition: AliEmcalJet.h:190
AliEmcalList * fOutput
!output list
AliJetKineParameters * fKineSim
! Particle level jet kinematics
Double_t fVertex[3]
!event vertex
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Container class for histograms.
Definition: THistManager.h:99
Double_t Pz() const
Definition: AliEmcalJet.h:108
AliSoftDropParameters MakeSoftDropParameters(const fastjet::PseudoJet &jet, const AliSoftdropDefinition &cut) const
AliJetTreeGlobalParameters * fGlobalTreeParams
! Global jet tree parameters (same for all jets in event)
void UserCreateOutputObjects()
Main initialization function on the worker.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:72
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Double_t NEF() const
Definition: AliEmcalJet.h:148
void SetDefaultClusterEnergy(Int_t d)
bool SelectJet(const AliEmcalJet &jet, const AliParticleContainer *particles) const
AliSoftDropParameters * fSoftDropTrue
! Data field for true soft drop parameters in jet tree
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
Double_t M() const
Definition: AliEmcalJet.h:120
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
Bool_t fUseTriggerSelectionForData
Use trigger selection on data (require trigger patch in addition to trigger selection string) ...
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
void SetClusHadCorrEnergyCut(Double_t cut)
static TString TrackContainerNameFactory(Bool_t isAOD)
Get name of the default track container.
static TString ClusterContainerNameFactory(Bool_t isAOD)
Get name of the default cluster container.