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