UserAnalysis.cxx

Go to the documentation of this file.

// $Id: UserAnalysis.cxx,v 1.29 2006/05/27 08:29:40 rhatcher Exp $
//
// A prototype user analysis module. This is meant to show how to make
// a simple analysis module to look at the data and serve as a boiler
// plate for more complicated analysis modules
//
// messier@huhepl.harvard.edu
#include "Demo/UserAnalysis.h"
#include <cstdio>
#include <cmath>
// ROOT includes
#include "TFile.h"
// MINOS includes
#include "MessageService/MsgService.h"     // MSG text output
#include "MinosObjectMap/MomNavigator.h"   // Data access
#include "CandData/CandRecord.h"           // DigitList handling
#include "CandDigit/CandDigitListHandle.h" // "                "
#include "CandDigit/CandDigitHandle.h"     // "                "
#include "Conventions/Munits.h"            // Unit conversions
#include "DataUtil/GetCandidate.h"         // Unpack data from Mom
#include <DataUtil/HistManModule.h>        // Easy access to our unique HistMan
#include "JobControl/JobCommand.h"         // JobCommand handling
#include "JobControl/JobCModuleRegistry.h" // JOBMODULE registration macro

CVSID("$Id: UserAnalysis.cxx,v 1.29 2006/05/27 08:29:40 rhatcher Exp $");
JOBMODULE(UserAnalysis,"UserAnalysis","Compute and record total charge\n");

//......................................................................

// Define the data structure for the ntuple
//  This structure defines what data is in the ntuple
typedef struct user_analysis_ntuple {
  Float_t qtot; // Total charge (p.e)
  Float_t nhit; // Total number of digitizations
} UserAnalysisNtuple_t;
//  This string sets the names of the ntuple variables
const char* gsChTags = "qtot:nhit";

//......................................................................

UserAnalysis::UserAnalysis() :
  // Classes for drawing
  fCanvas(std::tmpnam(0),"User Analysis",600,600),
  // Histogram convenience pointers
  fQtotHisto(0),
  fNhitHisto(0),
  fEvtQHisto(0),
  fEvtTHisto(0),
  // Ntuple
  fNtuple("fNtuple","Event Summary",gsChTags),
  // Module data
  fNevents(0), 
  fSumQ(0.0), 
  fSumQsqr(0.0)
{
}
void UserAnalysis::BeginJob()
{
    HistMan hm = HistManModule::GetHistMan(*this);

    // Histograms
    //         Name of object Histogram title # bins Low  High
    //         ============== =============== ====== ==== ======
    fQtotHisto = hm.Book<TH1F>("fQtotHisto",  "Total Charge", 50,  0.0, 5000.0);
    fNhitHisto = hm.Book<TH1F>("fNhitHisto",  "N?Hit!",       50,  0.0, 500.0);
    fEvtTHisto = hm.Book<TH1F>("fEvtQHisto",  "Event Charge", 50,  0.0, 50.0);

    // Check for user tweaked timing histogram parameters.  First
    // start by setting the defaults.
    int    nbin = 100;
    double xmin = 0.0, xmax = 200.0;

    // See if the user wants to change them.  Normally one might want
    // to check the return value of these Get()s to see if the key
    // exists, however, since we are guaranteed that the values are
    // unchanged in that case and we set defaults, we are safe.
    Registry& r = this->GetConfig();
    r.Get("NbinT",nbin);
    r.Get("Tmin", xmin);
    r.Get("Tmax", xmax);

    // Finally book the histogram with the given parameters
    fEvtQHisto = hm.Book<TH1F>("fEvtTHisto",  "Event Times",  nbin, xmin, xmax);

    // Alway set clear axis titles!!!
    fQtotHisto->SetXTitle("Total Charge (p.e.)");
    fQtotHisto->SetYTitle("Events");
    fNhitHisto->SetXTitle("N?Hit!");
    fNhitHisto->SetYTitle("Events");
    fEvtQHisto->SetXTitle("Charge (p.e)");
    fEvtQHisto->SetYTitle("N?stripends!");
    fEvtTHisto->SetXTitle(Form("Time (%s)",Munits::base_time_name));
    fEvtTHisto->SetYTitle("Charge (p.e)");

//======================================================================
// Set the initial state of your module. The code above initializes the
// variables listed to the values in parentheses. In this case that's
// all I need to do so the method itself is empty...
//======================================================================
    fCanvas.Divide(2,2); // Like PAW's zone command

}

