#include <Generator.h>

Public Member Functions
virtual const double	eventWeight () const

	Generator (const edm::ParameterSet &p)

virtual const HepMC::GenEvent *	genEvent () const

virtual const math::XYZTLorentzVector *	genVertex () const

void	HepMC2G4 (const HepMC::GenEvent g, G4Event e)

void	nonBeamEvent2G4 (const HepMC::GenEvent g, G4Event e)

void	setGenEvent (const HepMC::GenEvent *inpevt)

virtual	~Generator ()

Private Member Functions
bool	isExotic (int pdgcode) const

bool	isExoticNonDetectable (int pdgcode) const

bool	IsInTheFilterList (int pdgcode) const

void	particleAssignDaughters (G4PrimaryParticle p, HepMC::GenParticle hp, double length)

bool	particlePassesPrimaryCuts (const G4ThreeVector &p) const

void	setGenId (G4PrimaryParticle *p, int id) const

Private Attributes
HepMC::GenEvent *	evt_

bool	fEtaCuts

LumiMonitorFilter *	fLumiFilter

bool	fPCuts

bool	fPDGFilter

bool	fPhiCuts

bool	fPtransCut

std::vector< int >	pdgFilter

bool	pdgFilterSel

double	theDecLenCut

double	theDecRCut2

double	theEtaCutForHector

double	theMaxEtaCut

double	theMaxPCut

double	theMaxPhiCut

double	theMinEtaCut

double	theMinPCut

double	theMinPhiCut

double	theMinPtCut2

int	verbose

math::XYZTLorentzVector *	vtx_

double	weight_

double	Z_hector

double	Z_lmax

double	Z_lmin

Detailed Description

Definition at line 19 of file Generator.h.

Constructor & Destructor Documentation

Generator::Generator ( const edm::ParameterSet & p )

Definition at line 24 of file Generator.cc.

References LumiMonitorFilter::Describe(), edm::ParameterSet::exists(), JetChargeProducer_cfi::exp, fEtaCuts, fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, edm::ParameterSet::getParameter(), cuy::ii, csvLumiCalc::lumi, pdgFilter, pdgFilterSel, theDecLenCut, theDecRCut2, theEtaCutForHector, theMaxEtaCut, theMinEtaCut, theMinPCut, theMinPtCut2, Z_hector, Z_lmax, and Z_lmin.

                                            : 
   fPCuts(p.getParameter<bool>("ApplyPCuts")),
   fPtransCut(p.getParameter<bool>("ApplyPtransCut")),
   fEtaCuts(p.getParameter<bool>("ApplyEtaCuts")), 
   fPhiCuts(p.getParameter<bool>("ApplyPhiCuts")),
   theMinPhiCut(p.getParameter<double>("MinPhiCut")), // in radians (CMS standard)
   theMaxPhiCut(p.getParameter<double>("MaxPhiCut")),
   theMinEtaCut(p.getParameter<double>("MinEtaCut")),
   theMaxEtaCut(p.getParameter<double>("MaxEtaCut")),
   theMinPCut(p.getParameter<double>("MinPCut")),    // in GeV (CMS standard)
   theMaxPCut(p.getParameter<double>("MaxPCut")),   
   theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
   verbose(p.getUntrackedParameter<int>("Verbosity",0)),
   fLumiFilter(nullptr),
   evt_(nullptr),
   vtx_(nullptr),
   weight_(0),
   Z_lmin(0),
   Z_lmax(0),
   Z_hector(0),
   pdgFilterSel(false), 
   fPDGFilter(false)
 {
   bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
   if(lumi) { fLumiFilter = new LumiMonitorFilter(); }
 
   double theRDecLenCut = p.getParameter<double>("RDecLenCut")*cm;
   theDecRCut2 = theRDecLenCut*theRDecLenCut;
 
   theMinPtCut2 = theMinPCut*theMinPCut;
 
   double theDecLenCut = p.getParameter<double>("LDecLenCut")*cm;
 
   pdgFilter.resize(0);
   if ( p.exists("PDGselection") ) {
     pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).
       getParameter<bool>("PDGfilterSel");
     pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).
       getParameter<std::vector< int > >("PDGfilter");
     if(!pdgFilter.empty()) { 
       fPDGFilter = true; 
       std::stringstream ss;
       ss << "SimG4Core/Generator: ";
       if (pdgFilterSel) {
     ss << " Selecting only PDG ID = ";
       } else {
     ss << " Filtering out PDG ID = ";
       }
       for ( unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
     ss << pdgFilter[ii] << "  ";
       }
       edm::LogInfo("SimG4CoreGenerator") << ss.str();
     }
   }
 
   if(fEtaCuts){
     Z_lmax = theRDecLenCut*((1 - exp(-2*theMaxEtaCut) )/( 2*exp(-theMaxEtaCut)));
     Z_lmin = theRDecLenCut*((1 - exp(-2*theMinEtaCut) )/( 2*exp(-theMinEtaCut)));
   }
 
   Z_hector = theRDecLenCut*((1 - exp(-2*theEtaCutForHector)) 
                 / ( 2*exp(-theEtaCutForHector) ) );
 
   edm::LogInfo("SimG4CoreGenerator") 
     << " Rdecaycut= " << theRDecLenCut/cm 
     << " cm;  Zdecaycut= " << theDecLenCut/cm 
     << "Z_min= " << Z_lmin/cm << " cm; Z_max= " << Z_lmax 
     << " cm;  Z_hector = " << Z_hector << " cm\n"
     << "ApplyPCuts: " << fPCuts << "  ApplyPtransCut: " << fPtransCut
     << "  ApplyEtaCuts: " << fEtaCuts 
     << "  ApplyPhiCuts: " << fPhiCuts
     << "  ApplyLumiMonitorCuts: " << lumi;
   if(lumi) { fLumiFilter->Describe(); }
 }

