AliPhysics  aaf9c62 (aaf9c62)
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->CreateTH1("hTriggerClusterCounter", "Event counter separating into trigger clusters", 7, -1.5, 5.5);
179  fQAHistos->CreateTH2("hClusterConstE", "EMCAL cluster energy vs jet pt; p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
180  fQAHistos->CreateTH2("hClusterConstTime", "EMCAL cluster time vs. jet pt; p_{t, jet} (GeV/c); t_{cl} (ns)", jetptbinning, timebinning);
181  fQAHistos->CreateTH2("hClusterConstM02", "EMCAL cluster M02 vs. jet pt; p{t, jet} (GeV/c); M02", jetptbinning, m02binning);
182  fQAHistos->CreateTH2("hClusterConstNcell", "EMCAL cluster ncell vs. jet pt; p{t, jet} (GeV/c); Number of cells", jetptbinning, ncellbinning);
183 
184  // Test of constituent QA
185 #ifdef EXPERIMENTAL_JETCONSTITUENTS
186  fQAHistos->CreateTH2("hChargedConstituentPt", "charged constituent pt vs jet pt (via constituent map); p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
187  fQAHistos->CreateTH2("hChargedIndexPt", "charged constituent pt vs jet pt (via index map); p_{t, jet} (GeV/c); p_{t, ch} (GeV/c)", jetptbinning, clusterenergybinning);
188 
189  fQAHistos->CreateTH2("hClusterConstituentEDefault", "cluster constituent default energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
190  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);
191  fQAHistos->CreateTH2("hClusterConstituentEHC", "cluster constituent hadronic-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
192  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);
193  fQAHistos->CreateTH2("hClusterIndexEHC", "cluster constituent hadronic-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
194  fQAHistos->CreateTH2("hLeadingChargedConstituentPt", "Pt of the leading charged constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
195  fQAHistos->CreateTH2("hLeadingClusterConstituentPt", "Pt of the leading cluster constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
196 #endif
197  for(auto h : *(fQAHistos->GetListOfHistograms())) fOutput->Add(h);
198 
199  OpenFile(2);
200  std::string treename = this->GetOutputSlot(2)->GetContainer()->GetName();
201  fJetSubstructureTree = new TTree(treename.data(), "Tree with jet substructure information");
202 
204  fGlobalTreeParams->LinkJetTreeBranches(fJetSubstructureTree, fFillRho);
207  if(fFillPart) {
210  }
211  if(fFillSoftDrop) {
214  if(fFillPart) {
217  }
218  }
219  if(fFillNSub) {
222  if(fFillPart) {
225  }
226  }
227 
228  if(fFillStructGlob){
231  if(fFillPart){
234  }
235  }
236 
237  PostData(1, fOutput);
238  PostData(2, fJetSubstructureTree);
239 }
240 
244  }
245 }
246 
250  AliParticleContainer *particles = GetParticleContainer("mcparticles");
251 
252  AliJetContainer *mcjets = GetJetContainer("mcjets");
253  AliJetContainer *datajets = GetJetContainer("datajets");
254 
255  FillLuminosity(); // Makes only sense in data
256 
257  // for(auto e : *(fInputEvent->GetList())) std::cout << e->GetName() << std::endl;
258 
259  std::stringstream rhoTagData, rhoTagMC;
260  if(datajets) rhoTagData << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(datajets->GetJetRadius() * 10.);
261  if(mcjets) rhoTagMC << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(mcjets->GetJetRadius() * 10.);
262 
263  if(fFillRho){
264  std::string rhoSparseData = "RhoSparse_Full_" + rhoTagData.str(), rhoSparseMC = "RhoSparse_Full_" + rhoTagMC.str(),
265  rhoMassData = "RhoMassSparse_Full_" + rhoTagData.str(), rhoMassMC = "RhoMassSparse_Full_" + rhoTagMC.str();
266  AliRhoParameter *rhoPtRec = GetRhoFromEvent(rhoSparseData.data()),
267  *rhoMassRec = GetRhoFromEvent(rhoMassData.data()),
268  *rhoPtSim = GetRhoFromEvent(rhoSparseMC.data()),
269  *rhoMassSim = GetRhoFromEvent(rhoMassMC.data());
270  AliDebugStream(2) << "Found rho parameter for reconstructed pt: " << (rhoPtRec ? "yes" : "no") << ", value: " << (rhoPtRec ? rhoPtRec->GetVal() : 0.) << std::endl;
271  AliDebugStream(2) << "Found rho parameter for sim pt: " << (rhoPtSim ? "yes" : "no") << ", value: " << (rhoPtSim ? rhoPtSim->GetVal() : 0.) << std::endl;
272  AliDebugStream(2) << "Found rho parameter for reconstructed Mass: " << (rhoMassRec ? "yes" : "no") << ", value: " << (rhoMassRec ? rhoMassRec->GetVal() : 0.) << std::endl;
273  AliDebugStream(2) << "Found rho parameter for sim Mass: " << (rhoMassSim ? "yes" : "no") << ", value: " << (rhoMassSim ? rhoMassSim->GetVal() : 0.) << std::endl;
274  Double_t rhopars[4] = {
275  rhoPtRec ? rhoPtRec->GetVal() : 0.,
276  rhoPtSim ? rhoPtSim->GetVal() : 0.,
277  rhoMassRec ? rhoMassRec->GetVal() : 0.,
278  rhoMassSim ? rhoMassSim->GetVal() : 0.
279  };
280  memcpy(this->fGlobalTreeParams->fRhoParamters, rhopars, sizeof(Double_t) * 4);
281  }
282 
283  AliDebugStream(1) << "Inspecting jet radius " << (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius()) << std::endl;
284  this->fGlobalTreeParams->fJetRadius = (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius());
285  fGlobalTreeParams->fTriggerClusterIndex = -1; // Reset trigger cluster index
286 
287  if(datajets && !mcjets){
288  // decode trigger string in order to determine the trigger clusters
289  std::vector<std::string> clusternames;
290  auto triggerinfos = PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data());
291  for(auto t : triggerinfos) {
292  if(std::find(clusternames.begin(), clusternames.end(), t.Triggercluster()) == clusternames.end()) clusternames.emplace_back(t.Triggercluster());
293  }
294  bool isCENT = (std::find(clusternames.begin(), clusternames.end(), "CENT") != clusternames.end()),
295  isCENTNOTRD = (std::find(clusternames.begin(), clusternames.end(), "CENTNOTRD") != clusternames.end()),
296  isCALO = (std::find(clusternames.begin(), clusternames.end(), "CALO") != clusternames.end()),
297  isCALOFAST = (std::find(clusternames.begin(), clusternames.end(), "CALOFAST") != clusternames.end());
298  if(isCENT || isCENTNOTRD) {
299  if(isCENT && isCENTNOTRD) fGlobalTreeParams->fTriggerClusterIndex = 0; // CENTBOTH
300  else if(isCENT) fGlobalTreeParams->fTriggerClusterIndex = 1; // OnlyCENT
301  else fGlobalTreeParams->fTriggerClusterIndex = 2; // OnlyCENTNOTRD
302  }
303  if(isCALO || isCALOFAST){
304  // CALO(FAST) and CENT(NOTRD) clusters disjunct - CALO cluster not read out in parallel with CENT cluster
305  // Therefore mixing between CALO and CENT triggers doesn't need to be handled.
306  if(isCALO && isCALOFAST) fGlobalTreeParams->fTriggerClusterIndex = 3;
307  else if(isCALO) fGlobalTreeParams->fTriggerClusterIndex = 4;
308  else if(isCALOFAST) fGlobalTreeParams->fTriggerClusterIndex = 5;
309  }
310  }
311 
312  double weight = 1.;
314  AliDebugStream(2) << "Trigger selection string: " << fTriggerSelectionString << std::endl;
315  TString selectionString = (fTriggerSelectionBits & AliVEvent::kINT7) ? "INT7" : fTriggerSelectionString;
316  auto triggerstring = MatchTrigger(selectionString.Data());
317  AliDebugStream(2) << "Getting downscale correction factor for trigger string " << triggerstring << std::endl;
319  }
320  AliDebugStream(1) << "Using downscale weight " << weight << std::endl;
321  this->fGlobalTreeParams->fEventWeight = weight;
322 
323 
324  // Count events (for spectrum analysis)
325  fQAHistos->FillTH1("hEventCounter", 1);
326  fQAHistos->FillTH1("hTriggerClusterCounter", fGlobalTreeParams->fTriggerClusterIndex);
327 
328  AliSoftdropDefinition softdropSettings;
329  softdropSettings.fBeta = fSDBetaCut;
330  softdropSettings.fZ = fSDZCut;
331  switch(fReclusterizer) {
332  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
333  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
334  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
335  };
336 
337  AliNSubjettinessDefinition nsubjettinessSettings;
338  nsubjettinessSettings.fBeta = 1.;
339  nsubjettinessSettings.fRadius = 0.4;
340 
341  if(datajets) {
342  AliDebugStream(1) << "In data jets branch: found " << datajets->GetNJets() << " jets, " << datajets->GetNAcceptedJets() << " were accepted\n";
343  AliDebugStream(1) << "Having MC information: " << (mcjets ? TString::Format("yes, with %d jets", mcjets->GetNJets()) : "no") << std::endl;
344  if(mcjets) {
345  AliDebugStream(1) << "In MC jets branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
346  }
347  for(auto jet : datajets->accepted()) {
348  double pt = jet->Pt(), pz = jet->Pz(), E = jet->E(), M = TMath::Sqrt(E*E - pt*pt - pz*pz);
349  AliDebugStream(2) << "Next jet: pt:" << jet->Pt() << ", E: " << E << ", pz: " << pz << ", M(self): " << M << "M(fj)" << jet->M() << std::endl;
350  AliEmcalJet *associatedJet = jet->ClosestJet();
351 
352  if(mcjets) {
353  if(!associatedJet) {
354  AliDebugStream(2) << "Not found associated jet" << std::endl;
355  continue;
356  }
357  if(!(SelectJet(*jet, tracks) && SelectJet(*associatedJet, particles))) continue;
358  try {
359  DoConstituentQA(jet, tracks, clusters);
360  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, datajets->GetJetRadius() * 2., tracks, clusters, {softdropSettings, nsubjettinessSettings}),
361  structureMC = fFillPart ? MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings}) : AliJetSubstructureData();
363  if(fKineSim) *fKineSim = MakeJetKineParameters(*associatedJet);
364  if(fSoftDropMeasured) *fSoftDropMeasured = structureData.fSoftDrop;
365  if(fSoftDropTrue) *fSoftDropTrue = structureMC.fSoftDrop;
366  if(fNSubMeasured) *fNSubMeasured = structureData.fNsubjettiness;
367  if(fNSubTrue) *fNSubTrue = structureMC.fNsubjettiness;
368  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
369  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*associatedJet, particles, nullptr), MakePtD(*associatedJet, particles, nullptr)};
370  fJetSubstructureTree->Fill();
371  } catch(ReclusterizerException &e) {
372  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
373  } catch(SubstructureException &e) {
374  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
375  }
376  } else {
377  if(!SelectJet(*jet, tracks)) continue;
378  try {
379  DoConstituentQA(jet, tracks, clusters);
380  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
383  if(fNSubMeasured) *fNSubMeasured = structure.fNsubjettiness;
384  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
385  fJetSubstructureTree->Fill();
386  } catch(ReclusterizerException &e) {
387  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
388  } catch(SubstructureException &e) {
389  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
390  }
391  }
392  }
393  } else {
394  if(mcjets) {
395  // for MCgen train
396  AliDebugStream(1) << "In MC pure jet branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
397  for(auto j : mcjets->accepted()){
398  AliEmcalJet *mcjet = static_cast<AliEmcalJet *>(j);
399  try {
400  AliJetSubstructureData structure = MakeJetSubstructure(*mcjet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
401  if(this->fKineSim) *fKineSim = MakeJetKineParameters(*mcjet);
402  if(fSoftDropTrue) *fSoftDropTrue = structure.fSoftDrop;
403  if(fNSubTrue) *fNSubTrue = structure.fNsubjettiness;
404  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*mcjet, particles, nullptr), MakePtD(*mcjet, particles, nullptr)};
405  fJetSubstructureTree->Fill();
406  } catch (ReclusterizerException &e) {
407  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
408  } catch (SubstructureException &e) {
409  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
410  }
411  }
412  }
413  }
414 
415  return true;
416 }
417 
419  AliDebugStream(1) << "Trigger selection called\n";
420  if(!(fHasRecEvent || fHasTrueEvent)){
421  AliErrorStream() << "Impossible combination: Neither rec nor true event available. Rejecting ..." << std::endl;
422  return false;
423  }
424  // Run trigger selection (not on pure MCgen train - pure MCgen train has no rec event, ESD event is fake there)
425  if(fHasRecEvent){
426  if(!fHasTrueEvent){
427  // Pure data - do EMCAL trigger selection from selection string
428  AliDebugStream(1) << "Applying trigger selection for trigger bits " << std::bitset<sizeof(decltype(fTriggerSelectionBits)) * 8>(fTriggerSelectionBits) << "and trigger selection string " << fTriggerSelectionString << std::endl;
429  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
430  AliDebugStream(1) << "Passed trigger bit selection" << std::endl;
431  if(fTriggerSelectionString.Length()) {
432  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
433  AliDebugStream(1) << "Passed trigger string section" << std::endl;
435  auto trgselresult = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
436  AliDebugStream(1) << "Found trigger decision object: " << (trgselresult ? "yes" : "no") << std::endl;
437  if(!trgselresult){
438  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
439  return false;
440  }
441  if(!trgselresult->IsEventSelected(fTriggerSelectionString)) return false;
442  AliDebugStream(1) << "Data event selected" << std::endl;
443  }
444  }
445  } else {
446  // Simulation - do EMCAL trigger selection from trigger selection object
447  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false; // Require INT7 trigger - EMCAL triggers will be a subset
449  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
450  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
451  if(!mctrigger){
452  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
453  return false;
454  }
455  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
456  }
457  }
458  }
459  return true; // trigger selected or pure MCgen information
460 }
461 
463  AliCDBManager * cdb = AliCDBManager::Instance();
464  if(!fMCEvent && cdb){
465  // Get List of trigger clusters
466  AliCDBEntry *en = cdb->Get("GRP/CTP/Config");
467  AliTriggerConfiguration *trg = static_cast<AliTriggerConfiguration *>(en->GetObject());
468  std::vector<std::string> clusternames;
469  for(auto c : trg->GetClusters()) {
470  AliTriggerCluster *clust = static_cast<AliTriggerCluster *>(c);
471  std::string clustname = clust->GetName();
472  auto iscent = clustname.find("CENT") != std::string::npos, iscalo = clustname.find("CALO") != std::string::npos;
473  if(!(iscalo || iscent)) continue;
474  clusternames.emplace_back(clustname);
475  }
476 
477  // Set the x-axis of the luminosity monitor histogram
478  fLumiMonitor = new TH1F("hLumiMonitor", "Luminosity monitor", clusternames.size(), 0, clusternames.size());
479  int currentbin(1);
480  for(auto c : clusternames) {
481  fLumiMonitor->GetXaxis()->SetBinLabel(currentbin++, c.data());
482  }
483  fOutput->Add(fLumiMonitor);
484  }
485 }
486 
489  auto downscalefactors = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance();
490  if(fInputEvent->GetFiredTriggerClasses().Contains("INT7")) {
491  for(auto trigger : PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data())){
492  auto int7trigger = trigger.IsTriggerClass("INT7");
493  auto bunchcrossing = trigger.BunchCrossing() == "B";
494  auto nopf = trigger.PastFutureProtection() == "NOPF";
495  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;
496  if(int7trigger && bunchcrossing && nopf) {
497  double downscale = downscalefactors->GetDownscaleFactorForTriggerClass(trigger.ExpandClassName());
498  AliDebugStream(5) << "Using downscale " << downscale << std::endl;
499  fLumiMonitor->Fill(trigger.Triggercluster().data(), 1./downscale);
500  }
501  }
502  }
503  }
504 }
505 
507  AliJetKineParameters result;
508  result.fPt = TMath::Abs(jet.Pt());
509  result.fE = jet.E();
510  result.fEta = jet.Eta();
511  result.fPhi = jet.Phi();
512  result.fArea = jet.Area();
513  result.fMass = jet.M();
514  result.fNEF = jet.NEF();
515  result.fNCharged = jet.GetNumberOfTracks();
516  result.fNNeutral = jet.GetNumberOfClusters();
517  return result;
518 }
519 
521  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
522  std::vector<fastjet::PseudoJet> constituents;
523  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
524  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
525  if(tracks && (fUseChargedConstituents || isMC)){ // Neutral particles part of particle container in case of MC
526  AliDebugStream(1) << "Jet substructure: Using charged constituents" << std::endl;
527  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
528  auto track = jet.TrackAt(itrk, tracks->GetArray());
529  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
530  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
531  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
532  constituentTrack.set_user_index(jet.TrackAt(itrk));
533  constituents.push_back(constituentTrack);
534  }
535  }
536 
537  if(clusters && fUseNeutralConstituents){
538  AliDebugStream(1) << "Jet substructure: Using neutral constituents" << std::endl;
539  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
540  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
541  TLorentzVector clustervec;
542  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
543  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
544  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
545  constituents.push_back(constituentCluster);
546  }
547  }
548 
549  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
550  if(!constituents.size()){
551  AliErrorStream() << "Jet has 0 constituents." << std::endl;
552  throw ReclusterizerException();
553  }
554  // Redo jet finding on constituents with a
555  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
556  std::vector<fastjet::PseudoJet> outputjets;
557  try {
558  fastjet::ClusterSequence jetfinder(constituents, jetdef);
559  outputjets = jetfinder.inclusive_jets(0);
561  return result;
562  } catch (fastjet::Error &e) {
563  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
564  throw ReclusterizerException();
565  } catch (SoftDropException &e) {
566  AliErrorStream() << "Softdrop exception caught: " << e.what() << std::endl;
567  throw ReclusterizerException();
568  }
569 }
570 
572  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
573  softdropAlgorithm.set_verbose_structure(kTRUE);
574  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
575  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
576  AliDebugStream(4) << "Jet has " << jet.constituents().size() << " constituents" << std::endl;
577  auto groomed = softdropAlgorithm(jet);
578  try {
579  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
580 
581  AliSoftDropParameters result({softdropstruct.symmetry(),
582  groomed.m(),
583  softdropstruct.delta_R(),
584  groomed.perp(),
585  softdropstruct.delta_R(),
586  softdropstruct.mu(),
587  softdropstruct.dropped_count()});
588  return result;
589  } catch(std::bad_cast &e) {
590  throw SoftDropException();
591  }
592 }
593 
596  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet),
597  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet)
598  });
599  return result;
600 }
601 
603  if(!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
604  throw SubstructureException();
605  TVector3 jetvec(jet.Px(), jet.Py(), jet.Pz());
606  Double_t den(0.), num(0.);
607  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
608  if(tracks && (fUseChargedConstituents || isMC)){
609  AliDebugStream(1) << "Angularity: Using charged constituents" << std::endl;
610  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
611  auto track = jet.TrackAt(itrk, tracks->GetArray());
612  if(!track){
613  AliErrorStream() << "Associated constituent particle / track not found\n";
614  continue;
615  }
616  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
617  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
618  TVector3 trackvec(track->Px(), track->Py(), track->Pz());
619 
620  num += track->Pt() * trackvec.DrEtaPhi(jetvec);
621  den += +track->Pt();
622  }
623  }
624  if(clusters && fUseNeutralConstituents) {
625  AliDebugStream(1) << "Using neutral constituents" << std::endl;
626  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
627  auto clust = jet.ClusterAt(icl, clusters->GetArray());
628  if(!clust) {
629  AliErrorStream() << "Associated constituent cluster not found\n";
630  continue;
631  }
632  TLorentzVector clusterp;
633  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
634 
635  num += clusterp.Pt() * clusterp.Vect().DrEtaPhi(jetvec);
636  den += clusterp.Pt();
637  }
638  }
639  return num/den;
640 }
641 
643  if (!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
644  throw SubstructureException();
645  Double_t den(0.), num(0.);
646  bool isMC = dynamic_cast<const AliMCParticleContainer *>(particles);
647  if(particles && (fUseChargedConstituents || isMC)){
648  AliDebugStream(1) << "Using charged constituents" << std::endl;
649  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
650  auto trk = jet.TrackAt(itrk, particles->GetArray());
651  if(!trk){
652  AliErrorStream() << "Associated constituent particle / track not found\n";
653  continue;
654  }
655  if(!trk->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
656  if(trk->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
657  num += trk->Pt() * trk->Pt();
658  den += trk->Pt();
659  }
660  }
661  if(clusters && fUseNeutralConstituents){
662  AliDebugStream(1) << "Using neutral constituents" << std::endl;
663  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
664  auto clust = jet.ClusterAt(icl, clusters->GetArray());
665  if(!clust) {
666  AliErrorStream() << "Associated constituent cluster not found\n";
667  continue;
668  }
669  TLorentzVector clusterp;
670  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
671  num += clusterp.Pt() * clusterp.Pt();
672  den += clusterp.Pt();
673  }
674  }
675  return TMath::Sqrt(num)/den;
676 }
677 
679  for(int icl = 0; icl < jet->GetNumberOfClusters(); icl++){
680  auto clust = jet->ClusterAt(icl, clusters->GetArray());
681  AliDebugStream(3) << "cluster time " << clust->GetTOF() << std::endl;
682  fQAHistos->FillTH2("hClusterConstE", jet->Pt(),clust->GetUserDefEnergy(clusters->GetDefaultClusterEnergy()));
683  fQAHistos->FillTH2("hClusterConstTime", jet->Pt(), clust->GetTOF());
684  fQAHistos->FillTH2("hClusterConstM02", jet->Pt(), clust->GetM02());
685  fQAHistos->FillTH2("hClusterConstNcell", jet->Pt(), clust->GetNCells());
686 
687 #ifdef EXPERIMENTAL_JETCONSTITUENTS
688  fQAHistos->FillTH2("hClusterIndexENLC", jet->Pt(), clust->GetNonLinCorrEnergy());
689  fQAHistos->FillTH2("hClusterIndexEHC", jet->Pt(), clust->GetHadCorrEnergy());
690 #endif
691  }
692 
693 #ifdef EXPERIMENTAL_JETCONSTITUENTS
694  // Loop over charged particles - fill test histogram
695  for(int itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
696  auto part = jet->TrackAt(itrk, cont->GetArray());
697  fQAHistos->FillTH2("hChargedIndexPt", jet->Pt(), part->Pt());
698  }
699 
700  // Look over charged constituents
701  AliDebugStream(2) << "Jet: Number of particle constituents: " << jet->GetParticleConstituents().size() << std::endl;
702  for(auto part : jet->GetParticleConstituents()) {
703  //auto part = static_cast<PWG::JETFW::AliEmcalParticleJetConstituent *>(pconst);
704  AliDebugStream(3) << "Found particle constituent with pt " << part.Pt() << ", from VParticle " << part.GetParticle()->Pt() << std::endl;
705  fQAHistos->FillTH2("hChargedConstituentPt", jet->Pt(), part.Pt());
706  }
707 
708  // Loop over neutral constituents
709  AliDebugStream(2) << "Jet: Number of cluster constituents: " << jet->GetClusterConstituents().size() << std::endl;
710  for(auto clust : jet->GetClusterConstituents()){
711  //auto clust = static_cast<PWG::JETFW::AliEmcalClusterJetConstituent *>(cconst);
712  AliDebugStream(3) << "Found cluster constituent with energy " << clust.E() << " using energy definition " << static_cast<int>(clust.GetDefaultEnergyType()) << std::endl;
713  fQAHistos->FillTH2("hClusterConstituentEDefault", jet->Pt(), clust.E());
714  fQAHistos->FillTH2("hClusterConstituentENLC", jet->Pt(), clust.GetCluster()->GetNonLinCorrEnergy());
715  fQAHistos->FillTH2("hClusterConstituentEHC", jet->Pt(), clust.GetCluster()->GetHadCorrEnergy());
716  }
717 
718  // Fill global observables: Leading charged and cluster constituents
719  auto leadingcharged = jet->GetLeadingParticleConstituent();
720  auto leadingcluster = jet->GetLeadingClusterConstituent();
721  if(leadingcluster){
722  fQAHistos->FillTH1("hLeadingClusterConstituentPt", jet->Pt(), leadingcluster->GetCluster()->GetHadCorrEnergy());
723  }
724  if(leadingcharged) {
725  fQAHistos->FillTH1("hLeadingChargedConstituentPt", jet->Pt(), leadingcharged->GetParticle()->Pt());
726  }
727 #endif
728 }
729 
731  int ncharged = 0, nneutral = jet.GetNumberOfClusters();
732  if(particles) {
733  for(decltype(jet.GetNumberOfTracks()) ipart = 0; ipart < jet.GetNumberOfTracks(); ipart++){
734  auto part = jet.TrackAt(ipart, particles->GetArray());
735  if(!part) continue;
736  if(part->Charge()) ncharged++;
737  else nneutral++;
738  }
739  }
740  // check if the jet has at least one consituent for jet substructure
741  int nallowed = 0;
742  nallowed += fUseChargedConstituents ? ncharged : 0;
743  nallowed += fUseNeutralConstituents ? nneutral : 0;
744  return nallowed > 0;
745 }
746 
747 std::string AliAnalysisTaskEmcalJetSubstructureTree::MatchTrigger(const std::string &triggertoken) const {
748  std::vector<std::string> tokens;
749  std::string result;
750  std::stringstream decoder(fInputEvent->GetFiredTriggerClasses().Data());
751  while(std::getline(decoder, result, ' ')) tokens.emplace_back(result);
752  result.clear();
753  for(auto t : tokens) {
754  if(t.find(triggertoken) != std::string::npos) {
755  // take first occurrence - downscale factor should normally be the same
756  result = t;
757  break;
758  }
759  }
760  return result;
761 }
762 
763 bool AliAnalysisTaskEmcalJetSubstructureTree::IsSelectEmcalTriggers(const std::string &triggerstring) const {
764  const std::array<std::string, 8> kEMCALTriggers = {
765  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
766  };
767  bool isEMCAL = false;
768  for(auto emcaltrg : kEMCALTriggers) {
769  if(triggerstring.find(emcaltrg) != std::string::npos) {
770  isEMCAL = true;
771  break;
772  }
773  }
774  return isEMCAL;
775 }
776 
778  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
779 
780  Bool_t isAOD(kFALSE);
781  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
782  if(inputhandler) {
783  if(inputhandler->IsA() == AliAODInputHandler::Class()){
784  std::cout << "Analysing AOD events\n";
785  isAOD = kTRUE;
786  } else {
787  std::cout << "Analysing ESD events\n";
788  }
789  }
790 
791  std::stringstream taskname;
792  taskname << "JetSubstructureTreemaker_R" << std::setw(2) << std::setfill('0') << int(jetradius*10) << trigger;
794  mgr->AddTask(treemaker);
795  treemaker->SetMakeGeneralHistograms(kTRUE);
796  if(isMC) treemaker->SetHasTrueEvent(true);
797  if(isData) treemaker->SetHasRecEvent(true);
798 
799  // Adding containers
800  if(isMC) {
801  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
802  particles->SetMinPt(0.);
803 
804  AliJetContainer *mcjets = treemaker->AddJetContainer(
805  jettype,
807  recombinationScheme,
808  jetradius,
810  particles, nullptr);
811  mcjets->SetName("mcjets");
812  mcjets->SetJetPtCut(20.);
813  }
814 
815  if(isData) {
816  AliTrackContainer *tracks(nullptr);
817  if((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet)){
819  std::cout << "Track container name: " << tracks->GetName() << std::endl;
820  tracks->SetMinPt(0.15);
821  }
822  AliClusterContainer *clusters(nullptr);
823  if((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet)){
824  std::cout << "Using full or neutral jets ..." << std::endl;
826  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
827  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
828  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
829  } else {
830  std::cout << "Using charged jets ... " << std::endl;
831  }
832 
833  AliJetContainer *datajets = treemaker->AddJetContainer(
834  jettype,
836  recombinationScheme,
837  jetradius,
839  tracks, clusters);
840  datajets->SetName("datajets");
841  datajets->SetJetPtCut(20.);
842 
843  treemaker->SetUseAliAnaUtils(true, true);
844  treemaker->SetVzRange(-10., 10);
845 
846  // configure trigger selection
847  std::string triggerstring(trigger);
848  if(triggerstring.find("INT7") != std::string::npos) {
849  treemaker->SetTriggerBits(AliVEvent::kINT7);
850  } else if(triggerstring.find("EJ1") != std::string::npos) {
851  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
852  treemaker->SetTriggerString("EJ1");
853  } else if(triggerstring.find("EJ2") != std::string::npos) {
854  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
855  treemaker->SetTriggerString("EJ2");
856  } else if(triggerstring.find("EG1") != std::string::npos) {
857  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
858  treemaker->SetTriggerString("EG1");
859  } else if(triggerstring.find("EG2") != std::string::npos) {
860  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
861  treemaker->SetTriggerString("EG2");
862  }
863  }
864 
865  std::string jettypestring;
866  switch(jettype) {
867  case AliJetContainer::kFullJet: jettypestring = "FullJets"; break;
868  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; break;
869  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; break;
870  default: jettypestring = "Undef";
871  };
872 
873  // Connecting containers
874  std::stringstream outputfile, histname, treename;
875  outputfile << mgr->GetCommonFileName() << ":JetSubstructure_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
876  histname << "JetSubstructureHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
877  treename << "JetSubstructureTree_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
878  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
879  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
880  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer(treename.str().data(), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()));
881 
882  return treemaker;
883 }
884 
885 void AliSoftDropParameters::LinkJetTreeBranches(TTree *jettree, const char *tag) {
886  LinkBranch(jettree, &fZg, Form("Zg%s", tag), "D");
887  LinkBranch(jettree, &fRg, Form("Rg%s", tag), "D");
888  LinkBranch(jettree, &fMg, Form("Mg%s", tag), "D");
889  LinkBranch(jettree, &fPtg, Form("Ptg%s", tag), "D");
890  LinkBranch(jettree, &fMug, Form("Mug%s", tag), "D");
891  LinkBranch(jettree, &fDeltaR, Form("DeltaRg%s", tag), "D");
892  LinkBranch(jettree, &fNDropped, Form("NDropped%s", tag), "I");
893 };
894 
896  LinkBranch(jettree, &fOneSubjettiness, Form("OneSubjettiness%s", tag), "D");
897  LinkBranch(jettree, &fTwoSubjettiness, Form("TwoSubjettiness%s", tag), "D");
898 }
899 
901  LinkBranch(jettree, &fAngularity, Form("Angularity%s", tag), "D");
902  LinkBranch(jettree, &fPtD, Form("PtD%s", tag), "D");
903 }
904 
905 void AliJetKineParameters::LinkJetTreeBranches(TTree *jettree, const char *tag){
906  LinkBranch(jettree, &fPt, Form("PtJet%s", tag), "D");
907  LinkBranch(jettree, &fE, Form("EJet%s", tag), "D");
908  LinkBranch(jettree, &fEta, Form("Eta%s", tag), "D");
909  LinkBranch(jettree, &fPhi, Form("Phi%s", tag), "D");
910  LinkBranch(jettree, &fArea, Form("Area%s", tag), "D");
911  LinkBranch(jettree, &fMass, Form("Mass%s", tag), "D");
912  LinkBranch(jettree, &fNEF, Form("NEF%s", tag), "D");
913  LinkBranch(jettree, &fNCharged, Form("NCharged%s", tag), "I");
914  LinkBranch(jettree, &fNNeutral, Form("NNeutral%s", tag), "I");
915 }
916 
918  LinkBranch(jettree, &fJetRadius, "Radius", "D");
919  LinkBranch(jettree, &fEventWeight, "EventWeight", "D");
920  LinkBranch(jettree, &fTriggerClusterIndex, "TriggerClusterIndex", "I");
921  if(fillRho) {
922  std::string varnames[] = {"RhoPtRec", "RhoPtSim", "RhoMassRec", "RhoMassSim"};
923  for(int i = 0; i < 4; i++){
924  LinkBranch(jettree, fRhoParamters + i, varnames[i].data(), "D");
925  }
926  }
927 }
928 
929 void LinkBranch(TTree *jettree, void *data, const char *branchname, const char *type) {
930  jettree->Branch(branchname, data, Form("%s/%s", branchname, type));
931 }
932 } /* 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.