9.5.2012 Kamil Sedlak

Several smaller improvements, should not have impact on the simulation.

doc/musrSim.tex

@@ -953,6 +953,12 @@ in order to simulate optical photons:
 		of the optical detection device as e.g.\ G-APD, and 
 		``OPSA\_CFD\_\emph{eventID}\_\emph{detID}\_\emph{n}'', which shows
 		the signal from OPSAshape after the constant fraction discriminator.
+
+		There are two special cases:\\
+		\emph{eventID} = -1 ... OPSA histograms for all events will be summed into 
+		                        a signle histogram.\\
+		\emph{eventID} = -2 ... OPSA histograms will be created for all events
+		                        (all histograms saved individually).
 	    \item{\bf pulseShapeArray \emph{pulseShapeFileName}}\\
 	        \emph{pulseShapeFileName} is the name of the file that contains response 
 		function (pulse shape) of a single cell (single photon) detected by 
@@ -1067,6 +1073,9 @@ in order to simulate optical photons:
 	caused by the decay positron.  This command would be
 	harmful in most physics studies.
 
+\item{\bf /musr/command killAllElectrons true}\\
+        See ``/musr/command killAllPositrons true'' for the explanation.
+
 \item{\bf /musr/command killAllGammas true}\\
 	See ``/musr/command killAllPositrons true'' for the explanation.
 
@@ -1098,6 +1107,17 @@ in order to simulate optical photons:
 	possible to set different \emph{timeSeparation} for a slow
 	and fast scintillator detectors of the instrument.
 
+\item{\bf /musr/command maximumNrOfStepsPerTrack \emph{n}}\\
+        Maximum number of steps allowed for a single track.
+	If the number of steps per track exceeds \emph{n}, the track is killed,
+	and the simulation proceeds with the next track.
+
+\item{\bf /musr/command maximumTimePerEvent \emph{time}}\\
+        Maximum time (in seconds) allowed for the simulation of a single event.
+	If this time is exceeded, the event is killed and the simulation
+        of the next event starts.  The default value is 60 seconds.
+	The variable \emph{time} is integer!
+
 \item{\bf /musr/command maximumRunTimeAllowed \emph{timeMax}}\\
 	If a musrSim job is run on a pc farm with a time limit on the job execution,
 	and the job exceeds the time limit, the simulation will be killed.

include/musrParameters.hh

