Fri 2 Aug 22:46:21 BST 2019

7 files changed, 265 insertions, 43 deletions

diff --git a/src/#G4BeamTestSiHit.cxx# b/src/#G4BeamTestSiHit.cxx#
new file mode 100644
index 0000000..663a03f
--- /dev/null
+++ b/src/#G4BeamTestSiHit.cxx#
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <iomanip>
+
+#include "G4BeamTestSiHit.h"
+#include "G4UnitsTable.hh"
+#include "G4VVisManager.hh"
+#include "G4Circle.hh"
+#include "G4Colour.hh"
+#include "G4VisAttributes.hh"
+
+G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator=0;
+static std::fstream testnew("./testnew.txt", std::ofstream::out);
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSiHit::G4BeamTestSiHit()
+  : G4VHit(),
+    fTrackID(-1),
+  
+    fEdep(0.),
+              fTime(0.),
+    fPos(G4ThreeVector())
+{
+  /* G4cout << "opening file" << G4endl; */
+  /* testnew.open(testnew_out); */
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSiHit ::~G4BeamTestSiHit() {
+  /* G4cout << "closing file" << G4endl; */
+  /* testnew.close(); */
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4BeamTestSiHit::G4BeamTestSiHit(const G4BeamTestSiHit& right)
+  : G4VHit()
+{
+  fTrackID   = right.fTrackID;
+  // fChamberNb = right.fChamberNb;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+  fTime = right.fTime;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+const G4BeamTestSiHit& G4BeamTestSiHit::operator=(const G4BeamTestSiHit& right)
+{
+  fTrackID   = right.fTrackID;
+  //  fChamberNb = right.fChamberNb;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+  fTime = right.fTime;
+
+  return *this;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4int G4BeamTestSiHit::operator==(const G4BeamTestSiHit& right) const
+{
+  return ( this == &right ) ? 1 : 0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BeamTestSiHit::Draw()
+{
+  /* G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance(); */
+  /* if(pVVisManager) */
+  /*   { */
+  /*     G4Circle circle(fPos); */
+  /*     circle.SetScreenSize(4.); */
+  /*     circle.SetFillStyle(G4Circle::filled); */
+  /*     G4Colour colour(1.,0.,0.); */
+  /*     G4VisAttributes attribs(colour); */
+  /*     circle.SetVisAttributes(attribs); */
+  /*     pVVisManager->Draw(circle); */
+  /*   } */
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+void G4BeamTestSiHit::Print()
+{
+ G4cout
+    << " trackID: " << fTrackID
+    << "Edep: "
+    << std::setw(7) << G4BestUnit(fEdep,"Energy")
+    << " Position: "
+    << std::setw(7) << G4BestUnit(fPos ,"Length")
+    << "Time: "
+    << std::setw(7) << G4BestUnit(fTime,"Time")
+    << G4endl;
+}
+
+void G4BeamTestSiHit::Dataout()
+{
+testnew
+  << " trackID: " << fTrackID
+  << "Edep: "
+  << std::setw(7) << G4BestUnit(fEdep,"Energy")
+  << " Position: "
+  << std::setw(7) << G4BestUnit( fPos ,"Length")
+  << "Time: "
+  << std::setw(7) << G4BestUnit( fTime,"Time")
+  << G4endl;
+
+ }
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx
index ebbd9b8..495330a 100644
--- a/src/G4BeamTestDetectorConstruction.cxx
+++ b/src/G4BeamTestDetectorConstruction.cxx
@@ -3,6 +3,7 @@
 
 #include <G4LogicalVolume.hh>
 #include <G4PVPlacement.hh>
+#include <G4SDManager.hh>
 
 #include <G4Box.hh>
 
@@ -34,12 +35,17 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
   /* origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength); */
 
   // Determine World dimensions
-  G4double xWorld = 2.0 * CLHEP::m;
-  G4double yWorld = 2.0 * CLHEP::m;
-  G4double zWorld = 2.0 * CLHEP::m;
+  G4double xWorld = 4.0 * CLHEP::m;
+  G4double yWorld = 4.0 * CLHEP::m;
+  G4double zWorld = 4.0 * CLHEP::m;
+
+  // SC4 dimensions
+  G4double scinHeight_ = 1 * 2.54 * CLHEP::cm;
+  G4double scinWidth_ = 1 * 2.54 * CLHEP::cm;
+  G4double scinThickness_ = 1 * 2.54 * CLHEP::cm;
 
   // Create world volume
-  G4Box* world_box = new G4Box("solid_world", xWorld, yWorld, zWorld);
+  G4Box* world_box = new G4Box("solid_world", xWorld*0.5, yWorld*0.5, zWorld*0.5);
   G4LogicalVolume* worldLog =
       new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0);
   G4VPhysicalVolume* worldPhys =
@@ -56,15 +62,29 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct()
   tank_->InstallTank(worldPhys, origin_);
   // }
 
+  // Define SC4
+  G4Box* sc4_box = new G4Box("sc4",scinHeight_*0.5, scinWidth_*0.5, scinThickness_*0.5);
+  G4LogicalVolume* sc4Log =
+      new G4LogicalVolume(sc4_box, G4Material::GetMaterial("SC4"), "log_sc4", 0, 0, 0);
+  G4VPhysicalVolume* sc4Phys =
+      new G4PVPlacement(0, G4ThreeVector(1.2*CLHEP::m,0,0), sc4Log, "sc4", worldLog, false, 0);
+
   // User limits (energy cutoffs)
   // Do not create photons or electrons below cherenkov threshold
   // See also corresponding UserSpecialCuts in Physicslist !!!!
   G4UserLimits* energyLimit = new G4UserLimits();
   /* energyLimit->SetUserMinEkine(280.0 * CLHEP::keV);  // Cherenkov threshold of electrons in ice */
-  energyLimit->SetUserMinEkine(1.907 * CLHEP::eV);  // Lower threshold of PMT - 600nm
+  energyLimit->SetUserMinEkine(2.26 * CLHEP::eV);  // Lower threshold of PMT - 550nm
+  // energyLimit->SetUserMaxEkine(3.55 * CLHEP::eV); //upper threshold of PMT - 350nm
   worldLog->SetUserLimits(energyLimit);
   /* snowLog->SetUserLimits(energyLimit); */
+  sc4Log->SetUserLimits(energyLimit);
   
+  // G4SDManager* sdManager = G4SDManager::GetSDMpointer();
+  // sc4SD_ = new G4BeamTestSC4SD("sc4_SD_", "HitsCollection");
+  // sdManager->AddNewDetector(sc4SD_);
+  // sc4Log->SetSensitiveDetector(sc4SD_);
+
   return worldPhys;
 }
 
@@ -80,6 +100,7 @@ void G4BeamTestDetectorConstruction::CreateMaterials()
   /* CreatePerlite(); */
   CreateGlassSphere();
   CreateEffectiveDOMMaterial();
+  CreateSC4();
 
   //if(verboseLevel_>0) G4cout << *G4Material::GetMaterialTable() << G4endl;
 }
@@ -105,7 +126,8 @@ void G4BeamTestDetectorConstruction::CreateWater()
   // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste>
   // TODO(shivesh): add more properties?
   const G4int water_bins = 2;
-  G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
+  // G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
+  G4double water_ephot[water_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
   G4double water_refr[water_bins] = {1.33, 1.33};
   G4MaterialPropertiesTable *mpt_water = new G4MaterialPropertiesTable ();
   mpt_water->AddProperty ("RINDEX", water_ephot, water_refr, water_bins);
@@ -144,9 +166,40 @@ void G4BeamTestDetectorConstruction::CreateGlassSphere()
   // 20 lbs. weight = 9072 g
   // 6.5" outer radius & 0.5" thickness = 4024 cm3
   G4NistManager* nistManager = G4NistManager::Instance();
-  G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
-  glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
-  glass->AddElement(nistManager->FindOrBuildElement("O"),  2);
+  // G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
+  // glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
+  // glass->AddElement(nistManager->FindOrBuildElement("O"),  2);
+
+  // http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2013/03/11/12.39-85121-chDetectorConstruction.cc
+  // Define elements for all materials not found in the NIST database
+  G4Element* Si = nistManager->FindOrBuildElement("Si");
+  G4Element* B = nistManager->FindOrBuildElement("B");
+  G4Element* O = nistManager->FindOrBuildElement("O");
+  G4Element* Na = nistManager->FindOrBuildElement("Na");
+  G4Element* Al = nistManager->FindOrBuildElement("Al");
+  G4Element* K = nistManager->FindOrBuildElement("K");  
+
+  G4double density;
+  G4int ncomponents;
+  G4double fractionmass;
+  G4Material* glass = new G4Material("Glass", density= 2.23*CLHEP::g/CLHEP::cm3, ncomponents=6);
+  glass->AddElement(B, fractionmass=0.040064);
+  glass->AddElement(O, fractionmass=0.539562); 
+  glass->AddElement(Na, fractionmass=0.028191);
+  glass->AddElement(Al, fractionmass=0.011644);
+  glass->AddElement(Si, fractionmass=0.377220);
+  glass->AddElement(K, fractionmass=0.003321);
+
+  // pmt spectral response 300-650nm
+  // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste>
+  // TODO(shivesh): add more properties?
+  const G4int glass_bins = 2;
+  // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
+  G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
+  G4double glass_refr[glass_bins] = {1.47, 1.47};
+  G4MaterialPropertiesTable *mpt_glass = new G4MaterialPropertiesTable ();
+  mpt_glass->AddProperty ("RINDEX", glass_ephot, glass_refr, glass_bins);
+  glass->SetMaterialPropertiesTable(mpt_glass);
 }
 
 /*****************************************************************/
@@ -160,4 +213,23 @@ void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial()
   G4Material* glass = new G4Material("effectiveDOM", 0.2 * CLHEP::g / CLHEP::cm3, 2, kStateSolid);
   glass->AddElement(nistManager->FindOrBuildElement("Si"), 1);
   glass->AddElement(nistManager->FindOrBuildElement("O"),  2);
+
+  // const G4int glass_bins = 2;
+  // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV};
+  // glass 
+  // G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV};
+  // G4double glass_refr[glass_bins] = {1.47, 1.47};
+  // G4MaterialPropertiesTable *mpt_glass = new G4MaterialPropertiesTable ();
+  // mpt_glass->AddProperty ("RINDEX", glass_ephot, glass_refr, glass_bins);
+  // glass->SetMaterialPropertiesTable(mpt_glass);
+}
+
+void G4BeamTestDetectorConstruction::CreateSC4()
+{
+  G4NistManager* nistManager = G4NistManager::Instance();
+  // POM
+  G4Material* plastic = new G4Material("SC4", 1.425 * CLHEP::g / CLHEP::cm3, 3, kStateSolid);
+  plastic->AddElement(nistManager->FindOrBuildElement("H"), 2);
+  plastic->AddElement(nistManager->FindOrBuildElement("C"), 1);
+  plastic->AddElement(nistManager->FindOrBuildElement("O"), 1);
 }
diff --git a/src/G4BeamTestEventAction.cxx b/src/G4BeamTestEventAction.cxx
index 176991d..53402c9 100644
--- a/src/G4BeamTestEventAction.cxx
+++ b/src/G4BeamTestEventAction.cxx
@@ -33,13 +33,14 @@ void G4BeamTestEventAction::BeginOfEventAction(const G4Event* event )
  
   
   G4SDManager * SDman = G4SDManager::GetSDMpointer();
+  SDman->ListTree();
   if(SiCollID<0)
     {
       G4String colNam;
       SiCollID = SDman->GetCollectionID(colNam="G4BeamTestSiSDCollection");
      
     }
-  
+
 
 }
 
@@ -49,6 +50,9 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event)
 {   
 
   G4cout << ">>> Summary of Event " << event->GetEventID() << G4endl;
+  // testnew << ">>> Summary of Event " << event->GetEventID() << G4endl;
+
+  G4cout << SiCollID << G4endl;
 
   if(SiCollID<0) return;
 
@@ -63,17 +67,23 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event)
 
 
 
