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