update GtGunToll

Examples/options/tut_detsim_pan_matrix.py

@@ -62,20 +62,13 @@ from Configurables import HepMCRdr
 from Configurables import GenPrinter
 
 gun = GtGunTool("GtGunTool")
-gun.Particles = ["gamma"]
-#gun.Energies =  [0.5, 1] # GeV
-#gun.EnergyMin = [0.1] # GeV
-gun.EnergyMin = [10] # GeV
-gun.EnergyMax = [10] # GeV
-gun.ThetaMins = [90] # rad; 45deg
-gun.ThetaMaxs = [90] # rad; 45deg
-gun.PhiMins   = [0] # rad; 0deg
-gun.PhiMaxs   = [0] # rad; 360deg
-
-gun.Create_second = True
-gun.Energy1       = 10
-gun.dtheta        = 0
-gun.dphi          = 1
+gun.Particles = ["gamma","gamma"]
+gun.Energies =  [10, 10] # GeV
+gun.ThetaMins = [90, 90] # degree
+gun.ThetaMaxs = [90, 90] # degree
+gun.PhiMins   = [0,  1 ] # degree
+gun.PhiMaxs   = [0,  1 ] # degree
+
 
 stdheprdr = StdHepRdr("StdHepRdr")
 #stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"

Generator/src/GtGunTool.cpp

@@ -15,12 +15,12 @@ GtGunTool::initialize() {
         error() << "Please specify the list of particle names/pdgs" << endmsg;
         return StatusCode::FAILURE;
     }
-    /*
+    
     if (m_energies.value().size() != m_particles.value().size()) {
         error() << "Mismatched energies and particles." << endmsg;
         return StatusCode::FAILURE;
     }
-    */
+    
     // others should be empty or specify
     if (m_thetamins.value().size()
         && m_thetamins.value().size() != m_particles.value().size()) {
@@ -79,9 +79,7 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
             }
         }
 
-        //double energy = m_energies.value()[i];
-        double energy_min = m_energies_min.value()[0];
-        double energy_max = m_energies_max.value()[0];
+        double energy = m_energies.value()[i];
 
         // create the MC particle
         edm4hep::MCParticle mcp = event.m_mc_vec.create();
@@ -95,51 +93,15 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
         // mcp.setEndpoint();
 
         // assume energy is momentum
-        double p = energy_min==energy_max ? energy_max : CLHEP::RandFlat::shoot(energy_min, energy_max);
-        double p1 = m_p1.value() ;
+        double p = energy;
         
         // direction
         // by default, randomize the direction
-        double theta = m_thetamins.value()[0]==m_thetamaxs.value()[0] ? m_thetamins.value()[0] : CLHEP::RandFlat::shoot(m_thetamins.value()[0], m_thetamaxs.value()[0]);
-        double phi =   m_phimins  .value()[0]==m_phimaxs  .value()[0] ? m_phimins  .value()[0] : CLHEP::RandFlat::shoot(m_phimins  .value()[0], m_phimaxs  .value()[0]);
-        double theta1 = m_dtheta.value() + theta;
-        double phi1   = m_dphi.value()   + phi  ;
+        double theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
+        double phi =   m_phimins  .value()[i]==m_phimaxs  .value()[i] ? m_phimins  .value()[i] : CLHEP::RandFlat::shoot(m_phimins  .value()[i], m_phimaxs  .value()[i]);
         double costheta = cos(theta*acos(-1)/180);
-        double costheta1 = cos(theta1*acos(-1)/180);
         double phi_  = phi*acos(-1)/180;
-        double phi1_ = phi1*acos(-1)/180;
         double sintheta = sqrt(1.-costheta*costheta);
-        double sintheta1 = sqrt(1.-costheta1*costheta1);
-        // check if theta min/max is set
-        /*
-        if (i < m_thetamins.value().size() 
-            && i < m_thetamaxs.value().size()) {
-            double thetamin = m_thetamins.value()[i];
-            double thetamax = m_thetamaxs.value()[i];
-
-            if (thetamin == thetamax) { // fixed theta
-                costheta = cos(thetamin);
-                sintheta = sin(thetamin);
-                info() << "theta is fixed: " << thetamin << endmsg;
-            }
-        }
-
-        if (i < m_phimins.value().size()
-            && i < m_phimaxs.value().size()) {
-            double phimin = m_phimins.value()[i];
-            double phimax = m_phimaxs.value()[i];
-
-            if (phimin == phimax) { // fixed phi
-                phi = phimin;
-                info() << "phi is fixed: " << phimin << endmsg;
-            }
-        }
-
-        debug() << "Direction: "
-                << " cos(theta): " << costheta
-                << " phi: " << phi
-                << endmsg;
-        */ 
         double px = p*sintheta*cos(phi_);
         double py = p*sintheta*sin(phi_);
         double pz = p*costheta;
@@ -148,20 +110,6 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
         // mcp.setMomentumAtEndpoint();
         // mcp.setSpin();
         // mcp.setColorFlow();
-        if(m_create_second)
-        {
-            edm4hep::MCParticle mcp1 = event.m_mc_vec.create();
-            mcp1.setPDG(pdgcode);
-            mcp1.setGeneratorStatus(1);
-            mcp1.setSimulatorStatus(1);
-            mcp1.setTime(0.0);
-            mcp1.setMass(mass);
-            double px1 = p1*sintheta1*cos(phi1_);
-            double py1 = p1*sintheta1*sin(phi1_);
-            double pz1 = p1*costheta1;
-            std::cout<<"GenGt p1="<<p1<<", px1="<<px1<<",py1="<<py1<<",pz1="<<pz1<<",theta1="<<theta1<<",phi1="<<phi1<<std::endl;
-            mcp1.setMomentum(edm4hep::Vector3f(px1,py1,pz1));
-        }
 
     }
 

Generator/src/GtGunTool.h

@@ -35,8 +35,6 @@ private:
     Gaudi::Property<std::vector<std::string>> m_particles{this, "Particles"};
 
     Gaudi::Property<std::vector<double>> m_energies{this, "Energies"};
-    Gaudi::Property<std::vector<double>> m_energies_min{this, "EnergyMin"};
-    Gaudi::Property<std::vector<double>> m_energies_max{this, "EnergyMax"};
 
     Gaudi::Property<std::vector<double>> m_thetamins{this, "ThetaMins"};
     Gaudi::Property<std::vector<double>> m_thetamaxs{this, "ThetaMaxs"};
@@ -44,10 +42,6 @@ private:
     Gaudi::Property<std::vector<double>> m_phimins{this, "PhiMins"};
     Gaudi::Property<std::vector<double>> m_phimaxs{this, "PhiMaxs"};
 
-    Gaudi::Property<bool> m_create_second{this, "Create_second", false};
-    Gaudi::Property< double > m_p1{this, "Energy1",1};
-    Gaudi::Property< double > m_dtheta{this, "dtheta",1};
-    Gaudi::Property< double > m_dphi  {this, "dphi"  ,1};
 };