PhotosInterface.cc

Go to the documentation of this file.

#include <iostream>
#include "GeneratorInterface/PhotosInterface/interface/PhotosInterface.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"
 
#include "FWCore/ServiceRegistry/interface/Service.h"
 
#include "HepMC/GenEvent.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/HEPEVT_Wrapper.h"
 
using namespace gen;
using namespace edm;
using namespace std;
 
CLHEP::HepRandomEngine* PhotosInterface::fRandomEngine = nullptr;
 
extern "C" {
 
void phoini_(void);
void photos_(int&);
 
double phoran_(int* idummy) { return PhotosInterface::flat(); }
 
extern struct {
  // bool qedrad[NMXHEP];
  bool qedrad[4000];  // hardcoded for now...
} phoqed_;
}
 
PhotosInterface::PhotosInterface() : fOnlyPDG(-1) {
  fSpecialSettings.push_back("QED-brem-off:all");
  fAvoidTauLeptonicDecays = false;
  fIsInitialized = false;
}
 
PhotosInterface::PhotosInterface(const edm::ParameterSet&) : fOnlyPDG(-1) {
  fSpecialSettings.push_back("QED-brem-off:all");
  fIsInitialized = false;
}
 
PhotosInterface::~PhotosInterface() {}
 
void PhotosInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine) { fRandomEngine = decayRandomEngine; }
 
void PhotosInterface::configureOnlyFor(int ipdg) {
  fOnlyPDG = ipdg;
  //   std::ostringstream command;
  //   command << "QED-brem-off:" << fOnlyPDG ;
  fSpecialSettings.clear();
  //   fSpecialSettings.push_back( command.str() );
 
  return;
}
 
void PhotosInterface::init() {
  if (fIsInitialized)
    return;  // do init only once
 
  phoini_();
 
  fIsInitialized = true;
 
  return;
}
 
HepMC::GenEvent* PhotosInterface::apply(HepMC::GenEvent* evt) {
  if (!fIsInitialized)
    return evt;  // conv.read_next_event();
 
  // loop over HepMC::GenEvent, find vertices
 
  // for ( int ip=0; ip<evt->particles_size(); ip++ )
  for (int ip = 0; ip < 4000; ip++)  // 4000 is the max size of the array
  {
    phoqed_.qedrad[ip] = true;
  }
 
  //
  // now do actual job
  //
 
  for (int iv = 1; iv <= evt->vertices_size(); iv++) {
    bool legalVtx = false;
 
    fSecVtxStore.clear();
 
    HepMC::GenVertex* vtx = evt->barcode_to_vertex(-iv);
 
    if (vtx->particles_in_size() != 1)
      continue;  // more complex than we need
    if (vtx->particles_out_size() <= 1)
      continue;  // no outcoming particles
 
    if ((*(vtx->particles_in_const_begin()))->pdg_id() == 111)
      continue;  // pi0 decay vtx - no point to try
 
    if (fOnlyPDG != 1 && (*(vtx->particles_in_const_begin()))->pdg_id() != fOnlyPDG) {
      continue;
    } else {
      // requested for specific PDG ID only, typically tau (15)
      //
      // first check if a brem vertex, where outcoming are the same pdg id and a photon
      //
      bool same = false;
      for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children);
           pitr != vtx->particles_end(HepMC::children);
           ++pitr) {
        if ((*pitr)->pdg_id() == fOnlyPDG) {
          same = true;
          break;
        }
      }
      if (same)
        continue;
 
      // OK, we get here if incoming fOnlyPDG and something else outcoming
      // call it for the whole branch starting at vtx barcode iv, and go on
      // NOTE: theoretically, it has a danger of double counting in vertices
      // down the decay branch originating from fOnlyPDG, but in practice
      // it's unlikely that down the branchg there'll be more fOnlyPDG's
 
      // cross-check printout
      // vtx->print();
 
      // overprotection...
      //
      if (fOnlyPDG == 15 && fAvoidTauLeptonicDecays && isTauLeptonicDecay(vtx))
        continue;
 
      applyToBranch(evt, -iv);
      continue;
    }
 
    // configured for all types of particles
    //
    for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children);
         pitr != vtx->particles_end(HepMC::children);
         ++pitr) {
      // quark or gluon out of this vertex - no good !
      if (abs((*pitr)->pdg_id()) >= 1 && abs((*pitr)->pdg_id()) <= 8)
        break;
      if (abs((*pitr)->pdg_id()) == 21)
        break;
 
      if ((*pitr)->status() == 1 || (*pitr)->end_vertex()) {
        // OK, legal already !
        legalVtx = true;
        break;
      }
    }
 
    if (!legalVtx)
      continue;
 
    applyToVertex(evt, -iv);
 
  }  // end of master loop
 
  // restore event number in HEPEVT (safety measure, somehow needed by Hw6)
  HepMC::HEPEVT_Wrapper::set_event_number(evt->event_number());
 
  return evt;
}
 