-
   if(SiHC)
     {
       int n_hit = SiHC->entries();
+      // testnew << std::flush;
       G4cout << G4endl;
-      G4cout << "Si hits " <<
-      "--------------------------------------------------------------"
-	     << G4endl;
+      // G4cout << "Si hits " <<
+      // "--------------------------------------------------------------"
+      //    << G4endl;
       G4cout << n_hit << " hits are stored in G4BeamTestSiHitsCollection."
 	     << G4endl;
       /* G4cout << "List of hits in tracker" << G4endl; */
+      // testnew << G4endl;
+      // testnew << "Si hits " <<
+      // "--------------------------------------------------------------"
+      testnew << n_hit << " hits are stored in G4BeamTestSiHitsCollection."
+	     << G4endl;
+      /* testnew << "List of hits in tracker" << G4endl; */
       for(int i=0;i<n_hit;i++)
       {
           /* (*SiHC)[i]->Print(); */
@@ -81,7 +91,9 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event)
       }
 
       //      G4cout << "sid + " << SiCollID << G4endl;
-   
+      //      testnew << "sid + " << SiCollID << G4endl;
+      testnew << std::flush;
+
     } 
 
 }
diff --git a/src/G4BeamTestSiHit.cxx b/src/G4BeamTestSiHit.cxx
index 591274c..0c90623 100644
--- a/src/G4BeamTestSiHit.cxx
+++ b/src/G4BeamTestSiHit.cxx
@@ -9,7 +9,6 @@
 #include "G4VisAttributes.hh"
 
 G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator=0;