//......................................................................

double UserAnalysis::GetTotalCharge(const CandDigitListHandle* cdlh, 
                                    int* nHit) 
{
//======================================================================
// Calculate the total charge and number of hit strips in an event
//
// Inputs:  cdlh - Pointer to a candidate digit list handle
//
// Outputs: nHit - total number of digits in the list
//======================================================================
  int    nhit = 0;   // Number of hit strips
  double qtot = 0.0; // Total charge of this event in p.e.

  fEvtTHisto->Reset(); // Reset event histograms
  fEvtQHisto->Reset(); // Reset event histograms
  
  // The loop is done with a class that iterates over the list
  CandDigitHandleItr cdhItr(cdlh->GetDaughterIterator());
  for (; cdhItr.IsValid(); cdhItr.Next()) {
    double t; // Time of single digit
    double q; // Charge of single digit

    t = (*cdhItr)->GetTime(); // time 
    q = (*cdhItr)->GetCharge(CalDigitType::kPE); // charge in pe's
  
    // If this digit passes cuts then include it in tallies
    if (q > fQminCut) {
      fEvtTHisto->Fill(t/Munits::ns,q); // Fill histograms (time weighted by charge)
      fEvtQHisto->Fill(q);   // Fill histograms
      
      qtot += q; // Accumulate total charge
      ++nhit;    // Accumulate total number of hits
    }
  }
  // If we've been give a pointer for the nhit variable put nhit there
  if (nHit) *nHit = nhit;
  return qtot;
}

//......................................................................

void UserAnalysis::HandleCommand(JobCommand* cmd)
{
//======================================================================
// Handle a command send from the job controller
//
// ** This method should no longer be used. Use Config() method instead **
// 
// Inputs: cmd - The parsed job command
//======================================================================
  if (cmd->HaveCmd()) {                      // If we have a command...
    string sc = cmd->PopCmd();               // Get the command
    if (sc == "Draw") {                      // If the command is "Draw"
      fDrawHistos = true;                    // Default is on
      if (cmd->HaveOpt()) {                  // If we have options
        string sopt = cmd->PopOpt();         // Get the option
        // If the option is "Off" turn drawing off
        if (sopt=="Off" || sopt=="off" || sopt=="OFF") fDrawHistos = false;
      }
      return;
    }
    if (sc == "WriteHistos") {
      fWriteHistos = true;
      return;
    }
    if (sc == "Set") {                   // "Set" command
      if (cmd->HaveOpt()) {
        string sopt = cmd->PopOpt();
        if (sopt == "QminCut") {         // Set Qmin cut
          fQminCut = cmd->PopFloatOpt(); // Get cut value
          return;
        }
        if (sopt == "tBins") { // Set Time Range
          int    nbin = cmd->PopIntOpt();
          double xmin = cmd->PopFloatOpt();
          double xmax = cmd->PopFloatOpt();
          fEvtTHisto->SetBins(nbin,xmin,xmax);
          return;
        }
      }
    }
    // Don't understand command
    MSG("Demo",Msg::kWarning) << "Don't understand '" << sc.c_str() << "'\n";
  }
}

//......................................................................

void UserAnalysis::Help() 
{
//======================================================================
// ** No longer needed. Use Module::Report() to see list of parameters
// ** which can be set using Set("name=value").
//
// Print help for this module
//======================================================================
  MSG("Demo",Msg::kInfo) << 
    "Help For UserAnalysis Module:\n" <<
    "\n" <<
    " UserAnalysis accepts the following commands:\n" <<
    "\n" <<
    "  Draw [off]  -- Turn drawing histograms on (default) or off\n" <<
    "\n"
    "  WriteHistos -- Write histogram file at end of run\n"<<
    "                 Default is not to write histograms\n"<<
    "\n"
    "  Set/QminCut [charge]            -- Set minimum charge cuts\n" <<
    "     /tBins [nbins] [tmin] [tmax] -- Set binning of time histogram\n";
}

//......................................................................

