24.3.2010 Kamil Sedlak

New variable "timeToNextEvent" added in order to simulate time differences between subsequent events.
	1)  /gun/meanarrivaltime meanArrivalTime (defines the mean time difference betweent 
                                                  the subsequent events).
        2)  timeToNextEvent - new variable written out to the root tree.

   This new variable is needed for the pile-up studies

doc/musrSim.tex

@@ -468,6 +468,12 @@ Three special volumes ``Target, M0, M1 and M2''.
 	Set the primary particle type, if it is not positive muon.  For example, the negative muons
 	are specified by ``/gun/primaryparticle mu-''.
 
+\item{\bf /gun/meanarrivaltime \emph{meanArrivalTime}}\\
+        (default: /gun/meanarrivaltime 33.33333 microsecond)\\
+        Set mean arrival time difference between two subsequent muons (at continuos beam).
+	The output variable ``timeToNextEvent'' is subsequently generated using
+	the value of \emph{meanArrivalTime} and filled into the Root tree.
+
 \item{\bf /gun/vertex \emph{x0} \emph{y0} \emph{z0} \emph{unit}}\\
 	(default: /gun/vertex 0 0 -100 mm) \\
 	Set mean values of the $x$, $y$ and $z$ coordinates of the generated particles (muons).
@@ -827,6 +833,8 @@ The list of variables that can be stored in the Root tree:
 \begin{description}
 \item{\bf runID} (Int\_t) -- run ID number.        
 \item{\bf eventID} (Int\_t) -- event ID number.
+\item{\bf timeToNextEvent} (Double\_t) -- time difference between two subsequent events (at continuous beam facility) (in $\mu$s).
+        This time difference is generated randomly using the ``meanArrivalTime'' defined in ``/gun/meanarrivaltime meanArrivalTime''.
 \item{\bf weight} (Double\_t) -- event weight.
 \item{\bf BFieldAtDecay\_Bx, BFieldAtDecay\_By, BFieldAtDecay\_Bz, BFieldAtDecay\_B3, BFieldAtDecay\_B4, BFieldAtDecay\_B5} (Double\_t) --
 	value of the 6 coordinates of the electromagnetic field at the position and time where and when the muon decayed.

include/musrPrimaryGeneratorAction.hh

@@ -55,6 +55,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     void SetvertexRelativeR(G4double val)   {relativeRMaxAllowed=val;}
     void SetboxBoundary(G4ThreeVector v)    {xMaxSource=v[0]; yMaxSource=v[1]; zMaxSource=v[2];}          //P.B. 15 Dec 2009
     void SetboxBoundaryCentre(G4ThreeVector v) {xMaxSource0=v[0]; yMaxSource0=v[1]; zMaxSource0=v[2];}    //P.B. 15 Dec 2009
+    void SetMeanArrivalTime(G4double val) {meanArrivalTime=val;}
     void SetMuonTime(G4double val)          {t0=val;}     //P.B. 13 May 2009
     void SetMuonTimeSigma(G4double val)     {tSigma=val;} //P.B. 13 May 2009
     void SetKEnergy(G4double val);
@@ -101,6 +102,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     static G4String thePrimaryParticleName ;
 
     G4double x0, y0, z0, xSigma, ySigma, zSigma, rMaxAllowed, zMinAllowed, zMaxAllowed;
+    G4double meanArrivalTime;
     G4double t0, tSigma;                                 //P.B. 13 May 2009
     G4double relativeRMaxAllowed;
     G4double xMaxSource0, yMaxSource0, zMaxSource0;      //P.B. 15 Dec 2009

include/musrPrimaryGeneratorMessenger.hh

@@ -52,6 +52,7 @@ class musrPrimaryGeneratorMessenger: public G4UImessenger
     G4UIcmdWith3VectorAndUnit* setvertexSigmaCmd;
     G4UIcmdWith3VectorAndUnit* setvertexBoundaryCmd;
     G4UIcmdWithADoubleAndUnit* setvertexRelativeRCmd;