void PhotosInterface::applyToVertex(HepMC::GenEvent* evt, int vtxbcode) {
  HepMC::GenVertex* vtx = evt->barcode_to_vertex(vtxbcode);
 
  if (fAvoidTauLeptonicDecays && isTauLeptonicDecay(vtx))
    return;
 
  // cross-check printout
  //
  // vtx->print();
 
  // first, flush out HEPEVT & tmp barcode storage
  //
  HepMC::HEPEVT_Wrapper::zero_everything();
  fBarcodes.clear();
 
  // add incoming particle
  //
  int index = 1;
  HepMC::HEPEVT_Wrapper::set_id(index, (*(vtx->particles_in_const_begin()))->pdg_id());
  HepMC::FourVector vec4;
  vec4 = (*(vtx->particles_in_const_begin()))->momentum();
  HepMC::HEPEVT_Wrapper::set_momentum(index, vec4.x(), vec4.y(), vec4.z(), vec4.e());
  HepMC::HEPEVT_Wrapper::set_mass(index, (*(vtx->particles_in_const_begin()))->generated_mass());
  HepMC::HEPEVT_Wrapper::set_position(
      index, vtx->position().x(), vtx->position().y(), vtx->position().z(), vtx->position().t());
  HepMC::HEPEVT_Wrapper::set_status(index, (*(vtx->particles_in_const_begin()))->status());
  HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);
  fBarcodes.push_back((*(vtx->particles_in_const_begin()))->barcode());
 
  // add outcoming particles (decay products)
  //
  int lastDau = 1;
  for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children);
       pitr != vtx->particles_end(HepMC::children);
       ++pitr) {
    if ((*pitr)->status() == 1 || (*pitr)->end_vertex()) {
      index++;
      vec4 = (*pitr)->momentum();
      HepMC::HEPEVT_Wrapper::set_id(index, (*pitr)->pdg_id());
      HepMC::HEPEVT_Wrapper::set_momentum(index, vec4.x(), vec4.y(), vec4.z(), vec4.e());
      HepMC::HEPEVT_Wrapper::set_mass(index, (*pitr)->generated_mass());
      vec4 = (*pitr)->production_vertex()->position();
      HepMC::HEPEVT_Wrapper::set_position(index, vec4.x(), vec4.y(), vec4.z(), vec4.t());
      HepMC::HEPEVT_Wrapper::set_status(index, (*pitr)->status());
      HepMC::HEPEVT_Wrapper::set_parents(index, 1, 1);
      fBarcodes.push_back((*pitr)->barcode());
      lastDau++;
    }
    if ((*pitr)->end_vertex()) {
      fSecVtxStore.push_back((*pitr)->end_vertex()->barcode());
    }
  }
 
  // store, further to set NHEP in HEPEVT
  //
  int nentries = index;
 
  // reset master pointer to mother
  index = 1;
  HepMC::HEPEVT_Wrapper::set_children(index, 2, lastDau);  // FIXME: need to check
                                                           // if last daughter>=2 !!!
 
  // finally, set number of entries (NHEP) in HEPEVT
  //
  HepMC::HEPEVT_Wrapper::set_number_entries(nentries);
 
  // cross-check printout HEPEVT
  //  HepMC::HEPEVT_Wrapper::print_hepevt();
 
  // OK, 1-level vertex is formed - now, call PHOTOS
  //
  photos_(index);
 
  // another cross-check printout HEPEVT - after photos
  // HepMC::HEPEVT_Wrapper::print_hepevt();
 
  // now check if something has been generated
  // and make all adjustments to underlying vtx/parts
  //
  attachParticles(evt, vtx, nentries);
 
  // ugh, done with this vertex !
 
  return;
}
 
