ParticleReplacerClass.cc

Go to the documentation of this file.

#include "TauAnalysis/MCEmbeddingTools/interface/ParticleReplacerClass.h"

#include "GeneratorInterface/ExternalDecays/interface/DecayRandomEngine.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "HepMC/IO_HEPEVT.h"

#ifndef TXGIVE
#define TXGIVE txgive_
extern "C" {
  void TXGIVE(const char*,int length);
}
#endif

#ifndef TXGIVE_INIT
#define TXGIVE_INIT txgive_init_
extern "C" {
  void TXGIVE_INIT();
}
#endif

ParticleReplacerClass::ParticleReplacerClass(const edm::ParameterSet& pset, bool verbose):
  ParticleReplacerBase(pset),
  generatorMode_(pset.getParameter<std::string>("generatorMode")),
  tauola_(gen::TauolaInterface::getInstance()),
  printEvent_(verbose),
  outTree(0),
  maxNumberOfAttempts_(pset.getUntrackedParameter<int>("maxNumberOfAttempts", 1000))
{
    tauola_->setPSet(pset.getParameter< edm::ParameterSet>("TauolaOptions"));
//  using namespace reco;
    using namespace edm;
    using namespace std;

    //HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
    //HepMC::HEPEVT_Wrapper::set_sizeof_real(8);

    // transformationMode =
    //  0 - no transformation
    //  1 - mumu -> tautau
    transformationMode_ = pset.getUntrackedParameter<int>("transformationMode",1);
    switch (transformationMode_)
    {
        case 0:
        {   
            LogInfo("Replacer") << "won't do any transformation with the given mumu";
            break;
        }
        case 1:
        {
            LogInfo("Replacer") << "will transform mumu into tautau";
            break;
        }
        case 2:
    {
      LogInfo("Replacer") << "will transform mumu into taunu (as coming from a W boson)";
      break;
    }    
        case 3:
        {
            LogInfo("Replacer") << "Will transform  mu-nu into tau-nu. No mass correction will be made.";
            break;
        }
        default:
        {
            throw cms::Exception("ParticleReplacerClass")  << "Unknown transformation mode!\n";
            break;
        }
            
    }
    
    // If one wants to use two instances of this module in one
    // configuration file, there might occur some segmentation
    // faults due to the second initialisation of Tauola. This
    // can be prevented by setting noInitialisation to false.
    //          Caution: This option is not tested!
    noInitialisation_ = pset.getUntrackedParameter<bool>("noInitialisation",false);

    motherParticleID_ = pset.getUntrackedParameter<int>("motherParticleID",23);

    // requires the visible decay products of a tau to have a sum transverse momentum
    std::string minVisibleTransverseMomentumLine = pset.getUntrackedParameter<std::string>("minVisibleTransverseMomentum","");

  // fallback for backwards compatibility: If it's a single number then use this as a threshold for both particles
  const char* startptr = minVisibleTransverseMomentumLine.c_str();
  char* endptr;
  double d = strtod(startptr, &endptr);
  if(*endptr == '\0' && endptr != startptr)
  {
        MinVisPtCut cuts[2];
        cuts[0].type_ = cuts[1].type_ = MinVisPtCut::TAU;
    cuts[0].pt_ = cuts[1].pt_ = d;
        cuts[0].index_ = 0; cuts[1].index_ = 1;
        minVisPtCuts_.push_back(std::vector<MinVisPtCut>(cuts, cuts+2));
  }
  else
  {
      // string has new format: parse the minvistransversemomentum string
        for(std::string::size_type prev = 0, pos = 0; prev < minVisibleTransverseMomentumLine.length(); prev = pos + 1)
        {
            pos = minVisibleTransverseMomentumLine.find(';', prev);
            if(pos == std::string::npos) pos = minVisibleTransverseMomentumLine.length();

            std::string sub = minVisibleTransverseMomentumLine.substr(prev, pos - prev);
            std::vector<MinVisPtCut> cuts;
            const char* sub_c = sub.c_str();
            while(*sub_c != '\0')
            {
                const char* sep = std::strchr(sub_c, '_');
                if(sep == NULL) throw cms::Exception("Configuration") << "Minimum transverse parameter string must contain an underscore to separate type from pt threshold" << std::endl;
                std::string type(sub_c, sep);

                MinVisPtCut cut;
                if(type == "elec1") { cut.type_ = MinVisPtCut::ELEC; cut.index_ = 0; }
                else if(type == "mu1") { cut.type_ = MinVisPtCut::MU; cut.index_ = 0; }
                else if(type == "had1") { cut.type_ = MinVisPtCut::HAD; cut.index_ = 0; }
                else if(type == "tau1") { cut.type_ = MinVisPtCut::TAU; cut.index_ = 0; }
                else if(type == "elec2") { cut.type_ = MinVisPtCut::ELEC; cut.index_ = 1; }
                else if(type == "mu2") { cut.type_ = MinVisPtCut::MU; cut.index_ = 1; }
                else if(type == "had2") { cut.type_ = MinVisPtCut::HAD; cut.index_ = 1; }
                else if(type == "tau2") { cut.type_ = MinVisPtCut::TAU; cut.index_ = 1; }
                else throw cms::Exception("Configuration") << "'" << type << "' is not a valid type. Allowed values are elec1,mu1,had1,tau1,elec2,mu2,had2,tau2" << std::endl;

                char* endptr;
                cut.pt_ = strtod(sep + 1, &endptr);
                if(endptr == sep + 1) throw cms::Exception("Configuration") << "No pt threshold given" << std::endl;

                cuts.push_back(cut);
                sub_c = endptr;
            }
        minVisPtCuts_.push_back(cuts);
        }
  }

    edm::Service<TFileService> fileService_;
        if(fileService_.isAvailable()) {
          outTree = fileService_->make<TTree>( "event_generation","This tree stores information about the event generation");
          outTree->Branch("attempts",&attempts,"attempts/I");
        }

        edm::Service<RandomNumberGenerator> rng;
        if(!rng.isAvailable()) {
          throw cms::Exception("Configuration")
            << "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
            "which appears to be absent.  Please add that service to your configuration\n"
            "or remove the modules that require it." << std::endl;
        } 
        // this is a global variable defined in GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc
        decayRandomEngine = &rng->getEngine();

    edm::LogInfo("Replacer") << "generatorMode = "<< generatorMode_<< "\n";

    return;
}

