diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index 5feb581..a2293d5 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index 6a8e5ba..a21af0b 100644 --- a/doc/musrSim.tex +++ b/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. diff --git a/include/musrPrimaryGeneratorAction.hh b/include/musrPrimaryGeneratorAction.hh index 026d6a0..4965535 100644 --- a/include/musrPrimaryGeneratorAction.hh +++ b/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 diff --git a/include/musrPrimaryGeneratorMessenger.hh b/include/musrPrimaryGeneratorMessenger.hh index aa1d1c2..ae7af0e 100644 --- a/include/musrPrimaryGeneratorMessenger.hh +++ b/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 diff --git a/include/musrRootOutput.hh b/include/musrRootOutput.hh index c4e8d12..0d9320e 100644 --- a/include/musrRootOutput.hh +++ b/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; diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 722f1ef..1455ae8 100644 --- a/src/musrDetectorConstruction.cc +++ b/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="<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="<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 diff --git a/src/musrRootOutput.cc b/src/musrRootOutput.cc index 4c51647..4a33f99 100644 --- a/src/musrRootOutput.cc +++ b/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;