QwMockDataAnalysis.cc File Reference

Parity mock data analysis, test code. More...

#include <iostream>
#include <boost/random.hpp>
#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include "QwBeamLine.h"
#include "QwOptionsParity.h"
#include "QwEventBuffer.h"
#include "QwParityDB.h"
#include "QwHelicity.h"
#include "QwHelicityPattern.h"
#include "QwHistogramHelper.h"
#include "QwMainCerenkovDetector.h"
#include "QwLumi.h"
#include "QwSubsystemArrayParity.h"

Include dependency graph for QwMockDataAnalysis.cc:

Go to the source code of this file.

Macros
#define	NEVENTS   1000
#define	NVARS   4

Functions
int	main (int argc, char *argv[])

Variables
static const int	kMultiplet = 4
static bool	bDebug = false
static bool	bTree = true
static bool	bHisto = true
static bool	bHelicity = true
static bool	bDatabase = false

Detailed Description

Parity mock data analysis, test code.

Definition in file QwMockDataAnalysis.cc.

Macro Definition Documentation

#define NEVENTS   1000

Definition at line 35 of file QwMockDataAnalysis.cc.

#define NVARS   4

Definition at line 37 of file QwMockDataAnalysis.cc.

Function Documentation

int main	(	int	argc,
		char *	argv[]
	)

Fill the search paths for the parameter files; this sets a static variable within the QwParameterFile class which will be used by all instances. The "scratch" directory should be first.

Definition at line 52 of file QwMockDataAnalysis.cc.

References QwHelicityPattern::AccumulateRunningBurstSum(), QwSubsystemArrayParity::AccumulateRunningSum(), QwParameterFile::AppendToSearchPath(), bDatabase, bDebug, bHelicity, bHisto, bTree, QwHelicityPattern::CalculateAsymmetry(), QwHelicityPattern::CalculateBurstAverage(), QwSubsystemArrayParity::CalculateRunningAverage(), QwHelicityPattern::CalculateRunningAverage(), QwHelicityPattern::CalculateRunningBurstAverage(), QwHelicityPattern::ClearBurstSum(), QwHelicityPattern::ClearEventData(), QwEventBuffer::CloseDataFile(), QwSubsystemArrayParity::ConstructBranchAndVector(), QwHelicityPattern::ConstructBranchAndVector(), QwHelicityPattern::ConstructHistograms(), QwSubsystemArray::ConstructHistograms(), DefineOptionsParity(), QwLog::endl(), QwHelicityPattern::FillDB(), QwSubsystemArrayParity::FillHistograms(), QwHelicityPattern::FillHistograms(), QwEventBuffer::FillSubsystemConfigurationData(), QwEventBuffer::FillSubsystemData(), QwSubsystemArrayParity::FillTreeVector(), QwHelicityPattern::FillTreeVector(), QwParityDB::GetAnalysisID(), getenv_safe_string(), getenv_safe_TString(), QwEventBuffer::GetEvent(), QwEventBuffer::GetEventNumber(), QwHelicity::GetHelicityActual(), QwHelicity::GetHelicityDelayed(), QwOptions::GetIntValuePairFirst(), QwOptions::GetIntValuePairLast(), QwParityDB::GetLumiDetectorID(), QwParityDB::GetMainDetectorID(), QwParityDB::GetMonitorID(), QwHelicity::GetPhaseNumber(), QwHelicity::GetRandomSeedActual(), QwHelicity::GetRandomSeedDelayed(), QwParityDB::GetRunID(), QwEventBuffer::GetRunNumber(), QwSubsystemArrayParity::GetSubsystemByName(), QwDatabase::GetVersion(), gQwHists, gQwOptions, QwSubsystemArrayParity::IncrementErrorCounters(), QwHelicityPattern::IsCompletePattern(), QwEventBuffer::IsPhysicsEvent(), QwEventBuffer::IsROCConfigurationEvent(), Qw::kBoldMagenta, kMultiplet, Qw::kNormal, QwHelicityPattern::LoadEventData(), QwHistogramHelper::LoadHistParamsFromFile(), QwEventBuffer::OpenDataFile(), QwSubsystemArray::ProcessEvent(), QwHelicityPattern::ProcessOptions(), QwSubsystemArray::ProcessOptions(), QwMessage, MQwCodaControlEvent::ReportRunSummary(), QwOptions::SetCommandLine(), and QwOptions::SetConfigFile().