Generator::~Generator ( )

virtual

Definition at line 99 of file Generator.cc.

References fLumiFilter.

 {
   delete fLumiFilter;
 }

Member Function Documentation

virtual const double Generator::eventWeight ( ) const

inlinevirtual

Definition at line 31 of file Generator.h.

References GenParticle::GenParticle, AnalysisDataFormats_SUSYBSMObjects::hp, isExotic(), isExoticNonDetectable(), IsInTheFilterList(), particleAssignDaughters(), particlePassesPrimaryCuts(), and weight_.

Referenced by RunManagerMTWorker::produce(), and RunManager::produce().

31 { return weight_; }

Generator::weight_

double weight_

Definition: Generator.h:63

virtual const HepMC::GenEvent* Generator::genEvent ( ) const

inlinevirtual

Definition at line 29 of file Generator.h.

References evt_.

Referenced by RunManagerMTWorker::produce(), and RunManager::produce().

29 { return evt_; }

Generator::evt_

HepMC::GenEvent * evt_

Definition: Generator.h:61

virtual const math::XYZTLorentzVector* Generator::genVertex ( ) const

inlinevirtual

Definition at line 30 of file Generator.h.

References vtx_.

Referenced by RunManagerMTWorker::produce(), and RunManager::produce().

30 { return vtx_; }

Generator::vtx_

math::XYZTLorentzVector * vtx_

Definition: Generator.h:62

void Generator::HepMC2G4	(	const HepMC::GenEvent *	g,
		G4Event *	e
	)

Definition at line 104 of file Generator.cc.

References funct::abs(), ALCARECOTkAlJpsiMuMu_cff::charge, class-composition::children, fEtaCuts, fLumiFilter, fPCuts, fPDGFilter, fPhiCuts, fPtransCut, GenParticle::GenParticle, GeV, isExotic(), isExoticNonDetectable(), LumiMonitorFilter::isGoodForLumiMonitor(), IsInTheFilterList(), LogDebug, nullptr, AlCaHLTBitMon_ParallelJobs::p, particleAssignDaughters(), pdgFilterSel, phi, setGenId(), mathSSE::sqrt(), mps_update::status, theDecRCut2, theMaxPCut, theMaxPhiCut, theMinPtCut2, vtx_, weight_, globals_cff::x1, globals_cff::x2, hybridSuperClusters_cfi::xi, geometryCSVtoXML::xx, geometryCSVtoXML::yy, Z_hector, Z_lmax, and Z_lmin.

