QwCompton.cc File Reference

Compton data analysis. More...

#include <iostream>
#include <boost/random.hpp>
#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include "QwLog.h"
#include "QwRootFile.h"
#include "QwOptionsParity.h"
#include "QwEventBuffer.h"
#include "QwEPICSEvent.h"
#include "QwHelicity.h"
#include "QwHelicityPattern.h"
#include "QwHistogramHelper.h"
#include "QwSubsystemArrayParity.h"
#include "QwScaler.h"
#include "QwBeamLine.h"
#include "QwComptonPhotonDetector.h"
#include "QwComptonElectronDetector.h"

Include dependency graph for QwCompton.cc:

Go to the source code of this file.

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

Detailed Description

Compton data analysis.

Definition in file QwCompton.cc.

Function Documentation

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

Change some default settings for use by the Compton analyzer

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.

Set the command line arguments and the configuration filename

Define the command line options

Message logging facilities

Create the event buffer

Start loop over all runs

Begin processing for the first run

Set the current event number for parameter file lookup

Create an EPICS event

Load the detectors from file

Create the helicity pattern

Start loop over events

Definition at line 103 of file QwCompton.cc.

References QwHelicityPattern::AccumulateRunningBurstSum(), QwParameterFile::AppendToSearchPath(), QwSubsystemArrayParity::ApplySingleEventCuts(), QwHelicityPattern::CalculateAsymmetry(), QwHelicityPattern::CalculateBurstAverage(), QwHelicityPattern::CalculateRunningAverage(), QwHelicityPattern::CalculateRunningBurstAverage(), QwSubsystemArrayParity::CheckForEndOfBurst(), QwHelicityPattern::ClearBurstSum(), QwHelicityPattern::ClearEventData(), QwHelicityPattern::ClearRunningSum(), QwRootFile::Close(), QwEventBuffer::CloseStream(), QwRootFile::ConstructHistograms(), QwRootFile::ConstructIndices(), QwRootFile::ConstructTreeBranches(), DefineOptionsParity(), QwLog::endl(), QwEventBuffer::FillEPICSData(), QwRootFile::FillHistograms(), QwEventBuffer::FillSubsystemConfigurationData(), QwEventBuffer::FillSubsystemData(), QwRootFile::FillTree(), QwRootFile::FillTreeBranches(), QwHelicityPattern::GetBurstAsymmetry(), QwHelicityPattern::GetBurstDifference(), QwHelicityPattern::GetBurstYield(), getenv_safe_string(), QwEventBuffer::GetEventNumber(), QwEventBuffer::GetNextEvent(), QwEventBuffer::GetRunLabel(), QwEventBuffer::GetRunNumber(), gQwHists, gQwLog, gQwOptions, QwHelicityPattern::IsBurstSumEnabled(), QwHelicityPattern::IsCompletePattern(), QwHelicityPattern::IsEndOfBurst(), QwEventBuffer::IsEPICSEvent(), QwHelicityPattern::IsGoodAsymmetry(), QwEventBuffer::IsPhysicsEvent(), QwEventBuffer::IsROCConfigurationEvent(), QwHelicityPattern::IsRunningSumEnabled(), QwHelicityPattern::LoadEventData(), QwHistogramHelper::LoadHistParamsFromFile(), QwEventBuffer::OpenNextStream(), QwRootFile::PrintDirs(), QwEventBuffer::PrintRunTimes(), QwRootFile::PrintTrees(), QwSubsystemArray::ProcessEvent(), QwEPICSEvent::ProcessOptions(), QwHelicityPattern::ProcessOptions(), QwEventBuffer::ProcessOptions(), QwLog::ProcessOptions(), QwSubsystemArray::ProcessOptions(), QwError, QwMessage, MQwCodaControlEvent::ReportRunSummary(), QwOptions::SetCommandLine(), QwParameterFile::SetCurrentRunNumber(), QwEventBuffer::SetDefaultDataFileStem(), QwRootFile::SetDefaultRootFileStem(), QwRootFile::Write(), and QwRootFile::WriteParamFileList().