-static std::fstream testnew("./testnew.txt", std::ofstream::out);
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -86,28 +85,30 @@ void G4BeamTestSiHit::Draw()
 
 void G4BeamTestSiHit::Print()
 {
- G4cout
-    << " trackID: " << fTrackID
-    << "Edep: "
-    << std::setw(7) << G4BestUnit(fEdep,"Energy")
-    << " Position: "
-    << std::setw(7) << G4BestUnit(fPos ,"Length")
-    << "Time: "
-    << std::setw(7) << G4BestUnit(fTime,"Time")
-    << G4endl;
+ // G4cout
+ //    << " trackID: " << fTrackID
+ //    << "Edep: "
+ //    << std::setw(7) << G4BestUnit(fEdep,"Energy")
+ //    << " Position: "
+ //    << std::setw(7) << G4BestUnit(fPos ,"Length")
+ //    << "Time: "
+ //    << std::setw(7) << G4BestUnit(fTime,"Time")
+ //    << G4endl;
 }
 
 void G4BeamTestSiHit::Dataout()
 {
-testnew
-  << " trackID: " << fTrackID
-  << "Edep: "
-  << std::setw(7) << G4BestUnit(fEdep,"Energy")
-  << " Position: "
-  << std::setw(7) << G4BestUnit( fPos ,"Length")
-  << "Time: "
-  << std::setw(7) << G4BestUnit( fTime,"Time")
-  << G4endl;
+// testnew << std::flush;
+// testnew
+//   << " trackID: " << fTrackID
+//   << "Edep: "
+//   << std::setw(7) << G4BestUnit(fEdep,"Energy")
+//   << " Position: "
+//   << std::setw(7) << G4BestUnit( fPos ,"Length")
+//   << "Time: "
+//   << std::setw(7) << G4BestUnit( fTime,"Time")
+//   << G4endl;
+// testnew << std::flush;
 
  }
 