ParticleReplacerClass::~ParticleReplacerClass()
{
    // do anything here that needs to be done at desctruction time
    // (e.g. close files, deallocate resources etc.)
}

// ------------ method called to produce the data  ------------
std::auto_ptr<HepMC::GenEvent> ParticleReplacerClass::produce(const reco::MuonCollection& muons, const reco::Vertex *pvtx, const HepMC::GenEvent *genEvt)
{
    using namespace edm;
    using namespace std;
    using namespace HepMC;

        if(pvtx != 0)
          throw cms::Exception("Configuration") << "ParticleReplacerClass does NOT support using primary vertex as the origin for taus" << std::endl;

    HepMC::GenEvent * evt=0;

    GenVertex * zvtx = new GenVertex();

    reco::GenParticle * part1=0;
    reco::GenParticle * part2=0;

    std::vector<reco::Particle> particles;  
    switch (transformationMode_)
    {
        case 0: // mumu->mumu
        {
            if (muons.size()!=2)
            {
                LogError("Replacer") << "the decay mode Z->mumu requires exactly two muons, aborting processing";
                return std::auto_ptr<HepMC::GenEvent>(0);
            }
    
            targetParticleMass_  = 0.105658369;
            targetParticlePdgID_ = 13;
                
            reco::Muon muon1 = muons.at(0);
            reco::Muon muon2 = muons.at(1);
            reco::Particle tau1(muon1.charge(), muon1.p4(), muon1.vertex(), muon1.pdgId(), 0, true);
            reco::Particle tau2(muon2.charge(), muon2.p4(), muon2.vertex(), muon2.pdgId(), 0, true);
            particles.push_back(tau1);
            particles.push_back(tau2);
            break;
        } 
        case 1: // mumu->tautau
        {
            if (muons.size()!=2)
            {
                LogError("Replacer") << "the decay mode Z->tautau requires exactly two muons, aborting processing";
                return std::auto_ptr<HepMC::GenEvent>(0);
            }

            targetParticleMass_  = 1.77690;
            targetParticlePdgID_ = 15;
            
            reco::Muon muon1 = muons.at(0);
            reco::Muon muon2 = muons.at(1);
            reco::Particle tau1(muon1.charge(), muon1.p4(), muon1.vertex(), muon1.pdgId(), 0, true);
            reco::Particle tau2(muon2.charge(), muon2.p4(), muon2.vertex(), muon2.pdgId(), 0, true);
            transformMuMu2TauTau(&tau1, &tau2);
            particles.push_back(tau1);
            particles.push_back(tau2);          
            break;
        }
    case 2: // mumu->taunu (W boson)
    {
      if (muons.size()!=2)
      {
        LogError("Replacer") << "the decay mode Z->tautau requires exactly two muons, aborting processing";
        return std::auto_ptr<HepMC::GenEvent>(0);
      }

      targetParticleMass_  = 1.77690;
      targetParticlePdgID_ = 15;
      
      reco::Muon muon1 = muons.at(0);
      reco::Muon muon2 = muons.at(1);
      reco::Particle tau1(muon1.charge(), muon1.p4(), muon1.vertex(), muon1.pdgId(), 0, true);
      reco::Particle tau2(muon2.charge(), muon2.p4(), muon2.vertex(), muon2.pdgId(), 0, true);
      transformMuMu2TauNu(&tau1, &tau2);
      particles.push_back(tau1);
      particles.push_back(tau2);                      
      break;
    }  
    case 3: // mu-nu->tau-nu
    {
      if (muons.size()!=2)
      {
        LogError("Replacer") << "transformation mode mu-nu ->tau-nu - wrong input";
        return std::auto_ptr<HepMC::GenEvent>(0);
      }

      targetParticleMass_  = 1.77690;
      targetParticlePdgID_ = 15;
      int targetParticlePdgIDNu_ = 16;
      
      reco::Muon muon1 = muons.at(0);
      reco::Muon::LorentzVector l(muon1.px(), muon1.py(), muon1.pz(), 
                                sqrt(
                                muon1.px()*muon1.px()+
                                muon1.py()*muon1.py()+
                                muon1.pz()*muon1.pz()+targetParticleMass_*targetParticleMass_));

      reco::Particle tau1(muon1.charge(), l, muon1.vertex(), targetParticlePdgID_*std::abs(muon1.pdgId())/muon1.pdgId() 
                                , 0, true
                         );
      tau1.setStatus(1);
      particles.push_back(tau1);

      reco::Muon nu    = muons.at(1);
      reco::Particle nutau( 0, nu.p4(), nu.vertex(), -targetParticlePdgIDNu_*std::abs(muon1.pdgId())/muon1.pdgId(), 0, true);
      nutau.setStatus(1);
      particles.push_back(nutau);
 
      break;
    }  

    }
    
    if (particles.size()==0)
    {
        LogError("Replacer") << "the creation of the new particles failed somehow"; 
        return std::auto_ptr<HepMC::GenEvent>(0);
    }
    else
    {
        LogInfo("Replacer") << particles.size() << " particles found, continue processing";
    }

    if (genEvt)
    {
    
            evt = new HepMC::GenEvent(*genEvt);
    
            for ( GenEvent::vertex_iterator p = evt->vertices_begin(); p != evt->vertices_end(); p++ ) 
            {
                GenVertex * vtx=(*p);
                
                if ( vtx->particles_out_size()<=0 || vtx->particles_in_size()<=0)
                    continue;
                
                bool temp_muon1=false,temp_muon2=false,temp_z=false;
                for (GenVertex::particles_out_const_iterator it = vtx->particles_out_const_begin(); it!=vtx->particles_out_const_end(); it++)
                {
                    if ((*it)->pdg_id()==13) temp_muon1=true;
                    if ((*it)->pdg_id()==-13) temp_muon2=true;
                    if ((*it)->pdg_id()==23) temp_z=true;
                }
    
                int mother_pdg=(*vtx->particles_in_const_begin())->pdg_id();
    
                if ((vtx->particles_out_size()==2 && vtx->particles_in_size()>0 
                    && mother_pdg == 23  
                    && temp_muon1
                    && temp_muon2)
                || (vtx->particles_out_size()>2 && vtx->particles_in_size()>0 
                    && mother_pdg == 23
                    && temp_muon1 && temp_muon2 && temp_z ))
                {
                    zvtx=*p;
                }
            }

            cleanEvent(evt, zvtx);

            // prevent a decay of existing particles
            // this is due to a bug in the PythiaInterface that should be fixed in newer versions
            for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
                (*it)->set_status(0);

            for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
            {
                zvtx->add_particle_out(new HepMC::GenParticle((FourVector)it->p4(), it->pdgId(), 1, Flow(), Polarization(0,0)));
            }
    }
    // new product with tau decays
    else
    {
        reco::Particle::LorentzVector mother_particle_p4;
        for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
            mother_particle_p4+=it->p4();

        reco::Particle::Point production_point = particles.begin()->vertex();

                GenVertex* startVtx = new GenVertex(FourVector(production_point.x()*10,production_point.y()*10,production_point.z()*10,0));
                startVtx->add_particle_in( new GenParticle( FourVector(0,0,7000,7000), 2212, 3 ) );
                startVtx->add_particle_in( new GenParticle( FourVector(0,0,-7000,7000), 2212, 3 ) );

                GenVertex * decayvtx = new GenVertex(FourVector(production_point.x()*10,production_point.y()*10,production_point.z()*10,0));
        HepMC::GenParticle * mother_particle = new HepMC::GenParticle((FourVector)mother_particle_p4, motherParticleID_, (generatorMode_=="Pythia" ? 3 : 2), Flow(), Polarization(0,0));
                if (transformationMode_ == 3) {
                  //std::cout << "Overriding mother particle id\n" << std::endl;
                  int muPDG = particles.begin()->pdgId();
                  int id = -24*muPDG/std::abs(muPDG);
                  mother_particle->set_pdg_id(id);
                }

                startVtx->add_particle_out(mother_particle);
                decayvtx->add_particle_in(mother_particle);
        evt = new HepMC::GenEvent();
        for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
        {
                        //std::cout << "XXX" << it->p4().pt() << " " << it->pdgId() << std::endl;
            decayvtx->add_particle_out(new HepMC::GenParticle((FourVector)it->p4(), it->pdgId(), 1, Flow(), Polarization(0,0)));            
        }

        evt->add_vertex(startVtx);
        evt->add_vertex(decayvtx);
    }
    repairBarcodes(evt);
    
    HepMC::GenEvent * retevt = 0;
    HepMC::GenEvent * tempevt = 0;

    int nr_of_trials=0;

    unsigned int cntVisPt_all = 0;
    unsigned int cntVisPt_pass = 0;
    
        HepMC::IO_HEPEVT conv;
    for (int i = 0; i<maxNumberOfAttempts_; i++)
    {
        ++cntVisPt_all;
        if (generatorMode_ == "Pythia") // Pythia
        {
            LogError("Replacer") << "Pythia is currently not supported!";
            return std::auto_ptr<HepMC::GenEvent>(evt);
        }

        if (generatorMode_ == "Tauola") // TAUOLA
        {
            conv.write_event(evt);
            tempevt=tauola_->decay(evt);
        }

        if (testEvent(tempevt))
        {
            if (retevt==0) {
              retevt=tempevt;
                        } else {
                          delete tempevt; 
                        }
            ++cntVisPt_pass;
        } else {
                  delete tempevt;
                }
    }

    tried = cntVisPt_all;
    passed = cntVisPt_pass;

    std::cout << /*minVisibleTransverseMomentum_ <<*/ " " << cntVisPt_pass << "\t" << cntVisPt_all << "\n";
    if (!retevt)
    {
        LogError("Replacer") << "failed to create an event which satisfies the minimum visible transverse momentum cuts ";
        attempts=-1;
                if(outTree) outTree->Fill();
        return std::auto_ptr<HepMC::GenEvent>(0);
    }
    attempts=nr_of_trials;
    if(outTree) outTree->Fill();    

    // recover the status codes
    if (genEvt)
    {
        for (GenEvent::particle_iterator it=retevt->particles_begin();it!=retevt->particles_end();it++)
        {
            if ((*it)->end_vertex())
                (*it)->set_status(2);
            else
                (*it)->set_status(1);
        }
    }

        std::auto_ptr<HepMC::GenEvent> ret(retevt);

    if (printEvent_)
        retevt->print(std::cout);

    delete part1;
    delete part2;
    delete zvtx;
    delete evt;
    return ret;
}

