ATLAS Offline Software
SUSYCrossSection.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Class header
6 #include "SUSYTools/SUSYCrossSection.h"
7 
8 // Find the xsec files in datapath
9 #include "PathResolver/PathResolver.h"
10 
11 // For directory commands
12 #include <dirent.h>
13 
14 // Error messages
15 #include <iostream>
16 
17 // Text file i/o
18 #include <fstream>
19 #include <array>
20 
21 namespace{
22  enum Sparticles{
23  UNKNOWN = -1,
24  ngluino,
25  nsquark, // (up and down type without bottom/top)
26  nantisquark, // (up and down type without bottom/top)
27 
28  nsbottom,
29  nstop,
30  nsbottom2,
31  nstop2,
32  nantisbottom,
33  nantistop,
34  nantisbottom2,
35  nantistop2,
36 
37  nchi01,
38  nchi02,
39  nchi03,
40  nchi04,
41  nch1plus,
42  nch2plus,
43  nch1minus,
44  nch2minus,
45 
46  //sleptons
47  nsmuonRplus,
48  nsmuonRminus,
49  nselecRplus,
50  nselecRminus,
51 
52  nsmuonLplus,
53  nsmuonLminus,
54  nselecLplus,
55  nselecLminus,
56 
57  nstau1plus,
58  nstau1minus,
59  nstau2plus,
60  nstau2minus,
61 
62  //snutrinos
63  nselnuL,
64  nsmunuL,
65  nstaunuL,
66  NUM_SPARTICLES
67  };
68 
69  Sparticles
70  classifyOne(int pdgId){
71  if (abs(pdgId) == 1000022) return nchi01;
72  else if (abs(pdgId) == 1000023) return nchi02;
73  else if (abs(pdgId) == 1000025) return nchi03;
74  else if (abs(pdgId) == 1000035) return nchi04;
75  else if ( pdgId == 1000024) return nch1plus;
76  else if ( pdgId == -1000024) return nch1minus;
77  else if ( pdgId == 1000037) return nch2plus;
78  else if ( pdgId == -1000037) return nch2minus;
79  else if ( pdgId == 1000021) return ngluino;
80  else if ((abs(pdgId) > 1000000 && abs(pdgId) <= 1000004) || (abs(pdgId) > 2000000 && abs(pdgId) <= 2000004)) {
81  if (pdgId > 0) return nsquark;
82  else return nantisquark;
83  }
84  else if (pdgId == 1000005) return nsbottom;
85  else if (pdgId == 1000006) return nstop;
86  else if (pdgId == 2000005) return nsbottom2;
87  else if (pdgId == 2000006) return nstop2;
88  else if (pdgId == -1000005) return nantisbottom;
89  else if (pdgId == -1000006) return nantistop;
90  else if (pdgId == -2000005) return nantisbottom2;
91  else if (pdgId == -2000006) return nantistop2;
92  else if (pdgId == 2000011) return nselecRminus;
93  else if (pdgId == -2000011) return nselecRplus;
94  else if (pdgId == 1000011) return nselecLminus;
95  else if (pdgId == -1000011) return nselecLplus;
96  else if (abs(pdgId) == 1000012) return nselnuL;
97  else if (pdgId == 2000013) return nsmuonRminus;
98  else if (pdgId == -2000013) return nsmuonRplus;
99  else if (pdgId == 1000013) return nsmuonLminus;
100  else if (pdgId == -1000013) return nsmuonLplus;
101  else if (abs(pdgId) == 1000014) return nsmunuL;
102  else if (pdgId == 1000015) return nstau1minus;
103  else if (pdgId == -1000015) return nstau1plus;
104  else if (pdgId == 2000015) return nstau2minus;
105  else if (pdgId == -2000015) return nstau2plus;
106  else if (abs(pdgId) == 1000016) return nstaunuL;
107  return UNKNOWN;
108  }
109 
110 
111 }
112 
113 SUSY::CrossSectionDB::CrossSectionDB(const std::string& txtfilename, bool usePathResolver, bool isExtended, bool usePMGTool)
114  : m_pmgxs("")
115 {
116  setExtended(isExtended);
117  setUsePMGTool(usePMGTool);
118 
119  // configuring the PMG tool...
120  m_pmgxs.setTypeAndName("PMGTools::PMGCrossSectionTool/PMGCrossSectionTool");
121  m_pmgxs.retrieve().ignore(); // Ignore the status code
122  std::vector<std::string> inFiles={};
123  std::string fullPath = usePathResolver ? PathResolverFindCalibFile(txtfilename) : txtfilename;
124  std::ifstream in(fullPath.c_str());
125  if (!in) return;
126  inFiles.push_back(fullPath);
127  m_pmgxs->readInfosFromFiles( inFiles );
128 
129 }
130 
131 void SUSY::CrossSectionDB::loadFile(const std::string& txtfilename){
132 
133  std::string line;
134 
135  std::ifstream in(txtfilename.c_str());
136  if (!in) return;
137  while ( getline(in, line) ){
138  // skip leading blanks (in case there are some in front of a comment)
139  if ( !line.empty() ){
140  while ( line[0] == ' ' ) line.erase(0, 1);
141  }
142  // skip lines that do not start with a number, they are comments
143  if ( !line.empty() && isdigit(line[0]) ){
144  std::stringstream is(line);
145  int id;
146  std::string name;
147  float xsect, kfactor, efficiency, relunc;
148  float sumweight = -1, stat = -1;
149  is >> id >> name >> xsect >> kfactor >> efficiency >> relunc;
150  if (m_extended == true){
151  is >> sumweight >> stat;
152  }
153  m_xsectDB[Key(id, name)] = Process(id, name, xsect, kfactor, efficiency, relunc, sumweight, stat);
154  }
155  }
156 }
157 
158 // Convenient accessor for finding based on *only* a process ID
159 SUSY::CrossSectionDB::xsDB_t::iterator SUSY::CrossSectionDB::my_find( const int proc ) {
160  for (SUSY::CrossSectionDB::xsDB_t::iterator it = m_xsectDB.begin(); it!=m_xsectDB.end();++it){
161  if (it->second.ID()==proc) return it;
162  }
163  return m_xsectDB.end();
164 }
165 
166 // Extend the record based on information from a second file
167 void SUSY::CrossSectionDB::extend(const std::string& txtfilename){
168  // Just like the above function, but with more functionality
169  std::string line;
170 
171  std::ifstream in(txtfilename.c_str());
172  if (!in) return;
173  while ( getline(in, line) )
174  {
175  // skip leading blanks (in case there are some in front of a comment)
176  if ( !line.empty() ){
177  while ( line[0] == ' ' ) line.erase(0, 1);
178  }
179  // skip lines that do not start with a number, they are comments
180  if ( !line.empty() && isdigit(line[0]) ){
181  std::stringstream is(line);
182  int id;
183  std::string name;
184  float xsect, kfactor, efficiency, relunc;
185  float sumweight = -1, stat = -1;
186  is >> id;
187  auto my_it = my_find( id );
188  if (my_it==m_xsectDB.end()){
189  is >> name >> xsect >> kfactor >> efficiency >> relunc;
190  if (m_extended == true){
191  is >> sumweight >> stat;
192  }
193  m_xsectDB[Key(id, name)] = Process(id, name, xsect, kfactor, efficiency, relunc, sumweight, stat);
194  } else {
195  // Now we have extended records
196  if (!m_extended) m_extended=true;
197  is >> sumweight >> stat;
198  my_it->second.sumweight(sumweight);
199  my_it->second.stat(stat);
200  }
201  }
202  }
203 
204 }
205 
206 SUSY::CrossSectionDB::Process SUSY::CrossSectionDB::process(int id, int proc) const
207 {
208  // for background x-sections, use the PMG tool
209  if(proc==0 && m_usePMGTool) {
210  return Process( id, m_pmgxs->getSampleName(id), m_pmgxs->getAMIXsection(id), m_pmgxs->getKfactor(id), m_pmgxs->getFilterEff(id), m_pmgxs->getXsectionUncertainty(id), -1, -1 );
211  } else {
212  const Key k(id, proc);
213  xsDB_t::const_iterator pos = m_cache.find(k);
214  if (pos != m_cache.end()) {
215  return pos->second;
216  } else {
217  pos = m_xsectDB.find(k);
218  if (pos != m_xsectDB.end()) {
219  return pos->second;
220  }
221  }
222  }
223  return Process();
224 }
225 
226 
227 unsigned int SUSY::finalState(const int SUSY_Spart1_pdgId, const int SUSY_Spart2_pdgId)
228 {
229  std::array<int, NUM_SPARTICLES> n{};
230  //Classification of the event follows (gg, sq...):
231  int idx = classifyOne(SUSY_Spart1_pdgId);
232  if ( (idx>=0) and (idx<std::ssize(n) )){
233  ++n[idx];
234  }
235  int idx2 = classifyOne(SUSY_Spart2_pdgId);
236  if ( (idx2>=0) and (idx2<std::ssize(n) )){
237  ++n[idx2];
238  }
239 
241  // gluino/squark + X
242  auto nanysquark=[&n](int i)->bool {
243  return (n[nsquark] == i) or (n[nantisquark] == i);
244  };
245  if (n[ngluino] == 1 && nanysquark(1)) return 1;
246  else if (n[ngluino] == 2) return 2;
247  else if (nanysquark(2)) return 3;
248  else if (n[nsquark] == 1 && n[nantisquark] == 1) return 4;
249 
250  else if (n[nsbottom] == 1 && n[nantisbottom] == 1) return 51;
251  else if (n[nsbottom2] == 1 && n[nantisbottom2] == 1) return 52;
252  else if (n[nstop] == 1 && n[nantistop] == 1) return 61;
253  else if (n[nstop2] == 1 && n[nantistop2] == 1) return 62;
254 
255  else if (n[ngluino] == 1 && n[nchi01] == 1) return 71;
256  else if (n[ngluino] == 1 && n[nchi02] == 1) return 72;
257  else if (n[ngluino] == 1 && n[nchi03] == 1) return 73;
258  else if (n[ngluino] == 1 && n[nchi04] == 1) return 74;
259 
260  else if (n[ngluino] == 1 && n[nch1plus] == 1) return 75;
261  else if (n[ngluino] == 1 && n[nch2plus] == 1) return 76;
262  else if (n[ngluino] == 1 && n[nch1minus] == 1) return 77;
263  else if (n[ngluino] == 1 && n[nch2minus] == 1) return 78;
264 
265  else if (nanysquark(1) && n[nchi01] == 1) return 81;
266  else if (nanysquark(1) && n[nchi02] == 1) return 82;
267  else if (nanysquark(1) && n[nchi03] == 1) return 83;
268  else if (nanysquark(1) && n[nchi04] == 1) return 84;
269 
270  else if (nanysquark(1) && n[nch1plus] == 1) return 85;
271  else if (nanysquark(1) && n[nch2plus] == 1) return 86;
272  else if (nanysquark(1) && n[nch1minus] == 1) return 87;
273  else if (nanysquark(1) && n[nch2minus] == 1) return 88;
274 
275 
276  // Gaugino pair-production
277  // chi^{0}_1 + X
278  else if (n[nchi01] == 2) return 111;
279  else if (n[nchi01] == 1 && n[nchi02] == 1) return 112;
280  else if (n[nchi01] == 1 && n[nchi03] == 1) return 113;
281  else if (n[nchi01] == 1 && n[nchi04] == 1) return 114;
282  else if (n[nchi01] == 1 && n[nch1plus] == 1) return 115;
283  else if (n[nchi01] == 1 && n[nch2plus] == 1) return 116;
284  else if (n[nchi01] == 1 && n[nch1minus] == 1) return 117;
285  else if (n[nchi01] == 1 && n[nch2minus] == 1) return 118;
286 
287  // chi^{0}_2 + X
288  else if (n[nchi02] == 2) return 122;
289  else if (n[nchi02] == 1 && n[nchi03] == 1) return 123;
290  else if (n[nchi02] == 1 && n[nchi04] == 1) return 124;
291  else if (n[nchi02] == 1 && n[nch1plus] == 1) return 125;
292  else if (n[nchi02] == 1 && n[nch2plus] == 1) return 126;
293  else if (n[nchi02] == 1 && n[nch1minus] == 1) return 127;
294  else if (n[nchi02] == 1 && n[nch2minus] == 1) return 128;
295 
296  // chi^{0}_3 + X
297  else if (n[nchi03] == 2) return 133;
298  else if (n[nchi03] == 1 && n[nchi04] == 1) return 134;
299  else if (n[nchi03] == 1 && n[nch1plus] == 1) return 135;
300  else if (n[nchi03] == 1 && n[nch2plus] == 1) return 136;
301  else if (n[nchi03] == 1 && n[nch1minus] == 1) return 137;
302  else if (n[nchi03] == 1 && n[nch2minus] == 1) return 138;
303 
304  // chi^{0}_4 + X
305  else if (n[nchi04] == 2) return 144;
306  else if (n[nchi04] == 1 && n[nch1plus] == 1) return 145;
307  else if (n[nchi04] == 1 && n[nch2plus] == 1) return 146;
308  else if (n[nchi04] == 1 && n[nch1minus] == 1) return 147;
309  else if (n[nchi04] == 1 && n[nch2minus] == 1) return 148;
310 
311  // chi^{+}_1/2 + chi^{-}_1/2
312  else if (n[nch1plus] == 1 && n[nch1minus] == 1) return 157;
313  else if (n[nch1plus] == 1 && n[nch2minus] == 1) return 158;
314 
315  else if (n[nch2plus] == 1 && n[nch1minus] == 1) return 167;
316  else if (n[nch2plus] == 1 && n[nch2minus] == 1) return 168;
317 
318  // slepton
319  else if (n[nselecLplus] == 1 && n[nselecLminus] == 1) return 201; // sElectronLPair
320  else if (n[nselecRplus] == 1 && n[nselecRminus] == 1) return 202; // sElectronRPair
321  else if (n[nselnuL] == 2) return 203; // sElectron neutrino pair
322  else if (n[nselecLplus] == 1 && n[nselnuL] == 1) return 204; // sElectron+ sNutrino
323  else if (n[nselecLminus] == 1 && n[nselnuL] == 1) return 205; // sElectron- sNutrino
324  else if (n[nstau1plus] == 1 && n[nstau1minus] == 1) return 206;
325  else if (n[nstau2plus] == 1 && n[nstau2minus] == 1) return 207;
326  else if ((n[nstau1plus] == 1 || n[nstau1minus] == 1) && (n[nstau2plus] == 1 || n[nstau2minus] == 1)) return 208;
327  else if (n[nstaunuL] == 2) return 209; // sTau neutrino pair
328  else if (n[nstau1plus] == 1 && n[nstaunuL] == 1) return 210;
329  else if (n[nstau1minus] == 1 && n[nstaunuL] == 1) return 211;
330  else if (n[nstau2plus] == 1 && n[nstaunuL] == 1) return 212;
331  else if (n[nstau2minus] == 1 && n[nstaunuL] == 1) return 213;
332 
333  else if (n[nsmuonLplus] == 1 && n[nsmuonLminus] == 1) return 216; // sMuonPair
334  else if (n[nsmuonRplus] == 1 && n[nsmuonRminus] == 1) return 217; // sMuonPair
335  else if (n[nsmunuL] == 2) return 218; // sMuon neutrino pair
336  else if (n[nsmuonLplus] == 1 && n[nsmunuL] == 1) return 219; // sMuon+ sNutrino
337  else if (n[nsmuonLminus] == 1 && n[nsmunuL] == 1) return 220; // sMuon- sNutrino
338 
339  std::cerr << "ERROR. could not determine finalState for:" << std::endl;
340  std::cerr << " SUSY_Spart1_pdgId: " << SUSY_Spart1_pdgId << std::endl;
341  std::cerr << " SUSY_Spart2_pdgId: " << SUSY_Spart2_pdgId << std::endl;
342  std::cerr << "Returning 0" << std::endl;
343 
344  return 0;
345 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
python.StoreID.UNKNOWN
int UNKNOWN
Definition: StoreID.py:16
SUSY::finalState
unsigned int finalState(const int SUSY_Spart1_pdgId, const int SUSY_Spart2_pdgId)
Definition: SUSYCrossSection.cxx:227
checkFileSG.line
line
Definition: checkFileSG.py:75
asg::AnaToolHandle::retrieve
StatusCode retrieve()
initialize the tool
AthenaPoolExample_ReadWrite.Key
Key
Definition: AthenaPoolExample_ReadWrite.py:53
SUSY::CrossSectionDB::CrossSectionDB
CrossSectionDB(const std::string &txtfilename="dev/PMGTools/PMGxsecDB_mc16.txt", bool usePathResolver=true, bool isExtended=false, bool usePMGTool=true)
Definition: SUSYCrossSection.cxx:113
SUSY::CrossSectionDB::my_find
xsDB_t::iterator my_find(const int proc)
Definition: SUSYCrossSection.cxx:159
Make4DCorrelationMatrix.inFiles
list inFiles
Definition: Make4DCorrelationMatrix.py:59
CscCalibQuery.fullPath
string fullPath
Definition: CscCalibQuery.py:360
SUSY::CrossSectionDB::Key
Definition: SUSYCrossSection.h:69
asg::AnaToolHandle::setTypeAndName
void setTypeAndName(const std::string &val_typeAndName)
set the value of type and name
skel.it
it
Definition: skel.GENtoEVGEN.py:396
efficiency
void efficiency(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:128
SUSY::CrossSectionDB::setUsePMGTool
void setUsePMGTool(bool usePMGTool=true)
Definition: SUSYCrossSection.h:86
PMGTools::IPMGCrossSectionTool::readInfosFromFiles
virtual bool readInfosFromFiles(std::vector< std::string >)=0
read infos from file, store them in the structure and make a vector that keeps all of them
SUSY::CrossSectionDB::setExtended
void setExtended(bool isExtended=true)
Definition: SUSYCrossSection.h:85
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
SUSY::CrossSectionDB::Process
Definition: SUSYCrossSection.h:40
beamspotman.stat
stat
Definition: beamspotman.py:266
SUSYCrossSection.h
PathResolver.h
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:227
mc.proc
proc
Definition: mc.PhPy8EG_A14NNPDF23_gg4l_example.py:22
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
SUSY::CrossSectionDB::extend
void extend(const std::string &)
Definition: SUSYCrossSection.cxx:167
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:431
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
SUSY::CrossSectionDB::process
Process process(int id, int proc=0) const
Definition: SUSYCrossSection.cxx:206
fitman.k
k
Definition: fitman.py:528
SUSY::CrossSectionDB::loadFile
void loadFile(const std::string &)
Definition: SUSYCrossSection.cxx:131
SUSY::CrossSectionDB::m_pmgxs
asg::AnaToolHandle< PMGTools::IPMGCrossSectionTool > m_pmgxs
Definition: SUSYCrossSection.h:125
Generated on Sun Dec 22 2024 21:18:41 for ATLAS Offline Software by   doxygen 1.8.18