Referenced by RunManagerMTWorker::generateEvent(), RunManager::generateEvent(), and setGenEvent().

 {
 
   HepMC::GenEvent *evt=new HepMC::GenEvent(*evt_orig);
 
   if ( *(evt->vertices_begin()) == nullptr ) {
     std::stringstream ss;
     ss << "SimG4Core/Generator: in event " << g4evt->GetEventID()
        << " Corrupted Event - GenEvent with no vertex \n"; 
     throw SimG4Exception(ss.str());
   }  
   
   if (!evt->weights().empty()) {
 
     weight_ = evt->weights()[0] ;
     for (unsigned int iw=1; iw<evt->weights().size(); ++iw) {
 
       // terminate if the vector of weights contains a zero-weight
       if ( evt->weights()[iw] <= 0 ) break;
       weight_ *= evt->weights()[iw] ;
     }     
   }
   
   if (vtx_ != nullptr) { delete vtx_; }
   vtx_ = new math::XYZTLorentzVector((*(evt->vertices_begin()))->position().x(),
                                      (*(evt->vertices_begin()))->position().y(),
                                      (*(evt->vertices_begin()))->position().z(),
                                      (*(evt->vertices_begin()))->position().t());
 
   if(verbose > 0) {
     edm::LogInfo("SimG4CoreGenerator") << &evt;
     LogDebug("SimG4CoreGenerator") << "Primary Vertex = (" 
                    << vtx_->x() << "," 
                    << vtx_->y() << ","
                    << vtx_->z() << ")";
   }
     
   unsigned int ng4vtx = 0;
 
   for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
       vitr != evt->vertices_end(); ++vitr ) { 
 
     // loop for vertex, is it a real vertex?
     // Set qvtx to true for any particles that should be propagated by GEANT, i.e.,
     // status 1 particles or
     // status 2 particles that decay outside the beampipe.   
     G4bool qvtx=false;
     HepMC::GenVertex::particle_iterator pitr;
     for (pitr= (*vitr)->particles_begin(HepMC::children);
          pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
 
       // For purposes of this function, the status is defined as follows:
       // 1:  particles are not decayed by generator
       // 2:  particles are decayed by generator but need to be propagated by GEANT
       // 3:  particles are decayed by generator but do not need to be propagated by GEANT
       int status = (*pitr)->status();
       int pdg = (*pitr)->pdg_id();
       if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
         // In Pythia 8, there are many status codes besides 1, 2, 3.
         // By setting the status to 2 for exotic particles, they will be checked:
         // if its decay vertex is outside the beampipe, it will be propagated by GEANT.
     // Some Standard Model particles, e.g., K0, cannot be propagated by GEANT, 
     // so do not change their status code.  
         status = 2;
       }
 
       // Particles which are not decayed by generator
       if (status == 1) {
 
     // filter out unwanted particles and vertices
     if(fPDGFilter && !pdgFilterSel && IsInTheFilterList(pdg)) { 
       continue;
     }
 
         qvtx = true;
         if (verbose > 2) LogDebug("SimG4CoreGenerator") 
       << "GenVertex barcode = " << (*vitr)->barcode() 
       << " selected for GenParticle barcode = " << (*pitr)->barcode();
         break;
       }  
       // The selection is made considering if the partcile with status = 2 
       // have the end_vertex with a radius greater than the radius of beampipe 
       // cylinder (no requirement on the Z of the vertex is applyed).
       else if (status == 2) {
         if ( (*pitr)->end_vertex() != nullptr  ) { 
           double xx = (*pitr)->end_vertex()->position().x();
           double yy = (*pitr)->end_vertex()->position().y();
           double r_dd = xx*xx+yy*yy;
           if (r_dd > theDecRCut2) {
             qvtx = true;
             if (verbose > 2) LogDebug("SimG4CoreGenerator") 
           << "GenVertex barcode = " << (*vitr)->barcode() 
           << " selected for GenParticle barcode = " 
           << (*pitr)->barcode() << " radius = " << std::sqrt(r_dd);
             break;
           }
         } else {
       // particles with status 2 without end_vertex are 
       // equivalent to stable
       qvtx = true;
       break;
     }
       }
     }
 
     // if this vertex is inside fiductial volume inside the beam pipe
     // and has no long-lived secondary the vertex is not saved 
     if (!qvtx) {
       continue;
     }
     
     double x1 = (*vitr)->position().x()*mm;
     double y1 = (*vitr)->position().y()*mm;
     double z1 = (*vitr)->position().z()*mm;
     double t1 = (*vitr)->position().t()*mm/c_light;
 
     G4PrimaryVertex* g4vtx = new G4PrimaryVertex(x1, y1, z1, t1);
 
     for (pitr= (*vitr)->particles_begin(HepMC::children);
          pitr != (*vitr)->particles_end(HepMC::children); ++pitr){
 
       int status = (*pitr)->status();
       int pdg = (*pitr)->pdg_id();
       bool hasDecayVertex = (nullptr != (*pitr)->end_vertex()) ? true : false;
 
       // Filter on allowed particle species if required      
       if (fPDGFilter) { 
     bool isInTheList = IsInTheFilterList(pdg);
     if((!pdgFilterSel && isInTheList) || (pdgFilterSel && !isInTheList)) { 
       edm::LogInfo("SimG4CoreGenerator") 
         << " Skiped GenParticle barcode= " << (*pitr)->barcode() 
         << " PDGid= " << pdg
         << " status= " << status 
         << " isExotic: " << isExotic(pdg)
         << " isExoticNotDet: " << isExoticNonDetectable(pdg)
         << " isInTheList: " << isInTheList
         << " hasDecayVertex: " << hasDecayVertex;
           continue;
         }
       }
 
       edm::LogInfo("SimG4CoreGenerator") 
     << " pdg= " << pdg
     << " status= " << status 
     << " hasPreDefinedDecay: " << hasDecayVertex
     << " isExotic: " << isExotic(pdg)
     << " isExoticNotDet: " << isExoticNonDetectable(pdg)
     << " isInTheList: " << IsInTheFilterList(pdg);
 
       if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg))) ) {
     status = hasDecayVertex ? 2 : 1;
       } 
 
       // this particle has predefined decay but has no vertex
       if (2 == status && !hasDecayVertex) { 
     edm::LogWarning("SimG4CoreGenerator: in event ") 
       << g4evt->GetEventID()
       << " a particle " 
       << " pdgid= " << pdg
       << " has status=2 but has no decay vertex, so will be fully tracked by Geant4";
     status = 1; 
       }
 
       double x2 = x1;
       double y2 = y1;
       double z2 = z1;
       double decay_length = 0.0;
       if(2 == status) {
     x2 = (*pitr)->end_vertex()->position().x();
     y2 = (*pitr)->end_vertex()->position().y();
     z2 = (*pitr)->end_vertex()->position().z();
     decay_length = std::sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
       } 
 
       bool toBeAdded = false;
       double px = (*pitr)->momentum().px();
       double py = (*pitr)->momentum().py();
       double pz = (*pitr)->momentum().pz();
       double ptot = std::sqrt(px*px + py*py + pz*pz);
       math::XYZTLorentzVector p(px, py, pz, (*pitr)->momentum().e());
       
       double ximpact = x1;
       double yimpact = y1;
       double zimpact = z1;
 
       // protection against numerical problems for extremely low momenta
       // compute impact point at transition to Hector
       const double minTan = 1.e-20;
       if(std::abs(z1) < Z_hector && std::abs(pz) >= minTan*ptot) {
         if(pz > 0.0) { zimpact =  Z_hector; }
     else         { zimpact = -Z_hector; }
         double del = (zimpact - z1)/pz;
         ximpact += del*px;
         yimpact += del*py;
       }
       double rimpact2 = ximpact*ximpact + yimpact*yimpact;
       
       if (verbose > 2) LogDebug("SimG4CoreGenerator") 
     << "Processing GenParticle barcode= " << (*pitr)->barcode() 
     << " pdg= " << pdg
     << " status= " << (*pitr)->status() 
     << " st= " << status 
     << " rimpact(cm)= " << std::sqrt(rimpact2)/cm
     << " zimpact(cm)= " << zimpact/cm
     << " ptot(GeV)= " << ptot
     << " pz(GeV)= " << pz;
 
       // Particles of status 1 trasnported along the beam pipe for forward 
       // detectors (HECTOR) always pass to Geant4 without cuts 
       if( 1 == status && 
       std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
         toBeAdded = true;
         if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
       << "GenParticle barcode = " << (*pitr)->barcode() 
       << " passed case 3";
       } else {
 
     // Standard case: particles not decayed by the generator
     if(1 == status && 
            (std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
 
       // Ptot cut for all particles
       if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
             continue;
       }
       // phi cut is applied mainly for particle gun
       if (fPhiCuts) {
         double phi = p.phi();
         if(phi < theMinPhiCut || phi  > theMaxPhiCut) {
           continue;
         }
       }
       // eta cut is applied if position of the decay
       // is within vacuum chamber and limited in Z
       if(fEtaCuts) {
 
         // eta cut
         double xi = x1;
         double yi = y1;
         double zi = z1;
 
         // can be propagated along Z
         if(std::abs(pz) >= minTan*ptot) {
           if((zi >= Z_lmax) & (pz < 0.0)) {
         zi = Z_lmax;
           } else if((zi <= Z_lmin) & (pz > 0.0)) {
         zi = Z_lmin;
           } else {
         if(pz > 0) { zi = Z_lmax; }
         else { zi = Z_lmin; }
           }
           double del = (zi - z1)/pz;
           xi += del*px;
           yi += del*py;
         }
         // check eta cut
         if((zi >= Z_lmin) & (zi <= Z_lmax) 
            & (xi*xi + yi*yi < theDecRCut2)) {
           continue;
         }
       }
       if(fLumiFilter && !fLumiFilter->isGoodForLumiMonitor(*pitr)) 
         { continue; }
       toBeAdded = true;
       if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
         << "GenParticle barcode = " << (*pitr)->barcode() 
         << " passed case 1";
 
       // Decay chain outside the fiducial cylinder defined by theRDecLenCut
       // are used for Geant4 tracking with predefined decay channel
       // In the case of decay in vacuum particle is not tracked by Geant4
     } else if(2 == status &&
           x2*x2 + y2*y2 >= theDecRCut2 && std::abs(z2) < Z_hector) {
 
       toBeAdded = true;
       if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
         << "GenParticle barcode = " << (*pitr)->barcode() 
         << " passed case 2" 
         << " decay_length(cm)= " << decay_length/cm;
     }
       }
       if(toBeAdded){
         G4PrimaryParticle* g4prim= 
           new G4PrimaryParticle(pdg, px*GeV, py*GeV, pz*GeV);
         
         if ( g4prim->GetG4code() != nullptr ){ 
           g4prim->SetMass( g4prim->GetG4code()->GetPDGMass() );
           double charge = g4prim->GetG4code()->GetPDGCharge();
 
       // apply Pt cut
       if (fPtransCut && 
           0.0 != charge && px*px + py*py < theMinPtCut2) {
             delete g4prim;
             continue;
       }
           g4prim->SetCharge(charge);  
         }
 
     // V.I. do not use SetWeight but the same code
         // value of the code compute inside TrackWithHistory        
         //g4prim->SetWeight( 10000*(*vpitr)->barcode() ) ;
         setGenId( g4prim, (*pitr)->barcode() );
 
         if (2 == status) {
           particleAssignDaughters(g4prim,
                   (HepMC::GenParticle *) *pitr, decay_length);
     }
         if ( verbose > 1 ) g4prim->Print();
         g4vtx->SetPrimary(g4prim);
       }
     }
 
     if ( verbose > 1 ) g4vtx->Print();
     g4evt->AddPrimaryVertex(g4vtx);
     ++ng4vtx;
   }
   
   // Add a protection for completely empty events (produced by LHCTransport): 
   // add a dummy vertex with no particle attached to it
   if ( ng4vtx == 0 ) {
     G4PrimaryVertex* g4vtx = new G4PrimaryVertex(0.0, 0.0, 0.0, 0.0);
     if ( verbose > 1 ) g4vtx->Print();
     g4evt->AddPrimaryVertex(g4vtx);
   }
 
   delete evt;
 }

