AliPhysics  b0c77bb (b0c77bb)
AliAnalysisTaskEmcalJetSubstructureTree.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
27 #include <algorithm>
28 #include <array>
29 #include <bitset>
30 #include <iostream>
31 #include <string>
32 #include <set>
33 #include <sstream>
34 #include <vector>
35 
36 #include <fastjet/ClusterSequence.hh>
37 #include <fastjet/contrib/Nsubjettiness.hh>
38 #include <fastjet/contrib/SoftDrop.hh>
39 #include <fastjet/config.h>
40 #if FASJET_VERSION_NUMBER >= 30302
41 #include <fastjet/tools/Recluster.hh>
42 #else
43 #include <fastjet/contrib/Recluster.hh>
44 #endif
45 
46 #include <THistManager.h>
47 #include <TLinearBinning.h>
48 #include <TLorentzVector.h>
49 #include <TMath.h>
50 #include <TObjString.h>
51 #include <TString.h>
52 #include <TVector3.h>
53 
54 #include "AliAODEvent.h"
55 #include "AliAODInputHandler.h"
56 #include "AliAnalysisManager.h"
57 #include "AliAnalysisDataSlot.h"
58 #include "AliAnalysisDataContainer.h"
60 #include "AliCDBEntry.h"
61 #include "AliCDBManager.h"
62 #include "AliClusterContainer.h"
63 #include "AliJetContainer.h"
66 #include "AliEmcalJet.h"
67 #include "AliEmcalList.h"
70 #include "AliLog.h"
71 #include "AliParticleContainer.h"
72 #include "AliRhoParameter.h"
73 #include "AliTrackContainer.h"
74 #include "AliTriggerCluster.h"
75 #include "AliTriggerConfiguration.h"
76 #include "AliVCluster.h"
77 #include "AliVParticle.h"
78 
79 #ifdef EXPERIMENTAL_JETCONSTITUENTS
82 #endif
83 
84 
88 
89 namespace EmcalTriggerJets {
90 
93  fJetSubstructureTree(nullptr),
94  fGlobalTreeParams(nullptr),
95  fSoftDropMeasured(nullptr),
96  fSoftDropTrue(nullptr),
97  fNSubMeasured(nullptr),
98  fNSubTrue(nullptr),
99  fKineRec(nullptr),
100  fKineSim(nullptr),
101  fJetStructureMeasured(nullptr),
102  fJetStructureTrue(nullptr),
103  fQAHistos(nullptr),
104  fLumiMonitor(nullptr),
105  fSDZCut(0.1),
106  fSDBetaCut(0),
107  fReclusterizer(kCAAlgo),
108  fHasRecEvent(false),
109  fHasTrueEvent(false),
110  fTriggerSelectionBits(AliVEvent::kAny),
111  fTriggerSelectionString(""),
112  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
113  fUseTriggerSelectionForData(false),
114  fUseDownscaleWeight(false),
115  fUseChargedConstituents(true),
116  fUseNeutralConstituents(true),
117  fFillPart(true),
118  fFillRho(true),
119  fFillSoftDrop(true),
120  fFillNSub(true),
121  fFillStructGlob(true)
122 {
123 }
124 
126  AliAnalysisTaskEmcalJet(name, kTRUE),
133  fKineRec(nullptr),
134  fKineSim(nullptr),
139  fSDZCut(0.1),
140  fSDBetaCut(0),
142  fHasRecEvent(false),
143  fHasTrueEvent(false),
144  fTriggerSelectionBits(AliVEvent::kAny),
146  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
148  fUseDownscaleWeight(false),
151  fFillPart(true),
152  fFillRho(true),
153  fFillSoftDrop(true),
154  fFillNSub(true),
155  fFillStructGlob(true)
156 {
157  DefineOutput(2, TTree::Class());
158  SetUseAliAnaUtils(true);
159 }
160 
164  if(fSoftDropTrue) delete fSoftDropTrue;
165  if(fNSubMeasured) delete fNSubMeasured;
166  if(fNSubTrue) delete fNSubTrue;
167  if(fKineRec) delete fKineRec;
168  if(fKineSim) delete fKineSim;
171 }
172 
175 
176  // Make QA for constituent clusters
177  TLinearBinning jetptbinning(9, 20, 200),
178  clusterenergybinning(200, 0., 200),
179  cellenergybinning(1000, 0., 100),
180  timebinning(1000, -500., 500.),
181  m02binning(100, 0., 1.),
182  ncellbinning(101, -0.5, 100.5),
183  exoticsbinning(2, -0.5, 1.5);
184  fQAHistos = new THistManager("QAhistos");
185  fQAHistos->CreateTH1("hEventCounter", "Event counter", 1, 0.5, 1.5);
186  fQAHistos->CreateTH1("hTriggerClusterCounter", "Event counter separating into trigger clusters", 7, -1.5, 5.5);
187  fQAHistos->CreateTH2("hClusterConstE", "EMCAL cluster energy vs jet pt; p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
188  fQAHistos->CreateTH2("hClusterConstTime", "EMCAL cluster time vs. jet pt; p_{t, jet} (GeV/c); t_{cl} (ns)", jetptbinning, timebinning);
189  fQAHistos->CreateTH2("hClusterConstM02", "EMCAL cluster M02 vs. jet pt; p{t, jet} (GeV/c); M02", jetptbinning, m02binning);
190  fQAHistos->CreateTH2("hClusterConstNcell", "EMCAL cluster ncell vs. jet pt; p{t, jet} (GeV/c); Number of cells", jetptbinning, ncellbinning);
191  fQAHistos->CreateTH2("hClusterConstExotics", "EMCAL cluster exotics cut vs jet pt; p{t, jet} (GeV/c); Cluster exotics", jetptbinning, exoticsbinning);
192  fQAHistos->CreateTH2("hClusterConstMinCellEnergy", "EMCAL Cluster const min cell energy; p{t, jet} (GeV/c); E_{cell} (GeV/c)", jetptbinning, cellenergybinning);
193  fQAHistos->CreateTH2("hClusterConstMaxCellEnergy", "EMCAL Cluster const max (seed) cell energy; p{t, jet} (GeV/c); E_{cell} (GeV/c)", jetptbinning, cellenergybinning);
194 
195  // Test of constituent QA
196 #ifdef EXPERIMENTAL_JETCONSTITUENTS
197  fQAHistos->CreateTH2("hChargedConstituentPt", "charged constituent pt vs jet pt (via constituent map); p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
198  fQAHistos->CreateTH2("hChargedIndexPt", "charged constituent pt vs jet pt (via index map); p_{t, jet} (GeV/c); p_{t, ch} (GeV/c)", jetptbinning, clusterenergybinning);
199 
200  fQAHistos->CreateTH2("hClusterConstituentEDefault", "cluster constituent default energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
201  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);
202  fQAHistos->CreateTH2("hClusterConstituentEHC", "cluster constituent hadronic-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
203  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);
204  fQAHistos->CreateTH2("hClusterIndexEHC", "cluster constituent hadronic-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
205  fQAHistos->CreateTH2("hLeadingChargedConstituentPt", "Pt of the leading charged constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
206  fQAHistos->CreateTH2("hLeadingClusterConstituentPt", "Pt of the leading cluster constituent in jet; p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
207 #endif
208  for(auto h : *(fQAHistos->GetListOfHistograms())) fOutput->Add(h);
209 
210  OpenFile(2);
211  std::string treename = this->GetOutputSlot(2)->GetContainer()->GetName();
212  fJetSubstructureTree = new TTree(treename.data(), "Tree with jet substructure information");
213 
215  fGlobalTreeParams->LinkJetTreeBranches(fJetSubstructureTree, fFillRho);
218  if(fFillPart) {
221  }
222  if(fFillSoftDrop) {
225  if(fFillPart) {
228  }
229  }
230  if(fFillNSub) {
233  if(fFillPart) {
236  }
237  }
238 
239  if(fFillStructGlob){
242  if(fFillPart){
245  }
246  }
247 
248  PostData(1, fOutput);
249  PostData(2, fJetSubstructureTree);
250 }
251 
255  }
256 }
257 
261  AliParticleContainer *particles = GetParticleContainer("mcparticles");
262 
263  AliJetContainer *mcjets = GetJetContainer("mcjets");
264  AliJetContainer *datajets = GetJetContainer("datajets");
265 
266  FillLuminosity(); // Makes only sense in data
267 
268  // for(auto e : *(fInputEvent->GetList())) std::cout << e->GetName() << std::endl;
269 
270  std::stringstream rhoTagData, rhoTagMC;
271  if(datajets) rhoTagData << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(datajets->GetJetRadius() * 10.);
272  if(mcjets) rhoTagMC << "R" << std::setw(2) << std::setfill('0') << static_cast<Int_t>(mcjets->GetJetRadius() * 10.);
273 
274  if(fFillRho){
275  std::string rhoSparseData = "RhoSparse_Full_" + rhoTagData.str(), rhoSparseMC = "RhoSparse_Full_" + rhoTagMC.str(),
276  rhoMassData = "RhoMassSparse_Full_" + rhoTagData.str(), rhoMassMC = "RhoMassSparse_Full_" + rhoTagMC.str();
277  AliRhoParameter *rhoPtRec = GetRhoFromEvent(rhoSparseData.data()),
278  *rhoMassRec = GetRhoFromEvent(rhoMassData.data()),
279  *rhoPtSim = GetRhoFromEvent(rhoSparseMC.data()),
280  *rhoMassSim = GetRhoFromEvent(rhoMassMC.data());
281  AliDebugStream(2) << "Found rho parameter for reconstructed pt: " << (rhoPtRec ? "yes" : "no") << ", value: " << (rhoPtRec ? rhoPtRec->GetVal() : 0.) << std::endl;
282  AliDebugStream(2) << "Found rho parameter for sim pt: " << (rhoPtSim ? "yes" : "no") << ", value: " << (rhoPtSim ? rhoPtSim->GetVal() : 0.) << std::endl;
283  AliDebugStream(2) << "Found rho parameter for reconstructed Mass: " << (rhoMassRec ? "yes" : "no") << ", value: " << (rhoMassRec ? rhoMassRec->GetVal() : 0.) << std::endl;
284  AliDebugStream(2) << "Found rho parameter for sim Mass: " << (rhoMassSim ? "yes" : "no") << ", value: " << (rhoMassSim ? rhoMassSim->GetVal() : 0.) << std::endl;
285  Double_t rhopars[4] = {
286  rhoPtRec ? rhoPtRec->GetVal() : 0.,
287  rhoPtSim ? rhoPtSim->GetVal() : 0.,
288  rhoMassRec ? rhoMassRec->GetVal() : 0.,
289  rhoMassSim ? rhoMassSim->GetVal() : 0.
290  };
291  memcpy(this->fGlobalTreeParams->fRhoParamters, rhopars, sizeof(Double_t) * 4);
292  }
293 
294  AliDebugStream(1) << "Inspecting jet radius " << (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius()) << std::endl;
295  this->fGlobalTreeParams->fJetRadius = (datajets ? datajets->GetJetRadius() : mcjets->GetJetRadius());
296  fGlobalTreeParams->fTriggerClusterIndex = -1; // Reset trigger cluster index
297 
298  if(datajets && !mcjets){
299  // decode trigger string in order to determine the trigger clusters
300  std::vector<std::string> clusternames;
301  auto triggerinfos = PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data());
302  for(auto t : triggerinfos) {
303  if(std::find(clusternames.begin(), clusternames.end(), t.Triggercluster()) == clusternames.end()) clusternames.emplace_back(t.Triggercluster());
304  }
305  bool isCENT = (std::find(clusternames.begin(), clusternames.end(), "CENT") != clusternames.end()),
306  isCENTNOTRD = (std::find(clusternames.begin(), clusternames.end(), "CENTNOTRD") != clusternames.end()),
307  isCALO = (std::find(clusternames.begin(), clusternames.end(), "CALO") != clusternames.end()),
308  isCALOFAST = (std::find(clusternames.begin(), clusternames.end(), "CALOFAST") != clusternames.end());
309  if(isCENT || isCENTNOTRD) {
310  if(isCENT && isCENTNOTRD) fGlobalTreeParams->fTriggerClusterIndex = 0; // CENTBOTH
311  else if(isCENT) fGlobalTreeParams->fTriggerClusterIndex = 1; // OnlyCENT
312  else fGlobalTreeParams->fTriggerClusterIndex = 2; // OnlyCENTNOTRD
313  }
314  if(isCALO || isCALOFAST){
315  // CALO(FAST) and CENT(NOTRD) clusters disjunct - CALO cluster not read out in parallel with CENT cluster
316  // Therefore mixing between CALO and CENT triggers doesn't need to be handled.
317  if(isCALO && isCALOFAST) fGlobalTreeParams->fTriggerClusterIndex = 3;
318  else if(isCALO) fGlobalTreeParams->fTriggerClusterIndex = 4;
319  else if(isCALOFAST) fGlobalTreeParams->fTriggerClusterIndex = 5;
320  }
321  }
322 
323  double weight = 1.;
325  AliDebugStream(2) << "Trigger selection string: " << fTriggerSelectionString << std::endl;
326  TString selectionString = (fTriggerSelectionBits & AliVEvent::kINT7) ? "INT7" : fTriggerSelectionString;
327  auto triggerstring = MatchTrigger(selectionString.Data());
328  AliDebugStream(2) << "Getting downscale correction factor for trigger string " << triggerstring << std::endl;
330  }
331  AliDebugStream(1) << "Using downscale weight " << weight << std::endl;
332  this->fGlobalTreeParams->fEventWeight = weight;
333 
334 
335  // Count events (for spectrum analysis)
336  fQAHistos->FillTH1("hEventCounter", 1);
337  fQAHistos->FillTH1("hTriggerClusterCounter", fGlobalTreeParams->fTriggerClusterIndex);
338 
339  AliSoftdropDefinition softdropSettings;
340  softdropSettings.fBeta = fSDBetaCut;
341  softdropSettings.fZ = fSDZCut;
342  switch(fReclusterizer) {
343  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
344  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
345  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
346  };
347 
348  AliNSubjettinessDefinition nsubjettinessSettings;
349  nsubjettinessSettings.fBeta = 1.;
350  nsubjettinessSettings.fRadius = 0.4;
351 
352  if(datajets) {
353  AliDebugStream(1) << "In data jets branch: found " << datajets->GetNJets() << " jets, " << datajets->GetNAcceptedJets() << " were accepted\n";
354  AliDebugStream(1) << "Having MC information: " << (mcjets ? TString::Format("yes, with %d jets", mcjets->GetNJets()) : "no") << std::endl;
355  if(mcjets) {
356  AliDebugStream(1) << "In MC jets branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
357  }
358  for(auto jet : datajets->accepted()) {
359  double pt = jet->Pt(), pz = jet->Pz(), E = jet->E(), M = TMath::Sqrt(E*E - pt*pt - pz*pz);
360  AliDebugStream(2) << "Next jet: pt:" << jet->Pt() << ", E: " << E << ", pz: " << pz << ", M(self): " << M << "M(fj)" << jet->M() << std::endl;
361  AliEmcalJet *associatedJet = jet->ClosestJet();
362 
363  if(mcjets) {
364  if(!associatedJet) {
365  AliDebugStream(2) << "Not found associated jet" << std::endl;
366  continue;
367  }
368  if(!(SelectJet(*jet, tracks) && SelectJet(*associatedJet, particles))) continue;
369  try {
370  DoConstituentQA(jet, tracks, clusters);
371  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, datajets->GetJetRadius() * 2., tracks, clusters, {softdropSettings, nsubjettinessSettings}),
372  structureMC = fFillPart ? MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings}) : AliJetSubstructureData();
373  if(fKineRec) *fKineRec = MakeJetKineParameters(*jet, kDetLevel, tracks, clusters);
374  if(fKineSim) *fKineSim = MakeJetKineParameters(*associatedJet, kPartLevel, particles, nullptr);
375  if(fSoftDropMeasured) *fSoftDropMeasured = structureData.fSoftDrop;
376  if(fSoftDropTrue) *fSoftDropTrue = structureMC.fSoftDrop;
377  if(fNSubMeasured) *fNSubMeasured = structureData.fNsubjettiness;
378  if(fNSubTrue) *fNSubTrue = structureMC.fNsubjettiness;
379  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
380  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*associatedJet, particles, nullptr), MakePtD(*associatedJet, particles, nullptr)};
381  fJetSubstructureTree->Fill();
382  } catch(ReclusterizerException &e) {
383  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
384  } catch(SubstructureException &e) {
385  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
386  }
387  } else {
388  if(!SelectJet(*jet, tracks)) continue;
389  try {
390  DoConstituentQA(jet, tracks, clusters);
391  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
392  if(fKineRec) *fKineRec = MakeJetKineParameters(*jet, kDetLevel, tracks, clusters);
394  if(fNSubMeasured) *fNSubMeasured = structure.fNsubjettiness;
395  if(fJetStructureMeasured) *fJetStructureMeasured = {MakeAngularity(*jet, tracks, clusters), MakePtD(*jet, tracks, clusters)};
396  fJetSubstructureTree->Fill();
397  } catch(ReclusterizerException &e) {
398  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
399  } catch(SubstructureException &e) {
400  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
401  }
402  }
403  }
404  } else {
405  if(mcjets) {
406  // for MCgen train
407  AliDebugStream(1) << "In MC pure jet branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
408  for(auto j : mcjets->accepted()){
409  AliEmcalJet *mcjet = static_cast<AliEmcalJet *>(j);
410  try {
411  AliJetSubstructureData structure = MakeJetSubstructure(*mcjet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
412  if(this->fKineSim) *fKineSim = MakeJetKineParameters(*mcjet, kPartLevel, particles, nullptr);
413  if(fSoftDropTrue) *fSoftDropTrue = structure.fSoftDrop;
414  if(fNSubTrue) *fNSubTrue = structure.fNsubjettiness;
415  if(fJetStructureTrue) *fJetStructureTrue = {MakeAngularity(*mcjet, particles, nullptr), MakePtD(*mcjet, particles, nullptr)};
416  fJetSubstructureTree->Fill();
417  } catch (ReclusterizerException &e) {
418  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
419  } catch (SubstructureException &e) {
420  AliErrorStream() << "Error in substructure observable - skipping jet" << std::endl;
421  }
422  }
423  }
424  }
425 
426  return true;
427 }
428 
430  AliDebugStream(1) << "Trigger selection called\n";
431  if(!(fHasRecEvent || fHasTrueEvent)){
432  AliErrorStream() << "Impossible combination: Neither rec nor true event available. Rejecting ..." << std::endl;
433  return false;
434  }
435  // Run trigger selection (not on pure MCgen train - pure MCgen train has no rec event, ESD event is fake there)
436  if(fHasRecEvent){
437  if(!fHasTrueEvent){
438  // Pure data - do EMCAL trigger selection from selection string
439  AliDebugStream(1) << "Applying trigger selection for trigger bits " << std::bitset<sizeof(decltype(fTriggerSelectionBits)) * 8>(fTriggerSelectionBits) << "and trigger selection string " << fTriggerSelectionString << std::endl;
440  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
441  AliDebugStream(1) << "Passed trigger bit selection" << std::endl;
442  if(fTriggerSelectionString.Length()) {
443  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
444  AliDebugStream(1) << "Passed trigger string section" << std::endl;
446  auto trgselresult = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
447  AliDebugStream(1) << "Found trigger decision object: " << (trgselresult ? "yes" : "no") << std::endl;
448  if(!trgselresult){
449  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
450  return false;
451  }
452  if(!trgselresult->IsEventSelected(fTriggerSelectionString)) return false;
453  AliDebugStream(1) << "Data event selected" << std::endl;
454  }
455  }
456  } else {
457  // Simulation - do EMCAL trigger selection from trigger selection object
458  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false; // Require INT7 trigger - EMCAL triggers will be a subset
460  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
461  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
462  if(!mctrigger){
463  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
464  return false;
465  }
466  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
467  }
468  }
469  }
470  return true; // trigger selected or pure MCgen information
471 }
472 
474  AliCDBManager * cdb = AliCDBManager::Instance();
475  if(!fMCEvent && cdb){
476  // Get List of trigger clusters
477  AliCDBEntry *en = cdb->Get("GRP/CTP/Config");
478  AliTriggerConfiguration *trg = static_cast<AliTriggerConfiguration *>(en->GetObject());
479  std::vector<std::string> clusternames;
480  for(auto c : trg->GetClusters()) {
481  AliTriggerCluster *clust = static_cast<AliTriggerCluster *>(c);
482  std::string clustname = clust->GetName();
483  auto iscent = clustname.find("CENT") != std::string::npos, iscalo = clustname.find("CALO") != std::string::npos;
484  if(!(iscalo || iscent)) continue;
485  AliInfoStream() << "Adding trigger cluster " << clustname << " to cluster lumi monitor" << std::endl;
486  clusternames.emplace_back(clustname);
487  }
488 
489  // Set the x-axis of the luminosity monitor histogram
490  fLumiMonitor = new TH1F("hLumiMonitor", "Luminosity monitor", clusternames.size(), 0, clusternames.size());
491  int currentbin(1);
492  for(auto c : clusternames) {
493  fLumiMonitor->GetXaxis()->SetBinLabel(currentbin++, c.data());
494  }
495  fOutput->Add(fLumiMonitor);
496  }
497 }
498 
501  auto downscalefactors = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance();
502  if(fInputEvent->GetFiredTriggerClasses().Contains("INT7")) {
503  for(auto trigger : PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data())){
504  auto int7trigger = trigger.IsTriggerClass("INT7");
505  auto bunchcrossing = trigger.BunchCrossing() == "B";
506  auto nopf = trigger.PastFutureProtection() == "NOPF";
507  bool centcalo = (trigger.Triggercluster().find("CENT") != std::string::npos) || (trigger.Triggercluster().find("CALO") != std::string::npos);
508  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;
509  if(int7trigger && bunchcrossing && nopf && centcalo) {
510  double downscale = downscalefactors->GetDownscaleFactorForTriggerClass(trigger.ExpandClassName());
511  AliDebugStream(5) << "Using downscale " << downscale << std::endl;
512  fLumiMonitor->Fill(trigger.Triggercluster().data(), 1./downscale);
513  }
514  }
515  }
516  }
517 }
518 
520  AliJetKineParameters result;
521  result.fPt = TMath::Abs(jet.Pt());
522  result.fE = jet.E();
523  result.fEta = jet.Eta();
524  result.fPhi = jet.Phi();
525  result.fArea = jet.Area();
526  result.fMass = jet.M();
527  result.fNEF = jet.NEF();
528  result.fNCharged = jet.GetNumberOfTracks();
529  result.fNNeutral = jet.GetNumberOfClusters();
530  std::vector<double> zcharged, zneutral;
531  if(tracks) {
532  // Find the leading track
533  for(auto icharged = 0; icharged < jet.GetNumberOfTracks(); icharged++){
534  auto trk = jet.TrackAt(icharged, tracks->GetArray());
535  bool charged = true;
536  if(rectype == kPartLevel){
537  if(!trk->Charge()) charged = false;
538  }
539  auto z = jet.GetZ(trk);
540  if(charged) zcharged.push_back(z);
541  else zneutral.push_back(z);
542  }
543  }
544  if(clusters) {
545  for(auto iclust = 0; iclust < jet.GetNumberOfClusters(); iclust++){
546  auto clust = jet.ClusterAt(iclust, clusters->GetArray());
547  TLorentzVector clustervec;
548  clust->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
549  auto z = jet.GetZ(clustervec.Px(), clustervec.Py(), clustervec.Pz());
550  zneutral.push_back(z);
551  }
552  }
553  result.fZLeading = -1.;
554  result.fZLeadingCharged = -1.;
555  result.fZLeadingNeutral = -1;
556  if(zcharged.size()) {
557  std::sort(zcharged.begin(), zcharged.end(), std::greater<double>());
558  result.fZLeadingCharged = zcharged[0];
559  result.fZLeading = result.fZLeadingCharged;
560  }
561  if(zneutral.size()){
562  std::sort(zneutral.begin(), zneutral.end(), std::greater<double>());
563  result.fZLeadingNeutral = zneutral[0];
564  if(result.fZLeadingNeutral > result.fZLeading) result.fZLeading = result.fZLeadingNeutral;
565  }
566  return result;
567 }
568 
570  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
571  std::vector<fastjet::PseudoJet> constituents;
572  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
573  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
574  if(tracks && (fUseChargedConstituents || isMC)){ // Neutral particles part of particle container in case of MC
575  AliDebugStream(1) << "Jet substructure: Using charged constituents" << std::endl;
576  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
577  auto track = jet.TrackAt(itrk, tracks->GetArray());
578  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
579  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
580  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
581  constituentTrack.set_user_index(jet.TrackAt(itrk));
582  constituents.push_back(constituentTrack);
583  }
584  }
585 
586  if(clusters && fUseNeutralConstituents){
587  AliDebugStream(1) << "Jet substructure: Using neutral constituents" << std::endl;
588  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
589  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
590  TLorentzVector clustervec;
591  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
592  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
593  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
594  constituents.push_back(constituentCluster);
595  }
596  }
597 
598  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
599  if(!constituents.size()){
600  AliErrorStream() << "Jet has 0 constituents." << std::endl;
601  throw ReclusterizerException();
602  }
603  // Redo jet finding on constituents with a
604  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
605  std::vector<fastjet::PseudoJet> outputjets;
606  try {
607  fastjet::ClusterSequence jetfinder(constituents, jetdef);
608  outputjets = jetfinder.inclusive_jets(0);
610  return result;
611  } catch (fastjet::Error &e) {
612  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
613  throw ReclusterizerException();
614  } catch (SoftDropException &e) {
615  AliErrorStream() << "Softdrop exception caught: " << e.what() << std::endl;
616  throw ReclusterizerException();
617  }
618 }
619 
621  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
622  softdropAlgorithm.set_verbose_structure(kTRUE);
623 #if FASTJET_VERSION_NUMBER >= 30302
624  fastjet::Recluster reclusterizer(cutparameters.fRecluserAlgo, 1, fastjet::Recluster::keep_only_hardest);
625 #else
626  fastjet::contrib::Recluster reclusterizer(cutparameters.fRecluserAlgo, 1, true);
627 #endif
628  softdropAlgorithm.set_reclustering(kTRUE, &reclusterizer);
629  AliDebugStream(4) << "Jet has " << jet.constituents().size() << " constituents" << std::endl;
630  auto groomed = softdropAlgorithm(jet);
631  try {
632  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
633 
634  AliSoftDropParameters result({softdropstruct.symmetry(),
635  groomed.m(),
636  softdropstruct.delta_R(),
637  groomed.perp(),
638  softdropstruct.delta_R(),
639  softdropstruct.mu(),
640  softdropstruct.dropped_count()});
641  return result;
642  } catch(std::bad_cast &e) {
643  throw SoftDropException();
644  }
645 }
646 
649  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet),
650  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet)
651  });
652  return result;
653 }
654 
656  if(!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
657  throw SubstructureException();
658  TVector3 jetvec(jet.Px(), jet.Py(), jet.Pz());
659  Double_t den(0.), num(0.);
660  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
661  if(tracks && (fUseChargedConstituents || isMC)){
662  AliDebugStream(1) << "Angularity: Using charged constituents" << std::endl;
663  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
664  auto track = jet.TrackAt(itrk, tracks->GetArray());
665  if(!track){
666  AliErrorStream() << "Associated constituent particle / track not found\n";
667  continue;
668  }
669  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
670  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
671  TVector3 trackvec(track->Px(), track->Py(), track->Pz());
672 
673  num += track->Pt() * trackvec.DrEtaPhi(jetvec);
674  den += +track->Pt();
675  }
676  }
677  if(clusters && fUseNeutralConstituents) {
678  AliDebugStream(1) << "Using neutral constituents" << std::endl;
679  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
680  auto clust = jet.ClusterAt(icl, clusters->GetArray());
681  if(!clust) {
682  AliErrorStream() << "Associated constituent cluster not found\n";
683  continue;
684  }
685  TLorentzVector clusterp;
686  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
687 
688  num += clusterp.Pt() * clusterp.Vect().DrEtaPhi(jetvec);
689  den += clusterp.Pt();
690  }
691  }
692  return num/den;
693 }
694 
696  if (!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
697  throw SubstructureException();
698  Double_t den(0.), num(0.);
699  bool isMC = dynamic_cast<const AliMCParticleContainer *>(particles);
700  if(particles && (fUseChargedConstituents || isMC)){
701  AliDebugStream(1) << "Using charged constituents" << std::endl;
702  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
703  auto trk = jet.TrackAt(itrk, particles->GetArray());
704  if(!trk){
705  AliErrorStream() << "Associated constituent particle / track not found\n";
706  continue;
707  }
708  if(!trk->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
709  if(trk->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
710  num += trk->Pt() * trk->Pt();
711  den += trk->Pt();
712  }
713  }
714  if(clusters && fUseNeutralConstituents){
715  AliDebugStream(1) << "Using neutral constituents" << std::endl;
716  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
717  auto clust = jet.ClusterAt(icl, clusters->GetArray());
718  if(!clust) {
719  AliErrorStream() << "Associated constituent cluster not found\n";
720  continue;
721  }
722  TLorentzVector clusterp;
723  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
724  num += clusterp.Pt() * clusterp.Pt();
725  den += clusterp.Pt();
726  }
727  }
728  return TMath::Sqrt(num)/den;
729 }
730 
732  for(int icl = 0; icl < jet->GetNumberOfClusters(); icl++){
733  auto clust = jet->ClusterAt(icl, clusters->GetArray());
734  AliDebugStream(3) << "cluster time " << clust->GetTOF() << std::endl;
735  fQAHistos->FillTH2("hClusterConstE", jet->Pt(),clust->GetUserDefEnergy(clusters->GetDefaultClusterEnergy()));
736  fQAHistos->FillTH2("hClusterConstTime", jet->Pt(), clust->GetTOF()*1e9); // convert to nanoseconds
737  fQAHistos->FillTH2("hClusterConstM02", jet->Pt(), clust->GetM02());
738  fQAHistos->FillTH2("hClusterConstNcell", jet->Pt(), clust->GetNCells());
739  fQAHistos->FillTH2("hClusterConstExotics", jet->Pt(), clust->GetIsExotic() ? 1. : 0.);
740 
741  double mincell(100000.), maxcell(0.);
742  for(int icell = 0; icell < clust->GetNCells(); icell++){
743  double ecell = clust->E() * clust->GetCellAmplitudeFraction(icell);
744  if(ecell < mincell) mincell = ecell;
745  if(ecell > maxcell) maxcell = ecell;
746  }
747  fQAHistos->FillTH2("hClusterConstMinCellEnergy", jet->Pt(), mincell);
748  fQAHistos->FillTH2("hClusterConstMaxCellEnergy", jet->Pt(), maxcell);
749 
750 #ifdef EXPERIMENTAL_JETCONSTITUENTS
751  fQAHistos->FillTH2("hClusterIndexENLC", jet->Pt(), clust->GetNonLinCorrEnergy());
752  fQAHistos->FillTH2("hClusterIndexEHC", jet->Pt(), clust->GetHadCorrEnergy());
753 #endif
754  }
755 
756 #ifdef EXPERIMENTAL_JETCONSTITUENTS
757  // Loop over charged particles - fill test histogram
758  for(int itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
759  auto part = jet->TrackAt(itrk, cont->GetArray());
760  fQAHistos->FillTH2("hChargedIndexPt", jet->Pt(), part->Pt());
761  }
762 
763  // Look over charged constituents
764  AliDebugStream(2) << "Jet: Number of particle constituents: " << jet->GetParticleConstituents().size() << std::endl;
765  for(auto part : jet->GetParticleConstituents()) {
766  //auto part = static_cast<PWG::JETFW::AliEmcalParticleJetConstituent *>(pconst);
767  AliDebugStream(3) << "Found particle constituent with pt " << part.Pt() << ", from VParticle " << part.GetParticle()->Pt() << std::endl;
768  fQAHistos->FillTH2("hChargedConstituentPt", jet->Pt(), part.Pt());
769  }
770 
771  // Loop over neutral constituents
772  AliDebugStream(2) << "Jet: Number of cluster constituents: " << jet->GetClusterConstituents().size() << std::endl;
773  for(auto clust : jet->GetClusterConstituents()){
774  //auto clust = static_cast<PWG::JETFW::AliEmcalClusterJetConstituent *>(cconst);
775  AliDebugStream(3) << "Found cluster constituent with energy " << clust.E() << " using energy definition " << static_cast<int>(clust.GetDefaultEnergyType()) << std::endl;
776  fQAHistos->FillTH2("hClusterConstituentEDefault", jet->Pt(), clust.E());
777  fQAHistos->FillTH2("hClusterConstituentENLC", jet->Pt(), clust.GetCluster()->GetNonLinCorrEnergy());
778  fQAHistos->FillTH2("hClusterConstituentEHC", jet->Pt(), clust.GetCluster()->GetHadCorrEnergy());
779  }
780 
781  // Fill global observables: Leading charged and cluster constituents
782  auto leadingcharged = jet->GetLeadingParticleConstituent();
783  auto leadingcluster = jet->GetLeadingClusterConstituent();
784  if(leadingcluster){
785  fQAHistos->FillTH1("hLeadingClusterConstituentPt", jet->Pt(), leadingcluster->GetCluster()->GetHadCorrEnergy());
786  }
787  if(leadingcharged) {
788  fQAHistos->FillTH1("hLeadingChargedConstituentPt", jet->Pt(), leadingcharged->GetParticle()->Pt());
789  }
790 #endif
791 }
792 
794  int ncharged = 0, nneutral = jet.GetNumberOfClusters();
795  if(particles) {
796  for(decltype(jet.GetNumberOfTracks()) ipart = 0; ipart < jet.GetNumberOfTracks(); ipart++){
797  auto part = jet.TrackAt(ipart, particles->GetArray());
798  if(!part) continue;
799  if(part->Charge()) ncharged++;
800  else nneutral++;
801  }
802  }
803  // check if the jet has at least one consituent for jet substructure
804  int nallowed = 0;
805  nallowed += fUseChargedConstituents ? ncharged : 0;
806  nallowed += fUseNeutralConstituents ? nneutral : 0;
807  return nallowed > 0;
808 }
809 
810 std::string AliAnalysisTaskEmcalJetSubstructureTree::MatchTrigger(const std::string &triggertoken) const {
811  std::vector<std::string> tokens;
812  std::string result;
813  std::stringstream decoder(fInputEvent->GetFiredTriggerClasses().Data());
814  while(std::getline(decoder, result, ' ')) tokens.emplace_back(result);
815  result.clear();
816  for(auto t : tokens) {
817  if(t.find(triggertoken) != std::string::npos) {
818  // take first occurrence - downscale factor should normally be the same
819  result = t;
820  break;
821  }
822  }
823  return result;
824 }
825 
826 bool AliAnalysisTaskEmcalJetSubstructureTree::IsSelectEmcalTriggers(const std::string &triggerstring) const {
827  const std::array<std::string, 8> kEMCALTriggers = {
828  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
829  };
830  bool isEMCAL = false;
831  for(auto emcaltrg : kEMCALTriggers) {
832  if(triggerstring.find(emcaltrg) != std::string::npos) {
833  isEMCAL = true;
834  break;
835  }
836  }
837  return isEMCAL;
838 }
839 
841  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
842 
843  Bool_t isAOD(kFALSE);
844  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
845  if(inputhandler) {
846  if(inputhandler->IsA() == AliAODInputHandler::Class()){
847  std::cout << "Analysing AOD events\n";
848  isAOD = kTRUE;
849  } else {
850  std::cout << "Analysing ESD events\n";
851  }
852  }
853 
854  std::stringstream taskname;
855  taskname << "JetSubstructureTreemaker_R" << std::setw(2) << std::setfill('0') << int(jetradius*10) << trigger;
857  mgr->AddTask(treemaker);
858  treemaker->SetMakeGeneralHistograms(kTRUE);
859  if(isMC) treemaker->SetHasTrueEvent(true);
860  if(isData) treemaker->SetHasRecEvent(true);
861 
862  // Adding containers
863  if(isMC) {
864  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
865  particles->SetMinPt(0.);
866 
867  AliJetContainer *mcjets = treemaker->AddJetContainer(
868  jettype,
870  recombinationScheme,
871  jetradius,
873  particles, nullptr);
874  mcjets->SetName("mcjets");
875  mcjets->SetJetPtCut(20.);
876  }
877 
878  if(isData) {
879  AliTrackContainer *tracks(nullptr);
880  if((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet)){
882  std::cout << "Track container name: " << tracks->GetName() << std::endl;
883  tracks->SetMinPt(0.15);
884  }
885  AliClusterContainer *clusters(nullptr);
886  if((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet)){
887  std::cout << "Using full or neutral jets ..." << std::endl;
889  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
890  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
891  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
892  } else {
893  std::cout << "Using charged jets ... " << std::endl;
894  }
895 
896  AliJetContainer *datajets = treemaker->AddJetContainer(
897  jettype,
899  recombinationScheme,
900  jetradius,
902  tracks, clusters);
903  datajets->SetName("datajets");
904  datajets->SetJetPtCut(20.);
905 
906  treemaker->SetUseAliAnaUtils(true, true);
907  treemaker->SetVzRange(-10., 10);
908 
909  // configure trigger selection
910  std::string triggerstring(trigger);
911  if(triggerstring.find("INT7") != std::string::npos) {
912  treemaker->SetTriggerBits(AliVEvent::kINT7);
913  } else if(triggerstring.find("EJ1") != std::string::npos) {
914  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
915  treemaker->SetTriggerString("EJ1");
916  } else if(triggerstring.find("EJ2") != std::string::npos) {
917  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
918  treemaker->SetTriggerString("EJ2");
919  } else if(triggerstring.find("EG1") != std::string::npos) {
920  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
921  treemaker->SetTriggerString("EG1");
922  } else if(triggerstring.find("EG2") != std::string::npos) {
923  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
924  treemaker->SetTriggerString("EG2");
925  }
926  }
927 
928  std::string jettypestring;
929  switch(jettype) {
930  case AliJetContainer::kFullJet: jettypestring = "FullJets"; break;
931  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; break;
932  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; break;
933  default: jettypestring = "Undef";
934  };
935 
936  // Connecting containers
937  std::stringstream outputfile, histname, treename;
938  outputfile << mgr->GetCommonFileName() << ":JetSubstructure_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
939  histname << "JetSubstructureHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
940  treename << "JetSubstructureTree_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
941  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
942  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
943  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer(treename.str().data(), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()));
944 
945  return treemaker;
946 }
947 
948 void AliSoftDropParameters::LinkJetTreeBranches(TTree *jettree, const char *tag) {
949  LinkBranch(jettree, &fZg, Form("Zg%s", tag), "D");
950  LinkBranch(jettree, &fRg, Form("Rg%s", tag), "D");
951  LinkBranch(jettree, &fMg, Form("Mg%s", tag), "D");
952  LinkBranch(jettree, &fPtg, Form("Ptg%s", tag), "D");
953  LinkBranch(jettree, &fMug, Form("Mug%s", tag), "D");
954  LinkBranch(jettree, &fDeltaR, Form("DeltaRg%s", tag), "D");
955  LinkBranch(jettree, &fNDropped, Form("NDropped%s", tag), "I");
956 };
957 
959  LinkBranch(jettree, &fOneSubjettiness, Form("OneSubjettiness%s", tag), "D");
960  LinkBranch(jettree, &fTwoSubjettiness, Form("TwoSubjettiness%s", tag), "D");
961 }
962 
964  LinkBranch(jettree, &fAngularity, Form("Angularity%s", tag), "D");
965  LinkBranch(jettree, &fPtD, Form("PtD%s", tag), "D");
966 }
967 
968 void AliJetKineParameters::LinkJetTreeBranches(TTree *jettree, const char *tag){
969  LinkBranch(jettree, &fPt, Form("PtJet%s", tag), "D");
970  LinkBranch(jettree, &fE, Form("EJet%s", tag), "D");
971  LinkBranch(jettree, &fEta, Form("Eta%s", tag), "D");
972  LinkBranch(jettree, &fPhi, Form("Phi%s", tag), "D");
973  LinkBranch(jettree, &fArea, Form("Area%s", tag), "D");
974  LinkBranch(jettree, &fMass, Form("Mass%s", tag), "D");
975  LinkBranch(jettree, &fNEF, Form("NEF%s", tag), "D");
976  LinkBranch(jettree, &fNCharged, Form("NCharged%s", tag), "I");
977  LinkBranch(jettree, &fNNeutral, Form("NNeutral%s", tag), "I");
978  LinkBranch(jettree, &fZLeading, Form("ZLeading%s", tag), "D");
979  LinkBranch(jettree, &fZLeadingCharged, Form("ZLeadingCharged%s", tag), "D");
980  LinkBranch(jettree, &fZLeadingNeutral, Form("ZLeadingNeutral%s", tag), "D");
981 }
982 
984  LinkBranch(jettree, &fJetRadius, "Radius", "D");
985  LinkBranch(jettree, &fEventWeight, "EventWeight", "D");
986  LinkBranch(jettree, &fTriggerClusterIndex, "TriggerClusterIndex", "I");
987  if(fillRho) {
988  std::string varnames[] = {"RhoPtRec", "RhoPtSim", "RhoMassRec", "RhoMassSim"};
989  for(int i = 0; i < 4; i++){
990  LinkBranch(jettree, fRhoParamters + i, varnames[i].data(), "D");
991  }
992  }
993 }
994 
995 void LinkBranch(TTree *jettree, void *data, const char *branchname, const char *type) {
996  jettree->Branch(branchname, data, Form("%s/%s", branchname, type));
997 }
998 } /* namespace EmcalTriggerJets */
Double_t fZLeadingNeutral
z of the leading neutral constituent
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Bool_t fHasRecEvent
Has reconstructed event (for trigger selection)
Double_t fZLeadingCharged
z of the leading charged constituent
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:341
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)
Set pre-configured event cut object.
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
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.
Double_t GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz) const
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
AliJetKineParameters MakeJetKineParameters(const AliEmcalJet &jet, JetRecType_t rectype, const AliParticleContainer *const particles, const AliClusterContainer *const clusters) const
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.