{
  // First, we set the command line arguments and the configuration filename,
  // and we define the options that can be used in them (using QwOptions).
  gQwOptions.SetCommandLine(argc, argv);
  //  gQwOptions.SetConfigFile("qwmockdataanalysis.conf");
  gQwOptions.SetConfigFile("Parity/prminput/qwanalysis_beamline.conf");
  gQwOptions.SetConfigFile("Parity/prminput/qweak_mysql.conf");
  // Define the command line options
  DefineOptionsParity(gQwOptions);

  ///  Fill the search paths for the parameter files; this sets a static
  ///  variable within the QwParameterFile class which will be used by
  ///  all instances.
  ///  The "scratch" directory should be first.
  QwParameterFile::AppendToSearchPath(getenv_safe_string("QW_PRMINPUT"));
  QwParameterFile::AppendToSearchPath(getenv_safe_string("QWANALYSIS")+"/Parity/prminput");
  QwParameterFile::AppendToSearchPath(getenv_safe_string("QWANALYSIS") + "/Analysis/prminput");

  // Load histogram definitions
  gQwHists.LoadHistParamsFromFile("qweak_parity_hists.in");

  // Detector array
  QwSubsystemArrayParity detectors(gQwOptions);
  detectors.ProcessOptions(gQwOptions);

  // Helicity pattern
  QwHelicityPattern helicitypattern(detectors);
  helicitypattern.ProcessOptions(gQwOptions);

  // Running sum
  QwSubsystemArrayParity runningsum(detectors);

  // Get the helicity
  QwHelicity* helicity = dynamic_cast<QwHelicity*>(detectors.GetSubsystemByName("Helicity info"));

  // Event buffer
  QwEventBuffer eventbuffer;


  // Loop over all runs
  UInt_t runnumber_min = (UInt_t) gQwOptions.GetIntValuePairFirst("run");
  UInt_t runnumber_max = (UInt_t) gQwOptions.GetIntValuePairLast("run");
  for (UInt_t run  = runnumber_min;
              run <= runnumber_max;
              run++) {

    // Data file (input)
    TString datafilename = TString("QwMock_") + Form("%d.log",run);
    if (eventbuffer.OpenDataFile(datafilename,"R") != CODA_OK) {
      std::cout << "Error: could not open file!" << std::endl;
      return 0;
    }


    // ROOT file output (histograms)
    TString rootfilename = getenv_safe_TString("QW_ROOTFILES") + TString("/QwMock_") + Form("%d.root",run);
    TFile rootfile(rootfilename, "RECREATE", "QWeak ROOT file");
    if (bHisto) {
      rootfile.cd();
      detectors.ConstructHistograms(rootfile.mkdir("mps_histo"));
      if (bHelicity) {
        rootfile.cd();
        helicitypattern.ConstructHistograms(rootfile.mkdir("hel_histo"));
      }
    }

    // ROOT file output (trees)
    TTree *mpstree;
    TTree *heltree;
    Int_t eventnumber;
    std::vector <Double_t> mpsvector;
    std::vector <Double_t> helvector;
    if (bTree) {
      rootfile.cd();
      mpstree = new TTree("MPS_Tree","MPS event data tree");
      mpsvector.reserve(6000);
      mpstree->Branch("eventnumber",&eventnumber,"eventnumber/F");
      TString dummystr="";
      detectors.ConstructBranchAndVector(mpstree, dummystr, mpsvector);
      rootfile.cd();
      if (bHelicity) {
        rootfile.cd();
        heltree = new TTree("HEL_Tree","Helicity event data tree");
        helvector.reserve(6000);
        TString dummystr="";
        helicitypattern.ConstructBranchAndVector(heltree, dummystr, helvector);
      }
    }


    // Loop over events in this CODA file
    Int_t eventnumber_min = gQwOptions.GetIntValuePairFirst("event");
    Int_t eventnumber_max = gQwOptions.GetIntValuePairLast("event");
    while (eventbuffer.GetEvent() == CODA_OK) {

      // First, do processing of non-physics events...
      if (eventbuffer.IsROCConfigurationEvent()) {
        // Send ROC configuration event data to the subsystem objects.
        eventbuffer.FillSubsystemConfigurationData(detectors);
      }

      // Now, if this is not a physics event, go back and get a new event.
      if (! eventbuffer.IsPhysicsEvent()) continue;

      //  Check to see if we want to process this event.
      eventnumber = eventbuffer.GetEventNumber();
      if      (eventnumber < eventnumber_min) continue;
      else if (eventnumber > eventnumber_max) break;

      // Fill the subsystem objects with their respective data for this event.
      eventbuffer.FillSubsystemData(detectors);

      // Process this events
      detectors.ProcessEvent();
      detectors.IncrementErrorCounters();

      // Helicity pattern
      if (bHelicity)
        helicitypattern.LoadEventData(detectors);

      // Print the helicity information
      if (bHelicity && bDebug) {
        // - actual helicity
        if      (helicity->GetHelicityActual() == 0) std::cout << "-";
        else if (helicity->GetHelicityActual() == 1) std::cout << "+";
        else std::cout << "?";
        // - delayed helicity
        if      (helicity->GetHelicityDelayed() == 0) std::cout << "(-) ";
        else if (helicity->GetHelicityDelayed() == 1) std::cout << "(+) ";
        else std::cout << "(?) ";
        if (helicity->GetPhaseNumber() == kMultiplet) {
          std::cout << std::hex << helicity->GetRandomSeedActual() << std::dec << ",  \t";
          std::cout << std::hex << helicity->GetRandomSeedDelayed() << std::dec << std::endl;
        }
      }

      // Check for helicity validity (TODO I'd prefer to use kUndefinedHelicity)
      if (bHelicity && helicity->GetHelicityDelayed() == -9999) continue;

      // Fill the histograms
      if (bHisto) detectors.FillHistograms();

      // Accumulate the running sum to calculate the event based running average
      runningsum.AccumulateRunningSum(detectors);

      // Fill the MPS tree
      if (bTree) {
        eventnumber = eventbuffer.GetEventNumber();
        detectors.FillTreeVector(mpsvector);
        mpstree->Fill();
      }
      // Fill the helicity tree
      if (bHelicity && helicitypattern.IsCompletePattern()) {
        helicitypattern.CalculateAsymmetry();
        if (bHisto) helicitypattern.FillHistograms();
        if (bTree) {
          helicitypattern.FillTreeVector(helvector);
          heltree->Fill();
        }
        helicitypattern.ClearEventData();
      }

      // Periodically print event number
      if ((bDebug && eventbuffer.GetEventNumber() % 1000 == 0)
                  || eventbuffer.GetEventNumber() % 10000 == 0)
        std::cout << "Number of events processed so far: "
                  << eventbuffer.GetEventNumber() << std::endl;

      // TODO (wdc) QwEventBuffer should have Bool_t AtEndOfBurst()
      //if (QwEvt.AtEndOfBurst()){
      if (eventbuffer.GetEventNumber() % 1000 == 0) {
        {
          helicitypattern.AccumulateRunningBurstSum();
          helicitypattern.CalculateBurstAverage();
          helicitypattern.ClearBurstSum();
        }
      }

    } // end of loop over events

    // Calculate the running averages
    helicitypattern.CalculateRunningAverage();
    helicitypattern.CalculateRunningBurstAverage();
    runningsum.CalculateRunningAverage();

    // Close ROOT file
    rootfile.Write(0,TObject::kOverwrite);

    // Close data file and print run summary
    eventbuffer.CloseDataFile();
    eventbuffer.ReportRunSummary();

    // Write to database
    if (bDatabase) {
      QwParityDB* qweak_database  = new QwParityDB();
      QwMessage << "GetMonitorID(qwk_batext2) = " << qweak_database->GetMonitorID("qwk_batext2") << QwLog::endl;
      QwMessage << "GetMonitorID(phasemonitor) = " << qweak_database->GetMonitorID("phasemonitor") << QwLog::endl;
      QwMessage << "GetMonitorID(qwk_junk) = " << qweak_database->GetMonitorID("qwk_junk") << QwLog::endl;
      QwMessage << "GetMainDetectorID(md1neg) = " << qweak_database->GetMainDetectorID("md1neg") << QwLog::endl;
      QwMessage << "GetMainDetectorID(spare3) = " << qweak_database->GetMainDetectorID("spare3") << QwLog::endl;
      QwMessage << "GetMainDetectorID(combinationallmd) = " << qweak_database->GetMainDetectorID("combinationallmd") << QwLog::endl;
      QwMessage << "GetLumiDetectorID(dlumi8) = " << qweak_database->GetLumiDetectorID("dlumi8") << QwLog::endl;
      QwMessage << "GetVersion() = " << qweak_database->GetVersion() << QwLog::endl;
      // GetRunID() and GetAnalysisID have their own Connect() and Disconnect() functions.
      UInt_t run_id      = qweak_database->GetRunID(eventbuffer);
      UInt_t analysis_id = qweak_database->GetAnalysisID(eventbuffer);

      QwMessage << "QwMockDataAnalysis.cc::"
                << " Run Number "  << QwColor(Qw::kBoldMagenta) << eventbuffer.GetRunNumber() << QwColor(Qw::kNormal)
                << " Run ID "      << QwColor(Qw::kBoldMagenta) << run_id<< QwColor(Qw::kNormal)
                << " Analysis ID " << QwColor(Qw::kBoldMagenta) << analysis_id
                << QwLog::endl;

      // Each sussystem has its own Connect() and Disconnect() functions.
      helicitypattern.FillDB(qweak_database);

      delete qweak_database; qweak_database = NULL;

    } // end of database write


  } // end of loop over runs


}

Variable Documentation

bool bDatabase = false

static

Definition at line 50 of file QwMockDataAnalysis.cc.

Referenced by main().

bool bDebug = false

static

Definition at line 44 of file QwMockDataAnalysis.cc.

Referenced by main().

bool bHelicity = true

static

Definition at line 49 of file QwMockDataAnalysis.cc.

Referenced by main().

bool bHisto = true

static

Definition at line 48 of file QwMockDataAnalysis.cc.

Referenced by main().

bool bTree = true

static

Definition at line 47 of file QwMockDataAnalysis.cc.

Referenced by main().

const int kMultiplet = 4

static

Definition at line 41 of file QwMockDataAnalysis.cc.

Referenced by main().