bool Generator::isExotic ( int pdgcode ) const

private

Definition at line 542 of file Generator.cc.

References funct::abs(), and BPhysicsValidation_cfi::pdgid.

Referenced by eventWeight(), and HepMC2G4().

 {
   int pdgid = std::abs(pdgcode);  
   return ((pdgid >= 1000000 && pdgid <  4000000) || // SUSY, R-hadron, and technicolor particles
       pdgid == 17 || // 4th generation lepton 
       pdgid == 34 || // W-prime
       pdgid == 37)   // charged Higgs
     ? true : false;
 }

bool Generator::isExoticNonDetectable ( int pdgcode ) const

private

Definition at line 552 of file Generator.cc.

References funct::abs(), ALCARECOTkAlJpsiMuMu_cff::charge, source_particleGun_cfi::ParticleID, BPhysicsValidation_cfi::pdgid, and sysUtil::pid.

Referenced by eventWeight(), and HepMC2G4().

 {
   int pdgid = std::abs(pdgcode);  
   HepPDT::ParticleID pid(pdgcode);
   int charge = pid.threeCharge();
   return ((charge==0) && (pdgid >= 1000000 && pdgid <  1000040)) // SUSY 
     ? true : false;
 }

bool Generator::IsInTheFilterList ( int pdgcode ) const