+    G4UIcmdWithADoubleAndUnit* setMeanArrivalTimeCmd;
     G4UIcmdWithADoubleAndUnit* setStarttimeCmd;           //P.B. 13 May 2009
     G4UIcmdWithADoubleAndUnit* setStarttimeSigmaCmd;      //P.B. 13 May 2009
     G4UIcmdWith3VectorAndUnit* setboxBoundaryCmd;         //P.B. 15 Dec 2009

include/musrRootOutput.hh

@@ -60,6 +60,7 @@ class musrRootOutput  {
   // Setting variables common to the whole event:
     void SetRunID          (G4int id) {runID_t = id;};
     void SetEventID        (G4int id) {eventID_t = id;};
+    void SetTimeToNextEvent(G4double deltaT) {timeToNextEvent_t = deltaT/microsecond;}
     void SetDecayDetectorID (std::string detectorName) {muDecayDetID_t = SensDetectorMapping[detectorName];};
     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;};
@@ -128,6 +129,7 @@ class musrRootOutput  {
     static G4bool store_runID;
     static G4bool store_eventID;
     static G4bool store_weight;
+    static G4bool store_timeToNextEvent;
     static G4bool store_BFieldAtDecay;
     static G4bool store_muIniTime;
     static G4bool store_muIniPosX;
@@ -230,6 +232,7 @@ class musrRootOutput  {
     Int_t runID_t;
     Int_t eventID_t;
     Double_t weight_t;
+    Double_t timeToNextEvent_t;
     Double_t B_t[6];
     Double_t muIniTime_t;
     Double_t muIniPosX_t, muIniPosY_t, muIniPosZ_t;

src/musrDetectorConstruction.cc

@@ -161,7 +161,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 
       
       else if (strcmp(tmpString1,"construct")==0){
-	float x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0; 
+	float x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0,x8=0,x9=0,x10=0,x11,x12; 
 	float posx,posy,posz;
 	char name[100];
 	char mothersName[100];
@@ -333,6 +333,22 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg);
 	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetBox);
 	}
