added mass, pdgid, and charge to stopped particles

cms-sw · Mar 17, 2015 · ff8baa9 · ff8baa9
1 parent 5363a24
commit ff8baa9
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 9 deletions.
diff --git a/SimG4Core/CustomPhysics/plugins/RHStopDump.cc b/SimG4Core/CustomPhysics/plugins/RHStopDump.cc
@@ -18,14 +18,24 @@ RHStopDump::RHStopDump(edm::ParameterSet const & parameters)
    fEvent.getByLabel (mProducer, "StoppedParticlesY", ys);
    edm::Handle<std::vector<float> > zs;
    fEvent.getByLabel (mProducer, "StoppedParticlesZ", zs);
+   edm::Handle<std::vector<float> > ts;
+   fEvent.getByLabel (mProducer, "StoppedParticlesTime", ts);
+   edm::Handle<std::vector<int> > ids;
+   fEvent.getByLabel (mProducer, "StoppedParticlesPdgId", ids);
+   edm::Handle<std::vector<float> > masses;
+   fEvent.getByLabel (mProducer, "StoppedParticlesMass", masses);
+   edm::Handle<std::vector<float> > charges;
+   fEvent.getByLabel (mProducer, "StoppedParticlesCharge", charges);
+
    if (names->size() != xs->size() || xs->size() != ys->size() || ys->size() != zs->size()) {
      edm::LogError ("RHStopDump") << "mismatch array sizes name/x/y/z:"
 				  << names->size() << '/' << xs->size() << '/' << ys->size() << '/' << zs->size()
 				  << std::endl;
    }
    else {
      for (size_t i = 0; i < names->size(); ++i) {
-       mStream << (*names)[i] << ' ' << (*xs)[i] << ' ' << (*ys)[i] << ' ' << (*zs)[i] << std::endl;
+       mStream << (*names)[i] << ' ' << (*xs)[i] << ' ' << (*ys)[i] << ' ' << (*zs)[i] << ' ' << (*ts)[i] << std::endl;
+       mStream << (*ids)[i] << ' ' << (*masses)[i] << ' ' << (*charges)[i] << std::endl;
      }
    }
  }
diff --git a/SimG4Core/CustomPhysics/plugins/RHStopTracer.cc b/SimG4Core/CustomPhysics/plugins/RHStopTracer.cc
@@ -25,7 +25,10 @@ RHStopTracer::RHStopTracer(edm::ParameterSet const & p) {
   produces< std::vector<float> >("StoppedParticlesY");
   produces< std::vector<float> >("StoppedParticlesZ");
   produces< std::vector<float> >("StoppedParticlesTime");
-
+  produces< std::vector<int> >("StoppedParticlesPdgId");
+  produces< std::vector<float> >("StoppedParticlesMass");
+  produces< std::vector<float> >("StoppedParticlesCharge");
+
   LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::RHStopTracer->" 
 				     << mTraceParticleNameRegex << '/' << mTraceEnergy;
 }
