DAS  3.0
Das Analysis System
Classes | Enumerations | Functions | Variables
DAS::TPS Namespace Reference

Classes

struct  has_p4
 
struct  has_p4< T, std::void_t< decltype(std::declval< T >().p4.Eta())> >
 
struct  Hist
 
struct  MultiJets
 
class  TPSBranchFlattener
 
class  TPSTreeFlattener
 
class  TPSVariables
 

Enumerations

enum  RemovalStrategy { DeltaR, AntikT, None }
 
enum  TrainingMode { ThreeClasses = 0, TwoClasses = 1 }
 

Functions

RemovalStrategy parseStrategy (const string &strategyStr, TTree *tree)
 
template<class Jet >
vector< fjcore::PseudoJet > cluster (const vector< Jet > &jets, const double R)
 
template<class Jet >
bool overlappingJets (const vector< Jet > &jets)
 
map< Long64_t, vector< Long64_t > > buildRunToEntriesMap (TTree *tree, RecEvent *recEvent)
 
void applyEventMixing (const vector< fs::path > &inputs1, const vector< fs::path > &inputs2, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
float ShiftJetOrigin (const FourVector &p4_orig, const float &vertexZorig)
 
template<typename ObjT >
void ShiftP4EtaAndFillHist (std::vector< ObjT > &objects, float vertexZ, TH2D &hist)
 
void applyPVshift (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
template<class Jet >
std::vector< Jet > PhaseSpaceSelection (const std::vector< Jet > &jets, float minpt)
 
void getTPSobservables (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
 
void getTPSTrain (const vector< fs::path > &inputs_TPS, const vector< fs::path > &inputs_DPS1, const vector< fs::path > &inputs_DPS2, const vector< fs::path > &inputs_SPS, const fs::path &output, const pt::ptree &config, const int steering)
 

Variables

static const int N_PAIRS = 6
 
static const int N_S_4_PERM = 3
 

Enumeration Type Documentation

◆ RemovalStrategy

enum RemovalStrategy
strong
Enumerator
DeltaR 
AntikT 
None 
38 { DeltaR, AntikT, None };

◆ TrainingMode

enum TrainingMode
strong

Training mode enumeration.

Enumerator
ThreeClasses 
TwoClasses 
47  {
48  ThreeClasses = 0, // TPSvs(DPS1vsDPS2)vsSPS
49  TwoClasses = 1 // TPSvs(DPS1+DPS2+SPS)
50 };

Function Documentation

◆ applyEventMixing()

void DAS::TPS::applyEventMixing ( const vector< fs::path > &  inputs1,
const vector< fs::path > &  inputs2,
const fs::path &  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Mix the jets from second input to first input to simulate DPS or TPS events. Overlapping jets are properly treated. Only events with matching run numbers are mixed. The input must be taken from applyNjetSkim and must have the same isMC flag and pTcut.

Todo:
fix output logic
Todo:
Merge metainfos
Parameters
inputs1input ROOT files (n-tuples)
inputs2input ROOT files (n-tuples)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
105  {1,0}
106  )
107 {
108  cout << __func__ << ' ' << slice << " start" << endl;
109 
110  DT::Flow flow1(steering, inputs1);
111  DT::Flow flow2(steering, inputs2);
112  auto tIn1 = flow1.GetInputTree(slice);
113  auto tIn2 = flow2.GetInputTree(slice);
114  auto tOut2 = flow2.GetOutputTree(output);
115  auto tOut1 = flow1.GetOutputTree(output);
116 
117  DT::MetaInfo metainfo1(tOut1);
118  DT::MetaInfo metainfo2(tOut2);
120  metainfo1.Check(config);
121  metainfo2.Check(config);
122  auto isMC = metainfo1.Get<bool>("flags", "isMC");
123  auto R = metainfo1.Get<int>("flags","R");
124 
125  auto recEvent1 = flow1.GetBranchReadWrite<RecEvent>("recEvent");
126  auto recEvent2 = flow2.GetBranchReadOnly <RecEvent>("recEvent");
127  auto genEvent1 = isMC ? flow1.GetBranchReadWrite<GenEvent>("genEvent") : nullptr;
128  auto genEvent2 = isMC ? flow2.GetBranchReadOnly <GenEvent>("genEvent") : nullptr;
129  auto recJets1 = flow1.GetBranchReadWrite<vector<RecJet>>("recJets");
130  auto recJets2 = flow2.GetBranchReadOnly <vector<RecJet>>("recJets");
131  auto recMuons1 = flow1.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
132  auto recMuons2 = flow2.GetBranchReadOnly <vector<RecMuon>>("recMuons", DT::facultative);
133  auto genJets1 = isMC ? flow1.GetBranchReadWrite<vector<GenJet>>("genJets") : nullptr;
134  auto genJets2 = isMC ? flow2.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
135  auto genMuons1 = isMC ? flow1.GetBranchReadWrite<vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
136  auto genMuons2 = isMC ? flow2.GetBranchReadOnly <vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
137  auto recPVs1 = flow1.GetBranchReadOnly<vector<PrimaryVertex>>("recPVs");
138  auto recPVs2 = flow2.GetBranchReadOnly<vector<PrimaryVertex>>("recPVs");
139  auto genPVs1 = isMC ? flow1.GetBranchReadOnly<vector<PrimaryVertex>>("genPVs") : nullptr;
140  auto genPVs2 = isMC ? flow2.GetBranchReadOnly<vector<PrimaryVertex>>("genPVs") : nullptr;
141 
142  const auto minpt = config.get<float>("skims.jets.minpt");
143  const auto strategyStr = config.get<string>("corrections.event_mixing.overlap");
144  const auto strategy = parseStrategy(strategyStr, tIn1);
145 
146  // Build the map from run numbers to entry indices and track the last used index for each run in the 2nd tree
147  map<Long64_t, vector<Long64_t>> runToEntries = buildRunToEntriesMap(tIn2, recEvent2);
148  map<Long64_t, size_t> lastUsedIndex;
149 
150  Long64_t matchedEvents = 0;
151 
152  for (DT::Looper looper(tIn1); looper(); ++looper) {
153  [[maybe_unused]]
154  static auto& cout = (steering & DT::verbose) ? ::cout : DT::dev_null;
155 
156  int run1 = recEvent1->runNo;
157  // Skip if no matching run number
158  if (runToEntries.find(run1) == runToEntries.end())
159  continue;
160 
161  const auto& matchingEntries = runToEntries[run1];
162  size_t currentIndex = lastUsedIndex[run1];
163  if (currentIndex >= matchingEntries.size()) continue; // no more entries to match
164  Long64_t entry2 = matchingEntries[currentIndex];
165  tIn2->GetEntry(entry2);
166  lastUsedIndex[run1] = currentIndex + 1;
167 
168  // sanity check
169  if (recEvent1->runNo != recEvent2->runNo)
170  cerr << orange << "Run number mismatch: " << recEvent1->runNo << " != " << recEvent2->runNo << def << endl;
171  if (recPVs1->front().z != recPVs2->front().z)
172  cerr << orange << "Unmatched rec-level primary vertex: " << recPVs1->front().z << " != " << recPVs2->front().z << def << endl;
173 
174  if (recJets1 && recJets2) {
175  // Merge input2's jets into input1's jets
176  recJets1->insert(recJets1->end(), recJets2->begin(), recJets2->end());
177  sort(recJets1->begin(), recJets1->end(), pt_sort); // sort again the jets by pt
178 
179  switch (strategy) {
180  case RemovalStrategy::DeltaR:
181  if (overlappingJets(*recJets1)) continue;
182  break;
183  case RemovalStrategy::AntikT:{
184  // reclustering the jets with anti-kt algorithm
185  vector<fjcore::PseudoJet> inclusive_recJets = cluster(*recJets1, R/10.0);
186  // Jet multiplicity may vary due to the reclustering.
187  // We only leave the event with unchanged jet multiplicity.
188  int nJets1 = PhaseSpaceSelection(*recJets1, minpt).size();
189  int nJets2 = PhaseSpaceSelection(inclusive_recJets, minpt).size();
190  if (nJets1 != nJets2) {
191  cout << "Jet multiplicity changed: " << nJets1 << " -> " << nJets2 << endl;
192  cout << "recJets: ";
193  for (const auto& jet : *recJets1) cout << jet.CorrPt() << " ";
194  cout << "\nAnti-kT: ";
195  for (const auto& jet : inclusive_recJets) cout << jet.pt() << " ";
196  cout << endl;
197  continue;
198  }
199  break;
200  }
201  case RemovalStrategy::None:;
202  }
203  }
204  if (recMuons1 && recMuons2) {
205  // Merge input2's muons into input1's muons
206  recMuons1->insert(recMuons1->end(), recMuons2->begin(), recMuons2->end());
207  sort(recMuons1->begin(), recMuons1->end(), pt_sort);
208  }
209 
210  if ((recJets1 && recJets2) || (recMuons1 && recMuons2))
211  recEvent1->weights.front() *= recEvent2->weights.front();
212 
213  if (isMC) {
214  if (genJets1 && genJets2) {
215  // sanity check
216  if (genPVs1->front().z != genPVs2->front().z)
217  BOOST_THROW_EXCEPTION(DE::BadInput(("Unmatched gen-level primary vertex: " + to_string(genPVs1->front().z) + " != " + to_string(genPVs2->front().z)).c_str(), tIn1));
218  // Merge input2's jets into input1's jets
219  genJets1->insert(genJets1->end(), genJets2->begin(), genJets2->end());
220 
221  sort(genJets1->begin(), genJets1->end(), pt_sort); // sort again the jets by pt
222 
223  switch (strategy) {
224  case RemovalStrategy::DeltaR:
225  if (overlappingJets(*genJets1)) continue;
226  break;
227  case RemovalStrategy::AntikT:{
228  // reclustering the jets with anti-kt algorithm
229  vector<fjcore::PseudoJet> inclusive_genJets = cluster(*genJets1, R/10.0);
230  // Jet multiplicity may vary due to the reclustering.
231  // We only leave the event with unchanged jet multiplicity.
232  int nJets1 = PhaseSpaceSelection(*genJets1, minpt).size();
233  int nJets2 = PhaseSpaceSelection(inclusive_genJets, minpt).size();
234  if (nJets1 != nJets2) {
235  cout << "Jet multiplicity changed: " << nJets1 << " -> " << nJets2 << endl;
236  cout << "genJets: ";
237  for (const auto& jet : *genJets1) cout << jet.CorrPt() << " ";
238  cout << "\nAnti-kT: ";
239  for (const auto& jet : inclusive_genJets) cout << jet.pt() << " ";
240  cout << endl;
241  continue;
242  }
243  break;
244  }
245  case RemovalStrategy::None:;
246  }
247  }
248  if (genMuons1 && genMuons2) {
249  // Merge input2's muons into input1's muons
250  genMuons1->insert(genMuons1->end(), genMuons2->begin(), genMuons2->end());
251  sort(genMuons1->begin(), genMuons1->end(), pt_sort);
252  }
253 
254  if ((genJets1 && genJets2) || (genMuons1 && genMuons2))
255  genEvent1->weights.front() *= genEvent2->weights.front();
256  }
257  ++matchedEvents;
258  if (steering & DT::fill) tOut1->Fill();
259  }
260 
261  // Calculate final match rate
262  cout << "Processed " << tIn1->GetEntries() << " events, matched " << matchedEvents << " events ("
263  << fixed << setprecision(2) << matchedEvents * 100.0 / min(tIn1->GetEntries(), tIn2->GetEntries()) << "%)." << endl;
264 
265  metainfo1.Set<bool>("git", "complete", true);
266  cout << __func__ << ' ' << slice << " stop" << endl;
267 }

◆ applyPVshift()

void DAS::TPS::applyPVshift ( const vector< fs::path > &  inputs,
const fs::path &  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Translation of the primary vertex to the center of the detector. The jet four-momentum is adjusted to account for the vertex shift.

Todo:
Shift photons too
Parameters
inputsinput ROOT files (n-tuples)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
slicenumber and index of slice
92  {1,0}
93  )
94 {
95  cout << __func__ << ' ' << slice << " start" << endl;
96 
97  DT::Flow flow(steering, inputs);
98  auto tIn = flow.GetInputTree(slice);
99  auto [fOut, tOut] = flow.GetOutput(output);
100 
101  DT::MetaInfo metainfo(tOut);
102  metainfo.Check(config);
103  auto isMC = metainfo.Get<bool>("flags", "isMC");
104 
105  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
106  auto genJets = isMC ? flow.GetBranchReadWrite<vector<GenJet>>("genJets") : nullptr;
107  auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
108  auto genMuons = isMC ? flow.GetBranchReadWrite<vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
110 
111  auto recPVs = flow.GetBranchReadWrite<vector<PrimaryVertex>>("recPVs");
112  auto genPVs = isMC ? flow.GetBranchReadWrite<vector<PrimaryVertex>>("genPVs") : nullptr;
113 
114  TH2D rec_Deta_PVshift_jet("recDeta_PVshift_jet", "#Delta#eta vs V_{z};V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
115  TH2D gen_Deta_PVshift_jet("genDeta_PVshift_jet", "#Delta#eta vs V_{z};V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
116  TH2D rec_Deta_PVshift_muon("recDeta_PVshift_muon", "#Delta#eta vs V_{z} (muons);V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
117  TH2D gen_Deta_PVshift_muon("genDeta_PVshift_muon", "#Delta#eta vs V_{z} (gen muons);V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
118 
119  for (DT::Looper looper(tIn); looper(); ++looper) {
120  [[maybe_unused]]
121  static auto& cout = (steering & DT::verbose) ? ::cout : DT::dev_null;
122 
123  auto& vertexZ = recPVs->front().z;
124  ShiftP4EtaAndFillHist(*recJets, vertexZ, rec_Deta_PVshift_jet);
125  if (recMuons != nullptr)
126  ShiftP4EtaAndFillHist(*recMuons, vertexZ, rec_Deta_PVshift_muon);
127  vertexZ = 0;
128  if (isMC) {
129  auto& vertexZ = genPVs ->front().z;
130  ShiftP4EtaAndFillHist(*genJets, vertexZ, gen_Deta_PVshift_jet);
131  if (genMuons != nullptr)
132  ShiftP4EtaAndFillHist(*genMuons, vertexZ, gen_Deta_PVshift_muon);
133  vertexZ = 0;
134  }
135 
136  if (steering & DT::fill) tOut->Fill();
137  }
138  fOut->cd();
139  rec_Deta_PVshift_jet.Write();
140  gen_Deta_PVshift_jet.Write();
141  rec_Deta_PVshift_muon.Write();
142  gen_Deta_PVshift_muon.Write();
143  metainfo.Set<bool>("git", "complete", true);
144  cout << __func__ << ' ' << slice << " stop" << endl;
145 }

◆ buildRunToEntriesMap()

map<Long64_t, vector<Long64_t> > DAS::TPS::buildRunToEntriesMap ( TTree *  tree,
RecEvent *  recEvent 
)

Build a map that groups event entries by run number.

84  {
85  map<Long64_t, vector<Long64_t>> runToEntries;
86  cout << "Creating run-to-entries map..." << endl;
87  for (DT::Looper looper(tree); looper(); ++looper) {
88  runToEntries[recEvent->runNo].push_back(*looper);
89  }
90  cout << "Found " << runToEntries.size() << " unique run numbers in tree" << endl;
91  return runToEntries;
92 }

◆ cluster()

vector<fjcore::PseudoJet> DAS::TPS::cluster ( const vector< Jet > &  jets,
const double  R 
)

Anti-kt clustering algorithm to recluster the jets. The jets are sorted by pT and the reclustered jets are returned.

55  {
56  vector<fjcore::PseudoJet> input_jets;
57  input_jets.reserve(jets.size()+1);
58  for (const auto& jet: jets) {
59  const auto& p4 = jet.CorrP4();
60  input_jets.push_back(fjcore::PseudoJet(p4.px(),p4.py(),p4.pz(),p4.E()));
61  }
62  fjcore::JetDefinition jet_def(fjcore::antikt_algorithm, R);
63  fjcore::ClusterSequence clust_seq(input_jets, jet_def);
64  vector<fjcore::PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(0.0));
65  return inclusive_jets;
66 }

◆ getTPSobservables()

void DAS::TPS::getTPSobservables ( const vector< fs::path >  inputs,
const fs::path  output,
const pt::ptree &  config,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Calculate the observables that can be used to caracterise 6 jets events.

Parameters
inputsinput ROOT files (n-tuples)
outputname of output root file containing the histograms
configconfig file from `DTOptions`
steeringsteering parameters from `DTOptions`
sliceslices for running
121  {1,0}
122  )
123 {
124  cout << __func__ << ' ' << slice << " start" << endl;
125 
126  DT::Flow flow(steering, inputs);
127  auto tIn = flow.GetInputTree(slice);
128  auto [fOut, tOut] = flow.GetOutput(output);
129 
130  DT::MetaInfo metainfo(tOut);
131  metainfo.Check(config);
132  auto isMC = metainfo.Get<bool>("flags", "isMC");
133 
134  const auto minpt = config.get<float>("skims.jets.minpt");
135  metainfo.Set<float>("skims", "jets", "minpt", minpt);
136 
137  auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
138  auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
139  auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
140  auto genJets = isMC ? flow.GetBranchReadOnly<vector<GenJet>>("genJets") : nullptr;
141 
142  Hist genHist("gen"),
143  recHist("rec");
144 
145  for (DT::Looper looper(tIn); looper(); ++looper) {
146  [[ maybe_unused ]]
147  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
148 
149  MultiJets recMultijets(*recJets, minpt);
150  if (recMultijets) {
151  double evweight = rEv->weights.front();
152  if (isMC) evweight *= gEv->weights.front();
153  recHist.Fill(recMultijets, evweight);
154  }
155 
156  if (!isMC) continue;
157  MultiJets genMultijets(*genJets, minpt);
158  if (genMultijets) {
159  double evweight = gEv->weights.front();
160  genHist.Fill(genMultijets, evweight);
161  }
162  }
163 
164  recHist.Write(fOut);
165  if (isMC)
166  genHist.Write(fOut);
167  metainfo.Set<bool>("git", "complete", true);
168 
169  cout << __func__ << ' ' << slice << " end" << endl;
170 }

◆ getTPSTrain()

void DAS::TPS::getTPSTrain ( const vector< fs::path > &  inputs_TPS,
const vector< fs::path > &  inputs_DPS1,
const vector< fs::path > &  inputs_DPS2,
const vector< fs::path > &  inputs_SPS,
const fs::path &  output,
const pt::ptree &  config,
const int  steering 
)

Unified TPS training function (Multi-class and Classification)

Parameters
inputs_TPSTPS input ROOT files (n-tuples)
inputs_DPS1DPS1 input ROOT files (n-tuples)
inputs_DPS2DPS2 input ROOT files (n-tuples)
inputs_SPSSPS input ROOT files (n-tuples)
outputoutput ROOT file (n-tuple)
configconfig handled with `Darwin::Tools::options`
steeringparameters obtained from explicit options
63 {
64  DT::Flow flow_TPS(steering, inputs_TPS);
65  DT::Flow flow_DPS1(steering, inputs_DPS1);
66  DT::Flow flow_DPS2(steering, inputs_DPS2);
67  DT::Flow flow_SPS(steering, inputs_SPS);
68  auto tIn_TPS = flow_TPS.GetInputTree();
69  auto tIn_DPS1 = flow_DPS1.GetInputTree();
70  auto tIn_DPS2 = flow_DPS2.GetInputTree();
71  auto tIn_SPS = flow_SPS.GetInputTree();
72  auto [fOut, tOut] = flow_TPS.GetOutput(output);
73 
74  DT::MetaInfo metainfo(tOut, config);
75 
76  TrainingMode mode = static_cast<TrainingMode>(config.get<int>("training.mode"));
77 
78  int nEntries_TPS = tIn_TPS->GetEntries();
79  int nEntries_DPS1 = tIn_DPS1->GetEntries();
80  int nEntries_DPS2 = tIn_DPS2->GetEntries();
81  int nEntries_SPS = tIn_SPS->GetEntries();
82 
83  // Training sample size calculation depends on mode
84  int nTrain;
85  switch (mode) {
86  case TrainingMode::ThreeClasses:
87  {
88  int nEntries_DPS = nEntries_DPS1 + nEntries_DPS2;
89  nTrain = ranges::min({nEntries_TPS, nEntries_DPS, nEntries_SPS}) * 0.8;
90  cout << "Three-classes training mode selected" << endl;
91  }
92  break;
93  case TrainingMode::TwoClasses:
94  nTrain = ranges::min({nEntries_TPS, nEntries_DPS1, nEntries_DPS2, nEntries_SPS}) * 0.8;
95  cout << "Two-classes training mode selected" << endl;
96  break;
97  }
98 
99  // ------------ TMVA Training Setup ------------
100  TMVA::Tools::Instance();
101  unique_ptr<TMVA::Factory> factory;
102  switch (mode) {
103  case TrainingMode::ThreeClasses:
104  factory = make_unique<TMVA::Factory>(
105  "TMVAMulticlass", fOut,
106  "!V:!Silent:Color:DrawProgressBar:Transformations=I,D,P,G,N:AnalysisType=Multiclass");
107  break;
108  case TrainingMode::TwoClasses:
109  factory = make_unique<TMVA::Factory>(
110  "TMVAClassification", fOut,
111  "!V:!Silent:Color:DrawProgressBar:Transformations=I,D,P,G,N:AnalysisType=Classification");
112  break;
113  }
114 
115  auto dataloader = make_unique<TMVA::DataLoader>("dataset");
116 
117  DAS::TPS::TPSVariables tpsVars(dataloader.get());
118  tpsVars.printVariables();
119 
120  switch (mode) {
121  case TrainingMode::ThreeClasses:
122  dataloader->AddTree(tIn_TPS, "TPS", 1.0);
123  dataloader->AddTree(tIn_DPS1, "DPS", 1.0);
124  dataloader->AddTree(tIn_DPS2, "DPS", 1.0);
125  dataloader->AddTree(tIn_SPS, "SPS", 1.0);
126  dataloader->SetWeightExpression("EventWeight", "TPS");
127  dataloader->SetWeightExpression("EventWeight", "DPS");
128  dataloader->SetWeightExpression("EventWeight", "SPS");
129  break;
130  case TrainingMode::TwoClasses:
131  dataloader->AddSignalTree(tIn_TPS, 1.0);
132  dataloader->AddBackgroundTree(tIn_DPS1, 1.0);
133  dataloader->AddBackgroundTree(tIn_DPS2, 1.0);
134  dataloader->AddBackgroundTree(tIn_SPS, 1.0);
135  dataloader->SetWeightExpression("EventWeight", "Signal");
136  dataloader->SetWeightExpression("EventWeight", "Background");
137  break;
138  }
139 
140  cout << green << "Trees added successfully. About to prepare training and test trees..." << def << endl;
141 
142  TCut mycuts = "";
143  switch (mode) {
144  case TrainingMode::ThreeClasses:
145  {
146  dataloader->PrepareTrainingAndTestTree(
147  mycuts,
148  Form("SplitMode=Random:NormMode=NumEvents:"
149  "SplitSeed=100:nTrain_TPS=%d:nTrain_DPS=%d:"
150  "nTrain_SPS=%d:!V",
151  nTrain, nTrain, nTrain));
152 
153  TString BDToptions =
154  "!H:V:NTrees=800:MaxDepth=2:MinNodeSize=5%:"
155  "nCuts=100:BoostType=Grad:Shrinkage=0.03:"
156  "UseBaggedBoost:BaggedSampleFraction=0.7:"
157  "UseRandomisedTrees";
158  factory->BookMethod(dataloader.get(), TMVA::Types::kBDT, "BDT", BDToptions);
159 
160  TString DNNoptions =
161  "!H:V:ErrorStrategy=CROSSENTROPY"
162  ":WeightInitialization=XAVIER"
163  ":Layout=RELU|256,RELU|128,RELU|64,LINEAR"
164  ":TrainingStrategy=LearningRate=1e-4,Momentum=0.9,"
165  "Repetitions=1,ConvergenceSteps=20,BatchSize=1024,"
166  "TestRepetitions=1,MaxEpochs=100,WeightDecay=1e-3,"
167  "Regularization=L2,Optimizer=ADAM,"
168  "DropConfig=0.0+0.1+0.1+0.0";
169  factory->BookMethod(dataloader.get(), TMVA::Types::kDL, "DNN", DNNoptions);
170  }
171  break;
172  case TrainingMode::TwoClasses:
173  {
174  TCut mycutb = "";
175  dataloader->PrepareTrainingAndTestTree(
176  mycuts, mycutb,
177  Form("SplitMode=Random:NormMode=EqualNumEvents:"
178  "VerboseLevel=Debug:nTrain_Signal=%d:"
179  "nTrain_Background=%d:nTest_Signal=%d:"
180  "nTest_Background=%d",
181  nTrain, nTrain * 3, nTrain / 4, nTrain * 3 / 4));
182 
183  TString BDToptions =
184  "H:V:NTrees=800:MaxDepth=3:MinNodeSize=2%:"
185  "nCuts=100:BoostType=Grad:Shrinkage=0.05:"
186  "UseBaggedBoost:BaggedSampleFraction=0.8:"
187  "UseRandomisedTrees";
188  factory->BookMethod(dataloader.get(), TMVA::Types::kBDT, "BDT", BDToptions);
189 
190  TString DNNoptions =
191  "H:V:ErrorStrategy=CROSSENTROPY:"
192  "WeightInitialization=XAVIER"
193  ":Layout=RELU|512,RELU|256,RELU|128,RELU|64,LINEAR"
194  ":TrainingStrategy=LearningRate=5e-4,Momentum=0.9,"
195  "Repetitions=1,ConvergenceSteps=20,BatchSize=512,"
196  "TestRepetitions=1,MaxEpochs=100,WeightDecay=1e-5,"
197  "Regularization=L2,Optimizer=ADAM,"
198  "DropConfig=0.0+0.3+0.3+0.0";
199  factory->BookMethod(dataloader.get(), TMVA::Types::kDL, "DNN", DNNoptions);
200 
201  TString RFoptions =
202  "!H:!V:NTrees=200:MinNodeSize=1%:MaxDepth=4:"
203  "UseRandomisedTrees:UsePoissonNvars";
204  factory->BookMethod(dataloader.get(), TMVA::Types::kBDT, "RF", RFoptions);
205 
206  TString XGBoostOptions = "!H:!V:NTrees=1000:MaxDepth=4";
207  factory->BookMethod(dataloader.get(), TMVA::Types::kBDT, "XGBoost", XGBoostOptions);
208 
209  cout << "BDT options: " << BDToptions << endl;
210  cout << "DNN options: " << DNNoptions << endl;
211  cout << "RF options: " << RFoptions << endl;
212  }
213  break;
214  }
215 
216  factory->TrainAllMethods();
217  factory->TestAllMethods();
218  factory->EvaluateAllMethods();
219 
220  metainfo.Set<bool>("git", "complete", true);
221 
222 }

◆ overlappingJets()

bool DAS::TPS::overlappingJets ( const vector< Jet > &  jets)

Loop all possible jets combinations checking azimuthal difference to avioding overlapping jets.

72  {
73  for (size_t i = 0; i < jets.size(); ++i) {
74  for (size_t j = i + 1; j < jets.size(); ++j) {
75  if (DeltaR(jets[i].p4, jets[j].p4) < 0.4)
76  return true;
77  }
78  }
79  return false;
80 }

◆ parseStrategy()

RemovalStrategy DAS::TPS::parseStrategy ( const string &  strategyStr,
TTree *  tree 
)

Define the strategy to remove overlapping jets.

42  {
43  cout << "Removal strategy = " << strategyStr << endl;
44  if (strategyStr == "DeltaR") return RemovalStrategy::DeltaR;
45  else if (strategyStr == "antikt") return RemovalStrategy::AntikT;
46  else if (strategyStr == "none") return RemovalStrategy::None;
47 
48  BOOST_THROW_EXCEPTION(DE::BadInput(("Unrecognized strategy: " + strategyStr).c_str(), tree));
49 }

◆ PhaseSpaceSelection()

std::vector<Jet> DAS::TPS::PhaseSpaceSelection ( const std::vector< Jet > &  jets,
float  minpt 
)
24  {
25  std::vector<Jet> tempjets;
26  for (const auto& jet : jets) {
27  double eta;
28  if constexpr(has_p4<Jet>::value) {
29  eta = jet.p4.Eta();
30  } else {
31  eta = jet.eta();
32  }
33  if (jet.CorrPt() > minpt && eta < 5.0) {
34  tempjets.push_back(jet);
35  } else {
36  break;
37  }
38  }
39  return tempjets;
40 }

◆ ShiftJetOrigin()

float DAS::TPS::ShiftJetOrigin ( const FourVector &  p4_orig,
const float &  vertexZorig 
)

Adjust jet four-momentum to account for a vertex shift

  1. Jet direction: $\vec{d} = (\rho=1, \eta_0, \phi_0)$ (invariant)
  2. Hit point from original vertex $(0,0,z_v)$:
    • Barrel: $\vec{h}_1 = (R_{cal}\cos\phi, R_{cal}\sin\phi, s \cdot d_z)$ where $s = R_{cal}/\rho_d$
    • Endcap: $\vec{h}_1 = (s\rho_d\cos\phi, s\rho_d\sin\phi, Z_{cal})$ where $s = (Z_{cal})/d_z$
  3. Shifted hit point: $\vec{h}_2 = \vec{h}_1 - (0,0,z_v)$ (detector shift)
  4. Return: $(p_T, \eta_{new}, \phi_0, m)$ where $\eta_{new} = \vec{h}_2.\text{Eta()}$
41  {
42  // CMS calorimeter geometry (units: cm)
43  // Only use as reference, not for actual calculations
44  const float R_CAL_BARREL = 129.0f; // ECAL Barrel inner radius
45  const float Z_CAL_ENDCAP = 314.0f; // ECAL Endcap z-position
46  const float ETA_TRANSITION = 1.479f; // Barrel/endcap transition
47 
48  ROOT::Math::RhoEtaPhiVector dir(1.0, p4_orig.Eta(), p4_orig.Phi());
49  ROOT::Math::XYZPoint hitPoint;
50 
51  if (abs(p4_orig.Eta()) < ETA_TRANSITION) {
52  // Barrel: hit point on cylindrical surface at R_CAL_BARREL
53  float s = R_CAL_BARREL / dir.Rho();
54  float hit_z = s * dir.Z();
55  hitPoint.SetXYZ(R_CAL_BARREL * cos(p4_orig.Phi()), R_CAL_BARREL * sin(p4_orig.Phi()), hit_z - vertexZorig);
56  } else {
57  // Endcap: hit point on planar surface at target_z
58  float target_z = copysign(Z_CAL_ENDCAP, p4_orig.Eta());
59  float s = target_z / dir.Z();
60  float hit_rho = s * dir.Rho();
61  hitPoint.SetXYZ(hit_rho * cos(p4_orig.Phi()), hit_rho * sin(p4_orig.Phi()), target_z - vertexZorig);
62  }
63  return hitPoint.Eta();
64 }

◆ ShiftP4EtaAndFillHist()

void DAS::TPS::ShiftP4EtaAndFillHist ( std::vector< ObjT > &  objects,
float  vertexZ,
TH2D &  hist 
)

Adjusts objects' eta and fills histogram with vertex shift effect.

This function:

  1. Loops over objects and gets original eta from four-momentum p4.
  2. Computes shifted eta using ShiftJetOrigin with vertexZ.
  3. Updates object's p4 with new eta.
  4. Fills histogram with (-vertexZ) and eta difference (new_eta - orig_eta).
75  {
76  for (auto& obj: objects) {
77  float orig_eta = obj.p4.Eta();
78  float new_eta = ShiftJetOrigin(obj.p4, vertexZ);
79  obj.p4.SetEta(new_eta);
80  hist.Fill(-vertexZ, new_eta - orig_eta);
81  }
82 }

Variable Documentation

◆ N_PAIRS

const int N_PAIRS = 6
static

Number of possible jet pairs when considering 6 leading jets (N_PAIRS = 6) Number of possible permutations for 4-jet balance observables (N_S_4_PERM = 3) Each permutation corresponds to different ways of pairing 4 jets into 2 dijet systems

◆ N_S_4_PERM

const int N_S_4_PERM = 3
static
Ntupliser_cfg.cerr
cerr
Definition: Ntupliser_cfg.py:105
Darwin::Tools::fill
@ fill
activate -f to fill the tree
Definition: Options.h:32
Darwin::Tools::Flow
User-friendly handling of input and output n-tuples.
Definition: Flow.h:80
Step::verbose
static bool verbose
Definition: Step.h:40
Step::def
static const char * def
Definition: Step.h:36
DAS::TPS::TPSVariables
Definition: TPSVariables.h:16
DAS::TPS::TrainingMode
TrainingMode
Training mode enumeration.
Definition: getTPStraining.cc:47
DAS::TPS::parseStrategy
RemovalStrategy parseStrategy(const string &strategyStr, TTree *tree)
Define the strategy to remove overlapping jets.
Definition: applyEventMixing.cc:42
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:22
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
DAS::TPS::cluster
vector< fjcore::PseudoJet > cluster(const vector< Jet > &jets, const double R)
Definition: applyEventMixing.cc:55
DAS::TPS::TrainingMode::ThreeClasses
@ ThreeClasses
DAS::TPS::TrainingMode::TwoClasses
@ TwoClasses
DAS::pt_sort
bool pt_sort(const PhysicsObject &j1, const PhysicsObject &j2)
Definition: toolbox.h:18
DAS::MN::minpt
static const double minpt
Definition: getMNobservables.cc:39
sorted_by_pt
std::vector< PseudoJet > sorted_by_pt(const std::vector< PseudoJet > &jets)
Definition: fjcore.cc:4063
Ntupliser_cfg.jets
string jets
Definition: Ntupliser_cfg.py:41
antikt_algorithm
@ antikt_algorithm
Definition: fjcore.hh:1040
DAS::TPS::ShiftP4EtaAndFillHist
void ShiftP4EtaAndFillHist(std::vector< ObjT > &objects, float vertexZ, TH2D &hist)
Adjusts objects' eta and fills histogram with vertex shift effect.
Definition: applyPVshift.cc:75
Ntupliser_cfg.config
config
Definition: Ntupliser_cfg.py:333
DAS::TPS::overlappingJets
bool overlappingJets(const vector< Jet > &jets)
Definition: applyEventMixing.cc:72
orange
static const char * orange
Definition: colours.h:6
DAS::TPS::RemovalStrategy::AntikT
@ AntikT
Step::green
static const char * green
Definition: Step.h:33
DAS::TPS::PhaseSpaceSelection
std::vector< Jet > PhaseSpaceSelection(const std::vector< Jet > &jets, float minpt)
Definition: common.h:24
Step::None
@ None
iterate until max value is reached
Definition: Step.h:431
Ntupliser_cfg.isMC
dictionary isMC
Definition: Ntupliser_cfg.py:62
DAS::TPS::RemovalStrategy::DeltaR
@ DeltaR
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
DAS::TPS::ShiftJetOrigin
float ShiftJetOrigin(const FourVector &p4_orig, const float &vertexZorig)
Definition: applyPVshift.cc:41
Darwin::Exceptions::BadInput
Generic exception for ill-defined input (before the event loop).
Definition: exceptions.h:83
jmarExample.eta
eta
DeepAK8/ParticleNet tagging.
Definition: jmarExample.py:19
Darwin::Tools::facultative
@ facultative
mounting branch is facultative
Definition: Flow.h:33
DAS::TPS::buildRunToEntriesMap
map< Long64_t, vector< Long64_t > > buildRunToEntriesMap(TTree *tree, RecEvent *recEvent)
Build a map that groups event entries by run number.
Definition: applyEventMixing.cc:84