@@ -46,6 +46,8 @@ class musrParameters {
     static G4bool storeOnlyTheFirstTimeHit;  // if true, only the hit that happened first will be
                                              // stored, anything else will be ignored 
                                              // (usefull for some special studies, not for a serious simulation)
+    static G4int storeOnlyEventsWithMuonsDecayedInDetID;
+    static G4bool killAllElectrons;          // if true, all electron tracks will be deleted (usefull for the studies of the muon beam)
     static G4bool killAllPositrons;          // if true, all positron tracks will be deleted (usefull for the studies of the muon beam)
     static G4bool killAllGammas;             // 
     static G4bool killAllNeutrinos;          //
@@ -68,6 +70,9 @@ class musrParameters {
     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.
+    static G4int maximumTimePerEvent;        // maximum time (in seconds) allowed for an event simulation - if exceeded, kill the event and
+                                             // proceed to the next one.
+    static G4int maximumNrOfStepsPerTrack;   // Maximum number of steps per track - if exceeded, kill the track and proceed with a next one.
   private:
     static musrParameters* pointerToParameters;
     G4bool boolG4RegionRequested;  // variable used internally just to check that no new volume is defined after

include/musrRootOutput.hh

@@ -62,6 +62,7 @@ class musrRootOutput  {
     void SetEventID        (G4int id) {eventID = id;};
     void SetTimeToNextEvent(G4double deltaT) {timeToNextEvent = deltaT/microsecond;}
     void SetDecayDetectorID (std::string detectorName) {muDecayDetID = SensDetectorMapping[detectorName];};
+    Int_t GetDecayDetectorID() {return  muDecayDetID;};
     void SetBField       (G4double F[6]) {B_t[0]=F[0]/tesla; B_t[1]=F[1]/tesla; B_t[2]=F[2]/tesla; 
                                           B_t[3]=F[3]/tesla; B_t[4]=F[4]/tesla; B_t[5]=F[5]/tesla;};
     void SetDecayPolarisation (G4ThreeVector pol) {muDecayPolX=pol.x(); muDecayPolY=pol.y(); muDecayPolZ=pol.z();};

include/musrScintSD.hh

@@ -155,6 +155,7 @@ class musrScintSD : public G4VSensitiveDetector
     G4int    myStoreOnlyEventsWithHitInDetID;
     G4double mySignalSeparationTime;
     G4bool   myStoreOnlyTheFirstTimeHit;
+    G4int    myStoreOnlyEventsWithMuonsDecayedInDetID;
     G4bool   boolIsVvvInfoRequested;
     musrRootOutput* myRootOutput;
 

musrSimAna/musrAnalysis.cxx

@@ -439,6 +439,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
 	else if (strcmp(conditionNameTMP,"promptPeakD")==0)             conditionMap[iConditionTMP]=&promptPeakD;
 	else if (strcmp(conditionNameTMP,"promptPeakL")==0)             conditionMap[iConditionTMP]=&promptPeakL;
 	else if (strcmp(conditionNameTMP,"promptPeakR")==0)             conditionMap[iConditionTMP]=&promptPeakR;
+	else if (strcmp(conditionNameTMP,"doubleHit")==0)             conditionMap[iConditionTMP]=&doubleHit;
 	else {
 	  std::cout<<" !!!  ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx   S T O P !!!"<<std::endl;
 	  exit(1);
@@ -1052,6 +1053,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
   Long64_t timeBin0                    = -100000000;
   Long64_t timeBin1                    = -100000000;
   Long64_t timeBin2                    = -100000000;
+  Long64_t timeBinDoubleHit;
   Int_t    kEntry                      = 0;
   Int_t    posEntry                      = 0;
   Int_t    idetM                       = 0;
@@ -1095,18 +1097,32 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
     pileup_muDecayPosY         = -1000000000;
     pileup_muDecayPosZ         = -1000000000;
     pileup_muDecayPosR         = -1000000000;
+    det_time31                 = -1001.;
+    timeBinDoubleHit           = -100000000;
+    doubleHit                  = false;
+    pos_detID                  = -1.;
+    pos_detID_doubleHit        = -1.;
+    Int_t pos_detID_doubleHit_INT    = -1;
     if (mCounterHitExistsForThisEventID) {
       numberOfGoodMuons++;
       Long64_t dataBinMin = timeBin0+dataWindowBinMin;
       Long64_t dataBinMax = timeBin0+dataWindowBinMax;
       //      Long64_t positronBinMax = timeBin0+pileupWindowBinMax+dataWindowBinMin;
       Long64_t positronBinMax = dataBinMax-dataWindowBinMin;  // note that "dataWindowBinMin" is normally negative, i.e. positronBinMax > dataBinMax !!!
-      pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,positronBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep);
+      pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,positronBinMax,timeBin1,timeBin2,timeBinDoubleHit,posEntry,idetP,idetP_ID,idetP_edep,pos_detID_doubleHit_INT);
+      pos_detID = Double_t(idetP_ID);
+      pos_detID_doubleHit = Double_t(pos_detID_doubleHit_INT);
       if (debugEventMap[eventID]>2) {
 	if (pCounterHitExistsForThisEventID) {std::cout<<"FillHistograms: GOOD positron candidate found: timeBin1="<<timeBin1
 						       <<",  deltat+800="<<timeBin1-timeBin0+800 <<std::endl;}
 	else {std::cout<<"FillHistograms: NO positron candidate found"<<std::endl;}
       }
+      if ((!pCounterHitExistsForThisEventID)&&(pos_detID_doubleHit_INT!=-1)&&mCounterHitExistsForThisEventID) {
+	doubleHit = true;
+	det_time31 = (timeBinDoubleHit-timeBin1)*tdcresolution;
+	if (pos_detID_doubleHit==pos_detID) {std::cout<<"DOUBLEHIT: pos_detID_doubleHit="<<pos_detID_doubleHit<<"; pos_detID="<<pos_detID<<std::endl;}
+
+      }
       //cDEL      if (pCounterHitExistsForThisEventID) std::cout<<"        timeBin1-timeBin2 = "<<timeBin1<<"-"<<timeBin2<<"="<<timeBin1-timeBin2 <<std::endl;
       //      
       if (pCounterHitExistsForThisEventID&&(posEntry>0)) {
@@ -1405,7 +1421,7 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
 //}
 //
 //================================================================
-Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep) {
+Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Long64_t& tBinDoubleHit, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Int_t& idetP_ID_doubleHit) {
 
   if (bool_debugingRequired) {
     if (debugEventMap[eventID]>4) {std::cout<<"PositronCounterHit:  pCounterMap.size()="<<pCounterMap.size()<<std::endl;}
@@ -1415,6 +1431,8 @@ Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_
   //  Bool_t positronHitFound  = false;
   Bool_t goodPositronFound = false;
   Int_t positronQuality = 0;
+  //  Int_t prev_counterID = 0;
+  //  Long64_t prev_time = 0;
   //  std::cout<<"Debug 10------------------------------"<<std::endl;
   if (musrMode=='D') {
     // Loop over all positron counters
@@ -1428,9 +1446,15 @@ Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_
       if (positronQuality==2) {
 	if (goodPositronFound) {
 	  if (debugEventMap[eventID]>3) {std::cout<<"PositronCounterHit:  positron candidate killed - double hit in a different counter"<<std::endl;}
+	  //	  std::cout<<"("<<evID<<") Double hit in counters "<<idetP_ID<<", "<<prev_counterID<<";  time bins ="<<tBin1<<", "<<prev_time<<std::endl;
+	  //	  if (idetP_ID_doubleHit == idetP_ID) std::cout<<"DOUBLEHIT  idetP_ID_doubleHit = "<<idetP_ID_doubleHit<<";  idetP_ID ="<< idetP_ID<<std::endl;
 	  return false;  // double hit was found in a different counter
 	}
 	goodPositronFound = true;
+	tBinDoubleHit = tBin1;
+	idetP_ID_doubleHit = idetP_ID;
+	//	prev_counterID = idetP_ID;
+	//	prev_time = tBin1;
       }
     }
 
@@ -1440,6 +1464,7 @@ Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_
     if (goodPositronFound&&(tBin1<dataBinMax))    return true;
   }
   //  std::cout<<"Debug 110"<<std::endl;
+  idetP_ID_doubleHit = -1;
   return false;
 }
 

musrSimAna/musrAnalysis.hh

@@ -217,6 +217,8 @@ public :
    Double_t        pos_Theta_MINUS_muDecayPol_Theta360;
    Double_t        pos_Phi_MINUS_muDecayPol_Phi;
    Double_t        pos_Phi_MINUS_muDecayPol_Phi360;
+   Double_t        pos_detID;
+   Double_t        pos_detID_doubleHit;
   //   Double_t        det_time0;
   //   Double_t        get_time0;
   //   Double_t        det_time1;
@@ -225,6 +227,7 @@ public :
    Double_t        gen_time10;
    Double_t        det_time10_MINUS_gen_time10;
    Double_t        det_time20;
+   Double_t        det_time31;
    Double_t        det_time1_MINUS_muDecayTime;
    Double_t        detP_x;
    Double_t        detP_y;
@@ -279,6 +282,7 @@ public :
    Bool_t          promptPeakD;
    Bool_t          promptPeakL;
    Bool_t          promptPeakR;
+   Bool_t          doubleHit;
 
    musrAnalysis(TTree *tree=0);
    virtual ~musrAnalysis();
@@ -301,7 +305,7 @@ public :
    virtual void     PrintHitsInAllCounters();
    virtual void     InitialiseEvent();
    virtual Double_t PreprocessEvent(Long64_t iEn);
-   virtual Bool_t   PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep);
+   virtual Bool_t   PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Long64_t& tBinDoubleHit, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Int_t& idetP_ID_doubleHit);
   //   virtual Bool_t   MuonCounterHit(Int_t evID, Long64_t timeBinMin, Long64_t& timeBin0, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep);
    void CopySubstring(char* inputChar,int iStart,int iEnd,char* outputChar);
    void MyPrintTree();
@@ -510,6 +514,8 @@ musrAnalysis::musrAnalysis(TTree *tree)
   variableMap["pos_Theta_MINUS_muDecayPol_Theta360"]=&pos_Theta_MINUS_muDecayPol_Theta360;
   variableMap["pos_Phi_MINUS_muDecayPol_Phi"]=&pos_Phi_MINUS_muDecayPol_Phi;
   variableMap["pos_Phi_MINUS_muDecayPol_Phi360"]=&pos_Phi_MINUS_muDecayPol_Phi360;
+  variableMap["pos_detID"]=&pos_detID;
+  variableMap["pos_detID_doubleHit"]=&pos_detID_doubleHit;
   //  variableMap["det_time0"]=&det_time0;
   //  variableMap["gen_time0"]=&gen_time0;
   //  variableMap["det_time1"]=&det_time1;
@@ -537,6 +543,7 @@ musrAnalysis::musrAnalysis(TTree *tree)
   variableMap["pileup_muDecayPosZ"]=&pileup_muDecayPosZ;
   variableMap["pileup_muDecayPosR"]=&pileup_muDecayPosR;
   variableMap["det_time20"]=&det_time20;
+  variableMap["det_time31"]=&det_time31;
 
   testIVar1=0;
   humanDecayHistograms=NULL;

src/musrDetectorConstruction.cc

@@ -1201,6 +1201,19 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	}
       } 
 
