Skip to content

Bugfix issue 266 #277

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
paulmking wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Bugfix issue 266 #277

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
74e8480
Added a method to initialize the global randomization seed within the…
paulmking Feb 7, 2026
f08069a
Removed the "randomnessGenerator" variable from the mock data generat…
paulmking Feb 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion Analysis/include/MQwMockable.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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();
Expand Down
106 changes: 1 addition & 105 deletions Parity/main/QwMockDataGenerator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -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;

Expand All @@ -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
Expand Down Expand Up @@ -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();
Expand Down
Loading