private

Definition at line 561 of file Generator.cc.

References funct::abs(), pdgFilter, and BPhysicsValidation_cfi::pdgid.

Referenced by eventWeight(), and HepMC2G4().

 {
   int pdgid = std::abs(pdgcode);
   for(auto & pdg : pdgFilter) { if(std::abs(pdg) == pdgid) { return true; } }
   return false;
 }

void Generator::nonBeamEvent2G4	(	const HepMC::GenEvent *	g,
		G4Event *	e
	)

Definition at line 568 of file Generator.cc.

References GenParticle::GenParticle, GeV, runTauDisplay::gp, mps_fire::i, particlePassesPrimaryCuts(), setGenId(), and findQualityFiles::v.

Referenced by RunManagerMTWorker::generateEvent(), RunManager::generateEvent(), and setGenEvent().

 {
   int i = 0; 
   for(HepMC::GenEvent::particle_const_iterator it = evt->particles_begin(); 
       it != evt->particles_end(); ++it ) {
     ++i;
     HepMC::GenParticle * gp = (*it);    
     int g_status = gp->status();
     // storing only particle with status == 1   
     if (g_status == 1) {
       int g_id = gp->pdg_id();      
       G4PrimaryParticle * g4p = 
     new G4PrimaryParticle(g_id,gp->momentum().px()*GeV,
                   gp->momentum().py()*GeV,
                   gp->momentum().pz()*GeV);
       if (g4p->GetG4code() != nullptr) { 
     g4p->SetMass(g4p->GetG4code()->GetPDGMass());
     g4p->SetCharge(g4p->GetG4code()->GetPDGCharge());
       }
       setGenId(g4p,i);
       if (particlePassesPrimaryCuts(g4p->GetMomentum())) {
     G4PrimaryVertex * v = 
       new G4PrimaryVertex(gp->production_vertex()->position().x()*mm,
                   gp->production_vertex()->position().y()*mm,
                   gp->production_vertex()->position().z()*mm,
                   gp->production_vertex()->position().t()*mm/c_light);
     v->SetPrimary(g4p);
     g4evt->AddPrimaryVertex(v);
     if(verbose > 0) { v->Print(); }
 
       } else {
     delete g4p;
       }
     }
   } // end loop on HepMC particles
 }