// ------------ method called once each job just before starting event loop  ------------
void ParticleReplacerClass::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
{
    tauola_->init(iSetup);
}

// ------------ method called once each job just after ending the event loop  ------------
void 
ParticleReplacerClass::endJob()
{
    tauola_->statistics();
}

bool ParticleReplacerClass::testEvent(HepMC::GenEvent * evt)
{
    using namespace HepMC;
        using namespace edm;
    
    if (minVisPtCuts_.empty()) //ibleTransverseMomentum_<=0)
        return true;

  std::vector<double> mus;
  std::vector<double> elecs;
  std::vector<double> hads;
  std::vector<double> taus;

  for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
    {
        if (abs((*it)->pdg_id())==15 && (*it)->end_vertex())
    {
        FourVector vis_mom();
        math::PtEtaPhiMLorentzVector visible_momentum;
        std::queue<const GenParticle *> decaying_particles;
            decaying_particles.push(*it);
            int t=0;
      enum { ELEC, MU, HAD } type = HAD;
            while(!decaying_particles.empty() && (++t < 30))
            {
                const GenParticle * front = decaying_particles.front();
            decaying_particles.pop();

                if (!front->end_vertex())
                {
                    int pdgId=abs(front->pdg_id());
                    if (pdgId>10 && pdgId!=12 && pdgId!=14 && pdgId!=16)
                        visible_momentum+=(math::PtEtaPhiMLorentzVector)front->momentum();

          if(pdgId == 11) type = ELEC;
          if(pdgId == 13) type = MU;
                }
                else
                {
                    GenVertex * temp_vert = front->end_vertex();
                    for (GenVertex::particles_out_const_iterator it2=temp_vert->particles_out_const_begin();it2!=temp_vert->particles_out_const_end();it2++)
                        decaying_particles.push((*it2));
                }
        }

      double vis_pt = visible_momentum.pt();
      taus.push_back(vis_pt);
      if(type == MU) mus.push_back(vis_pt);
      if(type == ELEC) elecs.push_back(vis_pt);
      if(type == HAD) hads.push_back(vis_pt);
        }
    }

  std::sort(taus.begin(), taus.end(), std::greater<double>());
  std::sort(elecs.begin(), elecs.end(), std::greater<double>());
  std::sort(mus.begin(), mus.end(), std::greater<double>());
  std::sort(hads.begin(), hads.end(), std::greater<double>());

  for(std::vector<std::vector<MinVisPtCut> >::const_iterator iter = minVisPtCuts_.begin(); iter != minVisPtCuts_.end(); ++iter)
  {
    std::vector<MinVisPtCut>::const_iterator iter2;
    for(iter2 = iter->begin(); iter2 != iter->end(); ++iter2)
    {
      std::vector<double>* collection;
      switch(iter2->type_)
      {
      case MinVisPtCut::ELEC: collection = &elecs; break;
      case MinVisPtCut::MU: collection = &mus; break;
      case MinVisPtCut::HAD: collection = &hads; break;
      case MinVisPtCut::TAU: collection = &taus; break;
      default: assert(false); break;
      }

      // subcut fail
      if(iter2->index_ >= collection->size() || (*collection)[iter2->index_] < iter2->pt_)
        break;
    }

    // no subcut failed: This cut passed
    if(iter2 == iter->end())
      return true;
  }

  LogInfo("Replacer") << "refusing the event as the sum of the visible transverse momenta is too small\n";
    return false;
}

