23.2.2010 Kamil Sedlak

1) changed treatment for the odet_timeCFD determination
  2) changed definition of odet_timeC and odet_timeD
  3) corrected bug in execution time calculation for randomOption 3
  4) corrected bug in eventID setting when randomOption=3

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]}.
+	multiplied by \emph{odet\_nPhot[i]}. Similary for \emph{odet\_timeB[i]}.
-	The variables \emph{OPSA\_fracA}, ..., \emph{OPSA\_fracE}, are defined by ``/musr/command OPSA photonFractions'' command
+	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

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) 
+    void Set_OPSA_variables(G4double a, G4double b, G4double c, G4double d) 
-             {OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=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_C_threshold;
-      G4double OPSA_fracD;
+      G4double OPSA_D_threshold;
       typedef std::multimap<G4int,signalInfo*>  OPSA_signal_MapType;
       OPSA_signal_MapType OPSA_signal_Map;
   

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;

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];

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<<G4endl;

src/musrScintSD.cc

@@ -70,8 +70,8 @@ musrScintSD::musrScintSD(G4String name)
   OPSA_signalSeparationTime   = 10.;
   OPSA_fracA = 0.01;
   OPSA_fracB = 0.05;
-  OPSA_fracC = 0.10;
+  OPSA_C_threshold = 0.5;
-  OPSA_fracD = 0.20;
+  OPSA_D_threshold = 0.5;
   bool_multimapOfEventIDsForOPSAhistosEXISTS = false;
   bool_StoreThisOPSAhistSUMMED = false;
   bool_StoreThisOPSAhistALL = false;
@@ -593,8 +593,6 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
 	  G4double OPSA_f_nPhot = OPSA_nPhot;
 	  G4int NA = int (OPSA_fracA * OPSA_f_nPhot + 0.5);  if (NA<=0) NA=1;
 	  G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5);  if (NB<=0) NB=1;
-	  G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5);  if (NC<=0) NC=1;
-	  G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5);  if (ND<=0) ND=1;
 
 	  Int_t nP=0;
 	  // Define OPSA histograms if required for this event
@@ -641,10 +639,21 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
 	    if (nP==1)          OPSA_timeFirst = timePhot;
 	    if (nP==NA)         OPSA_timeA     = timePhot;
 	    if (nP==NB)         OPSA_timeB     = timePhot;
-	    if (nP==NC)         OPSA_timeC     = timePhot;
-	    if (nP==ND)         OPSA_timeD     = timePhot;
 	    if (nP==OPSA_nPhot) OPSA_timeLast  = timePhot;
-	    if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) OPSAhisto->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<OPSAhistoNbin-iBin; iBinNew++) { // loop over all bins of the shape histogram
+		  Int_t iPS = int( (double(iBinNew))*OPSAhistoBinWidth1000+rest*1000 );  // iPS is time in picoseconds 
+		                                                                // and also the index of the shape array
+		  if (iPS>=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
+	    //	    for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of photon histogram
-	      Double_t photonCount = OPSAhisto ->GetBinContent(iBin);
+	    //	      Double_t photonCount = OPSAhisto ->GetBinContent(iBin);
-	      if (photonCount==0.) break;  //nothing to be done for bins without any photons
+	    //	    //	      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
+	    // here was error - instead of continue I used break !!!  
-		Int_t iPS = int( (double(iBinNew)+0.5)*OPSAhistoBinWidth1000 );
+	    //        if (photonCount==0.) continue;  //nothing to be done for bins without any photons
-		if (iPS>=iPSmax-1) break;  // out of the pulse shape info ==> leave this loop
+	    //	      for (Int_t iBinNew=0; iBinNew<(OPSAhistoNbin-iBin); iBinNew++) { // loop over all bins from the actual one to the end
-		//		  OPSAshape->AddBinContent(iBin+iBinNew,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth);  // iPS corresponds to time in picoseconds
+	    //		Int_t iPS = int( (double(iBinNew)+0.5)*OPSAhistoBinWidth1000 );
-		OPSAshape->Fill((iBin+iBinNew-0.5)*OPSAhistoBinWidth,photonCount*pulseShapeArray[iPS]*OPSAhistoBinWidth);
+	    //		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="<<OPSA_CFD_ampl<<"     OPSA_timeFirst"<<OPSA_timeFirst<<G4endl;
+
+	    // Find the timeC and timeD from the shape histogram
+	    yValue=-1000;
+	    G4bool   OPSA_C_time_found = false;
+	    G4bool   OPSA_D_time_found = false;
+	    Double_t D_threshold = OPSA_D_threshold * (OPSAshape->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="<<OPSA_CFD_ampl
+	    //		  <<",  OPSA_timeFirst="<<OPSA_timeFirst<<",  OPSA_timeC="<<OPSA_timeC <<",  OPSA_timeD="<<OPSA_timeD<<G4endl;
 	  } // end if (pulseShapeExists)
 
 	  // Delete the histograms from memory if they were created