18.2.2011 Kamil Sedlak

1) Correction in rundom nr generator (option 3)
     - still there is a bug in the printing of time, for which the
       run is running (shown is a large negative number) - should be
       fixed.
  2) Change in physics list which allows to set rise time in the scintillators

include/musrParameters.hh

@@ -34,7 +34,8 @@ class musrParameters {
     static musrParameters* GetInstance();
 
     static G4String mySteeringFileName;      // name of the steering file (e.g. the *.mac file)
-    static G4String myStopFileName;      // name of the stop file (e.g. the *.stop file), which will stop the run if it is created
+    static G4String myStopFileName;          // name of the stop file (e.g. the *.stop file), which will stop the run if it is created
+    static G4String myRandomNumberFileName;  // name of the file with random numbers for RandomOption=3 in musrDetectorMessenger
     static G4bool storeOnlyEventsWithHits;   // variable specifying whether to store interesting 
                                              // or all events into the ROOT tree. (default = true)
     static G4int  storeOnlyEventsWithHitInDetID; // simillar to "storeOnlyEventsWithHits".  The event is stored
@@ -58,6 +59,9 @@ class musrParameters {
                                              // the SteppinAction and reset to "false" in the GetFieldValue.
                                              // It is being changed on step by step basis.
     static G4int nrOfEventsToBeGenerated;    // Nr of events to be simulated in this run (set by /run/beamOn command)
+    static G4bool finiteRiseTimeInScintillator; // if true, rise time will be simulated in the scintillator.  For some strange
+                                             // reason, Geant4 requires that this is specifically allowed.  We set it true
+                                             // as soon as "FASTSCINTILLATIONRISETIME" or "SLOWSCINTILLATIONRISETIME" is set.
   private:
     static musrParameters* pointerToParameters;
     G4bool boolG4RegionRequested;  // variable used internally just to check that no new volume is defined after

include/musrPrimaryGeneratorAction.hh

@@ -79,7 +79,8 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     void SetTurtleInputFileToEventNo(G4int lineNumberOfTurtleFile);
     void SetTurtleZ0(G4double val)          {z0_InitialTurtle=val;}
     void SetTurtleInterpretAxes(G4String interpretAxes){turtleInterpretAxes=interpretAxes;}
-    void SetOrReadTheRandomNumberSeeds(G4int eventID);
+  //    void SetOrReadTheRandomNumberSeeds(G4int eventID);
+    void SetOrReadTheRandomNumberSeeds(G4Event* anEvent);
     void SetTurtleMomentumBite (G4ThreeVector smearingParam)
                                {turtleMomentumBite=true; turtleMomentumP0=smearingParam[0]*MeV; turtleSmearingFactor=smearingParam[1]*0.01;}
     void SetPrimaryParticule(G4String particleName);
@@ -134,7 +135,7 @@ public:
     static G4bool setRandomNrSeedFromFile_RNDM;
     static G4int  nRndmEventToSaveSeeds;
     static std::vector<int> * GetPointerToSeedVector();
-
+    static G4int lastEventID_in_pointerToSeedVector;
     G4double  decaytime;
 
 private:

src/musrDetectorConstruction.cc

@@ -1364,6 +1364,11 @@ void musrDetectorConstruction::DefineMaterials()
 	      double value;
 	      sscanf(&line[0],"%*s %*s %*s %*s %*d %lf",&value);
 	      MPT_tmp->AddConstProperty(propertyName,value);
+	      if ( (strcmp(propertyName,"FASTSCINTILLATIONRISETIME")==0) || (strcmp(propertyName,"SLOWSCINTILLATIONRISETIME")==0) ) {
+		if (value!=0) {
+		  musrParameters::finiteRiseTimeInScintillator = true;
+		}
+	      }
 	    }
 	    else {              // AddProperty
 	      char* pch = strstr(line,propertyName)+strlen(propertyName);

src/musrDetectorMessenger.cc

@@ -36,6 +36,7 @@
 #include <time.h>            //cks   -----------------------------||-------------------------------
 #include "musrEventAction.hh" // cks needed for setting the variable "nHowOftenToPrintEvent"
 #include "musrPrimaryGeneratorAction.hh" // cks needed for the initialisation of the random nr. generator by event nr.
+#include "musrParameters.hh"
 #include <vector>
 #include "globals.hh"
 
@@ -155,20 +156,30 @@ void musrDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
 	std::ifstream indata;
 	int num;
 
-	indata.open("randomNum.dat"); // opens the file
+	//	indata.open("randomNum.dat"); // opens the file
+	indata.open(musrParameters::myRandomNumberFileName);
 	if(!indata) { // file couldn't be opened
 	  G4cout << "Error: file could not be opened" << G4endl;
 	  exit(1);
 	}
 	std::vector<int> * seedVector = musrPrimaryGeneratorAction::GetPointerToSeedVector();
+	std::vector<int> tmpSeedVector;
 	indata >> num;
 	while ( !indata.eof() ) { // keep reading until end-of-file
-	  G4cout << "The next number is " << num << G4endl;
+	  //	  G4cout << "The next number is " << num << G4endl;
-	  seedVector->push_back(num);
+	  tmpSeedVector.push_back(num);
+	  //	  musrPrimaryGeneratorAction::lastEventID_in_pointerToSeedVector = num;
 	  indata >> num; // sets EOF flag if no value found
 	}
 	indata.close();
-	G4cout << "End-of-file reached.." << seedVector->size()<<G4endl;
+	// save last eventID for further use in musrPrimaryGeneratorAction.cc
+	musrPrimaryGeneratorAction::lastEventID_in_pointerToSeedVector = tmpSeedVector.back();
+	// save the eventIDs in the reversed order so that they can be popped-back (using pop_back)
+	std::vector<int>::reverse_iterator rit;
+	for ( rit=tmpSeedVector.rbegin() ; rit < tmpSeedVector.rend(); ++rit ) {
+	  seedVector->push_back(*rit);
+	}
+	tmpSeedVector.clear();
       }
       else if (RandomOption == 4) {
 	G4cout << "*********************************************" << G4endl;

src/musrParameters.cc

@@ -30,12 +30,10 @@ musrParameters::musrParameters(G4String steeringFileName)
   pointerToParameters = this;
   boolG4RegionRequested = false;
   mySteeringFileName = steeringFileName;
-  myStopFileName = mySteeringFileName;  
+  myStopFileName = mySteeringFileName;
-  //  myStopFileName.replace(myStopFileName.end()-2,myStopFileName.end(),"stop");
-  //  Int_t lengthOfFileName=
-  //  myStopFileName.replace(3,5,"stop");
   myStopFileName.replace(myStopFileName.length()-3,myStopFileName.length()-1,"stop");
-
+  myRandomNumberFileName = mySteeringFileName;
+  myRandomNumberFileName.replace(myRandomNumberFileName.length()-3,myRandomNumberFileName.length()-1,"rndm");
 
   // Read in the parameters, which have to be known before the detector construction is run
   // (and therefore the parameters can not be read in in the musrDetectorConstruction.cc class).
@@ -106,6 +104,7 @@ musrParameters* musrParameters::GetInstance() {
 
 G4String musrParameters::mySteeringFileName="Unset";
 G4String musrParameters::myStopFileName="Unsetblablabla034tdk40928jfmfnakfh921djf02UNSET";
+G4String musrParameters::myRandomNumberFileName="Unsetblablabla034tdk40928jfmfnakfh921djf02UNSED";
 G4bool   musrParameters::storeOnlyEventsWithHits=true;
 G4int    musrParameters::storeOnlyEventsWithHitInDetID=0;
 G4double musrParameters::signalSeparationTime=100*nanosecond;
@@ -119,3 +118,4 @@ G4bool   musrParameters::boolG4OpticalPhotons=false;
 //cks G4bool   musrParameters::includeMuoniumProcesses =true; // TS
 //G4bool   musrParameters::boolG4GeneralParticleSource=true;
 G4int    musrParameters::nrOfEventsToBeGenerated=0;
+G4bool   musrParameters::finiteRiseTimeInScintillator=false;

src/musrPhysicsList.cc

@@ -368,7 +368,11 @@ void musrPhysicsList::ConstructEM()
 	  //	  else if (stringProcessName=="G4hIonisation")        pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
 	  //	  else if (stringProcessName=="G4hLowEnergyIonisation")      pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
 	  else if (stringProcessName=="musrMuFormation")              pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
-	  else if (stringProcessName=="G4Scintillation")              pManager->AddProcess(new G4Scintillation,nr1,nr2,nr3);
+	  else if (stringProcessName=="G4Scintillation")              {
+        	      G4Scintillation* myScintillation = new G4Scintillation();
+		      if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true);
+	              pManager->AddProcess(myScintillation,nr1,nr2,nr3);
+	  }
 	  // cks:  musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
 	  //       in his original "musrPhysicsList.cc about implementing musrMuScatter.
 	  //	  else if (stringProcessName=="musrMuScatter")       pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);

src/musrPrimaryGeneratorAction.cc

@@ -46,6 +46,7 @@ G4bool musrPrimaryGeneratorAction::setRandomNrSeedAccordingEventNr=0;
 G4bool musrPrimaryGeneratorAction::setRandomNrSeedFromFile=0;
 G4bool musrPrimaryGeneratorAction::setRandomNrSeedFromFile_RNDM=0;
 G4int  musrPrimaryGeneratorAction::nRndmEventToSaveSeeds=-2;
+G4int  musrPrimaryGeneratorAction::lastEventID_in_pointerToSeedVector=0;
 
 std::vector<int> * musrPrimaryGeneratorAction::pointerToSeedVector=NULL;
 std::vector<int> * musrPrimaryGeneratorAction::GetPointerToSeedVector() {
@@ -127,7 +128,7 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 
   // Set or read the seeds of random number generator 
   boolPrintInfoAboutGeneratedParticles=false;
-  SetOrReadTheRandomNumberSeeds(anEvent->GetEventID());
+  SetOrReadTheRandomNumberSeeds(anEvent);
 
   // If radioactive source is used, use  G4GeneralParticleSource :
   if (musrParameters::boolG4GeneralParticleSource) {
@@ -462,7 +463,8 @@ void musrPrimaryGeneratorAction::SetTurtleInputFileToEventNo(G4int lineNumberOfT
   }
 }
 //===============================================================================
-void musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds(G4int eventID) {
+void musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds(G4Event* anEvent) {
+  G4int eventID = anEvent->GetEventID();
   if (eventID == nRndmEventToSaveSeeds) {
     G4cout<<"musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds:  S A V I N G   R A N D O M    N O.   S E E D S"<<G4endl;
     G4cout<<"                                                                     (for even nr. "<<eventID<<")"<<G4endl;
@@ -477,13 +479,27 @@ void musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds(G4int eventID) {
     }
   }  
   if (setRandomNrSeedFromFile) {
-    //    G4cout<<"RandomNrInitialisers.size()="<<RandomNrInitialisers->size()<<G4endl;
+    //    //    G4cout<<"RandomNrInitialisers.size()="<<RandomNrInitialisers->size()<<G4endl;
-    if (eventID < int(pointerToSeedVector->size())) {
+    //    if (eventID < int(pointerToSeedVector->size())) {
-      G4cout <<"musrEventAction.cc: seed will be set to="<< pointerToSeedVector->at(eventID)<<G4endl;
+    //    //    if (eventID < lastEventID_in_pointerToSeedVector)
-      CLHEP::HepRandom::setTheSeed(pointerToSeedVector->at(eventID));
+    //      G4cout <<"musrEventAction.cc: seed will be set to="<< pointerToSeedVector->at(eventID)<<G4endl;
-      CLHEP::RandGauss::setFlag(false);
+    //      anEvent -> SetEventID(pointerToSeedVector->at(eventID));
-      boolPrintInfoAboutGeneratedParticles = true;
+    //      //      CLHEP::HepRandom::setTheSeed(pointerToSeedVector->at(eventID));
+    //      CLHEP::HepRandom::setTheSeed(eventID);
+    //      CLHEP::RandGauss::setFlag(false);
+    //      boolPrintInfoAboutGeneratedParticles = true;
+    //    }
+    if (!pointerToSeedVector->empty()) {
+      eventID = pointerToSeedVector->back();
+      pointerToSeedVector->pop_back();
+      anEvent -> SetEventID(eventID);
     }
+    else {
+      eventID = ++lastEventID_in_pointerToSeedVector;
+    }
+    CLHEP::HepRandom::setTheSeed(eventID);
+    CLHEP::RandGauss::setFlag(false);
+    //    G4cout <<"musrPrimaryGeneratorAction.cc: seed will be set to="<< eventID<<G4endl;
   }
   else if (setRandomNrSeedAccordingEventNr) {
     //    long seeds[2];

src/musrScintSD.cc

@@ -653,7 +653,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
 	  
 	  //	    OPSAhisto->Fit(poiss,"Q");
 	  if (boolStoreThisOPSAhist) OPSAhisto->Write();
-	  OPSA_timeMean = OPSAhisto->GetMean();
+	  OPSA_timeMean = OPSAhisto->GetMean() + OPSA_timeFirst;
 
 	  // if required, convert the histogram with photon times to histogram of electronic pulse shape
 	  if (bool_pulseShapeExists) {