void ParticleReplacerClass::cleanEvent(HepMC::GenEvent * evt, HepMC::GenVertex * vtx)
{
    using namespace HepMC;
    using namespace std;
    using namespace edm;
    using namespace reco;

    stack<HepMC::GenParticle *> deleteParticle;
    
    stack<GenVertex *> deleteVertex;
    stack<GenVertex *> queueVertex;

    if (vtx->particles_out_size()>0)
    {
        for (GenVertex::particles_out_const_iterator it=vtx->particles_out_const_begin();it!=vtx->particles_out_const_end();it++)
        {
            deleteParticle.push(*it);
            if ((*it)->end_vertex())
                queueVertex.push((*it)->end_vertex());
        }
    }

    while (!queueVertex.empty())
    {
        GenVertex * temp_vtx=queueVertex.top();
        if (temp_vtx->particles_out_size()>0)
        {
            for (GenVertex::particles_out_const_iterator it=temp_vtx->particles_out_const_begin();it!=temp_vtx->particles_out_const_end();it++)
            {
                if ((*it)->end_vertex())
                    queueVertex.push((*it)->end_vertex());
            }
            delete temp_vtx;
        }
        deleteVertex.push(queueVertex.top());
        queueVertex.pop();
    }
    
    while (!deleteVertex.empty())
    {
        evt->remove_vertex(deleteVertex.top());
        deleteVertex.pop();
    }

    while (!deleteParticle.empty())
    {
        delete vtx->remove_particle(deleteParticle.top());
        deleteParticle.pop();
    }

    while (!deleteVertex.empty())
        deleteVertex.pop();
    while (!queueVertex.empty())
        queueVertex.pop();

    repairBarcodes(evt);
}

