diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index b9d1461..8a4c3b9 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index 163a6df..5a4f289 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -1177,13 +1177,24 @@ The list of variables that can be stored in the Root tree: the photon enters (and not \emph{from} which it comes). This e.g.\ allows one to attach more APDs to one scintillator, in which case odet\_ID[] will correspond to the volume IDs of individual APDs. \item{\bf odet\_nPhot[odet\_n]} (array of Int\_t) -- number of photons detected in the given optical photon signal. -\item{\bf odet\_timeFirst[odet\_n], odet\_timeLast[odet\_n]} (array of Int\_t) +\item{\bf odet\_timeFirst[odet\_n], odet\_timeLast[odet\_n]} (array of Double\_t) -- times, when the the first and last photons contributing to the given signal were detected (odet\_timeFirst and odet\_timeLast). -\item{\bf odet\_timeA[odet\_n], odet\_timeB[odet\_n], odet\_timeC[odet\_n], odet\_timeD[odet\_n], odet\_timeE[odet\_n]} (array of Int\_t) +\item{\bf odet\_timeA[odet\_n], odet\_timeB[odet\_n]} (array of Double\_t) -- time, when the $n^{th}$ photon was detected, where $n$ for \emph{odet\_timeA[i]} is given as \emph{OPSA\_fracA} - multiplied by \emph{odet\_nPhot[i]}. Similary for \emph{odet\_timeB[i]}, ..., \emph{odet\_timeE[i]}. - The variables \emph{OPSA\_fracA}, ..., \emph{OPSA\_fracE}, are defined by ``/musr/command OPSA photonFractions'' command + multiplied by \emph{odet\_nPhot[i]}. Similary for \emph{odet\_timeB[i]}. + The variables \emph{OPSA\_fracA} and \emph{OPSA\_fracB} are defined by the ``/musr/command OPSA photonFractions'' command (see chapter~\ref{sec:opticalPhotons}). +\item{\bf odet\_timeC[odet\_n]} (array of Double\_t) + -- time, when the shape of the signal becomes larger than the threshold ``OPSA\_C\_threshold'' defined by + the ``/musr/command OPSA photonFractions'' command (see chapter~\ref{sec:opticalPhotons}). +\item{\bf odet\_timeD[odet\_n]} (array of Double\_t) + -- time, when the shape of the {\bf normalised} signal becomes larger than the threshold ``OPSA\_D\_threshold'' defined by + the ``/musr/command OPSA photonFractions'' command (see chapter~\ref{sec:opticalPhotons}). The normalisation is done + to 1. +\item{\bf odet\_timeCFD[odet\_n]} (array of Double\_t) + -- time as from the ``Constant Fraction Discriminator''. +\item{\bf odet\_timeMean[odet\_n]} (array of Double\_t) + -- mean time of the arriving photons (calculated from the OPSAhist defined by the command ``/musr/command OPSA OPSAhist ''. \item{\bf save\_n} (Int\_t) -- number of special kind of ``save'' volume that were hit in this event. The ``save volume'' is any volume whose name starts with letters ``save''. Their purpose in the simulation is usually to check positions and momenta of particles at some position of the detector, even if the particle does not deposit any energy in diff --git a/include/musrScintSD.hh b/include/musrScintSD.hh index b9973c3..8445b88 100644 --- a/include/musrScintSD.hh +++ b/include/musrScintSD.hh @@ -78,8 +78,8 @@ class musrScintSD : public G4VSensitiveDetector // Optical Photon Signal Analysis (OPSA) void Set_OPSA_minNrOfDetectedPhotons(G4int val) {OPSA_minNrOfDetectedPhotons=val;} void Set_OPSA_SignalSeparationTime(G4double val) {OPSA_signalSeparationTime=val;} - void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d) - {OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d;} + void Set_OPSA_variables(G4double a, G4double b, G4double c, G4double d) + {OPSA_fracA=a; OPSA_fracB=b; OPSA_C_threshold=c; OPSA_D_threshold=d;} void Set_OPSA_CFD(G4double a1, G4double delay, G4double timeShiftOffset) {OPSA_CFD_a1=a1; OPSA_CFD_delay=delay; OPSA_CFD_timeShiftOffset = timeShiftOffset;} void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID) { @@ -115,8 +115,8 @@ class musrScintSD : public G4VSensitiveDetector G4double OPSA_signalSeparationTime; G4double OPSA_fracA; G4double OPSA_fracB; - G4double OPSA_fracC; - G4double OPSA_fracD; + G4double OPSA_C_threshold; + G4double OPSA_D_threshold; typedef std::multimap OPSA_signal_MapType; OPSA_signal_MapType OPSA_signal_Map; diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index bad9019..5a1bafe 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -713,7 +713,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { else if (strcmp(varName,"photonFractions")==0) { double a, b, c, d; sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf",&a, &b, &c, &d); - myMusrScintSD -> Set_OPSA_frac(a,b,c,d); + myMusrScintSD -> Set_OPSA_variables(a,b,c,d); } else if (strcmp(varName,"eventsForOPSAhistos")==0) { int i_eventID, i_detectorID; diff --git a/src/musrEventAction.cc b/src/musrEventAction.cc index 28cfb21..6a1dce2 100644 --- a/src/musrEventAction.cc +++ b/src/musrEventAction.cc @@ -46,6 +46,7 @@ G4double musrEventAction::maximumRunTimeAllowed=85000; musrEventAction::musrEventAction() { pointer=this; + timeOfRunStart = -1000; } musrEventAction* musrEventAction::pointer=0; musrEventAction* musrEventAction::GetInstance() { @@ -97,7 +98,7 @@ void musrEventAction::EndOfEventAction(const G4Event* evt) { // periodic printing // // if (thisEventNr != 0 and thisEventNr%10000 == 0) { - if (thisEventNr == 0) timeOfRunStart=time(0); + if (timeOfRunStart == -1000) timeOfRunStart=time(0); // initiate the time when execution starts once at the beginning of run if ((time(0)-timeOfRunStart)>maximumRunTimeAllowed) { // Stop the execution of the run - the run took already too long time char eMessage[200]; diff --git a/src/musrPrimaryGeneratorAction.cc b/src/musrPrimaryGeneratorAction.cc index b8880dc..065bbd1 100644 --- a/src/musrPrimaryGeneratorAction.cc +++ b/src/musrPrimaryGeneratorAction.cc @@ -492,11 +492,11 @@ void musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds(G4Event* anEvent) if (!pointerToSeedVector->empty()) { eventID = pointerToSeedVector->back(); pointerToSeedVector->pop_back(); - anEvent -> SetEventID(eventID); } else { eventID = ++lastEventID_in_pointerToSeedVector; } + anEvent -> SetEventID(eventID); CLHEP::HepRandom::setTheSeed(eventID); CLHEP::RandGauss::setFlag(false); // G4cout <<"musrPrimaryGeneratorAction.cc: seed will be set to="<< eventID<Fill(timePhot-OPSA_timeFirst+0.00000000001); + if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) { + G4double tt = timePhot-OPSA_timeFirst+0.00000000001; + OPSAhisto->Fill(tt); + if (bool_pulseShapeExists) { // fill contribution of this detected photon to the overall signal shape + Int_t iBin = int(tt/OPSAhistoBinWidth); + Double_t rest = tt - iBin*OPSAhistoBinWidth; + for (Int_t iBinNew=0; iBinNew=iPSmax-1) break; // out of the pulse shape info ==> leave this loop + OPSAshape->Fill((iBin+iBinNew+0.5)*OPSAhistoBinWidth,pulseShapeArray[iPS]); + } + } + } if (bool_StoreThisOPSAhistSUMMED) { OPSAhistoSUM ->Fill(timePhot-OPSA_timeFirst+0.00000000001); OPSAhistoSUM0->Fill(timePhot+0.00000000001); @@ -657,16 +666,18 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() { // if required, convert the histogram with photon times to histogram of electronic pulse shape if (bool_pulseShapeExists) { - for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of photon histogram - Double_t photonCount = OPSAhisto ->GetBinContent(iBin); - if (photonCount==0.) break; //nothing to be done for bins without any photons - for (Int_t iBinNew=0; iBinNew<(OPSAhistoNbin-iBin); iBinNew++) { // loop over all bins from the actual one to the end - Int_t iPS = int( (double(iBinNew)+0.5)*OPSAhistoBinWidth1000 ); - if (iPS>=iPSmax-1) break; // out of the pulse shape info ==> leave this loop - // OPSAshape->AddBinContent(iBin+iBinNew,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth); // iPS corresponds to time in picoseconds - OPSAshape->Fill((iBin+iBinNew-0.5)*OPSAhistoBinWidth,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth); - } - } + // for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of photon histogram + // Double_t photonCount = OPSAhisto ->GetBinContent(iBin); + // // if (photonCount==0.) break; //nothing to be done for bins without any photons + // here was error - instead of continue I used break !!! + // if (photonCount==0.) continue; //nothing to be done for bins without any photons + // for (Int_t iBinNew=0; iBinNew<(OPSAhistoNbin-iBin); iBinNew++) { // loop over all bins from the actual one to the end + // Int_t iPS = int( (double(iBinNew)+0.5)*OPSAhistoBinWidth1000 ); + // if (iPS>=iPSmax-1) break; // out of the pulse shape info ==> leave this loop + // // OPSAshape->AddBinContent(iBin+iBinNew,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth); // iPS corresponds to time in picoseconds + // OPSAshape->Fill((iBin+iBinNew-0.5)*OPSAhistoBinWidth,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth); + // } + // } if (boolStoreThisOPSAhist) OPSAshape -> Write(); // Now fill the histogram with the CFD signal @@ -696,7 +707,34 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() { OPSA_CFD_time += OPSA_timeFirst - OPSA_CFD_delay + OPSA_CFD_timeShiftOffset; // correct for } } - // G4cout<<" OPSA_CFD_time = "<< OPSA_CFD_time<<" OPSA_CFD_ampl="<GetMaximum()); // covert relative "OPSA_D_threshold" into + // the absolute threshold using signal amplitude + for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of CFD histogram + Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth; + oldYvalue = yValue; + yValue = OPSAshape->GetBinContent(iBin); + if ( (yValue>=OPSA_C_threshold) && (!OPSA_C_time_found) ) { // signal just moved above the threshold + OPSA_C_time_found = true; + OPSA_timeC = xValue - (yValue-OPSA_C_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation + OPSA_timeC += OPSA_timeFirst; + } + + if ( (yValue>=D_threshold) && (!OPSA_D_time_found) ) { // signal just moved above the threshold + OPSA_D_time_found = true; + OPSA_timeD = xValue - (yValue-D_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation + OPSA_timeD += OPSA_timeFirst; + } + + if (OPSA_C_time_found && OPSA_D_time_found) break; // times C and D found, no need to continue looping. + } + + // G4cout<<" OPSA_CFD_time = "<< OPSA_CFD_time<<", OPSA_CFD_ampl="<