@@ -34,7 +37,7 @@ RHStopTracer::~RHStopTracer() {
 }
 
 void RHStopTracer::update (const BeginOfRun * fRun) {
-  LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the run " << (*fRun)()->GetRunID();
+  LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the run " << (*fRun)()->GetRunID(); 
 }
 
 void RHStopTracer::update (const BeginOfEvent * fEvent) {
@@ -53,7 +56,7 @@ void RHStopTracer::update (const BeginOfTrack * fTrack) {
 				       << track->GetPosition().y() << '/' <<  track->GetPosition().z()
 				       << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
 				       << "    px/py/pz/p=" << track->GetMomentum().x() << '/' 
-				       << track->GetMomentum().y() << '/' << track->GetMomentum().z() << '/'<< track->GetMomentum().mag();
+				       << track->GetMomentum().y() << '/' << track->GetMomentum().z() << '/'<< track->GetMomentum().mag(); 
   }
   if (mStopRegular && !matched (track->GetDefinition()->GetParticleName())) { // kill regular particles
     const_cast<G4Track*>(track)->SetTrackStatus(fStopAndKill);
@@ -72,13 +75,16 @@ void RHStopTracer::update (const EndOfTrack * fTrack) {
 				       << track->GetPosition().y() << '/' <<  track->GetPosition().z()
 				       << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
 				       << "    px/py/pz/p=" << track->GetMomentum().x() << '/' 
-				       << track->GetMomentum().y() << '/' << track->GetMomentum().z() << '/'<< track->GetMomentum().mag();
+				       << track->GetMomentum().y() << '/' << track->GetMomentum().z() << '/'<< track->GetMomentum().mag(); 
     if (track->GetMomentum().mag () < 0.001) {
       mStopPoints.push_back (StopPoint (track->GetDefinition()->GetParticleName(),
 					track->GetPosition().x(),
 					track->GetPosition().y(),
 					track->GetPosition().z(),
-					track->GetGlobalTime()));
+					track->GetGlobalTime(),
+					track->GetDefinition()->GetPDGEncoding(),
+                                        track->GetDefinition()->GetPDGMass()/GeV,
+                                        track->GetDefinition()->GetPDGCharge() ));
     }
   }
 }
@@ -95,6 +101,9 @@ bool RHStopTracer::matched (const std::string& fName) const {
    std::auto_ptr<std::vector<float> > ys (new std::vector<float>);
    std::auto_ptr<std::vector<float> > zs (new std::vector<float>);
    std::auto_ptr<std::vector<float> > ts (new std::vector<float>);
+   std::auto_ptr<std::vector<int> > ids (new std::vector<int>);
+   std::auto_ptr<std::vector<float> > masses (new std::vector<float>);
+   std::auto_ptr<std::vector<float> > charges (new std::vector<float>);
 
    std::vector <StopPoint>::const_iterator stopPoint = mStopPoints.begin ();
    for (;  stopPoint != mStopPoints.end(); ++stopPoint) {
@@ -103,11 +112,17 @@ bool RHStopTracer::matched (const std::string& fName) const {
      ys->push_back (stopPoint->y);
      zs->push_back (stopPoint->z);
      ts->push_back (stopPoint->t);
+     ids->push_back (stopPoint->id);
+     masses->push_back (stopPoint->mass);
+     charges->push_back (stopPoint->charge);
    }
    fEvent.put (names, "StoppedParticlesName");
    fEvent.put (xs, "StoppedParticlesX");
    fEvent.put (ys, "StoppedParticlesY");
    fEvent.put (zs, "StoppedParticlesZ");
    fEvent.put (ts, "StoppedParticlesTime");
+   fEvent.put (ids, "StoppedParticlesPdgId");
+   fEvent.put (masses, "StoppedParticlesMass");
+   fEvent.put (charges, "StoppedParticlesCharge");
    mStopPoints.clear ();
  }
diff --git a/SimG4Core/CustomPhysics/plugins/RHStopTracer.h b/SimG4Core/CustomPhysics/plugins/RHStopTracer.h
@@ -28,16 +28,18 @@ class RHStopTracer :  public SimProducer,
   void produce(edm::Event&, const edm::EventSetup&);
  private:
   struct StopPoint {
-    StopPoint (const std::string& fName, double fX, double fY, double fZ, double fT) 
-      : name(fName), x(fX), y(fY), z(fZ), t(fT) 
+    StopPoint (const std::string& fName, double fX, double fY, double fZ, double fT, int fId, double fMass, double fCharge) 
+    : name(fName), x(fX), y(fY), z(fZ), t(fT), id(fId), mass(fMass), charge(fCharge) 
     {}
     std::string name;
     double x;
     double y;
     double z;
     double t;
+    int id;
+    double mass;
+    double charge;
   };
-  bool mDebug;
   bool mStopRegular;
   double mTraceEnergy;
   boost::regex mTraceParticleNameRegex;