#include <ParticleReplacerClass.h>

Inheritance diagram for ParticleReplacerClass:

Classes
struct	MinVisPtCut

Public Member Functions
virtual void	beginRun (edm::Run &iRun, const edm::EventSetup &iSetup)

virtual void	endJob ()

	ParticleReplacerClass (const edm::ParameterSet &, bool)

virtual std::auto_ptr < HepMC::GenEvent >	produce (const reco::MuonCollection &muons, const reco::Vertex pvtx=0, const HepMC::GenEvent genEvt=0)

	~ParticleReplacerClass ()

Public Member Functions inherited from ParticleReplacerBase
virtual void	beginJob ()

virtual void	endRun ()

	ParticleReplacerBase (const edm::ParameterSet &iConfig)

virtual	~ParticleReplacerBase ()

Private Member Functions
void	cleanEvent (HepMC::GenEvent evt, HepMC::GenVertex vtx)

HepMC::GenEvent *	processEventWithPythia (HepMC::GenEvent *evt)

HepMC::GenEvent *	processEventWithTauola (HepMC::GenEvent *evt)

void	repairBarcodes (HepMC::GenEvent *evt)

bool	testEvent (HepMC::GenEvent *evt)

void	transformMuMu2TauNu (reco::Particle muon1, reco::Particle muon2)
	transform a muon pair into tau nu (as coming from a W boson) More...

void	transformMuMu2TauTau (reco::Particle muon1, reco::Particle muon2)
	transform a muon pair into a tau pair More...

Private Attributes
int	attempts

std::string	generatorMode_

int	maxNumberOfAttempts_

std::vector< std::vector < MinVisPtCut > >	minVisPtCuts_

int	motherParticleID_

bool	noInitialisation_

TTree *	outTree

bool	printEvent_

unsigned int	replacementMode_
	replace mode: More...

double	targetParticleMass_

int	targetParticlePdgID_

gen::TauolaInterfaceBase *	tauola_

unsigned int	transformationMode_

bool	useExternalGenerators_

bool	useTauola_

bool	useTauolaPolarization_

Additional Inherited Members
Public Attributes inherited from ParticleReplacerBase
double	eventWeight

Protected Attributes inherited from ParticleReplacerBase
const double	tauMass

Detailed Description

Definition at line 56 of file ParticleReplacerClass.h.

Constructor & Destructor Documentation

ParticleReplacerClass::ParticleReplacerClass	(	const edm::ParameterSet &	pset,
		bool	verbose
	)

explicit

Definition at line 25 of file ParticleReplacerClass.cc.

References attempts, GOODCOLL_filter_cfg::cut, hpstanc_transforms::cuts, decayRandomEngine, ParticleReplacerClass::MinVisPtCut::ELEC, edm::hlt::Exception, generatorMode_, reco::get(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), ParticleReplacerClass::MinVisPtCut::HAD, ParticleReplacerClass::MinVisPtCut::index_, minVisPtCuts_, motherParticleID_, ParticleReplacerClass::MinVisPtCut::MU, noInitialisation_, NULL, outTree, pos, ParticleReplacerClass::MinVisPtCut::pt_, replacementMode_, gen::TauolaInterfaceBase::SetDecayRandomEngine(), ParticleReplacerClass::MinVisPtCut::TAU, tauola_, transformationMode_, and ParticleReplacerClass::MinVisPtCut::type_.

                                                                                      :
   ParticleReplacerBase(pset),
   generatorMode_(pset.getParameter<std::string>("generatorMode")),
   printEvent_(verbose),
   outTree(0),
   maxNumberOfAttempts_(pset.getUntrackedParameter<int>("maxNumberOfAttempts", 1000))
 {
   tauola_ = (gen::TauolaInterfaceBase*)(TauolaFactory::get()->create("Tauola271215", 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
     CLHEP::HepRandomEngine* decayRandomEngine = &rng->getEngine();
     tauola_->SetDecayRandomEngine(decayRandomEngine);
 
     edm::LogInfo("Replacer") << "generatorMode = "<< generatorMode_<< "\n";
     edm::LogInfo("Replacer") << "replacementMode = "<< replacementMode_<< "\n";
 
     return;
 }

ParticleReplacerClass::~ParticleReplacerClass ( )

Definition at line 161 of file ParticleReplacerClass.cc.

References tauola_.

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

Member Function Documentation

void ParticleReplacerClass::beginRun	(	edm::Run &	iRun,
		const edm::EventSetup &	iSetup
	)

virtual

Reimplemented from ParticleReplacerBase.

Definition at line 451 of file ParticleReplacerClass.cc.

References gen::TauolaInterfaceBase::init(), and tauola_.

 {
     tauola_->init(iSetup);
 }

void ParticleReplacerClass::cleanEvent	(	HepMC::GenEvent *	evt,
		HepMC::GenVertex *	vtx
	)

private

Definition at line 550 of file ParticleReplacerClass.cc.

References dt_dqm_sourceclient_common_cff::reco, and repairBarcodes().

Referenced by produce().

 {
     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::endJob ( void )

virtual

Reimplemented from ParticleReplacerBase.

Definition at line 458 of file ParticleReplacerClass.cc.

References gen::TauolaInterfaceBase::statistics(), and tauola_.

 {
     tauola_->statistics();
 }

HepMC::GenEvent* ParticleReplacerClass::processEventWithPythia ( HepMC::GenEvent * evt )

private

HepMC::GenEvent* ParticleReplacerClass::processEventWithTauola ( HepMC::GenEvent * evt )

private

std::auto_ptr< HepMC::GenEvent > ParticleReplacerClass::produce	(	const reco::MuonCollection &	muons,
		const reco::Vertex *	pvtx = `0`,
		const HepMC::GenEvent *	genEvt = `0`
	)

virtual

2) transform the muons to the desired particles

3) prepare the event

3) process the event

