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