void ParticleReplacerClass::repairBarcodes(HepMC::GenEvent * evt)
{
    using namespace HepMC;

    // repair the barcodes
    int max_barc=0;
    for (GenEvent::vertex_iterator it=evt->vertices_begin(), next;it!=evt->vertices_end();it=next)
    {
        next=it;++next;
        while (!(*it)->suggest_barcode(-1*(++max_barc)))
            ;
    }

    max_barc=0;
    for (GenEvent::particle_iterator it=evt->particles_begin(), next;it!=evt->particles_end();it=next)
    {
        next=it;++next;
        while (!(*it)->suggest_barcode(++max_barc))
            ;
    }
}

void ParticleReplacerClass::transformMuMu2TauTau(reco::Particle * muon1, reco::Particle * muon2)
{
    using namespace edm;
    using namespace reco;
    using namespace std;
    
    reco::Particle::LorentzVector muon1_momentum = muon1->p4();
    reco::Particle::LorentzVector muon2_momentum =  muon2->p4();
    reco::Particle::LorentzVector z_momentum = muon1_momentum + muon2_momentum;

    ROOT::Math::Boost booster(z_momentum.BoostToCM());
    ROOT::Math::Boost invbooster(booster.Inverse());
    
    reco::Particle::LorentzVector Zb = booster(z_momentum);

    reco::Particle::LorentzVector muon1b = booster(muon1_momentum);
    reco::Particle::LorentzVector muon2b = booster(muon2_momentum);
    
    double tau_mass2 = targetParticleMass_*targetParticleMass_;

    double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
    double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;

    float scaling1 = sqrt(tauxb_mom2 / muonxb_mom2);
    float scaling2 = scaling1;

    float tauEnergy= Zb.t() / 2.;

    if (tauEnergy*tauEnergy<tau_mass2)
        return;
    
    reco::Particle::LorentzVector tau1b_mom = reco::Particle::LorentzVector(scaling1*muon1b.x(),scaling1*muon1b.y(),scaling1*muon1b.z(),tauEnergy);
    reco::Particle::LorentzVector tau2b_mom = reco::Particle::LorentzVector(scaling2*muon2b.x(),scaling2*muon2b.y(),scaling2*muon2b.z(),tauEnergy);

    reco::Particle::LorentzVector tau1_mom = (invbooster(tau1b_mom));
    reco::Particle::LorentzVector tau2_mom = (invbooster(tau2b_mom));
    
    // some additional checks
    // the following tests guarantee a deviation of less
    // than 0.1% for the following values of the original
    // muons and the placed taus
    //  invariant mass
    //  transverse momentum
    assert(std::abs((muon1_momentum+muon2_momentum).mass()-(tau1_mom+tau2_mom).mass())/(muon1_momentum+muon2_momentum).mass()<0.001);
    assert(std::abs((muon1_momentum+muon2_momentum).pt()-(tau1_mom+tau2_mom).pt())/(muon1_momentum+muon2_momentum).pt()<0.001);

    muon1->setP4(tau1_mom);
    muon2->setP4(tau2_mom);

    muon1->setPdgId(targetParticlePdgID_*muon1->pdgId()/abs(muon1->pdgId()));
    muon2->setPdgId(targetParticlePdgID_*muon2->pdgId()/abs(muon2->pdgId()));

    muon1->setStatus(1);
    muon2->setStatus(1);

    return;
}
void ParticleReplacerClass::transformMuMu2TauNu(reco::Particle * part1, reco::Particle * part2)
{
    using namespace edm;
    using namespace reco;
    using namespace std;

    reco::Particle::LorentzVector muon1_momentum = part1->p4();
    reco::Particle::LorentzVector muon2_momentum =  part2->p4();
    reco::Particle::LorentzVector z_momentum = muon1_momentum + muon2_momentum;

    ROOT::Math::Boost booster(z_momentum.BoostToCM());
    ROOT::Math::Boost invbooster(booster.Inverse());

    reco::Particle::LorentzVector Zb = booster(z_momentum);

    const double breitWignerWidth_Z = 2.4952;
    const double breitWignerWidth_W = 2.141;
    const double knownMass_W = 80.398;
    const double knownMass_Z = 91.1876;
              
    double Wb_mass = ( Zb.mass() - knownMass_Z ) * ( breitWignerWidth_W / breitWignerWidth_Z ) + knownMass_W;
    std::cout << "Wb_mass: " << Wb_mass << "\n";

    reco::Particle::LorentzVector muon1b = booster(muon1_momentum);
    reco::Particle::LorentzVector muon2b = booster(muon2_momentum);

    double tau_mass2 = targetParticleMass_*targetParticleMass_;

    double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
    double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;

    float scaling1 = sqrt(tauxb_mom2 / muonxb_mom2) * Wb_mass/Zb.mass();
    float scaling2 = scaling1;

    float tauEnergy= Zb.t() / 2.;

    if (tauEnergy*tauEnergy<tau_mass2)
              return;

    reco::Particle::LorentzVector tau1b_mom = reco::Particle::LorentzVector(scaling1*muon1b.x(),scaling1*muon1b.y(),scaling1*muon1b.z(),tauEnergy* Wb_mass/Zb.mass());
    reco::Particle::LorentzVector tau2b_mom = reco::Particle::LorentzVector(scaling2*muon2b.x(),scaling2*muon2b.y(),scaling2*muon2b.z(),tauEnergy* Wb_mass/Zb.mass());

    std::cout << "muon1b_momentum: " << muon1b << "\n";
    std::cout << "muon2b_momentum: " << muon2b << "\n";

    std::cout << "tau1b_momentum: " << tau1b_mom << "\n";
    std::cout << "tau2b_momentum: " << tau2b_mom << "\n";

    std::cout << "zb_momentum: " << Zb << "\n";
    std::cout << "wb_momentum: " << (tau1b_mom+tau2b_mom) << "\n";
                      
    // some checks
    // the following test guarantees a deviation
    // of less than 0.1% for phi and theta for the
    // original muons and the placed taus
    // (in the centre-of-mass system of the z boson)
    assert((muon1b.phi()-tau1b_mom.phi())/muon1b.phi()<0.001);
    assert((muon2b.phi()-tau2b_mom.phi())/muon2b.phi()<0.001);
    assert((muon1b.theta()-tau1b_mom.theta())/muon1b.theta()<0.001);
    assert((muon2b.theta()-tau2b_mom.theta())/muon2b.theta()<0.001);        

    reco::Particle::LorentzVector tau1_mom = (invbooster(tau1b_mom));
    reco::Particle::LorentzVector tau2_mom = (invbooster(tau2b_mom));

    // some additional checks
    // the following tests guarantee a deviation of less
    // than 0.1% for the following values of the original
    // muons and the placed taus
    //      invariant mass
    //      transverse momentum
    //assert(((muon1_momentum+muon1_momentum).mass()-(tau1_mom+tau2_mom).mass())/(muon1_momentum+muon1_momentum).mass()<0.001);
    //assert(((muon1_momentum+muon2_momentum).pt()-(tau1_mom+tau2_mom).pt())/(muon1_momentum+muon1_momentum).pt()<0.001);

    part1->setP4(tau1_mom);
    part2->setP4(tau2_mom);

    part1->setPdgId(15*part1->pdgId()/abs(part1->pdgId()));
    part2->setPdgId(16*part2->pdgId()/abs(part2->pdgId()));

    part1->setStatus(1);
    part2->setStatus(1);

    return;
}


//define this as a plug-in
//DEFINE_FWK_MODULE(Replacer);