Implements ParticleReplacerBase.

Definition at line 169 of file ParticleReplacerClass.cc.

References abs, attempts, reco::LeafCandidate::charge(), cleanEvent(), conv, gather_cfg::cout, gen::TauolaInterfaceBase::decay(), ParticleReplacerBase::eventWeight, edm::hlt::Exception, generatorMode_, configurableAnalysis::GenParticle, i, prof2calltree::l, maxNumberOfAttempts_, motherParticleID_, outTree, AlCaHLTBitMon_ParallelJobs::p, reco::LeafCandidate::p4(), reco::LeafCandidate::pdgId(), printEvent_, reco::LeafCandidate::px(), reco::LeafCandidate::py(), reco::LeafCandidate::pz(), repairBarcodes(), replacementMode_, run_regression::ret, reco::Particle::setStatus(), mathSSE::sqrt(), targetParticleMass_, targetParticlePdgID_, tauola_, testEvent(), transformationMode_, transformMuMu2TauNu(), transformMuMu2TauTau(), funct::true, and reco::LeafCandidate::vertex().

 {
     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;
                 }
     }
     eventWeight = (double)cntVisPt_pass / (double)cntVisPt_all;
     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;
         eventWeight=0;
                 if(outTree) outTree->Fill();
         return std::auto_ptr<HepMC::GenEvent>(0);
     }
     attempts=nr_of_trials;
     if(outTree) outTree->Fill();    
 
     // recover the status codes
     if (replacementMode_==0)
     {
         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;
 }

void ParticleReplacerClass::repairBarcodes ( HepMC::GenEvent * evt )

private

Definition at line 608 of file ParticleReplacerClass.cc.

Referenced by cleanEvent(), and produce().

 {
     using namespace HepMC;
 
     // repair the barcodes
     int max_barc=0;
     for (GenEvent::vertex_iterator it=evt->vertices_begin();it!=evt->vertices_end();it++)
         while (!(*it)->suggest_barcode(-1*(++max_barc)))
             ;
 
     max_barc=0;
     for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
         while (!(*it)->suggest_barcode(++max_barc))
             ;
 }

bool ParticleReplacerClass::testEvent ( HepMC::GenEvent * evt )

private

Definition at line 463 of file ParticleReplacerClass.cc.

References abs, runEdmFileComparison::collection, ParticleReplacerClass::MinVisPtCut::ELEC, prof2calltree::front, configurableAnalysis::GenParticle, ParticleReplacerClass::MinVisPtCut::HAD, minVisPtCuts_, ParticleReplacerClass::MinVisPtCut::MU, benchmark_cfg::pdgId, python.multivaluedict::sort(), lumiQTWidget::t, and ParticleReplacerClass::MinVisPtCut::TAU.

Referenced by produce().

 {
     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::transformMuMu2TauNu	(	reco::Particle *	muon1,
		reco::Particle *	muon2
	)

private

transform a muon pair into tau nu (as coming from a W boson)

Definition at line 683 of file ParticleReplacerClass.cc.

References abs, gather_cfg::cout, reco::Particle::p4(), reco::Particle::pdgId(), dt_dqm_sourceclient_common_cff::reco, reco::Particle::setP4(), reco::Particle::setPdgId(), reco::Particle::setStatus(), mathSSE::sqrt(), and targetParticleMass_.

Referenced by produce().

 {
     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;
 }

void ParticleReplacerClass::transformMuMu2TauTau	(	reco::Particle *	muon1,
		reco::Particle *	muon2
	)

private

transform a muon pair into a tau pair

Definition at line 625 of file ParticleReplacerClass.cc.

References abs, scaleCards::mass, reco::Particle::p4(), reco::Particle::pdgId(), dt_dqm_sourceclient_common_cff::reco, reco::Particle::setP4(), reco::Particle::setPdgId(), reco::Particle::setStatus(), mathSSE::sqrt(), targetParticleMass_, and targetParticlePdgID_.

Referenced by produce().

 {
     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;
 }