Kept checks for CADMesh. Added Global time to tree

include/DetectorConstruction.hh

@@ -37,6 +37,8 @@
 #include "G4Box.hh"
 #include "G4LogicalVolume.hh"
 #include "G4VUserDetectorConstruction.hh"
+#include "G4GDMLParser.hh"
+
 #include "Materials.hh"
 
 
@@ -68,10 +70,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction
     std::string GDMLoutput;
     G4double    fRoughness = 0;
 
-    #ifdef G4LIB_USE_GDML
     G4GDMLParser fParser;
-    #elif defined CADMesh
-    CADMesh* mesh;
 };
 
 #endif /*DetectorConstruction_h*/

include/PMTHit.hh

@@ -54,11 +54,13 @@ public:
   void          setHit        (G4ThreeVector xyz) { hit = xyz; }
   void          setMomentum   (G4ThreeVector mom) { momentum = mom; }
   void          setEnergy     (G4double e)        { energy = e; }
+  void          setTime       (G4double t)        { time = t; }
   G4int         getTrackID    ( )                 { return trackID; }
   G4ThreeVector getPos        ( )                 { return pos; }
   G4ThreeVector getHit        ( )                 { return hit; }
   G4double      getEnergy     ( )                 { return energy; }
   G4ThreeVector getMomentum   ( )                 { return momentum; }
+  G4double      getTime       ( )                 { return time; }
 
 private:
   G4int         trackID;
@@ -66,6 +68,7 @@ private:
   G4ThreeVector hit;
   G4double      energy;
   G4ThreeVector momentum;
+  G4double      time;
 };
 
 

lightGuide.cc

@@ -106,7 +106,6 @@ int main(int argc,char** argv)
 
 
 //----------------- Check options against dependencies -----------------//
-
   if(CADmodel != ""){
   // Check if the user has GDML support and if they are asking for it
     #ifndef G4LIB_USE_GDML
@@ -120,11 +119,12 @@ int main(int argc,char** argv)
           CADoutFile = "";
       }
     #endif
-    #ifndef CADMESH
-      if(filetype != "gdml"){
-        G4cout << "Cannot read non-gdml geometries without CADMESH installed!!! aborting" << G4endl;
-      }
-    #endif
+
+    //If you want to use CADMesh but don't have it installed
+    if(filetype == "stl" && !CADMESH){
+      G4cout << "Cannot read STL models without CADMESH installed!!! aborting" << G4endl;
+    }
+
   }
 
   // Instantiate G4UIExecutive if interactive mode

run1.mac

@@ -34,8 +34,8 @@
 /gps/pos/type Plane
 /gps/pos/shape Rectangle
 /gps/pos/centre 0. 0. 0. mm
-/gps/pos/halfx 44.875 mm
-/gps/pos/halfy 82 mm
+/gps/pos/halfx 41.875 mm
+/gps/pos/halfy 80 mm
 
 #
 # the incident surface is in the x-y plane
@@ -59,6 +59,6 @@
 ##############################################################
 ##############################################################
 # number of events
-/run/beamOn 1
+/run/beamOn 10000000
 ##############################################################
 ##############################################################

src/DetectorConstruction.cc

@@ -50,11 +50,8 @@
 #include "G4VisAttributes.hh"
 #include "G4Colour.hh"
 
-#ifdef G4LIB_USE_GDML
 #include "G4GDMLParser.hh"
-#endif
 
-// CADMESH //
 #ifdef CADMESH
 #include "CADMesh.hh"
 #endif
@@ -122,7 +119,9 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
 
   //----------------- Make the light guide -----------------//
   ////////////////////This got all sorts of fucked up. I'll probably drop GDML support or make it required