+	else if (strcmp(tmpString2,"GPDcollimator")==0){
+	  // Create a box, from which a box is cut out.  x1, x2, x3 = box half-widths; 
+	  //           x4,x5,x6,x7 define the tube, x8, x9 and x10 are the distances between the tube and box centres.
+	  //           x11, x12 are the half-withs of the ractangular opening in the collimator
+	  sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %g %g %g %g %g %g %s",
+		 name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,&x8,&x9,&x10,&x11,&x12,material);
+	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %g %g %g %s %s %s %d %s",&posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
+	  solidName+=name;
+	  G4Box* solidDetBox = new G4Box("SolidDetBox",x1*mm,x2*mm,x3*mm);
+	  G4Box* solidHole   = new G4Box("SolidDetBox",x11*mm,x2*mm+0.1,x12*mm);
+	  G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.,x4*mm,x5*mm,x6*deg,x7*deg);
+	  G4RotationMatrix* yRot = new G4RotationMatrix();
+	  G4ThreeVector zTrans(-x8*mm,-x9*mm,-x10*mm);
+	  G4SubtractionSolid* solidPart1 =  new G4SubtractionSolid("solidPart1", solidDetBox, solidDetTube, yRot, zTrans);
+	  solid =  new G4SubtractionSolid(solidName,solidPart1,solidHole);
+	}
 	else if (strcmp(tmpString2,"tubsboxsegm")==0){
 	  // Create a volume that looks like an intersection of tube and box.
 	  char orientation[100];
@@ -359,17 +375,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	else if (strcmp(tmpString2,"GPDsampleHolderA")==0){
 	  // First part of the GPD sample holder, where posx, posy, posz = centre of the whole (long) tube 
 	  //                                                      (=111.25mm below the centre of the holes)
-	  // Create a tube, from which center a box is cut out.  x1=box half-width; x2,x3,x4,x5 define the tube.
-	  //	  sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
-	  //		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
-	  //	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
-	  //	  solidName+=name;
-	  //	  G4double roundingErr=0.01*mm;  // to avoid some displaying problems of the subtracted volumes
-	  //	  G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*mm,x1*mm,x3*mm+roundingErr);
-	  //	  G4Tubs* solidOuterDetTube = new G4Tubs("SolidOuterDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg);
-	  //	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetBox);
-
-
 	  sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
 		 name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
 	  sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
@@ -391,10 +396,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 	  G4SubtractionSolid* solidA123 =  new G4SubtractionSolid("solidA123", solidA12 , solidGPDBoxA2, yRotA12, zTransA13);
 	  G4ThreeVector zTransA5(0,0,111.25*mm);
 	  solid =  new G4SubtractionSolid(solidName, solidA123, solidGPDBoxA5, yRotA12, zTransA5);
-	  //G4cout<<"solidName="<<solidName<<"  solidA123="<<solidA123<<"   solidGPDBoxA5="<<solidGPDBoxA5<<"   yRotA12="<<yRotA12<<"  zTransA5="<<zTransA5<<G4endl;
-	  //solid =  new G4SubtractionSolid(solidName, solidGPDTubeA1, solidGPDBoxA2, yRotA12, zTransA12);
-	  
-	  G4cout<<"Debug 80"<<G4endl;
 	}
 
 	else ReportGeometryProblem(line);

src/musrPrimaryGeneratorAction.cc

@@ -58,7 +58,7 @@ musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
                                             musrDetectorConstruction* musrDC)
   :musrDetector(musrDC), x0(0), y0(0), z0(-10*cm), xSigma(0), ySigma(0), zSigma(0), 
    rMaxAllowed(1e10*mm), zMinAllowed(-1e10*mm), zMaxAllowed(1e10*mm),
-   t0(0), tSigma(0),
+   meanArrivalTime(1./30000.*second), t0(0), tSigma(0),
    relativeRMaxAllowed(1e10*mm), 
    xMaxSource0(0), yMaxSource0(0), zMaxSource0(0),
    xMaxSource(1e10*mm), yMaxSource(1e10*mm), zMaxSource(1e10*mm),
@@ -343,9 +343,14 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
     generatedMuon->SetProperTime(decaytime);
   }
 
+  // Set the variable "timeToNextEvent", which is the time difference between this event and the next one
+  // at a continuous muon beam.
+  G4double timeIntervalBetweenTwoEvents = meanArrivalTime * CLHEP::RandExponential::shoot(1);
+
   // Save variables into ROOT output file:
   myRootOutput->SetInitialMuonParameters(x,y,z,px,py,pz,xpolaris,ypolaris,zpolaris,ParticleTime);
   myRootOutput->StoreGeantParameter(7,float(numberOfGeneratedEvents));
