Bugfix issue 266 (#277)

* Added a method to initialize the global randomization seed within the MQwMockable
mixin class for each run, and called it from the mock data generator main routine.

* Removed the "randomnessGenerator" variable from the mock data generator main routine.
This had been previously used as the randomization source, but now all signal and device
randomization is doen using an internal variable in the MQwMockable mixin class.

Also removed the old randomization and correlation lines from the mock data generator main
routine, which are now handled through the correlator class.

Cleaned up some out of date comment lines from mock data generator main routine.

Author: paulmking

Analysis/include/MQwMockable.h

@@ -39,7 +39,7 @@ class MQwMockable {
 public:
   MQwMockable(): fUseExternalRandomVariable(false),
                  fCalcMockDataAsDiff(false),
-		 fMockAsymmetry(0.0), fMockGaussianMean(0.0),
+                 fMockAsymmetry(0.0), fMockGaussianMean(0.0),
                  fMockGaussianSigma(0.0)
   {
     // Mock drifts
@@ -68,6 +68,11 @@ class MQwMockable {
   /// Set the helicity asymmetry
   void  SetRandomEventAsymmetry(Double_t asymmetry);
 
+  ///  Seed the internal Random Variable
+  static void Seed(uint seedval){
+    fRandomnessGenerator.seed(seedval);
+  }
+
   /// Return a random value generated either from the internal or
   /// external Random Variable.
   Double_t GetRandomValue();

Parity/main/QwMockDataGenerator.cc

@@ -80,26 +80,10 @@ int main(int argc, char* argv[])
   QwHelicity* helicity = dynamic_cast<QwHelicity*>(detectors.GetSubsystemByName("Helicity Info"));
   if (! helicity) QwWarning << "No helicity subsystem defined!" << QwLog::endl;
 
-  // Possible scenarios:
-  // - everything is random, no correlations at all, no asymmetries at all
-  // - variations in the mean of position, current, yield over the course of
-  //   a run (linearly with run number, change mean/sigma as function of event
-  //   number)
-  // - one parameter has a helicity-correlated asymmetry, every other parameter
-  //   is random and independently distributed
-  // - two parameters have independent helicity-correlated asymmetries
-  // - two parameters have correlated helicity-correlated asymmetries
-  // - beam modulation
 
   // Get the beamline channels we want to correlate
   detectors.LoadMockDataParameters("mock_parameters_list.map");
 
-//-----------------------------------------------------------------------------------------------
-  // Get the main detector channels we want to correlate
-//  QwDetectorArray* maindetector =
-//    dynamic_cast<QwDetectorArray*>(detectors.GetSubsystemByName("Main Detector"));
-//  if (! maindetector) QwWarning << "No main detector subsystem defined!" << QwLog::endl;
-
   // new vectors for GetSubsystemByType
   std::vector <QwDetectorArray*> detchannels;
   std::vector <VQwSubsystem*> tempvector = detectors.GetSubsystemByType("QwDetectorArray");
@@ -111,52 +95,6 @@ int main(int argc, char* argv[])
   //return detchannels;
   }
 
-/*
-if(1==2){
-  Double_t bar_mean = 2.0e4;
-  Double_t bar_sigma = 3.0e2;
-  Double_t bar_asym = 4.0e-4;
-  maindetector->SetRandomEventParameters(bar_mean, bar_sigma);
-  maindetector->SetRandomEventAsymmetry(bar_asym);
-  // Specific values
-  // disabled, wdc, 2010-07-23
-  maindetector->GetIntegrationPMT("MD2Neg")->SetRandomEventAsymmetry(2.0e-2);
-  maindetector->GetIntegrationPMT("MD2Pos")->SetRandomEventAsymmetry(2.0e-2);
-  maindetector->GetIntegrationPMT("MD3Neg")->SetRandomEventAsymmetry(3.0e-3);
-  maindetector->GetIntegrationPMT("MD3Pos")->SetRandomEventAsymmetry(3.0e-3);
-  maindetector->GetIntegrationPMT("MD4Neg")->SetRandomEventAsymmetry(4.0e-4);
-  maindetector->GetIntegrationPMT("MD4Pos")->SetRandomEventAsymmetry(4.0e-4);
-  maindetector->GetIntegrationPMT("MD5Neg")->SetRandomEventAsymmetry(5.0e-5);
-  maindetector->GetIntegrationPMT("MD5Pos")->SetRandomEventAsymmetry(5.0e-5);
-  maindetector->GetIntegrationPMT("MD6Neg")->SetRandomEventAsymmetry(6.0e-6);
-  maindetector->GetIntegrationPMT("MD6Pos")->SetRandomEventAsymmetry(6.0e-6);
-  maindetector->GetIntegrationPMT("MD7Neg")->SetRandomEventAsymmetry(7.0e-7);
-  maindetector->GetIntegrationPMT("MD7Pos")->SetRandomEventAsymmetry(7.0e-7);
-  maindetector->GetIntegrationPMT("MD8Neg")->SetRandomEventAsymmetry(8.0e-8);
-  maindetector->GetIntegrationPMT("MD8Pos")->SetRandomEventAsymmetry(8.0e-8);
-
-  // Set a asymmetric helicity asymmetry on one of the bars
-  maindetector->GetIntegrationPMT("MD1Neg")->SetRandomEventAsymmetry(5.0e-5);
-  maindetector->GetIntegrationPMT("MD1Pos")->SetRandomEventAsymmetry(-5.0e-5);
-
-  // Set a drift component (amplitude, phase, frequency)
-  maindetector->GetIntegrationPMT("MD3Neg")->AddRandomEventDriftParameters(3.0e6, 0, 60*Qw::Hz);
-  maindetector->GetIntegrationPMT("MD3Neg")->AddRandomEventDriftParameters(6.0e5, 0, 120*Qw::Hz);
-  maindetector->GetIntegrationPMT("MD3Neg")->AddRandomEventDriftParameters(4.5e5, 0, 240*Qw::Hz);
-
-} //end if(1==2)
-*/
-
-
-  // Initialize randomness provider and distribution
-  std::mt19937 randomnessGenerator(999); // Mersenne twister with seed (see below)
-  std::normal_distribution<double> normalDistribution;
-  auto normal = [&]() -> double { return normalDistribution(randomnessGenerator); };
-  // WARNING: This variate_generator will return the SAME random values as the
-  // variate_generator in QwVQWK_Channel when used with the same default seed!
-  // Therefore you really should give an explicitly different seed for the
-  // mt19937 randomness generator.
-
   // Initialize the stopwatch
   TStopwatch stopwatch;
 
@@ -168,7 +106,7 @@ if(1==2){
               run++) {
 
     // Set the random seed for this run
-    randomnessGenerator.seed(run);
+    MQwMockable::Seed(run);
     QwCombinedBCM<QwVQWK_Channel>::SetTripSeed(0x56781234 ^ (run*run));
 
     // Open new output file
@@ -249,48 +187,6 @@ if(1==2){
       int myhelicity = helicity->GetHelicityActual() ? +1 : -1;
       //std::cout << myhelicity << std::endl;
 
-      // Secondly introduce correlations between variables
-      //
-      // N-dimensional correlated normal random variables:
-      //   X = C' * Z
-      // with
-      //   X correlated and normally distributed,
-      //   Z independent and normally distributed,
-      //   C the Cholesky decomposition of the positive-definite covariance matrix
-      //     (C should probably be calculated offline)
-      //
-      /* Sigma =
-           1.00000   0.50000   0.50000
-           0.50000   2.00000   0.30000
-           0.50000   0.30000   1.50000
-
-         C =
-           1.00000   0.50000   0.50000
-           0.00000   1.32288   0.03780
-           0.00000   0.00000   1.11739
-
-         Sigma = C' * C
-       */
-      double z[NVARS], x[NVARS];
-      double C[NVARS][NVARS];
-      for (int var = 0; var < NVARS; var++) {
-        x[var] = 0.0;
-        z[var] = normal();
-        C[var][var] = 1.0;
-      }
-      C[0][0] = 1.0; C[0][1] = 0.5;     C[0][2] = 0.5;
-      C[1][0] = 0.0; C[1][1] = 1.32288; C[1][2] = 0.03780;
-      C[2][0] = 0.0; C[2][1] = 0.0;     C[2][2] = 1.11739;
-      for (int i = 0; i < NVARS; i++)
-        for (int j = 0; j < NVARS; j++)
-          x[i] += C[j][i] * z[j];
-
-      // Assign to data elements
-      //maindetector->GetChannel("MD2Neg")->SetExternalRandomVariable(x[0]);
-      //lumidetector->GetChannel("dlumi1")->SetExternalRandomVariable(x[1]);
-      //beamline->GetBCM("qwk_bcm0l07")->SetExternalRandomVariable(x[2]);
-
-
       // Randomize data for this event
       detectors.RandomizeEventData(myhelicity, time);
 //      detectors.ProcessEvent();