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