+  myRootOutput->SetTimeToNextEvent(timeIntervalBetweenTwoEvents);
   if (boolPrintInfoAboutGeneratedParticles) {
     G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries:  x="<<x<<", y="<<y<<", z="<<z<<G4endl;
     G4cout<<"      px="<<px<<", py="<<py<<", pz="<<pz<<", xpolaris="<<xpolaris<<", ypolaris="<<ypolaris<<", zpolaris="<<zpolaris<<G4endl;

src/musrPrimaryGeneratorMessenger.cc

@@ -60,6 +60,11 @@ musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGenerato
   setvertexRelativeRCmd->SetParameterName("mes_relativeRMaxAllowed",true);
   setvertexRelativeRCmd->SetDefaultUnit("mm"); 
 
+  setMeanArrivalTimeCmd = new G4UIcmdWithADoubleAndUnit("/gun/meanarrivaltime",this);
+  setMeanArrivalTimeCmd->SetGuidance(" Set mean arrival time difference between two subsequent muons (at continuos beam)");
+  setMeanArrivalTimeCmd->SetParameterName("mes_meanarrivaltime",true);
+  setMeanArrivalTimeCmd->SetDefaultUnit("ns");
+
   setStarttimeCmd = new G4UIcmdWithADoubleAndUnit("/gun/starttime",this);                //P.B. 13 May 2009
   setStarttimeCmd->SetGuidance(" Set start time t of the generated muons (with unit)");  //P.B. 13 May 2009
   setStarttimeCmd->SetParameterName("mes_t0",true);                                      //P.B. 13 May 2009
@@ -178,6 +183,7 @@ musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
   delete setvertexSigmaCmd;
   delete setvertexBoundaryCmd;
   delete setvertexRelativeRCmd;
+  delete setMeanArrivalTimeCmd;
   delete setStarttimeCmd;         //P.B. 13 May 2009
   delete setStarttimeSigmaCmd;    //P.B. 13 May 2009
   delete setboxBoundaryCentreCmd; //P.B. 15 Dec 2009
@@ -213,6 +219,8 @@ void musrPrimaryGeneratorMessenger::SetNewValue(G4UIcommand * command,G4String n
     { musrAction->SetvertexBoundary(setvertexBoundaryCmd->GetNew3VectorValue(newValue));}
   if( command == setvertexRelativeRCmd)
     { musrAction->SetvertexRelativeR(setvertexRelativeRCmd->GetNewDoubleValue(newValue));}
+  if( command == setMeanArrivalTimeCmd)
+    { musrAction->SetMeanArrivalTime(setMeanArrivalTimeCmd->GetNewDoubleValue(newValue));}
   if( command == setStarttimeCmd)                                                       //P.B. 13 May 2009
     { musrAction->SetMuonTime(setStarttimeCmd->GetNewDoubleValue(newValue));}           //P.B. 13 May 2009
   if( command == setStarttimeSigmaCmd)                                                  //P.B. 13 May 2009

src/musrRootOutput.cc

@@ -68,6 +68,7 @@ musrRootOutput* musrRootOutput::GetRootInstance() {
 G4bool musrRootOutput::store_runID = true;
 G4bool musrRootOutput::store_eventID = true;
 G4bool musrRootOutput::store_weight = true;
+G4bool musrRootOutput::store_timeToNextEvent = true;
 G4bool musrRootOutput::store_BFieldAtDecay = true;
 G4bool musrRootOutput::store_muIniTime = true;
 G4bool musrRootOutput::store_muIniPosX = true;
@@ -161,6 +162,7 @@ void musrRootOutput::BeginOfRunAction() {
   if (store_runID)        {rootTree->Branch("runID",&runID_t,"runID/I");}
   if (store_eventID)      {rootTree->Branch("eventID",&eventID_t,"eventID/I");}
   if (store_weight)       {rootTree->Branch("weight",&weight_t,"weight/D");}
+  if (store_timeToNextEvent){rootTree->Branch("timeToNextEvent",&timeToNextEvent_t,"timeToNextEvent/D");}
   if (store_BFieldAtDecay) {rootTree->Branch("BFieldAtDecay",&B_t,"Bx/D:By:Bz:B3:B4:B5");}
   if (store_muIniTime)    {rootTree->Branch("muIniTime",&muIniTime_t,"muIniTime/D");}
   if (store_muIniPosX)    {rootTree->Branch("muIniPosX",&muIniPosX_t,"muIniPosX/D");}
@@ -327,6 +329,7 @@ void musrRootOutput::ClearAllRootVariables() {
   runID_t=-1000;
   eventID_t=-1000;
   weight_t=1.;
+  timeToNextEvent_t = -1000;
   B_t[0]=-1000.;B_t[1]=-1000.;B_t[2]=-1000.;B_t[3]=-1000.;B_t[4]=-1000.;B_t[5]=-1000.;
   muIniTime_t=-1000;
   muIniPosX_t=-1000; muIniPosY_t=-1000; muIniPosZ_t=-1000;