AliPhysics  abafffd (abafffd)
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();
367  if(fKineRec) *fKineRec = MakeJetKineParameters(*jet, kDetLevel, tracks, clusters);
368  if(fKineSim) *fKineSim = MakeJetKineParameters(*associatedJet, kPartLevel, particles, nullptr);
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});
386  if(fKineRec) *fKineRec = MakeJetKineParameters(*jet, kDetLevel, tracks, clusters);
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, kPartLevel, particles, nullptr);
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  AliInfoStream() << "Adding trigger cluster " << clustname << " to cluster lumi monitor" << std::endl;
480  clusternames.emplace_back(clustname);
481  }
482 
483  // Set the x-axis of the luminosity monitor histogram
484  fLumiMonitor = new TH1F("hLumiMonitor", "Luminosity monitor", clusternames.size(), 0, clusternames.size());
485  int currentbin(1);
486  for(auto c : clusternames) {
487  fLumiMonitor->GetXaxis()->SetBinLabel(currentbin++, c.data());
488  }
489  fOutput->Add(fLumiMonitor);
490  }
491 }
492 
495  auto downscalefactors = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance();
496  if(fInputEvent->GetFiredTriggerClasses().Contains("INT7")) {
497  for(auto trigger : PWG::EMCAL::Triggerinfo::DecodeTriggerString(fInputEvent->GetFiredTriggerClasses().Data())){
498  auto int7trigger = trigger.IsTriggerClass("INT7");
499  auto bunchcrossing = trigger.BunchCrossing() == "B";
500  auto nopf = trigger.PastFutureProtection() == "NOPF";
501  bool centcalo = (trigger.Triggercluster().find("CENT") != std::string::npos) || (trigger.Triggercluster().find("CALO") != std::string::npos);
502  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;
503  if(int7trigger && bunchcrossing && nopf && centcalo) {
504  double downscale = downscalefactors->GetDownscaleFactorForTriggerClass(trigger.ExpandClassName());
505  AliDebugStream(5) << "Using downscale " << downscale << std::endl;
506  fLumiMonitor->Fill(trigger.Triggercluster().data(), 1./downscale);
507  }
508  }
509  }
510  }
511 }
512 
514  AliJetKineParameters result;
515  result.fPt = TMath::Abs(jet.Pt());
516  result.fE = jet.E();
517  result.fEta = jet.Eta();
518  result.fPhi = jet.Phi();
519  result.fArea = jet.Area();
520  result.fMass = jet.M();
521  result.fNEF = jet.NEF();
522  result.fNCharged = jet.GetNumberOfTracks();
523  result.fNNeutral = jet.GetNumberOfClusters();
524  std::vector<double> zcharged, zneutral;
525  if(tracks) {
526  // Find the leading track
527  for(auto icharged = 0; icharged < jet.GetNumberOfTracks(); icharged++){
528  auto trk = jet.TrackAt(icharged, tracks->GetArray());
529  bool charged = true;
530  if(rectype == kPartLevel){
531  if(!trk->Charge()) charged = false;
532  }
533  auto z = jet.GetZ(trk);
534  if(charged) zcharged.push_back(z);
535  else zneutral.push_back(z);
536  }
537  }
538  if(clusters) {
539  for(auto iclust = 0; iclust < jet.GetNumberOfClusters(); iclust++){
540  auto clust = jet.ClusterAt(iclust, clusters->GetArray());
541  TLorentzVector clustervec;
542  clust->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
543  auto z = jet.GetZ(clustervec.Px(), clustervec.Py(), clustervec.Pz());
544  zneutral.push_back(z);
545  }
546  }
547  result.fZLeading = -1.;
548  result.fZLeadingCharged = -1.;
549  result.fZLeadingNeutral = -1;
550  if(zcharged.size()) {
551  std::sort(zcharged.begin(), zcharged.end(), std::greater<double>());
552  result.fZLeadingCharged = zcharged[0];
553  result.fZLeading = result.fZLeadingCharged;
554  }
555  if(zneutral.size()){
556  std::sort(zneutral.begin(), zneutral.end(), std::greater<double>());
557  result.fZLeadingNeutral = zneutral[0];
558  if(result.fZLeadingNeutral > result.fZLeading) result.fZLeading = result.fZLeadingNeutral;
559  }
560  return result;
561 }
562 
564  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
565  std::vector<fastjet::PseudoJet> constituents;
566  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
567  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
568  if(tracks && (fUseChargedConstituents || isMC)){ // Neutral particles part of particle container in case of MC
569  AliDebugStream(1) << "Jet substructure: Using charged constituents" << std::endl;
570  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
571  auto track = jet.TrackAt(itrk, tracks->GetArray());
572  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
573  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
574  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
575  constituentTrack.set_user_index(jet.TrackAt(itrk));
576  constituents.push_back(constituentTrack);
577  }
578  }
579 
580  if(clusters && fUseNeutralConstituents){
581  AliDebugStream(1) << "Jet substructure: Using neutral constituents" << std::endl;
582  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
583  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
584  TLorentzVector clustervec;
585  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
586  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
587  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
588  constituents.push_back(constituentCluster);
589  }
590  }
591 
592  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
593  if(!constituents.size()){
594  AliErrorStream() << "Jet has 0 constituents." << std::endl;
595  throw ReclusterizerException();
596  }
597  // Redo jet finding on constituents with a
598  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
599  std::vector<fastjet::PseudoJet> outputjets;
600  try {
601  fastjet::ClusterSequence jetfinder(constituents, jetdef);
602  outputjets = jetfinder.inclusive_jets(0);
604  return result;
605  } catch (fastjet::Error &e) {
606  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
607  throw ReclusterizerException();
608  } catch (SoftDropException &e) {
609  AliErrorStream() << "Softdrop exception caught: " << e.what() << std::endl;
610  throw ReclusterizerException();
611  }
612 }
613 
615  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
616  softdropAlgorithm.set_verbose_structure(kTRUE);
617  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
618  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
619  AliDebugStream(4) << "Jet has " << jet.constituents().size() << " constituents" << std::endl;
620  auto groomed = softdropAlgorithm(jet);
621  try {
622  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
623 
624  AliSoftDropParameters result({softdropstruct.symmetry(),
625  groomed.m(),
626  softdropstruct.delta_R(),
627  groomed.perp(),
628  softdropstruct.delta_R(),
629  softdropstruct.mu(),
630  softdropstruct.dropped_count()});
631  return result;
632  } catch(std::bad_cast &e) {
633  throw SoftDropException();
634  }
635 }
636 
639  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet),
640  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet)
641  });
642  return result;
643 }
644 
646  if(!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
647  throw SubstructureException();
648  TVector3 jetvec(jet.Px(), jet.Py(), jet.Pz());
649  Double_t den(0.), num(0.);
650  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
651  if(tracks && (fUseChargedConstituents || isMC)){
652  AliDebugStream(1) << "Angularity: Using charged constituents" << std::endl;
653  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
654  auto track = jet.TrackAt(itrk, tracks->GetArray());
655  if(!track){
656  AliErrorStream() << "Associated constituent particle / track not found\n";
657  continue;
658  }
659  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
660  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
661  TVector3 trackvec(track->Px(), track->Py(), track->Pz());
662 
663  num += track->Pt() * trackvec.DrEtaPhi(jetvec);
664  den += +track->Pt();
665  }
666  }
667  if(clusters && fUseNeutralConstituents) {
668  AliDebugStream(1) << "Using neutral constituents" << std::endl;
669  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
670  auto clust = jet.ClusterAt(icl, clusters->GetArray());
671  if(!clust) {
672  AliErrorStream() << "Associated constituent cluster not found\n";
673  continue;
674  }
675  TLorentzVector clusterp;
676  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
677 
678  num += clusterp.Pt() * clusterp.Vect().DrEtaPhi(jetvec);
679  den += clusterp.Pt();
680  }
681  }
682  return num/den;
683 }
684 
686  if (!(jet.GetNumberOfTracks() || jet.GetNumberOfClusters()))
687  throw SubstructureException();
688  Double_t den(0.), num(0.);
689  bool isMC = dynamic_cast<const AliMCParticleContainer *>(particles);
690  if(particles && (fUseChargedConstituents || isMC)){
691  AliDebugStream(1) << "Using charged constituents" << std::endl;
692  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
693  auto trk = jet.TrackAt(itrk, particles->GetArray());
694  if(!trk){
695  AliErrorStream() << "Associated constituent particle / track not found\n";
696  continue;
697  }
698  if(!trk->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
699  if(trk->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
700  num += trk->Pt() * trk->Pt();
701  den += trk->Pt();
702  }
703  }
704  if(clusters && fUseNeutralConstituents){
705  AliDebugStream(1) << "Using neutral constituents" << std::endl;
706  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
707  auto clust = jet.ClusterAt(icl, clusters->GetArray());
708  if(!clust) {
709  AliErrorStream() << "Associated constituent cluster not found\n";
710  continue;
711  }
712  TLorentzVector clusterp;
713  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
714  num += clusterp.Pt() * clusterp.Pt();
715  den += clusterp.Pt();
716  }
717  }
718  return TMath::Sqrt(num)/den;
719 }
720 
722  for(int icl = 0; icl < jet->GetNumberOfClusters(); icl++){
723  auto clust = jet->ClusterAt(icl, clusters->GetArray());
724  AliDebugStream(3) << "cluster time " << clust->GetTOF() << std::endl;
725  fQAHistos->FillTH2("hClusterConstE", jet->Pt(),clust->GetUserDefEnergy(clusters->GetDefaultClusterEnergy()));
726  fQAHistos->FillTH2("hClusterConstTime", jet->Pt(), clust->GetTOF()*1e9); // convert to nanoseconds
727  fQAHistos->FillTH2("hClusterConstM02", jet->Pt(), clust->GetM02());
728  fQAHistos->FillTH2("hClusterConstNcell", jet->Pt(), clust->GetNCells());
729  fQAHistos->FillTH2("hClusterConstExotics", jet->Pt(), clust->GetIsExotic() ? 1. : 0.);
730 
731  double mincell(100000.), maxcell(0.);
732  for(int icell = 0; icell < clust->GetNCells(); icell++){
733  double ecell = clust->E() * clust->GetCellAmplitudeFraction(icell);
734  if(ecell < mincell) mincell = ecell;
735  if(ecell > maxcell) maxcell = ecell;
736  }
737  fQAHistos->FillTH2("hClusterConstMinCellEnergy", jet->Pt(), mincell);
738  fQAHistos->FillTH2("hClusterConstMaxCellEnergy", jet->Pt(), maxcell);
739 
740 #ifdef EXPERIMENTAL_JETCONSTITUENTS
741  fQAHistos->FillTH2("hClusterIndexENLC", jet->Pt(), clust->GetNonLinCorrEnergy());
742  fQAHistos->FillTH2("hClusterIndexEHC", jet->Pt(), clust->GetHadCorrEnergy());
743 #endif
744  }
745 
746 #ifdef EXPERIMENTAL_JETCONSTITUENTS
747  // Loop over charged particles - fill test histogram
748  for(int itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
749  auto part = jet->TrackAt(itrk, cont->GetArray());
750  fQAHistos->FillTH2("hChargedIndexPt", jet->Pt(), part->Pt());
751  }
752 
753  // Look over charged constituents
754  AliDebugStream(2) << "Jet: Number of particle constituents: " << jet->GetParticleConstituents().size() << std::endl;
755  for(auto part : jet->GetParticleConstituents()) {
756  //auto part = static_cast<PWG::JETFW::AliEmcalParticleJetConstituent *>(pconst);
757  AliDebugStream(3) << "Found particle constituent with pt " << part.Pt() << ", from VParticle " << part.GetParticle()->Pt() << std::endl;
758  fQAHistos->FillTH2("hChargedConstituentPt", jet->Pt(), part.Pt());
759  }
760 
761  // Loop over neutral constituents
762  AliDebugStream(2) << "Jet: Number of cluster constituents: " << jet->GetClusterConstituents().size() << std::endl;
763  for(auto clust : jet->GetClusterConstituents()){
764  //auto clust = static_cast<PWG::JETFW::AliEmcalClusterJetConstituent *>(cconst);
765  AliDebugStream(3) << "Found cluster constituent with energy " << clust.E() << " using energy definition " << static_cast<int>(clust.GetDefaultEnergyType()) << std::endl;
766  fQAHistos->FillTH2("hClusterConstituentEDefault", jet->Pt(), clust.E());
767  fQAHistos->FillTH2("hClusterConstituentENLC", jet->Pt(), clust.GetCluster()->GetNonLinCorrEnergy());
768  fQAHistos->FillTH2("hClusterConstituentEHC", jet->Pt(), clust.GetCluster()->GetHadCorrEnergy());
769  }
770 
771  // Fill global observables: Leading charged and cluster constituents
772  auto leadingcharged = jet->GetLeadingParticleConstituent();
773  auto leadingcluster = jet->GetLeadingClusterConstituent();
774  if(leadingcluster){
775  fQAHistos->FillTH1("hLeadingClusterConstituentPt", jet->Pt(), leadingcluster->GetCluster()->GetHadCorrEnergy());
776  }
777  if(leadingcharged) {
778  fQAHistos->FillTH1("hLeadingChargedConstituentPt", jet->Pt(), leadingcharged->GetParticle()->Pt());
779  }
780 #endif
781 }
782 
784  int ncharged = 0, nneutral = jet.GetNumberOfClusters();
785  if(particles) {
786  for(decltype(jet.GetNumberOfTracks()) ipart = 0; ipart < jet.GetNumberOfTracks(); ipart++){
787  auto part = jet.TrackAt(ipart, particles->GetArray());
788  if(!part) continue;
789  if(part->Charge()) ncharged++;
790  else nneutral++;
791  }
792  }
793  // check if the jet has at least one consituent for jet substructure
794  int nallowed = 0;
795  nallowed += fUseChargedConstituents ? ncharged : 0;
796  nallowed += fUseNeutralConstituents ? nneutral : 0;
797  return nallowed > 0;
798 }
799 
800 std::string AliAnalysisTaskEmcalJetSubstructureTree::MatchTrigger(const std::string &triggertoken) const {
801  std::vector<std::string> tokens;
802  std::string result;
803  std::stringstream decoder(fInputEvent->GetFiredTriggerClasses().Data());
804  while(std::getline(decoder, result, ' ')) tokens.emplace_back(result);
805  result.clear();
806  for(auto t : tokens) {
807  if(t.find(triggertoken) != std::string::npos) {
808  // take first occurrence - downscale factor should normally be the same
809  result = t;
810  break;
811  }
812  }
813  return result;
814 }
815 
816 bool AliAnalysisTaskEmcalJetSubstructureTree::IsSelectEmcalTriggers(const std::string &triggerstring) const {
817  const std::array<std::string, 8> kEMCALTriggers = {
818  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
819  };
820  bool isEMCAL = false;
821  for(auto emcaltrg : kEMCALTriggers) {
822  if(triggerstring.find(emcaltrg) != std::string::npos) {
823  isEMCAL = true;
824  break;
825  }
826  }
827  return isEMCAL;
828 }
829 
831  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
832 
833  Bool_t isAOD(kFALSE);
834  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
835  if(inputhandler) {
836  if(inputhandler->IsA() == AliAODInputHandler::Class()){
837  std::cout << "Analysing AOD events\n";
838  isAOD = kTRUE;
839  } else {
840  std::cout << "Analysing ESD events\n";
841  }
842  }
843 
844  std::stringstream taskname;
845  taskname << "JetSubstructureTreemaker_R" << std::setw(2) << std::setfill('0') << int(jetradius*10) << trigger;
847  mgr->AddTask(treemaker);
848  treemaker->SetMakeGeneralHistograms(kTRUE);
849  if(isMC) treemaker->SetHasTrueEvent(true);
850  if(isData) treemaker->SetHasRecEvent(true);
851 
852  // Adding containers
853  if(isMC) {
854  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
855  particles->SetMinPt(0.);
856 
857  AliJetContainer *mcjets = treemaker->AddJetContainer(
858  jettype,
860  recombinationScheme,
861  jetradius,
863  particles, nullptr);
864  mcjets->SetName("mcjets");
865  mcjets->SetJetPtCut(20.);
866  }
867 
868  if(isData) {
869  AliTrackContainer *tracks(nullptr);
870  if((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet)){
872  std::cout << "Track container name: " << tracks->GetName() << std::endl;
873  tracks->SetMinPt(0.15);
874  }
875  AliClusterContainer *clusters(nullptr);
876  if((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet)){
877  std::cout << "Using full or neutral jets ..." << std::endl;
879  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
880  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
881  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
882  } else {
883  std::cout << "Using charged jets ... " << std::endl;
884  }
885 
886  AliJetContainer *datajets = treemaker->AddJetContainer(
887  jettype,
889  recombinationScheme,
890  jetradius,
892  tracks, clusters);
893  datajets->SetName("datajets");
894  datajets->SetJetPtCut(20.);
895 
896  treemaker->SetUseAliAnaUtils(true, true);
897  treemaker->SetVzRange(-10., 10);
898 
899  // configure trigger selection
900  std::string triggerstring(trigger);
901  if(triggerstring.find("INT7") != std::string::npos) {
902  treemaker->SetTriggerBits(AliVEvent::kINT7);
903  } else if(triggerstring.find("EJ1") != std::string::npos) {
904  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
905  treemaker->SetTriggerString("EJ1");
906  } else if(triggerstring.find("EJ2") != std::string::npos) {
907  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
908  treemaker->SetTriggerString("EJ2");
909  } else if(triggerstring.find("EG1") != std::string::npos) {
910  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
911  treemaker->SetTriggerString("EG1");
912  } else if(triggerstring.find("EG2") != std::string::npos) {
913  treemaker->SetTriggerBits(AliVEvent::kEMCEGA);
914  treemaker->SetTriggerString("EG2");
915  }
916  }
917 
918  std::string jettypestring;
919  switch(jettype) {
920  case AliJetContainer::kFullJet: jettypestring = "FullJets"; break;
921  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; break;
922  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; break;
923  default: jettypestring = "Undef";
924  };
925 
926  // Connecting containers
927  std::stringstream outputfile, histname, treename;
928  outputfile << mgr->GetCommonFileName() << ":JetSubstructure_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
929  histname << "JetSubstructureHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
930  treename << "JetSubstructureTree_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
931  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
932  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
933  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer(treename.str().data(), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()));
934 
935  return treemaker;
936 }
937 
938 void AliSoftDropParameters::LinkJetTreeBranches(TTree *jettree, const char *tag) {
939  LinkBranch(jettree, &fZg, Form("Zg%s", tag), "D");
940  LinkBranch(jettree, &fRg, Form("Rg%s", tag), "D");
941  LinkBranch(jettree, &fMg, Form("Mg%s", tag), "D");
942  LinkBranch(jettree, &fPtg, Form("Ptg%s", tag), "D");
943  LinkBranch(jettree, &fMug, Form("Mug%s", tag), "D");
944  LinkBranch(jettree, &fDeltaR, Form("DeltaRg%s", tag), "D");
945  LinkBranch(jettree, &fNDropped, Form("NDropped%s", tag), "I");
946 };
947 
949  LinkBranch(jettree, &fOneSubjettiness, Form("OneSubjettiness%s", tag), "D");
950  LinkBranch(jettree, &fTwoSubjettiness, Form("TwoSubjettiness%s", tag), "D");
951 }
952 
954  LinkBranch(jettree, &fAngularity, Form("Angularity%s", tag), "D");
955  LinkBranch(jettree, &fPtD, Form("PtD%s", tag), "D");
956 }
957 
958 void AliJetKineParameters::LinkJetTreeBranches(TTree *jettree, const char *tag){
959  LinkBranch(jettree, &fPt, Form("PtJet%s", tag), "D");
960  LinkBranch(jettree, &fE, Form("EJet%s", tag), "D");
961  LinkBranch(jettree, &fEta, Form("Eta%s", tag), "D");
962  LinkBranch(jettree, &fPhi, Form("Phi%s", tag), "D");
963  LinkBranch(jettree, &fArea, Form("Area%s", tag), "D");
964  LinkBranch(jettree, &fMass, Form("Mass%s", tag), "D");
965  LinkBranch(jettree, &fNEF, Form("NEF%s", tag), "D");
966  LinkBranch(jettree, &fNCharged, Form("NCharged%s", tag), "I");
967  LinkBranch(jettree, &fNNeutral, Form("NNeutral%s", tag), "I");
968  LinkBranch(jettree, &fZLeading, Form("ZLeading%s", tag), "D");
969  LinkBranch(jettree, &fZLeadingCharged, Form("ZLeadingCharged%s", tag), "D");
970  LinkBranch(jettree, &fZLeadingNeutral, Form("ZLeadingNeutral%s", tag), "D");
971 }
972 
974  LinkBranch(jettree, &fJetRadius, "Radius", "D");
975  LinkBranch(jettree, &fEventWeight, "EventWeight", "D");
976  LinkBranch(jettree, &fTriggerClusterIndex, "TriggerClusterIndex", "I");
977  if(fillRho) {
978  std::string varnames[] = {"RhoPtRec", "RhoPtSim", "RhoMassRec", "RhoMassSim"};
979  for(int i = 0; i < 4; i++){
980  LinkBranch(jettree, fRhoParamters + i, varnames[i].data(), "D");
981  }
982  }
983 }
984 
985 void LinkBranch(TTree *jettree, void *data, const char *branchname, const char *type) {
986  jettree->Branch(branchname, data, Form("%s/%s", branchname, type));
987 }
988 } /* 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: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
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.