plugins/TauSpinner/read_particles_from_HepMC.cc

0001 #include "GeneratorInterface/TauolaInterface/interface/read_particles_from_HepMC.h"
0002 using namespace TauSpinner;
0003 /*******************************************************************************
0004   Get daughters of HepMC::GenParticle
0005
0006   Recursively searches for final-state daughters of 'x'
0007 *******************************************************************************/
0008 inline std::vector<SimpleParticle> *getDaughters(HepMC::GenParticle *x) {
0009   std::vector<SimpleParticle> *daughters = new std::vector<SimpleParticle>();
0010   if (!x->end_vertex())
0011     return daughters;
0012
0013   // Check decay products of 'x'
0014   for (HepMC::GenVertex::particles_out_const_iterator p = x->end_vertex()->particles_out_const_begin();
0015        p != x->end_vertex()->particles_out_const_end();
0016        ++p) {
0017     HepMC::GenParticle *pp = *p;
0018     HepMC::FourVector mm = pp->momentum();
0019
0020     // If the daughter of 'x' has its end vertex - recursively read
0021     // all of its daughters.
0022     if (pp->end_vertex() && pp->pdg_id() != 111) {
0023       std::vector<SimpleParticle> *sub_daughters = getDaughters(pp);
0024       daughters->insert(daughters->end(), sub_daughters->begin(), sub_daughters->end());
0025
0026       delete sub_daughters;
0027     }
0028     // Otherwise - add this particle to the list of daughters.
0029     else if (pp->pdg_id() != x->pdg_id()) {
0030       SimpleParticle tp(mm.px(), mm.py(), mm.pz(), mm.e(), pp->pdg_id());
0031       daughters->push_back(tp);
0032     }
0033   }
0034
0035   return daughters;
0036 }
0037
0038 /*******************************************************************************
0039   Find last self
0040
0041   Recursively finds the last particle with the same PDG ID
0042   on the list of its decay products
0043 *******************************************************************************/
0044 inline HepMC::GenParticle *findLastSelf(HepMC::GenParticle *x) {
0045   if (!x->end_vertex())
0046     return x;
0047
0048   for (HepMC::GenVertex::particle_iterator p = x->end_vertex()->particles_begin(HepMC::children);
0049        p != x->end_vertex()->particles_end(HepMC::children);
0050        ++p) {
0051     if ((*p)->pdg_id() == x->pdg_id())
0052       return findLastSelf(*p);
0053   }
0054
0055   return x;
0056 }
0057
0058 inline bool isFirst(HepMC::GenParticle *x) {
0059   for (HepMC::GenVertex::particle_iterator p = x->production_vertex()->particles_begin(HepMC::parents);
0060        p != x->production_vertex()->particles_end(HepMC::parents);
0061        ++p) {
0062     if (x->pdg_id() == (*p)->pdg_id())
0063       return false;
0064   }
0065   return true;
0066 }
0067
0068 /*******************************************************************************
0069   Read HepMC::GenEvent.
0070
0071   Read HepMC event from data file
0072   and return particles filtered out from the event.
0073
0074   This routine is prepared for use with files generated by Pythia8.
0075   Fills:
0076
0077   'X'              - Heavy particle (W+/-, H+/-, H, Z)
0078   'tau'            - first tau
0079   'tau2'           - second tau or nu_tau, if 'X' decays to one tau only
0080   'tau_daughters'  - daughters of 'tau'
0081   'tau2_daughters' - daughters of 'tau2' or empty list, if 'tau2' is nu_tau.
0082
0083   Returns:
0084   0 - no more events to read               (finished processing the file)
0085   1 - no decay found in the event          (finished processing current event)
0086   2 - decay found and processed correctly.
0087       Event will continue to be processed
0088       with next function call.
0089 *******************************************************************************/
0090 int readParticlesFromHepMC(const HepMC::GenEvent *event,
0091                            SimpleParticle &X,
0092                            SimpleParticle &tau,
0093                            SimpleParticle &tau2,
0094                            std::vector<SimpleParticle> &tau_daughters,
0095                            std::vector<SimpleParticle> &tau2_daughters) {
0096   if (event == nullptr)
0097     return 1;
0098
0099   // Exctract particles from event
0100   HepMC::GenParticle *hX = nullptr, *hTau = nullptr, *hTau2 = nullptr;
0101
0102   for (HepMC::GenEvent::particle_const_iterator it = event->particles_begin(); it != event->particles_end(); ++it) {
0103     int pdgid = (*it)->pdg_id();
0104     if ((abs(pdgid) == 37 || pdgid == 25 || abs(pdgid) == 24 || pdgid == 23) && (*it)->end_vertex() &&
0105         (*it)->end_vertex()->particles_out_size() >= 2) {
0106       hX = *it;
0107       HepMC::GenParticle *LhX = hX;
0108       if (!isFirst(hX))
0109         continue;
0110       findLastSelf(LhX);
0111       hTau = hTau2 = nullptr;
0112
0113       for (HepMC::GenVertex::particle_iterator it2 = LhX->end_vertex()->particles_begin(HepMC::children);
0114            it2 != LhX->end_vertex()->particles_end(HepMC::children);
0115            ++it2) {
0116         if (abs((*it2)->pdg_id()) == 15 && (*it2)->end_vertex()) {
0117           if (!hTau)
0118             hTau = *it2;
0119           else if (!hTau2)
0120             hTau2 = *it2;
0121           else {
0122             std::cout << "TauSpiner: three taus in one decay" << std::endl;
0123             return 1;
0124           }
0125         }
0126         if (abs((*it2)->pdg_id()) == 16) {
0127           if (!hTau2)
0128             hTau2 = *it2;
0129           else {
0130             std::cout << "TauSpiner: two neutrinos or two taus and neutrino in one decay" << std::endl;
0131             return 1;
0132           }
0133         }
0134       }
0135       if (hTau && hTau2)
0136         break;
0137     }
0138   }
0139
0140   if (!hX)
0141     return 1;
0142
0143   if (!hTau || !hTau2) {
0144     //std::cout<<"TauSpiner: boson found but no proper tau pair or tau + neutrino found."<<std::endl;
0145     return 1;
0146   }
0147
0148   // Check for self-decaying taus
0149   hTau = findLastSelf(hTau);
0150   hTau2 = findLastSelf(hTau2);
0151
0152   // Fill SimpleParticles from HepMC particles
0153   X.setPx(hX->momentum().px());
0154   X.setPy(hX->momentum().py());
0155   X.setPz(hX->momentum().pz());
0156   X.setE(hX->momentum().e());
0157   X.setPdgid(hX->pdg_id());
0158
0159   tau.setPx(hTau->momentum().px());
0160   tau.setPy(hTau->momentum().py());
0161   tau.setPz(hTau->momentum().pz());
0162   tau.setE(hTau->momentum().e());
0163   tau.setPdgid(hTau->pdg_id());
0164
0165   tau2.setPx(hTau2->momentum().px());
0166   tau2.setPy(hTau2->momentum().py());
0167   tau2.setPz(hTau2->momentum().pz());
0168   tau2.setE(hTau2->momentum().e());
0169   tau2.setPdgid(hTau2->pdg_id());
0170
0171   // Create list of tau daughters
0172   std::vector<SimpleParticle> *buf = getDaughters(hTau);
0173   tau_daughters.clear();
0174   tau_daughters.insert(tau_daughters.end(), buf->begin(), buf->end());
0175
0176   delete buf;
0177
0178   // Second particle can be 2nd tau. In that case - read its daughters.
0179   // Otherwise it is nu_tau~
0180   buf = getDaughters(hTau2);
0181   tau2_daughters.clear();
0182   tau2_daughters.insert(tau2_daughters.end(), buf->begin(), buf->end());
0183
0184   delete buf;
0185
0186   return 0;
0187 }