+      else if (strcmp(tmpString1,"storeOnlyEventsWithMuonsDecayedInDetID")==0){
+	G4int variable;
+	sscanf(&line[0],"%*s %*s %d",&variable);
+	if (variable!=0){ 
+	  musrParameters::storeOnlyEventsWithMuonsDecayedInDetID = variable;
+	  char eMessage[200];
+	  sprintf(eMessage,
+		  "musrDetectorConstruction.cc:: Only the events, in which muon stop in the detector ID=%d, are stored",
+		  variable);
+	  musrErrorMessage::GetInstance()->musrError(INFO,eMessage,false);
+	}
+      } 
+
       else if (strcmp(tmpString1,"storeOnlyEventsWithHits")==0){
 	if (strcmp(tmpString2,"false")==0){ musrParameters::storeOnlyEventsWithHits = false; }
       } 
@@ -1209,6 +1222,11 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	if (strcmp(tmpString2,"true")==0){ musrParameters::storeOnlyTheFirstTimeHit = true; }
       } 
 
+      else if (strcmp(tmpString1,"killAllElectrons")==0){
+	if (strcmp(tmpString2,"true")==0){ musrParameters::killAllElectrons = true; }
+	else                             { musrParameters::killAllElectrons = false; }
+      } 
+      
       else if (strcmp(tmpString1,"killAllPositrons")==0){
 	if (strcmp(tmpString2,"true")==0){ musrParameters::killAllPositrons = true; }
 	else                             { musrParameters::killAllPositrons = false; }
@@ -1224,6 +1242,22 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	else                             { musrParameters::killAllNeutrinos = false; }
       } 
       
