CalcTopPartonHistory.cxx 29.8 KB
Newer Older
1
/*
Tomas Dado's avatar
Tomas Dado committed
2
3
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
4
5
6
7
8
9

// $Id: CalcTopPartonHistory.cxx 800464 2017-03-13 18:06:24Z tpelzer $
#include "TopPartons/CalcTopPartonHistory.h"
#include "TopPartons/PartonHistory.h"
#include "TopConfiguration/TopConfig.h"
#include "TopPartons/CalcTtbarPartonHistory.h"
Tomas Dado's avatar
Tomas Dado committed
10
#include "TopPartons/PartonHistoryUtils.h"
11
12
#include "xAODTruth/TruthVertex.h"

Tomas Dado's avatar
Tomas Dado committed
13
14
15
16
17
namespace top {
  CalcTopPartonHistory::CalcTopPartonHistory(const std::string& name) :
    asg::AsgTool(name),
    m_config(nullptr) {
    declareProperty("config", m_config);
18
19
  }
  
20
  StatusCode CalcTopPartonHistory::buildContainerFromMultipleCollections(const std::vector<std::string> &collections, const std::string& out_contName)
21
  {
22
    ConstDataVector<DataVector<xAOD::TruthParticle_v1> > *out_cont = new ConstDataVector<DataVector<xAOD::TruthParticle_v1> > (SG::VIEW_ELEMENTS);
23
    
Marco Vanadia's avatar
Marco Vanadia committed
24
    for(const std::string& collection : collections)
25
26
27
    {
      const xAOD::TruthParticleContainer* cont=nullptr;
      ATH_CHECK(evtStore()->retrieve(cont,collection));
28
      for(const xAOD::TruthParticle* p : *cont) out_cont->push_back(p);
29
30
    }
    
31
32
33
34
    //we give control of the container to the store, because in this way we are able to retrieve it as a const data vector, see https://twiki.cern.ch/twiki/bin/view/AtlasComputing/DataVector#ConstDataVector
    xAOD::TReturnCode save = evtStore()->tds()->record(out_cont,out_contName);
    if (!save) return StatusCode::FAILURE;

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    return StatusCode::SUCCESS;
  }
  
  StatusCode CalcTopPartonHistory::linkBosonCollections() 
  {
    return decorateCollectionWithLinksToAnotherCollection("TruthBoson","TruthBosonsWithDecayParticles","AT_linkToTruthBosonsWithDecayParticles");
  }
  StatusCode CalcTopPartonHistory::decorateCollectionWithLinksToAnotherCollection(const std::string &collectionToDecorate, const std::string &collectionToLink, const std::string &nameOfDecoration)
  {
    const xAOD::TruthParticleContainer* cont1(nullptr);
    const xAOD::TruthParticleContainer* cont2(nullptr);
    ATH_CHECK(evtStore()->retrieve(cont1,collectionToDecorate));
    ATH_CHECK(evtStore()->retrieve(cont2,collectionToLink));

    for(const xAOD::TruthParticle *p : *cont1)
    {
      
      const xAOD::TruthParticle* link =0;
      for(const xAOD::TruthParticle *p2 : *cont2)
      {
        if(p->pdgId()==p2->pdgId() && p->barcode()==p2->barcode())
        {
          link=p2;
          break;
        }
      } 
      p->auxdecor<const xAOD::TruthParticle*>(nameOfDecoration)=link;
      
    }
    return StatusCode::SUCCESS;
  }
  
67
  const xAOD::TruthParticle* CalcTopPartonHistory::getTruthParticleLinkedFromDecoration(const xAOD::TruthParticle* part, const std::string &decorationName)
68
  {
69
    if(!part->isAvailable<const xAOD::TruthParticle*>(decorationName)) return part;
70
71
  
    const xAOD::TruthParticle* link=part->auxdecor<const xAOD::TruthParticle*>(decorationName);
72
    if(link) return link;
73
    
74
    return part;
75
  }
Tomas Dado's avatar
Tomas Dado committed
76

77
78
  ///Store the four-momentum of the post-FSR top or anti-top found using statusCodes
  ///This would only work if there is at most one "true" top of each charge (i.e. won't work for SS tops or 4 tops)
Tomas Dado's avatar
Tomas Dado committed
79
80
81
82
  ///This code was adapted from the 7TeV parton-level differential ttbar routine:
  // https://svnweb.cern.ch/trac/atlasphys-top/browser/Physics/Top/Software/MCvalidation/Rivet/Rivet2.X/trunk/routines/ATLAS_2014_I1304289/ATLAS_2014_I1304289.cc
  bool CalcTopPartonHistory::topAfterFSR_SC(const xAOD::TruthParticleContainer* truthParticles, int start,
                                            TLorentzVector& top_afterFSR_SC_p4) {
83
84
85
86
87
88
89
90
91
    /// Step1: create vectors of particles of each status codes
    // Vectors to hold any status=3 (anti-)top quarks (Pythia 6)
    std::vector<const xAOD::TruthParticle*> v_status3_top;
    // Vectors to hold any status=155 (anti-)top quarks (Herwig 6)
    std::vector<const xAOD::TruthParticle*> v_status155_top;
    // Vectors to hold any status=11 (anti-)top quarks for Herwig++
    std::vector<const xAOD::TruthParticle*> v_status11_top;
    // Vectors to hold any status=22 (anti-)top quarks
    std::vector<const xAOD::TruthParticle*> v_statusOther_top;
Tomas Dado's avatar
Tomas Dado committed
92

93
94
95
    /// Step2: loop on the container of particles and fill the above created vectors
    for (const xAOD::TruthParticle* particle : *truthParticles) {
      if (particle->pdgId() != start) continue; // only keep particles of a given pdgID (e.g. 6 or -6)
Tomas Dado's avatar
Tomas Dado committed
96
97

      if (particle->status() == 3) {
98
        v_status3_top.push_back(particle);
Tomas Dado's avatar
Tomas Dado committed
99
      } else if (particle->status() == 155) {
100
        v_status155_top.push_back(particle);
Tomas Dado's avatar
Tomas Dado committed
101
      } else if (particle->status() == 11) {// for Herwig++: take only the tops that decay into Wb!!!
102
103
104
105
        if (!particle->hasDecayVtx()) continue;
        const xAOD::TruthVertex* vertex = particle->decayVtx();
        if (vertex == nullptr) continue;
        if (vertex->nOutgoingParticles() == 2) v_status11_top.push_back(particle);
Tomas Dado's avatar
Tomas Dado committed
106
      } else {
107
108
109
        v_statusOther_top.push_back(particle);
      }
    }
Tomas Dado's avatar
Tomas Dado committed
110

111
112
    /// Step3: for some of the statuscodes, keep only the last of the vector
    // If there are more than 1 status 3 tops or anti-tops, only keep the last one put into the vector
Tomas Dado's avatar
Tomas Dado committed
113
    if (v_status3_top.size() > 1) {
114
115
116
      v_status3_top = std::vector<const xAOD::TruthParticle*>(v_status3_top.end() - 1, v_status3_top.end());
    }
    // If there are more than 1 status 11 tops or anti-tops, only keep the last one put into the vector
Tomas Dado's avatar
Tomas Dado committed
117
    if (v_status11_top.size() > 1) {
118
119
120
121
      v_status11_top = std::vector<const xAOD::TruthParticle*>(v_status11_top.end() - 1, v_status11_top.end());
    }
    // Rach: check for Pythia 8 as well
    // If there are more than 1 status 3 tops or anti-tops, only keep the last one put into the vector
Tomas Dado's avatar
Tomas Dado committed
122
    if (v_statusOther_top.size() > 1) {
123
124
      v_statusOther_top = std::vector<const xAOD::TruthParticle*>(v_statusOther_top.end() - 1, v_statusOther_top.end());
    }
Tomas Dado's avatar
Tomas Dado committed
125

126
127
128
    /// Step4: chose which statuscode to take according to what is found in the event
    const xAOD::TruthParticle* top = nullptr;
    // If there are status 3 tops and no status 155 tops this is probably a Pythia event, so used the status 3s.
Tomas Dado's avatar
Tomas Dado committed
129
    if (v_status3_top.size() == 1 && v_status155_top.size() == 0) {
130
131
132
      top = v_status3_top[0];
    }
    // If there are status 155 tops this must be a Herwig event, so use the status 155s.
Tomas Dado's avatar
Tomas Dado committed
133
    if (v_status155_top.size() == 1 && v_status3_top.size() == 0) {
134
135
      top = v_status155_top[0];
    }
Tomas Dado's avatar
Tomas Dado committed
136
137
    // If there are tops with other status this must be a Pythia8 event, so use them.
    if (v_statusOther_top.size() == 1 && v_status155_top.size() == 0 && v_status3_top.size() == 0) {
138
139
140
      top = v_statusOther_top[0];
    }
    // If there are status 155 tops this must be a Herwig event, so use the status 155s.
Tomas Dado's avatar
Tomas Dado committed
141
    if (v_status11_top.size() == 1 && v_status3_top.size() == 0) {
142
143
      top = v_status11_top[0];
    }
Tomas Dado's avatar
Tomas Dado committed
144

145
146
147
148
149
150
151
    /// Step5: if everything worked, set the 4-vector to its value and return true
    if (top != nullptr) {
      top_afterFSR_SC_p4 = top->p4();
      return true;
    }
    return false;
  }
Tomas Dado's avatar
Tomas Dado committed
152
153
154
155
156

  // for b coming from W'->tb
  bool CalcTopPartonHistory::b(const xAOD::TruthParticleContainer* truthParticles,
                               TLorentzVector& b_beforeFSR, TLorentzVector& b_afterFSR) {
    for (const xAOD::TruthParticle* particle : *truthParticles) {
157
      if (std::abs(particle->pdgId()) != 5) continue;
Tomas Dado's avatar
Tomas Dado committed
158
159
160
161

      bool skipit(false);
      for (size_t i = 0; i < particle->nParents(); i++) {
        const xAOD::TruthParticle* parent = particle->parent(i);
162
        if (parent && (parent->isTop() || std::abs(parent->pdgId()) == 5)) {
Tomas Dado's avatar
Tomas Dado committed
163
164
165
166
167
168
169
170
          skipit = true;
          break;
        }//if
      }//for

      if (skipit) continue;

      b_beforeFSR = particle->p4();
Tomas Dado's avatar
Tomas Dado committed
171
      b_afterFSR = PartonHistoryUtils::findAfterFSR(particle)->p4();
Tomas Dado's avatar
Tomas Dado committed
172
173
174
175
176
177
178
179

      return true;
    }


    return false;
  }

180
  // The function topWb has been overloaded to include aditional information when the W from t decays to a tau.
181
182
183
184
185
186
187
188
189
190
191
  bool CalcTopPartonHistory::topWb(const xAOD::TruthParticleContainer* truthParticles,
                                   int start, TLorentzVector& t_beforeFSR_p4, TLorentzVector& t_afterFSR_p4,
                                   TLorentzVector& W_p4,
                                   TLorentzVector& b_p4, TLorentzVector& Wdecay1_p4,
                                   int& Wdecay1_pdgId, TLorentzVector& Wdecay2_p4, int& Wdecay2_pdgId,
                                   TLorentzVector& tau_decay_from_W_p4, int& tau_decay_from_W_isHadronic) {
    bool hasT = false;
    bool hasW = false;
    bool hasB = false;
    bool hasWdecayProd1 = false;
    bool hasWdecayProd2 = false;
192
193
194
    
    bool store_tau_info = true;
    if (tau_decay_from_W_isHadronic == -9999) { bool store_tau_info = false; }
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

    for (const xAOD::TruthParticle* particle : *truthParticles) {
      if (particle->pdgId() != start) continue;

      if (PartonHistoryUtils::hasParticleIdenticalParent(particle)) continue; // kepping only top before FSR

      t_beforeFSR_p4 = particle->p4(); // top before FSR
      hasT = true;

      // demanding the last tops after FSR
      particle = PartonHistoryUtils::findAfterFSR(particle);
      t_afterFSR_p4 = particle->p4(); // top after FSR
      
      for (size_t k = 0; k < particle->nChildren(); k++) {
        const xAOD::TruthParticle* topChildren = particle->child(k);

211
        if (std::abs(topChildren->pdgId()) == 24) {
212
213
214
215
          W_p4 = topChildren->p4();  // W boson after FSR
          hasW = true;
          
          // demanding the last W after FSR
216
          topChildren = PartonHistoryUtils::findAfterFSR(topChildren);          
217
218
          for (size_t q = 0; q < topChildren->nChildren(); ++q) {
	    const xAOD::TruthParticle* WChildren = topChildren->child(q);
219
            if (std::abs(WChildren->pdgId()) < 17 && store_tau_info == true){
220
	      // When W decays leptonically, Wdecay1 stores the lepton and Wdecay2 the neutrino
221
	      if (std::abs(WChildren->pdgId()) == 11 || std::abs(WChildren->pdgId()) == 13 || std::abs(WChildren->pdgId()) == 15){
222
223
		Wdecay1_p4 = WChildren->p4();                                  
                Wdecay1_pdgId = WChildren->pdgId();
224
		tau_decay_from_W_isHadronic = PartonHistoryUtils::TauIsHadronic(WChildren, tau_decay_from_W_p4);
225
226
                hasWdecayProd1 = true;
	      }
227
	      if (std::abs(WChildren->pdgId()) == 12 || std::abs(WChildren->pdgId()) == 14 || std::abs(WChildren->pdgId()) == 16){
228
229
230
231
		Wdecay2_p4 = WChildren->p4();
                Wdecay2_pdgId = WChildren->pdgId();
                hasWdecayProd2 = true;
	      }
232
	      if (std::abs(WChildren->pdgId()) < 11){ // W does not decay leptonically
233
234
235
236
		if (WChildren->pdgId() > 0) {
		  Wdecay1_p4 = WChildren->p4();
		  Wdecay1_pdgId = WChildren->pdgId();
		  tau_decay_from_W_isHadronic = -99;
237
		  tau_decay_from_W_p4 = WChildren->p4();
238
239
240
241
242
243
244
		  hasWdecayProd1 = true;
		}else{
		  Wdecay2_p4 = WChildren->p4();
		  Wdecay2_pdgId = WChildren->pdgId();
		  hasWdecayProd2 = true;
		}//else
	      }// if not leptonic decay
245
246
247
248
249
250
251
252
253
254
255
256
257
	    }// end of if store_tau_info ==  true
	    if (std::abs(WChildren->pdgId()) < 17 && store_tau_info == false){
	      if (WChildren->pdgId() > 0) {
                Wdecay1_p4 = WChildren->p4();
                Wdecay1_pdgId = WChildren->pdgId();
                hasWdecayProd1 = true;
              } else {
                Wdecay2_p4 = WChildren->p4();
                Wdecay2_pdgId = WChildren->pdgId();
                hasWdecayProd2 = true;
              }
	    } // end of if store_tau_info == false
          }//end of for
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
        } else if (abs(topChildren->pdgId()) == 5) {
          b_p4 = topChildren->p4();
          hasB = true;
        } //else if
      } //for (size_t k=0; k < particle->nChildren(); k++)

      if (hasT && hasW && hasB && hasWdecayProd1 && hasWdecayProd2) return true;
      
    } //for (const xAOD::TruthParticle* particle : *truthParticles)
    
    return false;
  }



Tomas Dado's avatar
Tomas Dado committed
273
274
275
276
277
  bool CalcTopPartonHistory::topWb(const xAOD::TruthParticleContainer* truthParticles,
                                   int start, TLorentzVector& t_beforeFSR_p4, TLorentzVector& t_afterFSR_p4,
                                   TLorentzVector& W_p4,
                                   TLorentzVector& b_p4, TLorentzVector& Wdecay1_p4,
                                   int& Wdecay1_pdgId, TLorentzVector& Wdecay2_p4, int& Wdecay2_pdgId) {
278
279
280
    TLorentzVector tau_decay_from_W_p4;
    int tau_decay_from_W_isHadronic = -9999;
    return topWb(truthParticles, start, t_beforeFSR_p4, t_afterFSR_p4, W_p4, b_p4, Wdecay1_p4,Wdecay1_pdgId, Wdecay2_p4, Wdecay2_pdgId, tau_decay_from_W_p4, tau_decay_from_W_isHadronic);
Tomas Dado's avatar
Tomas Dado committed
281
282

    return false;
283
  }
Tomas Dado's avatar
Tomas Dado committed
284

285
 
Tomas Dado's avatar
Tomas Dado committed
286
287
288
289
290
291
292
293
294
295
296
  bool CalcTopPartonHistory::topWq(const xAOD::TruthParticleContainer* truthParticles,
                                   int start, TLorentzVector& t_beforeFSR_p4, TLorentzVector& t_afterFSR_p4,
                                   TLorentzVector& W_p4,
                                   TLorentzVector& q_p4, int& q_pdgId, TLorentzVector& Wdecay1_p4,
                                   int& Wdecay1_pdgId, TLorentzVector& Wdecay2_p4, int& Wdecay2_pdgId) {
    bool hasT = false;
    bool hasW = false;
    bool hasQ = false;
    bool hasWdecayProd1 = false;
    bool hasWdecayProd2 = false;

297
    for (const xAOD::TruthParticle* particle : *truthParticles) {
Tomas Dado's avatar
Tomas Dado committed
298
299
      if (particle->pdgId() != start) continue;

Tomas Dado's avatar
Tomas Dado committed
300
      if (PartonHistoryUtils::hasParticleIdenticalParent(particle)) continue; // kepping only top before FSR
Tomas Dado's avatar
Tomas Dado committed
301
302
303
304
305

      t_beforeFSR_p4 = particle->p4(); // top before FSR
      hasT = true;

      // demanding the last tops after FSR
Tomas Dado's avatar
Tomas Dado committed
306
      particle = PartonHistoryUtils::findAfterFSR(particle);
Tomas Dado's avatar
Tomas Dado committed
307
308
309
310
311
      t_afterFSR_p4 = particle->p4(); // top after FSR

      for (size_t k = 0; k < particle->nChildren(); k++) {
        const xAOD::TruthParticle* topChildren = particle->child(k);

312
        if (std::abs(topChildren->pdgId()) == 24) {
Tomas Dado's avatar
Tomas Dado committed
313
314
315
316
          W_p4 = topChildren->p4();  // W boson after FSR
          hasW = true;

          // demanding the last W after FSR
Tomas Dado's avatar
Tomas Dado committed
317
          topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
318
319
320

          for (size_t q = 0; q < topChildren->nChildren(); ++q) {
            const xAOD::TruthParticle* WChildren = topChildren->child(q);
321
            if (std::abs(WChildren->pdgId()) < 17) {
Tomas Dado's avatar
Tomas Dado committed
322
323
324
325
326
327
328
329
330
331
332
              if (WChildren->pdgId() > 0) {
                Wdecay1_p4 = WChildren->p4();
                Wdecay1_pdgId = WChildren->pdgId();
                hasWdecayProd1 = true;
              } else {
                Wdecay2_p4 = WChildren->p4();
                Wdecay2_pdgId = WChildren->pdgId();
                hasWdecayProd2 = true;
              }//else
            }//if
          }//for
333
        } else if (std::abs(topChildren->pdgId()) == 5 || std::abs(topChildren->pdgId()) == 3 || std::abs(topChildren->pdgId()) == 1) {
Tomas Dado's avatar
Tomas Dado committed
334
335
336
337
338
339
          q_p4 = topChildren->p4();
          q_pdgId = topChildren->pdgId();
          hasQ = true;
        } //else if
      } //for (size_t k=0; k < particle->nChildren(); k++)
      if (hasT && hasW && hasQ && hasWdecayProd1 && hasWdecayProd2) return true;
340
    } //for (const xAOD::TruthParticle* particle : *truthParticles)
Tomas Dado's avatar
Tomas Dado committed
341
342

    return false;
343
344
  }

Tomas Dado's avatar
Tomas Dado committed
345
346
347
348
349
350
351
352
  bool CalcTopPartonHistory::Wlv(const xAOD::TruthParticleContainer* truthParticles,
                                 TLorentzVector& W_p4,
                                 TLorentzVector& Wdecay1_p4, int& Wdecay1_pdgId,
                                 TLorentzVector& Wdecay2_p4, int& Wdecay2_pdgId) {
    bool hasW = false;
    bool hasWdecayProd1 = false;
    bool hasWdecayProd2 = false;

353
    for (const xAOD::TruthParticle* particle : *truthParticles) {
354
      if (std::abs(particle->pdgId()) != 24) continue;
Tomas Dado's avatar
Tomas Dado committed
355
356
357
      //std::cout << "PDGID: " << particle->pdgId() << std::endl;

      // demanding the last W after FSR
Tomas Dado's avatar
Tomas Dado committed
358
      particle = PartonHistoryUtils::findAfterFSR(particle);
Tomas Dado's avatar
Tomas Dado committed
359
      W_p4 = particle->p4();  // W boson after FSR
360
      hasW = true;
Tomas Dado's avatar
Tomas Dado committed
361
362
363

      for (size_t k = 0; k < particle->nChildren(); k++) {
        const xAOD::TruthParticle* WChildren = particle->child(k);
364
        if (std::abs(WChildren->pdgId()) < 17) {
Tomas Dado's avatar
Tomas Dado committed
365
          if (WChildren->pdgId() % 2 == 1) { // charged lepton in the Wlv case
366
367
368
            Wdecay1_p4 = WChildren->p4();
            Wdecay1_pdgId = WChildren->pdgId();
            hasWdecayProd1 = true;
Tomas Dado's avatar
Tomas Dado committed
369
          } else {// neutral lepton in the Wlv case
370
371
372
373
            Wdecay2_p4 = WChildren->p4();
            Wdecay2_pdgId = WChildren->pdgId();
            hasWdecayProd2 = true;
          }//else
Tomas Dado's avatar
Tomas Dado committed
374
375
        }//if
      } //for (size_t k=0; k < particle->nChildren(); k++)
376

Tomas Dado's avatar
Tomas Dado committed
377
      if (hasW && hasWdecayProd1 && hasWdecayProd2) return true;
378
379
    } //for (const xAOD::TruthParticle* particle : *truthParticles)

380

Tomas Dado's avatar
Tomas Dado committed
381
    return false;
382
383
  }

384
385
  // for Wt ST events, find W that is not from top
  bool CalcTopPartonHistory::Wt_W(const xAOD::TruthParticleContainer* truthParticles,
Tomas Dado's avatar
Tomas Dado committed
386
387
                                  TLorentzVector& W_p4, int& W_pdgId, TLorentzVector& Wdecay1_p4,
                                  int& Wdecay1_pdgId, TLorentzVector& Wdecay2_p4, int& Wdecay2_pdgId) {
388
389
390
391
392
393
    bool hasW = false;
    bool hasWdecayProd1 = false;
    bool hasWdecayProd2 = false;

    for (const xAOD::TruthParticle* particle : *truthParticles) {
      if (particle == nullptr) continue;
394
      if (std::abs(particle->pdgId()) != 24) continue; // W boson
395
396

      // need to check if the W is from top
Tomas Dado's avatar
Tomas Dado committed
397
      // identify the first in chain and check
398
      // if that particle has top as parent
Tomas Dado's avatar
Tomas Dado committed
399
      if (PartonHistoryUtils::hasParticleIdenticalParent(particle)) continue; // kepping only W before FSR
400
401
402

      bool isFromTop = false;
      // now we should have only the first W in chain
Tomas Dado's avatar
Tomas Dado committed
403
      for (size_t iparent = 0; iparent < particle->nParents(); ++iparent) {
404
        if (particle->parent(iparent) == nullptr) continue;
405
        if (std::abs(particle->parent(iparent)->pdgId()) == 6) { // has top as parent
406
407
          isFromTop = true;
          break;
Tomas Dado's avatar
Tomas Dado committed
408
409
        }
      }
410
411
412

      if (isFromTop) continue;
      else {
Tomas Dado's avatar
Tomas Dado committed
413
        particle = PartonHistoryUtils::findAfterFSR(particle);
414
415
416
417
        W_p4 = particle->p4();
        W_pdgId = particle->pdgId();
        hasW = true;
      }
Tomas Dado's avatar
Tomas Dado committed
418

419
420
      // check the decay products of the W
      for (size_t q = 0; q < particle->nChildren(); ++q) {
Tomas Dado's avatar
Tomas Dado committed
421
        const xAOD::TruthParticle* WChildren = particle->child(q);
422
        if (WChildren == nullptr) continue;
423
        if (std::abs(WChildren->pdgId()) < 17) {
Tomas Dado's avatar
Tomas Dado committed
424
          if (WChildren->pdgId() > 0) {
425
            Wdecay1_p4 = WChildren->p4();
Tomas Dado's avatar
Tomas Dado committed
426
427
428
429
430
431
            Wdecay1_pdgId = WChildren->pdgId();
            hasWdecayProd1 = true;
          } else {
            Wdecay2_p4 = WChildren->p4();
            Wdecay2_pdgId = WChildren->pdgId();
            hasWdecayProd2 = true;
432
          }//else
Tomas Dado's avatar
Tomas Dado committed
433
        }//if
434
435
      }//for

Tomas Dado's avatar
Tomas Dado committed
436
      if (hasW && hasWdecayProd1 && hasWdecayProd2) return true;
437
438
439
440
    } // loop over truth particles

    return false;
  }
Tomas Dado's avatar
Tomas Dado committed
441

442
443
444
  // for Wt ST events, find b that is not from top
  bool CalcTopPartonHistory::Wt_b(const xAOD::TruthParticleContainer* truthParticles,
                                  TLorentzVector& b_beforeFSR, TLorentzVector& b_afterFSR,
Tomas Dado's avatar
Tomas Dado committed
445
                                  int& b_pdgId) {
446
447
448
449
450
    bool hasB = false;

    // identify "other" b quark that is not from radiation but from ME (Wtb)
    // logic is simple: search for b quark that doesn't have top, proton, or
    // nullptr as parent
Tomas Dado's avatar
Tomas Dado committed
451

452
453
    for (const xAOD::TruthParticle* particle : *truthParticles) {
      if (particle == nullptr) continue;
454
      if (std::abs(particle->pdgId()) != 5) continue;
455

Tomas Dado's avatar
Tomas Dado committed
456
      for (size_t iparent = 0; iparent < particle->nParents(); ++iparent) {
457
458
459
        if (particle->parent(iparent) == nullptr) continue;

        // we dont want b-quarks that have b as parent
460
        if (std::abs(particle->parent(iparent)->pdgId()) == 5) continue;
461
462

        // we dont want b-quarks that come from top
463
        if (std::abs(particle->parent(iparent)->pdgId()) == 6) continue;
Tomas Dado's avatar
Tomas Dado committed
464

465
        // we dont want b-quarks that come from W
466
        if (std::abs(particle->parent(iparent)->pdgId()) == 24) continue;
Tomas Dado's avatar
Tomas Dado committed
467

468
        // we dont want b-quarks that come from proton
469
        if (std::abs(particle->parent(iparent)->pdgId()) == 2212) continue;
470
471
472
473
474
475

        hasB = true;
        b_beforeFSR = particle->p4();
        b_pdgId = particle->pdgId();

        // find after FSR
Tomas Dado's avatar
Tomas Dado committed
476
        particle = PartonHistoryUtils::findAfterFSR(particle);
477
478
479
480
481
482
        b_afterFSR = particle->p4();
      }
    }


    if (hasB) return true;
Tomas Dado's avatar
Tomas Dado committed
483

484
485
486
    return false;
  }

487
  // for ttbar + photon events
Tomas Dado's avatar
Tomas Dado committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
  bool CalcTopPartonHistory::topPhWb(const xAOD::TruthParticleContainer* truthParticles, int topId,
                                     TLorentzVector& t_beforeFSR_p4, TLorentzVector& t_afterFSR_p4,
                                     TLorentzVector& Ph_p4, TLorentzVector& W_p4, TLorentzVector& b_p4,
                                     TLorentzVector& Wdecay1_p4, int& Wdecay1_pdgId, TLorentzVector& Wdecay2_p4,
                                     int& Wdecay2_pdgId, bool& has_ph, int& BranchType, int& IniPartonType,
                                     bool& missingTop) {
    bool hasT = false;
    bool hasW = false;
    bool hasAbsentW = false;
    bool hasB = false;

    has_ph = false;
    bool ph_W = false;
    bool ph_Top = false;
    bool ph_ISR = false;
    bool ph_b = false;
    bool hasWdecayProd1 = false;
    bool hasWdecayProd2 = false;
    missingTop = false;
507
508

    for (const xAOD::TruthParticle* particle : *truthParticles) {
Tomas Dado's avatar
Tomas Dado committed
509
510
      if (particle->pdgId() != topId) continue;

Tomas Dado's avatar
Tomas Dado committed
511
      if (PartonHistoryUtils::hasParticleIdenticalParent(particle)) continue; // kepping only top before FSR
Tomas Dado's avatar
Tomas Dado committed
512
513
514
515
516
      BranchType = -1;// 10(50): leptonic(hadronic), 12(52):topRad, 14(54):Wrad, 15(55):ISR, 18(58):b
      IniPartonType = -1;

      // finding siblings
      for (size_t iparent = 0; iparent < particle->nParents(); iparent++) {
517
        if (std::abs(particle->parent(iparent)->pdgId()) == 21) {
Tomas Dado's avatar
Tomas Dado committed
518
519
          IniPartonType = 1;
        } // gg fusion
520
        else if (std::abs(particle->parent(iparent)->pdgId()) < 6) {
Tomas Dado's avatar
Tomas Dado committed
521
522
523
524
525
          IniPartonType = 2;
        } //qq annihilation

        for (size_t ichild = 0; ichild < particle->parent(iparent)->nChildren(); ichild++) {
          if (particle->parent(iparent)->child(ichild)->pdgId() == 22) {
Tomas Dado's avatar
Tomas Dado committed
526
            const xAOD::TruthParticle* photon = PartonHistoryUtils::findAfterFSR(particle->parent(iparent)->child(ichild));
Tomas Dado's avatar
Tomas Dado committed
527
528
529
530
531
            Ph_p4 = photon->p4();
            has_ph = true;
            ph_ISR = true;
          }
          if (!missingTop &&
532
533
              (std::abs(particle->parent(iparent)->child(ichild)->pdgId()) == 5 ||
               std::abs(particle->parent(iparent)->child(ichild)->pdgId()) == 24)) {
Tomas Dado's avatar
Tomas Dado committed
534
535
536
537
            missingTop = true;
          }
        }
      }
538

Tomas Dado's avatar
Tomas Dado committed
539
540
541
      t_beforeFSR_p4 = particle->p4(); // top before FSR
      hasT = true;
      // demanding the last tops after FSR
Tomas Dado's avatar
Tomas Dado committed
542
      particle = PartonHistoryUtils::findAfterFSR(particle);
Tomas Dado's avatar
Tomas Dado committed
543
544
545
546
547
      t_afterFSR_p4 = particle->p4(); // top after FSR

      for (size_t k = 0; k < particle->nChildren(); k++) {// top children
        const xAOD::TruthParticle* topChildren = particle->child(k);

548
        if (std::abs(topChildren->pdgId()) == 24) {
Tomas Dado's avatar
Tomas Dado committed
549
550
551
552
          W_p4 = topChildren->p4();  // W boson before FSR
          hasW = true;

          // demanding the last W after FSR
Tomas Dado's avatar
Tomas Dado committed
553
          topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
554
555
556

          for (size_t q = 0; q < topChildren->nChildren(); q++) {// W children
            const xAOD::TruthParticle* WChildren = topChildren->child(q);
557
558
            if (std::abs(WChildren->pdgId()) > 0 && std::abs(WChildren->pdgId()) < 17) {
              if (std::abs(WChildren->pdgId()) < 7) {
Tomas Dado's avatar
Tomas Dado committed
559
560
                BranchType = 50;
              }// hadronic
561
              else if (std::abs(WChildren->pdgId()) > 10 && std::abs(WChildren->pdgId()) < 17) {
Tomas Dado's avatar
Tomas Dado committed
562
563
564
                BranchType = 10;
              }// leptonic
              if (WChildren->pdgId() > 0) {
Tomas Dado's avatar
Tomas Dado committed
565
                WChildren = PartonHistoryUtils::findAfterFSR(WChildren);
Tomas Dado's avatar
Tomas Dado committed
566
567
568
569
                Wdecay1_p4 = WChildren->p4();
                Wdecay1_pdgId = WChildren->pdgId();
                hasWdecayProd1 = true;
              } else {
Tomas Dado's avatar
Tomas Dado committed
570
                WChildren = PartonHistoryUtils::findAfterFSR(WChildren);
Tomas Dado's avatar
Tomas Dado committed
571
572
573
574
                Wdecay2_p4 = WChildren->p4();
                Wdecay2_pdgId = WChildren->pdgId();
                hasWdecayProd2 = true;
              }//else
575
            } else if (std::abs(WChildren->pdgId()) == 22) {// photon
Tomas Dado's avatar
Tomas Dado committed
576
577
578
579
580
581
582
583
584
585
              // JUST FOR EXTRA SAFETY (not necessary)
              // check if there exists a photon already
              // if it does, check the photon's Pt
              // if found harder then consider, else do nothing
              if (has_ph) {
                if (WChildren->p4().Pt() > Ph_p4.Pt()) {
                  ph_W = true;
                  ph_ISR = false;
                  ph_Top = false;
                  ph_b = false;
Tomas Dado's avatar
Tomas Dado committed
586
                  WChildren = PartonHistoryUtils::findAfterFSR(WChildren);
Tomas Dado's avatar
Tomas Dado committed
587
588
589
590
591
                  Ph_p4 = WChildren->p4();
                }
              } else {
                has_ph = true;
                ph_W = true;
Tomas Dado's avatar
Tomas Dado committed
592
                WChildren = PartonHistoryUtils::findAfterFSR(WChildren);
Tomas Dado's avatar
Tomas Dado committed
593
594
595
596
                Ph_p4 = WChildren->p4();
              }
            }
          }// W children
597
        } else if (std::abs(topChildren->pdgId()) == 5) { // b
Tomas Dado's avatar
Tomas Dado committed
598
          hasB = true;
Tomas Dado's avatar
Tomas Dado committed
599
          topChildren = PartonHistoryUtils::findAfterFSR(topChildren);// b After FSR
Tomas Dado's avatar
Tomas Dado committed
600
601
602
603
604
605
606
607
608
609
610
611
          b_p4 = topChildren->p4();
          // In MG5 generation of ttgamma it is not expected to have any b radiation 'recorded'
          for (size_t b = 0; b < topChildren->nChildren(); b++) {// b Children
            const xAOD::TruthParticle* bChildren = topChildren->child(b);
            if (bChildren && bChildren->pdgId() == 22) {
              // JUST FOR EXTRA SAFETY (not necessary)
              if (has_ph) {
                if (bChildren->p4().Pt() > Ph_p4.Pt()) {
                  ph_b = true;
                  ph_ISR = false;
                  ph_Top = false;
                  ph_W = false;
Tomas Dado's avatar
Tomas Dado committed
612
                  bChildren = PartonHistoryUtils::findAfterFSR(bChildren);
Tomas Dado's avatar
Tomas Dado committed
613
614
615
616
617
                  Ph_p4 = bChildren->p4();
                }
              } else {
                has_ph = true;
                ph_b = true;
Tomas Dado's avatar
Tomas Dado committed
618
                bChildren = PartonHistoryUtils::findAfterFSR(bChildren);
Tomas Dado's avatar
Tomas Dado committed
619
620
621
622
                Ph_p4 = bChildren->p4();
              }
            }
          }
623
        } else if (std::abs(topChildren->pdgId()) == 22) {
Tomas Dado's avatar
Tomas Dado committed
624
625
626
          // JUST FOR EXTRA SAFETY (not necessary)
          if (has_ph) {
            if (topChildren->p4().Pt() > Ph_p4.Pt()) {
Tomas Dado's avatar
Tomas Dado committed
627
              topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
628
629
630
631
              Ph_p4 = topChildren->p4();
              ph_Top = true;
            }
          } else {
Tomas Dado's avatar
Tomas Dado committed
632
            topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
633
634
635
636
637
638
639
640
641
            Ph_p4 = topChildren->p4();
            has_ph = true;
            ph_Top = true;
            ph_W = false;
            ph_ISR = false;
            ph_b = false;
          }
        }
        // sometimes the W is not recorded and the W products are recorded as top products
642
        else if (std::abs(topChildren->pdgId()) <= 4 || (std::abs(topChildren->pdgId()) > 10 && std::abs(topChildren->pdgId()) < 17)) {
Tomas Dado's avatar
Tomas Dado committed
643
644
645
646
647
          hasW = true;
          hasAbsentW = true;
          if (abs(topChildren->pdgId()) < 7) {
            BranchType = 50;
          }// hadronic
648
          else if (std::abs(topChildren->pdgId()) > 10 && std::abs(topChildren->pdgId()) < 17) {
Tomas Dado's avatar
Tomas Dado committed
649
650
651
            BranchType = 10;
          }// leptonic
          if (topChildren->pdgId() > 0) {
Tomas Dado's avatar
Tomas Dado committed
652
            topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
653
654
655
656
            Wdecay1_p4 = topChildren->p4();
            Wdecay1_pdgId = topChildren->pdgId();
            hasWdecayProd1 = true;
          } else {
Tomas Dado's avatar
Tomas Dado committed
657
            topChildren = PartonHistoryUtils::findAfterFSR(topChildren);
Tomas Dado's avatar
Tomas Dado committed
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
            Wdecay2_p4 = topChildren->p4();
            Wdecay2_pdgId = topChildren->pdgId();
            hasWdecayProd2 = true;
          }//else
          W_p4 = W_p4 + topChildren->p4();
        }// if top children
      } // for top children

      // BranchType Determination if there is a photon
      if (hasAbsentW && (ph_Top || ph_W)) {
        BranchType = -1;
      }// if the W is not recorded and still the photon is from the top, the source of the photon is then ambiguous
       // among top and W. BranchType would be +1. Category would be 0.
      if (has_ph && ph_Top) {
        BranchType = BranchType + 2;
      }
      if (has_ph && ph_W) {
        BranchType = BranchType + 4;
      }
      if (has_ph && ph_ISR) {
        BranchType = BranchType + 5;
      }
      if (has_ph && ph_b) {
        BranchType = BranchType + 8;
      }

      if (hasT && hasW && hasB && hasWdecayProd1 && hasWdecayProd2 && BranchType != -1) return true;
    }// particle
686

Tomas Dado's avatar
Tomas Dado committed
687
688
    return false;
  }
689

Tomas Dado's avatar
Tomas Dado committed
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
  StatusCode CalcTopPartonHistory::execute() {
    // Get the Truth Particles
    const xAOD::TruthParticleContainer* truthParticles(nullptr);

    ATH_CHECK(evtStore()->retrieve(truthParticles, m_config->sgKeyMCParticle()));

    // Create the partonHistory xAOD object
    xAOD::PartonHistoryAuxContainer* partonAuxCont = new xAOD::PartonHistoryAuxContainer {};
    xAOD::PartonHistoryContainer* partonCont = new xAOD::PartonHistoryContainer {};
    partonCont->setStore(partonAuxCont);

    xAOD::PartonHistory* partonHistory = new xAOD::PartonHistory {};
    partonCont->push_back(partonHistory);

    // Save to StoreGate / TStore
    std::string outputSGKey = m_config->sgKeyTopPartonHistory();
    std::string outputSGKeyAux = outputSGKey + "Aux.";

    xAOD::TReturnCode save = evtStore()->tds()->record(partonCont, outputSGKey);
    xAOD::TReturnCode saveAux = evtStore()->tds()->record(partonAuxCont, outputSGKeyAux);
    if (!save || !saveAux) {
      return StatusCode::FAILURE;
    }

    return StatusCode::SUCCESS;
715
  }
Tomas Dado's avatar
Tomas Dado committed
716
717
718
719
720
721
722

  void CalcTopPartonHistory::fillEtaBranch(xAOD::PartonHistory* partonHistory, std::string branchName,
                                           TLorentzVector& tlv) {
    if (tlv.CosTheta() == 1.) partonHistory->auxdecor< float >(branchName) = 1000.;
    else if (tlv.CosTheta() == -1.) partonHistory->auxdecor< float >(branchName) = 1000.;
    else partonHistory->auxdecor< float >(branchName) = tlv.Eta();
    return;
723
724
  }
}