void Generator::particleAssignDaughters	(	G4PrimaryParticle *	p,
		HepMC::GenParticle *	hp,
		double	length
	)

private

Definition at line 432 of file Generator.cc.

References class-composition::children, createTree::dd, GeV, LogDebug, nullptr, AlCaHLTBitMon_ParallelJobs::p, setGenId(), mathSSE::sqrt(), globals_cff::x1, and globals_cff::x2.

Referenced by eventWeight(), and HepMC2G4().

 {
   if ( verbose > 1 ) {
     LogDebug("SimG4CoreGenerator") 
       << "Special case of long decay length \n" 
       << "Assign daughters with to mother with decaylength=" 
       << decaylength/cm << " cm";
   }
   math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(), 
                 vp->momentum().pz(), vp->momentum().e());
 
   // defined preassigned decay time
   double proper_time = decaylength/(p.Beta()*p.Gamma()*c_light);
   g4p->SetProperTime(proper_time); 
 
   if( verbose > 2 ) {
     LogDebug("SimG4CoreGenerator") <<" px= "<< p.px()
                    <<" py= "<< p.py()
                    <<" pz= "<< p.pz()
                    <<" e= " <<  p.e()
                    <<" beta= "<< p.Beta()
                    <<" gamma= " << p.Gamma()
                    <<" Proper time= " <<proper_time/ns <<" ns";
   }
 
   // the particle will decay after the same length if it 
   // has not interacted before
   double x1 = vp->end_vertex()->position().x();
   double y1 = vp->end_vertex()->position().y();
   double z1 = vp->end_vertex()->position().z();
 
   for (HepMC::GenVertex::particle_iterator 
          vpdec= vp->end_vertex()->particles_begin(HepMC::children);
        vpdec != vp->end_vertex()->particles_end(HepMC::children); ++vpdec) {
     
     //transform decay products such that in the rest frame of mother
     math::XYZTLorentzVector pdec((*vpdec)->momentum().px(),
                                  (*vpdec)->momentum().py(),
                                  (*vpdec)->momentum().pz(),
                                  (*vpdec)->momentum().e());
 
     // children should only be taken into account once
     G4PrimaryParticle * g4daught= 
       new G4PrimaryParticle((*vpdec)->pdg_id(), pdec.x()*GeV, 
                 pdec.y()*GeV, pdec.z()*GeV);
 
     if ( g4daught->GetG4code() != nullptr )
       { 
         g4daught->SetMass( g4daught->GetG4code()->GetPDGMass() ) ;
         g4daught->SetCharge( g4daught->GetG4code()->GetPDGCharge() ) ;  
       }
 
     // V.I. do not use SetWeight but the same code
     // value of the code compute inside TrackWithHistory        
     setGenId( g4daught, (*vpdec)->barcode() );
 
     if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
       <<"Assigning a "<< (*vpdec)->pdg_id()
       <<" as daughter of a " << vp->pdg_id();
     if ( (*vpdec)->status() == 2 && (*vpdec)->end_vertex() != nullptr) 
       {
         double x2 = (*vpdec)->end_vertex()->position().x();
         double y2 = (*vpdec)->end_vertex()->position().y();
         double z2 = (*vpdec)->end_vertex()->position().z();
         double dd = std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
         particleAssignDaughters(g4daught,*vpdec,dd);
       }
     (*vpdec)->set_status(1000+(*vpdec)->status()); 
     g4p->SetDaughter(g4daught);
 
     if ( verbose > 1 ) g4daught->Print();
   }
 }