+      else if (strcmp(tmpString1,"maximumTimePerEvent")==0){
+	int variable;
+	sscanf(&line[0],"%*s %*s %d",&variable);
+	if (variable>0){ 
+	  musrParameters::maximumTimePerEvent = variable;
+	}
+      }
+
+      else if (strcmp(tmpString1,"maximumNrOfStepsPerTrack")==0){
+	int variable;
+	sscanf(&line[0],"%*s %*s %d",&variable);
+	if (variable>0){ 
+	  musrParameters::maximumNrOfStepsPerTrack = variable;
+	}
+      }
+
       else if (strcmp(tmpString1,"getDetectorMass")==0){
 	G4LogicalVolume* massVol = FindLogicalVolume(tmpString2);
 	if (massVol==NULL) {
@@ -1581,7 +1615,26 @@ void musrDetectorConstruction::DefineMaterials()
   
   new G4Material("ArgonGas", z= 18., a= 39.95*g/mole, density= 0.00000000001*mg/cm3);
 
-  new G4Material("HeliumGas5mbar",z=2., a=4.002602*g/mole, density= 0.00000088132*g/cm3);
+  new G4Material("HeliumGas5mbar", z=2., a=4.002602*g/mole, density= 0.00000088132*g/cm3);
+  new G4Material("HeliumGas6mbar", z=2., a=4.002602*g/mole, density= 0.000001057584*g/cm3);
+  new G4Material("HeliumGas7mbar", z=2., a=4.002602*g/mole, density= 0.000001233848*g/cm3);
+  new G4Material("HeliumGas8mbar", z=2., a=4.002602*g/mole, density= 0.000001410112*g/cm3);
+  new G4Material("HeliumGas9mbar", z=2., a=4.002602*g/mole, density= 0.000001586376*g/cm3);
+  new G4Material("HeliumGas10mbar",z=2., a=4.002602*g/mole, density= 0.00000176264*g/cm3);
+  new G4Material("HeliumGas11mbar",z=2., a=4.002602*g/mole, density= 0.000001938904*g/cm3);
+  new G4Material("HeliumGas12mbar",z=2., a=4.002602*g/mole, density= 0.000002115168*g/cm3);
+  new G4Material("HeliumGas13mbar",z=2., a=4.002602*g/mole, density= 0.000002291432*g/cm3);
+  new G4Material("HeliumGas14mbar",z=2., a=4.002602*g/mole, density= 0.000002467696*g/cm3);
+  new G4Material("HeliumGas15mbar",z=2., a=4.002602*g/mole, density= 0.00000264396*g/cm3);
+  new G4Material("HeliumGas20mbar",z=2., a=4.002602*g/mole, density= 0.00000352528*g/cm3);
+  new G4Material("HeliumGas30mbar",z=2., a=4.002602*g/mole, density= 0.00000528792*g/cm3);
+  new G4Material("HeliumGas40mbar",z=2., a=4.002602*g/mole, density= 0.00000705056*g/cm3);
+  new G4Material("HeliumGas50mbar",z=2., a=4.002602*g/mole, density= 0.00000881320*g/cm3);
+  new G4Material("HeliumGas60mbar",z=2., a=4.002602*g/mole, density= 0.00001057584*g/cm3);
+  new G4Material("HeliumGas70mbar",z=2., a=4.002602*g/mole, density= 0.00001233848*g/cm3);
+  new G4Material("HeliumGas80mbar",z=2., a=4.002602*g/mole, density= 0.00001410112*g/cm3);
+  new G4Material("HeliumGas90mbar",z=2., a=4.002602*g/mole, density= 0.00001586376*g/cm3);
+  new G4Material("HeliumGas100mbar",z=2.,a=4.002602*g/mole, density= 0.00001762640*g/cm3);
 
   if (musrParameters::boolG4OpticalPhotons) {
 

src/musrParameters.cc

@@ -112,7 +112,9 @@ G4bool   musrParameters::storeOnlyEventsWithHits=true;
 G4int    musrParameters::storeOnlyEventsWithHitInDetID=0;
 G4double musrParameters::signalSeparationTime=100*nanosecond;
 G4bool   musrParameters::storeOnlyTheFirstTimeHit=false;
+G4int    musrParameters::storeOnlyEventsWithMuonsDecayedInDetID=0;
 G4bool   musrParameters::field_DecayWithSpin=false;
+G4bool   musrParameters::killAllElectrons=false;
 G4bool   musrParameters::killAllPositrons=false;
 G4bool   musrParameters::killAllGammas=false;
 G4bool   musrParameters::killAllNeutrinos=true;
@@ -123,3 +125,5 @@ G4bool   musrParameters::boolG4OpticalPhotonsUnprocess=false;
 //G4bool   musrParameters::boolG4GeneralParticleSource=true;
 G4int    musrParameters::nrOfEventsToBeGenerated=0;
 G4bool   musrParameters::finiteRiseTimeInScintillator=false;
+G4int    musrParameters::maximumTimePerEvent=60;
+G4int    musrParameters::maximumNrOfStepsPerTrack=100000;

src/musrScintSD.cc

@@ -158,6 +158,7 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
   mySignalSeparationTime    = musrParameters::signalSeparationTime;
   myStoreOnlyTheFirstTimeHit= musrParameters::storeOnlyTheFirstTimeHit;
   myStoreOnlyEventsWithHitInDetID = musrParameters::storeOnlyEventsWithHitInDetID;
+  myStoreOnlyEventsWithMuonsDecayedInDetID = musrParameters::storeOnlyEventsWithMuonsDecayedInDetID;
   musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
   boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
   myRootOutput = musrRootOutput::GetRootInstance();
@@ -342,6 +343,10 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
     }
   }
 
+  if (myStoreOnlyEventsWithMuonsDecayedInDetID) {
+    if ((myRootOutput->GetDecayDetectorID())!=myStoreOnlyEventsWithMuonsDecayedInDetID) return;
+  }
+
   //  Sort out hits and fill them into root
   if (NbHits>0) {
     const G4int det_IDmax = musrRootOutput::det_nMax;
@@ -1074,13 +1079,18 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void musrScintSD::EndOfRun() {
-  for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsumHistograms.begin(); it!=mapOfOPSAsumHistograms.end(); it++) {
+  if (musrParameters::boolG4OpticalPhotons) {
-    (it->second) -> Write();
+    if (!mapOfOPSAsumHistograms.empty()) {
+      for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsumHistograms.begin(); it!=mapOfOPSAsumHistograms.end(); it++) {
+	(it->second) -> Write();
+      }
+    }
+    if (!mapOfOPSAsum0Histograms.empty()) {
+      for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsum0Histograms.begin(); it!=mapOfOPSAsum0Histograms.end(); it++) {
+	(it->second) -> Write();
+      }
+    }
   }
-  for (mapOfOPSAsumHistograms_Type::const_iterator it = mapOfOPSAsum0Histograms.begin(); it!=mapOfOPSAsum0Histograms.end(); it++) {
-    (it->second) -> Write();
-  }
-
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

src/musrSteppingAction.cc

@@ -89,18 +89,27 @@ void musrSteppingAction::DoAtTheBeginningOfEvent() {
 void musrSteppingAction::UserSteppingAction(const G4Step* aStep)  { 
 
   G4Track* aTrack = aStep->GetTrack();
+
+  // kill the track, if required by user:
+  if (aTrack->GetDefinition()) {
+    G4String p_name = aTrack->GetDynamicParticle()->GetDefinition()->GetParticleName();
+    if      ((musrParameters::killAllPositrons)&&(p_name == "e+")) {aTrack->SetTrackStatus(fStopAndKill); return;}   // suspend the track
+    else if ((musrParameters::killAllElectrons)&&(p_name == "e-")) {aTrack->SetTrackStatus(fStopAndKill); return;}
+    else if ((musrParameters::killAllGammas)&&(p_name == "gamma")) {aTrack->SetTrackStatus(fStopAndKill); return;}
+    else if ((musrParameters::killAllNeutrinos)&&((p_name == "nu_mu")||(p_name == "anti_nu_mu")
+		  ||(p_name == "nu_e")||(p_name == "anti_nu_e")))  {aTrack->SetTrackStatus(fStopAndKill); return;}
+  }
   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
   G4ThreeVector preStepPosition = preStepPoint->GetPosition();
   G4ThreeVector postStepPosition = postStepPoint->GetPosition();
 
   //  suspend the track if too many steps has already happened (relevant at high field)
-  if (aTrack->GetCurrentStepNumber()>100000) {
+  if (aTrack->GetCurrentStepNumber()>musrParameters::maximumNrOfStepsPerTrack) {
-    musrErrorMessage::GetInstance()->musrError(WARNING,
+    char eMessage[200]; 
-    	  "musrSteppingAction: Current number of steps for the track > 100000 ==> TRACK KILLED",true);
+    sprintf(eMessage,"musrSteppingAction: Current number of steps for the track > %d ==> TRACK KILLED",
-    //   G4double x=aStep->GetPostStepPoint()->GetPosition().x()/mm;
+	    musrParameters::maximumNrOfStepsPerTrack);
-    //   G4double y=aStep->GetPostStepPoint()->GetPosition().y()/mm;
+    musrErrorMessage::GetInstance()->musrError(WARNING,eMessage,true);
-    //   G4double z=aStep->GetPostStepPoint()->GetPosition().z()/mm;
     G4double x = postStepPosition.x()/mm;
     G4double y = postStepPosition.y()/mm;
     G4double z = postStepPosition.z()/mm;
@@ -112,13 +121,14 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep)  {
   }
 
   //  abort the event if it takes too long to finish (e.g. more than 60 seconds)
-  if ((time(0) - realTimeWhenThisEventStarted)>60) {
+  if ((time(0) - realTimeWhenThisEventStarted)>musrParameters::maximumTimePerEvent) {
     G4RunManager* fRunManager = G4RunManager::GetRunManager();
     G4cout<<"musrSteppingAction: event "<<fRunManager->GetCurrentEvent()->GetEventID()
-  	  <<" aborted because calculation took already 60 seconds."<<G4endl;
+  	  <<" aborted because calculation took already "<<musrParameters::maximumTimePerEvent<<" seconds."<<G4endl;
-    musrErrorMessage::GetInstance()->musrError(WARNING,
+    char eMessage[200]; 
-              "musrSteppingAction:  event aborted because its calculation takes more than 60 seconds.",true);
+    sprintf(eMessage,"musrSteppingAction:  event aborted because its calculation takes more than %d seconds.",
-    //delete    musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
+	    musrParameters::maximumTimePerEvent);
+    musrErrorMessage::GetInstance()->musrError(WARNING,eMessage,true);
     myRootOutput->SetEventWeight(0);
     fRunManager->AbortEvent();
   }
@@ -372,11 +382,6 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep)  {
       // There is an example how to delete the track in example/novice/N04.
       // It is done in a different way here, because the example/novice/N04 was not doing
       // exactly what I wanted.
-      if ( ((musrParameters::killAllPositrons)&&(p_name == "e+")) ||
-	   ((musrParameters::killAllGammas)&&(p_name == "gamma"))    ||
-	   ((musrParameters::killAllNeutrinos)&&((p_name == "nu_mu")||(p_name == "anti_nu_mu")||(p_name == "nu_e")||(p_name == "anti_nu_e"))) ){
-	aTrack->SetTrackStatus(fStopAndKill);   // suspend the track
-      }
       if((actualVolume(0,10)=="log_shield")||(actualVolume(0,10)=="log_Shield")) { 
 	aTrack->SetTrackStatus(fStopAndKill);   // suspend the track
       }

src/musrSteppingVerbose.cc

@@ -65,7 +65,8 @@ void musrSteppingVerbose::StepInfo()
 	<< std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length")
 	<< std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length")
 	<< std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy")
-	<< std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
+      //	<< std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
+	<< std::setw(6) << G4BestUnit(fStep->GetDeltaEnergy(),"Energy")
 	<< std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length")
 	<< std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length")
 	<< "  ";

src/musrTabulatedElementField.cc

@@ -518,17 +518,17 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
     //cks  The following check is necessary - even though xindex and zindex should never be out of range,
     //     it may happen (due to some rounding error ?).  It is better to leave the check here.
     if ((xindex<0)||(xindex>(nx-2))) {
-      std::cout<<"SERIOUS PROBLEM:  xindex out of range!  xindex="<<xindex<<"   x="<<x<<"  xfraction="<<xfraction<<std::endl;
+      std::cout<<"SERIOUS PROBLEM:  xindex out of range!  xindex="<<xindex<<"   x="<<x<<"  xfraction="<<xfraction<<"   maximumx="<<maximumx<<std::endl;
       if (xindex<0) xindex=0;
       else xindex=nx-2;
     }
     if ((yindex<0)||(yindex>(ny-2))) {
-      std::cout<<"SERIOUS PROBLEM:  yindex out of range!  yindex="<<yindex<<"   y="<<y<<"  yfraction="<<yfraction<<std::endl;
+      std::cout<<"SERIOUS PROBLEM:  yindex out of range!  yindex="<<yindex<<"   y="<<y<<"  yfraction="<<yfraction<<"   maximumy="<<maximumy<<std::endl;
       if (yindex<0) yindex=0;
       else yindex=ny-2;
     }
     if ((zindex<0)||(zindex>(nz-2))) {
-      std::cout<<"SERIOUS PROBLEM:  zindex out of range!  zindex="<<zindex<<"   z="<<z<<"  zfraction="<<zfraction<<std::endl;
+      std::cout<<"SERIOUS PROBLEM:  zindex out of range!  zindex="<<zindex<<"   z="<<z<<"  zfraction="<<zfraction<<"   maximumz="<<maximumz<<std::endl;
       if (zindex<0) zindex=0;
       else zindex=nz-2;
     }