{
  /// Change some default settings for use by the Compton analyzer
  QwEventBuffer::SetDefaultDataFileStem("Compton_");
  QwRootFile::SetDefaultRootFileStem("Compton_");


  /// 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");


  /// Set the command line arguments and the configuration filename
  gQwOptions.SetCommandLine(argc, argv);
  /// Define the command line options
  DefineOptionsParity(gQwOptions);

  /// Message logging facilities
  gQwLog.ProcessOptions(&gQwOptions);


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


  ///  Create the event buffer
  QwEventBuffer eventbuffer;
  eventbuffer.ProcessOptions(gQwOptions);

  ///  Start loop over all runs
  while (eventbuffer.OpenNextStream() == CODA_OK) {

    ///  Begin processing for the first run


    ///  Set the current event number for parameter file lookup
    QwParameterFile::SetCurrentRunNumber(eventbuffer.GetRunNumber());


    ///  Create an EPICS event
    QwEPICSEvent epicsevent;
    epicsevent.ProcessOptions(gQwOptions);
    epicsevent.LoadChannelMap("compton_epics_table.map");


    ///  Load the detectors from file
    QwSubsystemArrayParity detectors(gQwOptions);
    detectors.ProcessOptions(gQwOptions);

    ///  Create the helicity pattern
    QwHelicityPattern helicitypattern(detectors);
    helicitypattern.ProcessOptions(gQwOptions);

    // Create a configuration options output
    QwSubsystemArrayParity configoutput(detectors);


    //  Open the ROOT file
    QwRootFile* rootfile = new QwRootFile(eventbuffer.GetRunLabel());
    if (! rootfile) QwError << "QwAnalysis made a boo boo!" << QwLog::endl;
    rootfile->WriteParamFileList("mapfiles", detectors);
    //  Construct histograms
    rootfile->ConstructHistograms("mps_histo", detectors);
    rootfile->ConstructHistograms("hel_histo", helicitypattern);

    //  Construct tree branches
    rootfile->ConstructTreeBranches("Mps_Tree", "MPS event data tree", detectors);
    rootfile->ConstructTreeBranches("Hel_Tree", "Helicity event data tree", helicitypattern);
    rootfile->ConstructTreeBranches("Burst_Tree", "Burst level data tree", helicitypattern.GetBurstYield(),"yield_");
    rootfile->ConstructTreeBranches("Burst_Tree", "Burst level data tree", helicitypattern.GetBurstAsymmetry(),"asym_");
    rootfile->ConstructTreeBranches("Burst_Tree", "Burst level data tree", helicitypattern.GetBurstDifference(),"diff_");
    rootfile->ConstructTreeBranches("Slow_Tree", "EPICS and slow control tree", epicsevent);
    rootfile->ConstructTreeBranches("Config_Tree","Configuration Tree",configoutput,"conf");
    // Construct indices to get from one tree to the other
    rootfile->ConstructIndices("Mps_Tree","Slow_Tree");
    rootfile->ConstructIndices("Hel_Tree","Slow_Tree");
    rootfile->ConstructIndices("Mps_Tree","Config_Tree");
    rootfile->ConstructIndices("Hel_Tree","Config_Tree");

    // Summarize the ROOT file structure
    rootfile->PrintTrees();
    rootfile->PrintDirs();


    //  Clear the running sums
    helicitypattern.ClearRunningSum();
    helicitypattern.ClearBurstSum();


    ///  Start loop over events
    while (eventbuffer.GetNextEvent() == CODA_OK) {

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

      //  Secondly, process EPICS events
      if (eventbuffer.IsEPICSEvent()) {
        eventbuffer.FillEPICSData(epicsevent);
        if (epicsevent.HasDataLoaded()) {
          epicsevent.CalculateRunningValues();

          rootfile->FillTreeBranches(epicsevent);
          rootfile->FillTree("Slow_Tree");
        }
      }

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



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

      // Process this events
      detectors.ProcessEvent();

      // The event pass the event cut constraints
      if (detectors.ApplySingleEventCuts()) {

        // Fill the histograms
        rootfile->FillHistograms(detectors);

        // Fill the tree
        rootfile->FillTreeBranches(detectors);
        rootfile->FillTree("Mps_Tree");


        //  Somehow test to see if the laser state has changed.  If it has,
        //  then calculate and print the running sums so far.
        //  Maybe instead of here, it should be done as soon as we've picked up an
        //  event, instead of waiting to build a complete pattern.

        if (detectors.CheckForEndOfBurst()){
          //  Do the burst versions of
          helicitypattern.CalculateRunningAverage();
          // Maybe we need to clear the burst sums also
        }


        // Helicity pattern
        helicitypattern.LoadEventData(detectors);


        // TODO We need another check here to test for pattern validity.  Right
        // now the first 24 cycles are also added to the histograms.
        if (helicitypattern.IsCompletePattern()) {

          // Calculate the asymmetry
          helicitypattern.CalculateAsymmetry();
          if (helicitypattern.IsGoodAsymmetry()) {

            // Fill histograms
            rootfile->FillHistograms(helicitypattern);

            // Fill tree branches
            rootfile->FillTreeBranches(helicitypattern);
            rootfile->FillTree("Hel_Tree");

            // Clear the data
            helicitypattern.ClearEventData();
          }


          //  Somehow test to see if the laser state has changed.  If it has,
          //  then calculate and print the running sums so far.
          //  Maybe instead of here, it should be done as soon as we've picked up an
          //  event, instead of waiting to build a complete pattern.

        }

      } // detectors.ApplySingleEventCuts()


      // Burst mode
      if (helicitypattern.IsEndOfBurst()) {
        helicitypattern.CalculateBurstAverage();
        helicitypattern.AccumulateRunningBurstSum();

        // Fill tree branches
        rootfile->FillTreeBranches(helicitypattern.GetBurstYield());
        rootfile->FillTreeBranches(helicitypattern.GetBurstAsymmetry());
        rootfile->FillTreeBranches(helicitypattern.GetBurstDifference());
        rootfile->FillTree("Burst_Tree");

        // Clear the data
        helicitypattern.ClearBurstSum();
      }

    } // end of loop over events
    rootfile->FillTreeBranches(configoutput);
    rootfile->FillTree("Config_Tree");

    QwMessage << "Last event processed: "
              << eventbuffer.GetEventNumber() << QwLog::endl;

    // Calculate running averages over helicity patterns
    if (helicitypattern.IsRunningSumEnabled()) {
      helicitypattern.CalculateRunningAverage();
      if (helicitypattern.IsBurstSumEnabled()) {
        helicitypattern.CalculateRunningBurstAverage();
      }
    }


    /*  Write to the root file, being sure to delete the old cycles  *
     *  which were written by Autosave.                              *
     *  Doing this will remove the multiple copies of the ntuples    *
     *  from the root file.                                          *
     *                                                               *
     *  Then, we need to delete the histograms here.                 *
     *  If we wait until the subsystem destructors, we get a         *
     *  segfault; but in additional to that we should delete them     *
     *  here, in case we run over multiple runs at a time.           */
    rootfile->Write(0,TObject::kOverwrite);

    // Close ROOT file
    rootfile->Close();
    delete rootfile; rootfile = 0;

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

    //  Report run summary
    eventbuffer.ReportRunSummary();
    eventbuffer.PrintRunTimes();

  } // end of loop over runs

  QwMessage << "I have done everything I can do..." << QwLog::endl;

  return 0;
}