bool Generator::particlePassesPrimaryCuts ( const G4ThreeVector & p ) const

private

Definition at line 509 of file Generator.cc.

References MillePedeFileConverter_cfg::e, PVValHelper::eta, fEtaCuts, RemoveAddSevLevel::flag, fPCuts, fPhiCuts, GeV, LogDebug, phi, theMaxEtaCut, theMaxPCut, and theMaxPhiCut.

Referenced by eventWeight(), and nonBeamEvent2G4().

 {
   bool flag = true;
   double ptot = p.mag();
   if (fPCuts && (ptot < theMinPCut*GeV || ptot > theMaxPCut*GeV)) {
     flag = false;
   } 
   if (fEtaCuts && flag) {
     double pz = p.z();
     if(ptot < pz + 1.e-10) { 
       flag = false;
 
     } else {
       double eta = 0.5*G4Log((ptot + pz)/(ptot - pz));
       if(eta < theMinEtaCut || eta > theMaxEtaCut) {
     flag = false;
       }
     }
   }
   if (fPhiCuts && flag) {
     double phi = p.phi();
     if(phi < theMinPhiCut || phi > theMaxPhiCut) {
       flag = false;
     }
   }
   
   if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
     << "Generator ptot(GeV)= " << ptot/GeV << " eta= " << p.eta() 
     << "  phi= " << p.phi() << " Flag= " << flag;
 
   return flag;
 }