void PhotosInterface::applyToBranch(HepMC::GenEvent* evt, int vtxbcode) {
  fSecVtxStore.clear();
 
  // 1st level vertex
  //
  applyToVertex(evt, vtxbcode);
 
  // now look down the branch for more vertices, if any
  //
  // Note: fSecVtxStore gets filled up in applyToVertex, if necessary
  //
  unsigned int vcounter = 0;
 
  while (vcounter < fSecVtxStore.size()) {
    applyToVertex(evt, fSecVtxStore[vcounter]);
    vcounter++;
  }
 
  fSecVtxStore.clear();
 
  return;
}
 
void PhotosInterface::attachParticles(HepMC::GenEvent* evt, HepMC::GenVertex* vtx, int nentries) {
  if (HepMC::HEPEVT_Wrapper::number_entries() > nentries) {
    // yes, need all corrections and adjustments -
    // figure out how many photons and what particles in
    // the decay branch have changes;
    // also, follow up each one and correct accordingly;
    // at the same time, add photon(s) to the GenVertex
    //
 
    // vtx->print();
 
    int largestBarcode = -1;
    int Nbcodes = fBarcodes.size();
 
    for (int ip = 1; ip < Nbcodes; ip++) {
      int bcode = fBarcodes[ip];
      HepMC::GenParticle* prt = evt->barcode_to_particle(bcode);
      if (bcode > largestBarcode)
        largestBarcode = bcode;
      double px = HepMC::HEPEVT_Wrapper::px(ip + 1);
      double py = HepMC::HEPEVT_Wrapper::py(ip + 1);
      double pz = HepMC::HEPEVT_Wrapper::pz(ip + 1);
      double e = HepMC::HEPEVT_Wrapper::e(ip + 1);
      double m = HepMC::HEPEVT_Wrapper::m(ip + 1);
 
      if (prt->end_vertex()) {
        HepMC::GenVertex* endVtx = prt->end_vertex();
 
        std::vector<int> secVtxStorage;
        secVtxStorage.clear();
 
        secVtxStorage.push_back(endVtx->barcode());
 
        const HepMC::FourVector& mom4 = prt->momentum();
 
        // now rescale all descendants
        double bet1[3], bet2[3], gam1, gam2, pb;
        double mass = mom4.m();
        bet1[0] = -(mom4.px() / mass);
        bet1[1] = -(mom4.py() / mass);
        bet1[2] = -(mom4.pz() / mass);
        bet2[0] = px / m;
        bet2[1] = py / m;
        bet2[2] = pz / m;
        gam1 = mom4.e() / mass;
        gam2 = e / m;
 
        unsigned int vcounter = 0;
 
        while (vcounter < secVtxStorage.size()) {
          HepMC::GenVertex* theVtx = evt->barcode_to_vertex(secVtxStorage[vcounter]);
 
          for (HepMC::GenVertex::particle_iterator pitr = theVtx->particles_begin(HepMC::children);
               pitr != theVtx->particles_end(HepMC::children);
               ++pitr) {
            if ((*pitr)->end_vertex()) {
              secVtxStorage.push_back((*pitr)->end_vertex()->barcode());
            }
 
            if (theVtx->particles_out_size() == 1 && (*pitr)->pdg_id() == prt->pdg_id()) {
              // carbon copy
              (*pitr)->set_momentum(HepMC::FourVector(px, py, pz, e));
              continue;
            }
 
            HepMC::FourVector dmom4 = (*pitr)->momentum();
 
            // Boost vector to parent rest frame...
            pb = bet1[0] * dmom4.px() + bet1[1] * dmom4.py() + bet1[2] * dmom4.pz();
            double dpx = dmom4.px() + bet1[0] * (dmom4.e() + pb / (gam1 + 1.));
            double dpy = dmom4.py() + bet1[1] * (dmom4.e() + pb / (gam1 + 1.));
            double dpz = dmom4.pz() + bet1[2] * (dmom4.e() + pb / (gam1 + 1.));
            double de = gam1 * dmom4.e() + pb;
            // ...and boost back to modified parent frame
            pb = bet2[0] * dpx + bet2[1] * dpy + bet2[2] * dpz;
            dpx += bet2[0] * (de + pb / (gam2 + 1.));
            dpy += bet2[1] * (de + pb / (gam2 + 1.));
            dpz += bet2[2] * (de + pb / (gam2 + 1.));
            de *= gam2;
            de += pb;
 
            (*pitr)->set_momentum(HepMC::FourVector(dpx, dpy, dpz, de));
          }
          vcounter++;
        }
 
        secVtxStorage.clear();
      }
 
      prt->set_momentum(HepMC::FourVector(px, py, pz, e));
 
    }  // ok, all affected particles update, but the photon(s) still not inserted
 
    int newlyGen = HepMC::HEPEVT_Wrapper::number_entries() - nentries;
 
    if (largestBarcode < evt->particles_size()) {
      // need to adjust barcodes down from the affected vertex/particles
      // such that we "free up" barcodes for newly generated photons
      // in the middle of the event record
      //
      for (int ipp = evt->particles_size(); ipp > largestBarcode; ipp--) {
        (evt->barcode_to_particle(ipp))->suggest_barcode(ipp + newlyGen);
      }
    }
 
    // now attach new generated photons to the vertex
    //
    for (int ipnw = 1; ipnw <= newlyGen; ipnw++) {
      int nbcode = largestBarcode + ipnw;
      int pdg_id = HepMC::HEPEVT_Wrapper::id(nentries + ipnw);
      int status = HepMC::HEPEVT_Wrapper::status(nentries + ipnw);
      double px = HepMC::HEPEVT_Wrapper::px(nentries + ipnw);
      double py = HepMC::HEPEVT_Wrapper::py(nentries + ipnw);
      double pz = HepMC::HEPEVT_Wrapper::pz(nentries + ipnw);
      double e = HepMC::HEPEVT_Wrapper::e(nentries + ipnw);
      double m = HepMC::HEPEVT_Wrapper::m(nentries + ipnw);
 
      HepMC::GenParticle* NewPart = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, e), pdg_id, status);
      NewPart->set_generated_mass(m);
      NewPart->suggest_barcode(nbcode);
      vtx->add_particle_out(NewPart);
    }
 
    //vtx->print();
    //std::cout << " leaving attachParticles() " << std::endl;
 
  }  // end of global if-statement
 
  return;
}
 
bool PhotosInterface::isTauLeptonicDecay(HepMC::GenVertex* vtx) {
  if (abs((*(vtx->particles_in_const_begin()))->pdg_id()) != 15)
    return false;
 
  for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children);
       pitr != vtx->particles_end(HepMC::children);
       ++pitr) {
    if (abs((*pitr)->pdg_id()) == 11 || abs((*pitr)->pdg_id()) == 13) {
      return true;
    }
  }
 
  return false;
}
 
double PhotosInterface::flat() {
  if (!fRandomEngine) {
    throw cms::Exception("LogicError")
        << "PhotosInterface::flat: Attempt to generate random number when engine pointer is null\n"
        << "This might mean that the code was modified to generate a random number outside the\n"
        << "event and beginLuminosityBlock methods, which is not allowed.\n";
  }
  return fRandomEngine->flat();
}
 
DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosInterface, "Photos2155");