diff --git a/src/G4BeamTestSiSD.cxx b/src/G4BeamTestSiSD.cxx
index b67092b..afa6bf8 100644
--- a/src/G4BeamTestSiSD.cxx
+++ b/src/G4BeamTestSiSD.cxx
@@ -45,6 +45,7 @@ void G4BeamTestSiSD::Initialize(G4HCofThisEvent* hce)
 
   G4int hcID 
     = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
+  G4cout << "hcID " << hcID << G4endl;
   hce->AddHitsCollection( hcID, fHitsCollection ); 
 }
 
@@ -61,11 +62,25 @@ G4bool G4BeamTestSiSD::ProcessHits(G4Step* aStep,
   /* G4cout << " Particle_name = " <<  name  << G4endl;   */
 
   if (name == "opticalphoton" || name == "gamma") {
-/*  G4cout << " Particle_name = " <<  name  << G4endl; */
+  // G4cout << " Particle_name = " <<  name  << G4endl;
 
+// total energy
+  G4double etot = aStep->GetTrack()->GetTotalEnergy();
 // energy deposit
   G4double edep = aStep->GetTotalEnergyDeposit();
 
+  if (etot < 2.26 * CLHEP::eV) { // Lower threshold of PMT - 550nm
+  // if (etot < 2.48 * CLHEP::eV) { // Lower threshold of PMT - 500nm
+      // G4cout << "particle " << name << " under threshold with energy " << etot << G4endl;
+      return true;
+  }
+  if (etot > 3.55 * CLHEP::eV) { // Upper threshold of PMT - 350nm
+  // if (etot > 3.10 * CLHEP::eV) { // Upper threshold of PMT - 400nm
+      // G4cout << "particle " << name << " over threshold with energy " << etot << G4endl;
+      return true;
+  }
+  // G4cout << "inserting particle " << name << " with energy " << etot << " into record" << G4endl;
+
 // if (edep==0.) return false;
 
   /* G4cout << " Particle_name_after_edep = " <<  name  << G4endl;   */
diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx
index 56570a4..02f6751 100644
--- a/src/G4BeamTestTank.cxx
+++ b/src/G4BeamTestTank.cxx
@@ -28,7 +28,8 @@
 G4BeamTestTank::G4BeamTestTank()
 {
   // Get tank dimensions
-  tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness
+  // tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness
+  tankThickness_ = 0.44 * 2.54 *CLHEP::cm; // TODO(shivesh) : check thickness
   tankHeight_  = 76.83 * 2.54 * CLHEP::cm;
   innerRadius_ = 32 * 2.54 * CLHEP::cm;
   outerRadius_ = innerRadius_ + tankThickness_;
@@ -56,7 +57,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
   // See also corresponding UserSpecialCuts in Physicslist !!!!
   // TODO(shivesh): Maybe do all of this as stepping action ??????
   G4UserLimits* energyLimit = new G4UserLimits();
-  energyLimit->SetUserMinEkine(1.907 * CLHEP::eV);  // Lower threshold of PMT - 600nm
+  energyLimit->SetUserMinEkine(2.26 * CLHEP::eV);  // Lower threshold of PMT - 550nm
 
   // std::string tankName=boost::lexical_cast<std::string>(tankKey_);
   std::string tankName = "BTT";
@@ -140,7 +141,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
       new G4LogicalVolume(upperglasssphere, glass,
                           ("log_dom_up_" + omName).c_str(), 0, 0, 0);
   G4LogicalVolume* logLowerGlass =
-      new G4LogicalVolume(lowerglasssphere, water,
+      new G4LogicalVolume(lowerglasssphere, glass,
                           ("log_dom_lo_" + omName).c_str(), 0, 0, 0);
   G4LogicalVolume* logUpperDOM =
       new G4LogicalVolume(upperdomsphere, effectiveDOM,
@@ -172,7 +173,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const
   G4SDManager* sdManager = G4SDManager::GetSDMpointer();
   iceSD_ = new G4BeamTestSiSD(("ice_SD_" + tankName).c_str(), "HitsCollection");
   sdManager->AddNewDetector(iceSD_);
-  /* logLowerDOM->SetSensitiveDetector(iceSD_); */
+  // logLowerDOM->SetSensitiveDetector(iceSD_);
   logLowerGlass->SetSensitiveDetector(iceSD_);
 
   // Instantiation of a set of visualization attributes with red colour
diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx
index f7c053b..e41ccc3 100644
--- a/src/G4BeamTestUserTrackingAction.cxx
+++ b/src/G4BeamTestUserTrackingAction.cxx
@@ -15,6 +15,7 @@ void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track)
   G4UserLimits *limit = volume->GetUserLimits();
   if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl;
   G4double threshold = limit->GetUserMinEkine(*track);
+  G4double max_threshold = 3.54;
   G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
   if(secondaries)
   {
@@ -29,10 +30,15 @@ void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track)
         {
           //check if particle energy is below threshold; if true, kill the particle
           G4double energy = (*secondaries)[i]->GetTotalEnergy();
-          if(energy < threshold){
-              G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
-              (*secondaries)[i]->SetTrackStatus(fStopAndKill);
-          }
+          // if(energy < threshold){
+          //     G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl;
+          //     (*secondaries)[i]->SetTrackStatus(fStopAndKill);
+          // }
+	  // if (energy >  max_threshold * CLHEP::eV){
+	      // G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " > " << max_threshold << G4endl;
+	      // (*secondaries)[i]->SetTrackStatus(fStopAndKill); 
+	    
+	  // }
         }
       }
     }