ATLAS Offline Software
SUSYSignalTagger.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "DerivationFrameworkSUSY/SUSYSignalTagger.h"
6 
7 #include "xAODEventInfo/EventInfo.h"
8 
9 #include "TruthUtils/MagicNumbers.h"
10 
11 namespace DerivationFramework {
12 
13  static const SG::AuxElement::Decorator<int> dec_procID("SUSY_procID");
14  static const SG::AuxElement::Decorator<int> dec_pdgId1("SUSY_pid1");
15  static const SG::AuxElement::Decorator<int> dec_pdgId2("SUSY_pid2");
16 
17  SUSYSignalTagger::SUSYSignalTagger(const std::string& t, const std::string& n, const IInterface* p):
18  AthAlgTool(t,n,p)
19  {
20 
21  declareInterface<DerivationFramework::IAugmentationTool>(this);
22 
23  declareProperty("EventInfoName",m_eventInfoName="EventInfo");
24  declareProperty("MCCollectionName",m_mcName="TruthParticles");
25  }
26 
27 
28 
29  SUSYSignalTagger::~SUSYSignalTagger(){}
30 
31 
32 
33  StatusCode SUSYSignalTagger::initialize(){
34 
35  ATH_MSG_INFO("Initialize " );
36 
37  return StatusCode::SUCCESS;
38 
39  }
40 
41 
42 
43  StatusCode SUSYSignalTagger::finalize(){
44 
45  return StatusCode::SUCCESS;
46 
47  }
48 
49 
50 
51  StatusCode SUSYSignalTagger::addBranches() const{
52 
53  const xAOD::EventInfo* eventInfo;
54  if (evtStore()->retrieve(eventInfo,m_eventInfoName).isFailure()) {
55  ATH_MSG_ERROR("could not retrieve event info " <<m_eventInfoName);
56  return StatusCode::FAILURE;
57  }
58 
59  const xAOD::TruthParticleContainer* truthPC = 0;
60  if (evtStore()->retrieve(truthPC,m_mcName).isFailure()) {
61  ATH_MSG_DEBUG("WARNING could not retrieve TruthParticleContainer " <<m_mcName);
62  return StatusCode::FAILURE;
63  }
64 
65  //Identify SUSY hard proc
66  int pdgId1(0);
67  int pdgId2(0);
68 
69  bool found = FindSusyHardProc( truthPC, pdgId1, pdgId2);
70  if (!found) {
71  ATH_MSG_WARNING("could not identify SUSY process! ");
72 
73  dec_procID(*eventInfo) = 0;
74  dec_pdgId1(*eventInfo) = -99;
75  dec_pdgId2(*eventInfo) = -99;
76 
77  return StatusCode::SUCCESS;
78  }
79 
80 
81  //Get SUSY proc ID
82  unsigned int procID = finalStateID(pdgId1, pdgId2);
83 
84  dec_procID(*eventInfo) = procID;
85  dec_pdgId1(*eventInfo) = pdgId1;
86  dec_pdgId2(*eventInfo) = pdgId2;
87 
88  return StatusCode::SUCCESS;
89  }
90 
91 
92  bool SUSYSignalTagger::FindSusyHardProc(const xAOD::TruthParticleContainer* truthP, int& pdgid1, int& pdgid2) const
93  {
94  pdgid1 = 0;
95  pdgid2 = 0;
96 
97  const xAOD::TruthParticle* firstsp(0);
98  const xAOD::TruthParticle* secondsp(0);
99 
100  if (!truthP || truthP->empty()) {
101  return false;
102  }
103  for (const auto tp : *truthP) {
104 
105  //check ifSUSY particle
106  if ((std::abs(tp->pdgId()) > 1000000 && std::abs(tp->pdgId()) < 1000007) || // squarkL
107  (std::abs(tp->pdgId()) > 1000010 && std::abs(tp->pdgId()) < 1000017) || // sleptonL
108  (std::abs(tp->pdgId()) > 2000000 && std::abs(tp->pdgId()) < 2000007) || // squarkR
109  (std::abs(tp->pdgId()) > 2000010 && std::abs(tp->pdgId()) < 2000017) || // sleptonR
110  (std::abs(tp->pdgId()) > 1000020 && std::abs(tp->pdgId()) < 1000040)) { // gauginos
111 
112  if (tp->nParents() != 0) {
113  if ( tp->parent(0)->absPdgId() < 1000000) {
114  if (!firstsp) {
115  firstsp = tp;
116  } else if (!secondsp) {
117  secondsp = tp;
118  } else {
119  if (firstsp->nChildren() != 0 && HepMC::uniqueID(tp) == HepMC::uniqueID(firstsp->child(0))) {
120  firstsp = tp;
121  }
122  else if (secondsp->nChildren() != 0 && HepMC::uniqueID(tp) == HepMC::uniqueID(secondsp->child(0))) {
123  secondsp = tp;
124  }
125  else if (firstsp->nChildren() != 0 && HepMC::uniqueID(firstsp->child(0)) == HepMC::uniqueID(secondsp)) {
126  firstsp = secondsp;
127  secondsp = tp;
128  }
129  else if (secondsp->nChildren() != 0 && HepMC::uniqueID(secondsp->child(0)) == HepMC::uniqueID(firstsp)) {
130  secondsp = firstsp;
131  firstsp = tp;
132  }
133  }
134  }
135  }
136  }
137  }
138 
139  // quit if no sparticles found
140  if (!firstsp && !secondsp) return false; // should find none or two
141 
142  if (firstsp && firstsp->nChildren() == 1) {
143  for (const auto tp : *truthP) {
144  if (HepMC::uniqueID(tp) == HepMC::uniqueID(firstsp->child(0)) && tp->pdgId() != firstsp->pdgId()) {
145  firstsp = tp;
146  break;
147  }
148  }
149  }
150 
151  if (secondsp && secondsp->nChildren() == 1) {
152  for (const auto tp : *truthP) {
153  if (HepMC::uniqueID(tp) == HepMC::uniqueID(secondsp->child(0)) && tp->pdgId() != secondsp->pdgId()) {
154  secondsp = tp;
155  break;
156  }
157  }
158  }
159 
160  if (firstsp && abs(firstsp->pdgId()) > 1000000) pdgid1 = firstsp->pdgId();
161  if (secondsp && abs(secondsp->pdgId()) > 1000000) pdgid2 = secondsp->pdgId();
162 
163  // Return gracefully:
164  return true;
165  }
166 
167 
168  unsigned int SUSYSignalTagger::finalStateID(const int SUSY_Spart1_pdgId, const int SUSY_Spart2_pdgId) const
169  {
170  int ngluino=0;
171  int nsquark=0; // (up and down type without bottom/top)
172  int nantisquark=0; // (up and down type without bottom/top)
173 
174  int nsbottom=0;
175  int nstop=0;
176  int nsbottom2=0;
177  int nstop2=0;
178  int nantisbottom=0;
179  int nantistop=0;
180  int nantisbottom2=0;
181  int nantistop2=0;
182 
183  int nchi01=0;
184  int nchi02=0;
185  int nchi03=0;
186  int nchi04=0;
187  int nch1plus=0;
188  int nch2plus=0;
189  int nch1minus=0;
190  int nch2minus=0;
191 
192  // sleptons
193  int nsmuonRplus=0;
194  int nsmuonRminus=0;
195  int nselecRplus=0;
196  int nselecRminus=0;
197 
198  int nsmuonLplus=0;
199  int nsmuonLminus=0;
200  int nselecLplus=0;
201  int nselecLminus=0;
202 
203  int nstau1plus=0;
204  int nstau1minus=0;
205  int nstau2plus=0;
206  int nstau2minus=0;
207 
208  // snutrinos
209  int nselnuL=0;
210  int nsmunuL=0;
211  int nstaunuL=0;
212 
213  //Classification of the event follows (gg, sq...):
214 
215  if (std::abs(SUSY_Spart1_pdgId)== 1000022) nchi01++;
216  else if (std::abs(SUSY_Spart1_pdgId)== 1000023) nchi02++;
217  else if (std::abs(SUSY_Spart1_pdgId)== 1000025) nchi03++;
218  else if (std::abs(SUSY_Spart1_pdgId)== 1000035) nchi04++;
219  else if (SUSY_Spart1_pdgId == 1000024) nch1plus++;
220  else if (SUSY_Spart1_pdgId ==-1000024) nch1minus++;
221  else if (SUSY_Spart1_pdgId == 1000037) nch2plus++;
222  else if (SUSY_Spart1_pdgId ==-1000037) nch2minus++;
223  else if (SUSY_Spart1_pdgId == 1000021) ngluino++;
224  else if ((std::abs(SUSY_Spart1_pdgId)>1000000 && std::abs(SUSY_Spart1_pdgId)<= 1000004) || (std::abs(SUSY_Spart1_pdgId)>2000000 && std::abs(SUSY_Spart1_pdgId)<=2000004)) {
225  if (SUSY_Spart1_pdgId>0) nsquark++;
226  else nantisquark++;
227  }
228  else if (SUSY_Spart1_pdgId==1000005) nsbottom++;
229  else if (SUSY_Spart1_pdgId==1000006) nstop++;
230  else if (SUSY_Spart1_pdgId==2000005) nsbottom2++;
231  else if (SUSY_Spart1_pdgId==2000006) nstop2++;
232  else if (SUSY_Spart1_pdgId==-1000005) nantisbottom++;
233  else if (SUSY_Spart1_pdgId==-1000006) nantistop++;
234  else if (SUSY_Spart1_pdgId==-2000005) nantisbottom2++;
235  else if (SUSY_Spart1_pdgId==-2000006) nantistop2++;
236  else if (SUSY_Spart1_pdgId==2000011) nselecRminus++;
237  else if (SUSY_Spart1_pdgId==-2000011) nselecRplus++;
238  else if (SUSY_Spart1_pdgId==1000011) nselecLminus++;
239  else if (SUSY_Spart1_pdgId==-1000011) nselecLplus++;
240  else if (std::abs(SUSY_Spart1_pdgId)==1000012) nselnuL++;
241  else if (SUSY_Spart1_pdgId==2000013) nsmuonRminus++;
242  else if (SUSY_Spart1_pdgId==-2000013) nsmuonRplus++;
243  else if (SUSY_Spart1_pdgId==1000013) nsmuonLminus++;
244  else if (SUSY_Spart1_pdgId==-1000013) nsmuonLplus++;
245  else if (std::abs(SUSY_Spart1_pdgId)==1000014) nsmunuL++;
246  else if (SUSY_Spart1_pdgId==1000015) nstau1minus++;
247  else if (SUSY_Spart1_pdgId==-1000015) nstau1plus++;
248  else if (SUSY_Spart1_pdgId==2000015) nstau2minus++;
249  else if (SUSY_Spart1_pdgId==-2000015) nstau2plus++;
250  else if (std::abs(SUSY_Spart1_pdgId)==1000016) nstaunuL++;
251 
252 
253 
254 
255  if (std::abs(SUSY_Spart2_pdgId)==1000022) nchi01++;
256  else if (std::abs(SUSY_Spart2_pdgId)==1000023) nchi02++;
257  else if (std::abs(SUSY_Spart2_pdgId)==1000025) nchi03++;
258  else if (std::abs(SUSY_Spart2_pdgId)==1000035) nchi04++;
259  else if (SUSY_Spart2_pdgId==1000024) nch1plus++;
260  else if (SUSY_Spart2_pdgId==-1000024) nch1minus++;
261  else if (SUSY_Spart2_pdgId==1000037) nch2plus++;
262  else if (SUSY_Spart2_pdgId==-1000037) nch2minus++;
263 
264  else if (SUSY_Spart2_pdgId==1000021) ngluino++;
265  else if ((std::abs(SUSY_Spart2_pdgId)>1000000 && std::abs(SUSY_Spart2_pdgId)<=1000004) || (std::abs(SUSY_Spart2_pdgId)>2000000 && std::abs(SUSY_Spart2_pdgId)<=2000004)) {
266  if (SUSY_Spart2_pdgId>0) nsquark++;
267  else nantisquark++;
268  }
269  else if (SUSY_Spart2_pdgId==1000005) nsbottom++;
270  else if (SUSY_Spart2_pdgId==1000006) nstop++;
271  else if (SUSY_Spart2_pdgId==2000005) nsbottom2++;
272  else if (SUSY_Spart2_pdgId==2000006) nstop2++;
273  else if (SUSY_Spart2_pdgId==-1000005) nantisbottom++;
274  else if (SUSY_Spart2_pdgId==-1000006) nantistop++;
275  else if (SUSY_Spart2_pdgId==-2000005) nantisbottom2++;
276  else if (SUSY_Spart2_pdgId==-2000006) nantistop2++;
277 
278  else if (SUSY_Spart2_pdgId==2000011) nselecRminus++;
279  else if (SUSY_Spart2_pdgId==-2000011) nselecRplus++;
280  else if (SUSY_Spart2_pdgId==1000011) nselecLminus++;
281  else if (SUSY_Spart2_pdgId==-1000011) nselecLplus++;
282  else if (std::abs(SUSY_Spart2_pdgId)==1000012) nselnuL++;
283  else if (SUSY_Spart2_pdgId==2000013) nsmuonRminus++;
284  else if (SUSY_Spart2_pdgId==-2000013) nsmuonRplus++;
285  else if (SUSY_Spart2_pdgId==1000013) nsmuonLminus++;
286  else if (SUSY_Spart2_pdgId==-1000013) nsmuonLplus++;
287  else if (std::abs(SUSY_Spart2_pdgId)==1000014) nsmunuL++;
288  else if (SUSY_Spart2_pdgId==1000015) nstau1minus++;
289  else if (SUSY_Spart2_pdgId==-1000015) nstau1plus++;
290  else if (SUSY_Spart2_pdgId==2000015) nstau2minus++;
291  else if (SUSY_Spart2_pdgId==-2000015) nstau2plus++;
292  else if (std::abs(SUSY_Spart2_pdgId)==1000016) nstaunuL++;
293 
294 
296  // gluino/squark + X
297  if (ngluino==1 && (nsquark==1 || nantisquark ==1)) return 1;
298  else if (ngluino==2) return 2;
299  else if (nsquark==2 || nantisquark==2) return 3;
300  else if (nsquark==1 && nantisquark==1) return 4;
301 
302  else if (nsbottom==1 && nantisbottom==1) return 51;
303  else if (nsbottom2==1 && nantisbottom2==1) return 52;
304  else if (nstop==1 && nantistop==1) return 61;
305  else if (nstop2==1 && nantistop2==1) return 62;
306 
307  else if (ngluino==1 && nchi01==1) return 71;
308  else if (ngluino==1 && nchi02==1) return 72;
309  else if (ngluino==1 && nchi03==1) return 73;
310  else if (ngluino==1 && nchi04==1) return 74;
311 
312  else if (ngluino==1 && nch1plus==1) return 75;
313  else if (ngluino==1 && nch2plus==1) return 76;
314  else if (ngluino==1 && nch1minus==1) return 77;
315  else if (ngluino==1 && nch2minus==1) return 78;
316 
317  else if ((nsquark==1 || nantisquark ==1) && nchi01==1) return 81;
318  else if ((nsquark==1 || nantisquark ==1) && nchi02==1) return 82;
319  else if ((nsquark==1 || nantisquark ==1) && nchi03==1) return 83;
320  else if ((nsquark==1 || nantisquark ==1) && nchi04==1) return 84;
321 
322  else if ((nsquark==1 || nantisquark ==1) && nch1plus==1) return 85;
323  else if ((nsquark==1 || nantisquark ==1) && nch2plus==1) return 86;
324  else if ((nsquark==1 || nantisquark ==1) && nch1minus==1) return 87;
325  else if ((nsquark==1 || nantisquark ==1) && nch2minus==1) return 88;
326 
327 
328  // Gaugino pair-production
329  // chi^{0}_1 + X
330  else if (nchi01==2) return 111;
331  else if (nchi01==1 && nchi02==1) return 112;
332  else if (nchi01==1 && nchi03==1) return 113;
333  else if (nchi01==1 && nchi04==1) return 114;
334  else if (nchi01==1 && nch1plus==1) return 115;
335  else if (nchi01==1 && nch2plus==1) return 116;
336  else if (nchi01==1 && nch1minus==1) return 117;
337  else if (nchi01==1 && nch2minus==1) return 118;
338 
339  // chi^{0}_2 + X
340  else if (nchi02==2) return 122;
341  else if (nchi02==1 && nchi03==1) return 123;
342  else if (nchi02==1 && nchi04==1) return 124;
343  else if (nchi02==1 && nch1plus==1) return 125;
344  else if (nchi02==1 && nch2plus==1) return 126;
345  else if (nchi02==1 && nch1minus==1) return 127;
346  else if (nchi02==1 && nch2minus==1) return 128;
347 
348  // chi^{0}_3 + X
349  else if (nchi03==2) return 133;
350  else if (nchi03==1 && nchi04==1) return 134;
351  else if (nchi03==1 && nch1plus==1) return 135;
352  else if (nchi03==1 && nch2plus==1) return 136;
353  else if (nchi03==1 && nch1minus==1) return 137;
354  else if (nchi03==1 && nch2minus==1) return 138;
355 
356  // chi^{0}_4 + X
357  else if (nchi04==2) return 144;
358  else if (nchi04==1 && nch1plus==1) return 145;
359  else if (nchi04==1 && nch2plus==1) return 146;
360  else if (nchi04==1 && nch1minus==1) return 147;
361  else if (nchi04==1 && nch2minus==1) return 148;
362 
363  // chi^{+}_1/2 + chi^{-}_1/2
364  else if (nch1plus==1 && nch1minus==1) return 157;
365  else if (nch1plus==1 && nch2minus==1) return 158;
366 
367  else if (nch2plus==1 && nch1minus==1) return 167;
368  else if (nch2plus==1 && nch2minus==1) return 168;
369 
370  // slepton
371  else if (nselecLplus==1 && nselecLminus==1) return 201; // sElectronLPair
372  else if (nselecRplus==1 && nselecRminus==1) return 202; // sElectronRPair
373  else if (nselnuL==2) return 203; // sElectron neutrino pair
374  else if (nselecLplus==1 && nselnuL==1) return 204; // sElectron+ sNutrino
375  else if (nselecLminus==1 && nselnuL==1) return 205; // sElectron- sNutrino
376  else if (nstau1plus==1 && nstau1minus==1) return 206;
377  else if (nstau2plus==1 && nstau2minus==1) return 207;
378  else if ((nstau1plus==1 || nstau1minus==1) && (nstau2plus==1 || nstau2minus==1)) return 208;
379  else if (nstaunuL==2) return 209; // sTau neutrino pair
380  else if (nstau1plus==1 && nstaunuL==1) return 210;
381  else if (nstau1minus==1 && nstaunuL==1) return 211;
382  else if (nstau2plus==1 && nstaunuL==1) return 212;
383  else if (nstau2minus==1 && nstaunuL==1) return 213;
384 
385  else if (nsmuonLplus==1 && nsmuonLminus==1) return 216; // sMuonPair
386  else if (nsmuonRplus==1 && nsmuonRminus==1) return 217; // sMuonPair
387  else if (nsmunuL==2) return 218; // sMuon neutrino pair
388  else if (nsmuonLplus==1 && nsmunuL==1) return 219; // sMuon+ sNutrino
389  else if (nsmuonLminus==1 && nsmunuL==1) return 220; // sMuon- sNutrino
390 
391  ATH_MSG_WARNING("could not determine finalState for : (" << SUSY_Spart1_pdgId << " , " << SUSY_Spart2_pdgId << "). Returning 0.");
392 
393  return 0;
394  }
395 
396 
397 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
DerivationFramework::SUSYSignalTagger::FindSusyHardProc
bool FindSusyHardProc(const xAOD::TruthParticleContainer *truthP, int &pdgid1, int &pdgid2) const
Definition: SUSYSignalTagger.cxx:92
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
DerivationFramework::SUSYSignalTagger::SUSYSignalTagger
SUSYSignalTagger(const std::string &t, const std::string &n, const IInterface *p)
Definition: SUSYSignalTagger.cxx:17
ParticleTest.tp
tp
Definition: ParticleTest.py:25
DerivationFramework::SUSYSignalTagger::initialize
StatusCode initialize()
Definition: SUSYSignalTagger.cxx:33
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
SG::Decorator
Helper class to provide type-safe access to aux data.
Definition: Decorator.h:58
DerivationFramework::SUSYSignalTagger::finalStateID
unsigned int finalStateID(const int SUSY_Spart1_pdgId, const int SUSY_Spart2_pdgId) const
Definition: SUSYSignalTagger.cxx:168
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
DerivationFramework::SUSYSignalTagger::finalize
StatusCode finalize()
Definition: SUSYSignalTagger.cxx:43
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:41
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:113
DerivationFramework::SUSYSignalTagger::~SUSYSignalTagger
~SUSYSignalTagger()
Definition: SUSYSignalTagger.cxx:29
DerivationFramework::SUSYSignalTagger::addBranches
virtual StatusCode addBranches() const
Pass the thinning service
Definition: SUSYSignalTagger.cxx:51
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
xAOD::TruthParticle_v1::nChildren
size_t nChildren() const
Number of children of this particle.
Definition: TruthParticle_v1.cxx:140
MagicNumbers.h
SUSYSignalTagger.h
tool to decorate EventInfo with the SUSY signal process information
EventInfo.h
xAOD::EventInfo_v1
Class describing the basic event information.
Definition: EventInfo_v1.h:43
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::TruthParticle_v1::child
const TruthParticle_v1 * child(size_t i=0) const
Retrieve the i-th mother (TruthParticle) of this TruthParticle.
Definition: TruthParticle_v1.cxx:149
DerivationFramework::SUSYSignalTagger::m_eventInfoName
std::string m_eventInfoName
Definition: SUSYSignalTagger.h:38
DerivationFramework::SUSYSignalTagger::m_mcName
std::string m_mcName
Definition: SUSYSignalTagger.h:39
AthAlgTool
Definition: AthAlgTool.h:26
xAOD::TruthParticle_v1::pdgId
int pdgId() const
PDG ID code.
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
Generated on Thu Jul 4 2024 21:27:36 for ATLAS Offline Software by   doxygen 1.8.18