JobCResult UserAnalysis::Ana(const MomNavigator *mom) 
{
//======================================================================
// Given the object to which the data is attached perform some
// analysis of the data
//======================================================================
  // Access to the event is given through handles. In this case I want
  // the list of candidate digits
  const CandDigitListHandle* canddigit = 
    DataUtil::GetCandidate<CandDigitListHandle>(mom,
                                                "CandDigitListHandle", 
                                                "canddigitlist");
  
  // If the GetDigitList is successful, loop over the elements in the
  // list to add up the charge
  if (canddigit) {
    double qtot;                 // Total charge
    int    nhit;                 // Total number of hit strips
    UserAnalysisNtuple_t ntuple; // The data structure for the ntuple data
    
    // Get the total charge and Nhit and put it into the ntuple data
    // structure
    qtot = this->GetTotalCharge(canddigit,&nhit);
    
    // Accumulate statistics
    fQtotHisto->Fill(qtot);
    fNhitHisto->Fill(nhit);
    fSumQ    += qtot;
    fSumQsqr += qtot*qtot;
    ++fNevents;

    // Fill the ntuple
    ntuple.qtot = qtot;
    ntuple.nhit = nhit;
    fNtuple.Fill(&ntuple.qtot);

    // If graphics are enabled, draw the histograms
    if (fDrawHistos) this->DrawHistograms();

    // Make a filter decision on this event.
    if (ntuple.qtot>=fQtotCut) return JobCResult::kPassed; // Passed cut
    else                       return JobCResult::kFailed; // Failed cut
  }
  else {
    // Couldn't find a digit list...
    MSG("User",Msg::kError) << "No CandDigitList\n";
//    return (JobCResult::kError+JobCResult::;
    return (JobCResult::kError);
  }
  // If we get here somehow, don't make any pass/fail decision on the event
  return JobCResult::kAOK;
}

//......................................................................

const Registry& UserAnalysis::DefaultConfig() const
{
//======================================================================
// Create a registry which holds the default configuration and return it
//======================================================================
  static Registry r;

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values of config
  r.UnLockValues();
  r.Set("QminCut",    0.0);
  r.Set("QtotCut",    400.0);
  r.Set("DrawHistos", true);
  r.Set("WriteHistos",true);
  r.Set("NbinT",      200);
  r.Set("Tmin",       0.0);
  r.Set("Tmax",       200.0);
  r.LockValues();

  return r;
}

//......................................................................

void UserAnalysis::Config(const Registry& r)
{
//======================================================================
// Configure the module given the registry r
//======================================================================
//  char   tmpb;
  double tmpd;
  int    tmpi;

  if (r.Get("DrawHistos", tmpi)) { fDrawHistos  = tmpi; }
  if (r.Get("WriteHistos",tmpi)) { fWriteHistos = tmpi; }
  if (r.Get("QminCut",    tmpd)) { fQminCut     = tmpd; }
  if (r.Get("QtotCut",    tmpd)) { fQtotCut     = tmpd; }

  // Let parent class handle stashing any other config parameters.
  this->JobCModule::Config(r);

}

//......................................................................

void UserAnalysis::DrawHistograms() 
{
//======================================================================
// Draw the histograms
//======================================================================
  fCanvas.cd(1); gPad->SetLogy(1); fNhitHisto->DrawCopy();
  fCanvas.cd(2); gPad->SetLogy(1); fQtotHisto->DrawCopy();
  fCanvas.cd(3); gPad->SetLogy(0); fEvtQHisto->DrawCopy();
  fCanvas.cd(4); gPad->SetLogy(0); fEvtTHisto->DrawCopy();
  fCanvas.Update();
}

//......................................................................

void UserAnalysis::EndJob() 
{
//======================================================================
// At the end of the job print some stuff out and save the histogram
// to a file
//======================================================================
  double aveQ   = fNevents ? fSumQ/fNevents : 0;
  double sigmaQ = fNevents ? sqrt(fabs(fSumQsqr/fNevents - aveQ*aveQ)) : 0;
  MSG("User",Msg::kInfo) << "Total Events: " << fNevents << "\n";
  MSG("User",Msg::kInfo) << " Average Q=" << aveQ 
                         << " SigmaQ=" << sigmaQ << "\n";
  
  // Create a file to write the histogram to. By default it becomes
  // the current directory. Write uses current directory by default
  if (fWriteHistos) {
    TFile userFile("useranalysis.root","RECREATE");
    fQtotHisto->Write();        // also can be written by using
    fNhitHisto->Write();        // HistManModule in path
    fNtuple.Write();
  }
}

---