-  materials->AlSurface->SetSigmaAlpha(fRoughness);
+  //materials->AlSurface->SetFinish(ground);
+  //materials->AlSurface->SetPolish(0.1);
+  //materials->AlSurface->SetSigmaAlpha(2.0);
 
   //If we have defined a CAD file, use it
   if(filename != ""){
@@ -132,9 +131,8 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
     rot->rotateZ(90*deg);
     rot->rotateY(-90*deg);
 
-    //Import the cad model and place it
     #ifdef CADMESH
-    if(filetype != "gdml"){
+    if(filetype == "stl"){
       CADMesh* mesh = new CADMesh((char*) filename.c_str());
       mesh->SetScale(mm);
       mesh->SetOffset( G4ThreeVector(-20*cm, 0, 0) );
@@ -145,10 +143,10 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
                             Al,                      //material
                             "LightGuide");           //name
     }
-    #elif defined G4LIB_USE_GDML
+    #endif
+
     if(filetype == "gdml"){
-      fParser.Read(fReadFile);
-      break;
+      fParser.Read(filename);
     }
     if(filetype == "step"){
       cad_logical = fParser.ParseST(filename,Air,Al);
@@ -156,7 +154,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
 
     if(cad_logical !=0 ){
       cad_physical =
-        new G4PVPlacement(0, //rot,
+        new G4PVPlacement(rot,
                           G4ThreeVector(0,0,lgHeight-(3*cm)),
                           cad_logical,
                           "cad_physical",
@@ -165,13 +163,11 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
                           0);
     }
 
-    //If output to GDML if requested and supported
-    #ifdef G4LIB_USE_GDML
       if(GDMLoutput != ""){
         G4GDMLParser* gdml = new G4GDMLParser();
         gdml->Write("GDMLoutput",cad_physical);
       }
-    #endif
+
 
     //Make the surface optically reflective
     G4LogicalSkinSurface* alumLSS =

src/EventAction.cc

@@ -66,7 +66,6 @@ void EventAction::BeginOfEventAction(const G4Event* event){
  *
  */
 void EventAction::EndOfEventAction(const G4Event* event){
-
   hitsCollID = G4SDManager::GetSDMpointer()->GetCollectionID("MyPMT");
   G4HCofThisEvent* HCE = event->GetHCofThisEvent();
   if(HCE){
@@ -79,8 +78,8 @@ void EventAction::EndOfEventAction(const G4Event* event){
       analysisManager->FillNtupleDColumn(1, aHit->getPos().y() );
       analysisManager->FillNtupleDColumn(2, aHit->getHit().x() );
       analysisManager->FillNtupleDColumn(3, aHit->getHit().y() );
-      analysisManager->FillNtupleDColumn(4, aHit->getEnergy()  );
-      //Add global timeanalysisManager->FillNtupleDColumn(5, aHit->getTime()  );
+      //analysisManager->FillNtupleDColumn(4, aHit->getEnergy()  );
+      analysisManager->FillNtupleDColumn(4, aHit->getTime()  );
       analysisManager->AddNtupleRow();
     }
   }

src/Materials.cc

@@ -68,7 +68,7 @@ void Materials::DefineOpticalProperties(void){
   G4double quartz_RIND[nEntriesWLS], quartz_ABSL[nEntriesWLS],quartz_RFLT[nEntriesWLS],quartz_EFIC[nEntriesWLS], PhotonEnergy[nEntriesWLS];
      for(int i = 0; i < nEntriesWLS; i++){
          PhotonEnergy[i] = 2.00*eV + i*0.03*eV;
-         quartz_RIND[i] =tile_RI; //Refractive Index - constants
+         quartz_RIND[i] = 1.44221 + 5.9e-3*PhotonEnergy[i] + 6.71e-4*pow(PhotonEnergy[i],2) + 2.3e-5*pow(PhotonEnergy[i],3) + 2.745e-5*pow(PhotonEnergy[i],4); //Refractive Index
          quartz_ABSL[i] = 300.00*cm + i*20*cm; //Attenuation length
          quartz_RFLT[i] = 0.5;
          quartz_EFIC[i] = 0.5;

src/PMTHit.cc

@@ -54,6 +54,7 @@ PMTHit::PMTHit(const PMTHit& right)
   hit       = right.hit;
   energy    = right.energy;
   momentum  = right.momentum;
+  time      = right.time;
 }
 
 /*
@@ -66,6 +67,7 @@ const PMTHit& PMTHit::operator=(const PMTHit& right)
   hit       = right.hit;
   energy    = right.energy;
   momentum  = right.momentum;
+  time      = right.time;
   return *this;
 }
 

src/PMTSD.cc

@@ -86,11 +86,12 @@ G4bool PMTSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
   G4Track* theTrack = aStep->GetTrack();
   PMTHit* newHit = new PMTHit();
 
-  newHit->setTrackID  (aStep->GetTrack()->GetTrackID() );
-  newHit->setPos      (aStep->GetTrack()->GetVertexPosition() );
-  newHit->setHit      (aStep->GetTrack()->GetPosition() );
+  newHit->setTrackID  (theTrack->GetTrackID() );
+  newHit->setPos      (theTrack->GetVertexPosition() );
+  newHit->setHit      (theTrack->GetPosition() );
   newHit->setEnergy   (aStep->GetPreStepPoint()->GetTotalEnergy());
   newHit->setMomentum (aStep->GetPreStepPoint()->GetMomentum());
+  newHit->setTime     (theTrack->GetGlobalTime());
 
   hitCollection->insert( newHit );
 

src/RunAction.cc

@@ -65,7 +65,8 @@ RunAction::RunAction(G4String fileName)
   analysisManager->CreateNtupleDColumn("Y");
   analysisManager->CreateNtupleDColumn("hitX");
   analysisManager->CreateNtupleDColumn("hitY");
-  analysisManager->CreateNtupleDColumn("energy");
+  //analysisManager->CreateNtupleDColumn("energy");
+  analysisManager->CreateNtupleDColumn("time");
   analysisManager->FinishNtuple();
 }