sc_ana_module.cxx

Go to the documentation of this file.

/********************************************************************\

  Name:         sc_ana_module.cxx
  Created by:   Thomas Prokscha
  Date:         26-Apr-2006

  Contents:     Analyzer module for slow control histograms

\********************************************************************/

/*-- Include files -------------------------------------------------*/

/* standard includes */
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <time.h>
#include <string.h>
#include <math.h>

#include <vector>

/* midas includes */
#include "midas.h"
#include "experim.h"
#include "nemu_experim.h"

#ifdef HAVE_HBOOK
#include <cfortran.h>
#include <hbook.h>
#else
/* root includes */
#include <TFolder.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TTree.h>
#endif

#define SC_MOD_CF           0
#define SC_MOD_XTC          1
#define SC_MOD_HV           2
#define SC_SAMPLE_CF        3
#define SC_PRESSURE_SC_GJ   4
#define SC_SAMPLE_HV        5
#define SC_SAMPLE_ZERO_FLUX 6

/*-- ODB Parameters ------------------------------------------------*/
extern INFO             info;
extern RUNINFO          runinfo;
extern TRIGGER_SETTINGS trigger_settings;
MEAN_BANK mean;
SCANAMODULE_PARAM       sc_ana_param;
SCANAMODULE_PARAM_STR(sc_ana_param_str);
extern float get_magnetic_field(float current);

/*-- Module declaration --------------------------------------------*/
INT sc_init(void);
INT sc_exit();
INT sc_ana(EVENT_HEADER*,void*);
INT sc_bor(INT run_number);
INT sc_eor(INT run_number);

ANA_MODULE sc_ana_module = {
  "SCAnaModule",              /* module name           */
  "T. Prokscha, A. Suter",    /* author                */
  sc_ana,                     /* event routine         */
  sc_bor,                     /* BOR routine           */
  sc_eor,                     /* EOR routine           */
  sc_init,                    /* init routine          */
  sc_exit,                    /* exit routine          */
  &sc_ana_param,              /* parameter structure   */
  sizeof(sc_ana_param),       /* structure size        */
  (const char **)sc_ana_param_str, /* initial parameters    */
};

/*-- other global variables ----------------------------------------*/

std::vector< std::vector<double> > sc_histo_vec;

#define N_SC_HIST  8
#ifndef HAVE_HBOOK
static TH1F* hSC_Hist[N_SC_HIST];
//the following objects are from mana.c
extern TFolder   *gManaHistosFolder;
extern TObjArray *gHistoFolderStack;
#endif

static INT offset;

//-- INIT routine --------------------------------------------------
INT sc_init(void)
{
  sc_histo_vec.clear(); // clean up - paranoia
  sc_histo_vec.resize(N_SC_HIST); // resize it to the proper number of SC histograms

  return SUCCESS;
}

//-- EXIT routine ---------------------------------------------------

INT sc_exit()
{
  // clean up
  sc_histo_vec.clear();

  return SUCCESS;
}

//-- BOR routine ---------------------------------------------------

INT sc_bor(INT run_number)
{
  HNDLE hDB, hKey;
  INT status, size, ival, success;
  char str[32];

  sc_histo_vec.clear(); // clean up - paranoia
  sc_histo_vec.resize(N_SC_HIST); // resize it to the proper number of SC histograms

  // check the ctrl channel of the sample cryostat and set the offset for the analyzer module
  cm_get_experiment_database(&hDB, NULL);
  status = db_find_key(hDB, 0, "/Equipment/SampleCryo/Settings/Devices/Lake336_Sample_0/DD/Loop1/CTRL_CH", &hKey);
  if (status != DB_SUCCESS) {
    cm_msg(MINFO, "sc_bor", "analyzer/sc_ana_module/sc_bor: couldn't obtain sample cryo ctrl channel key.");
  } else {
    size = sizeof(str);
    success = db_get_data(hDB, hKey, str, &size, TID_STRING);
    if (status != DB_SUCCESS) {
      cm_msg(MINFO, "sc_bor", "analyzer/sc_ana_module/sc_bor: couldn't obtain sample cryo ctrl channel value");
    } else {
      success = db_find_key(hDB, 0, "/Analyzer/Parameters/SCAnaModule/sample_CF_CTRL_channel", &hKey);
      if (strstr(str, "A")) {
        ival = 0;
        size = sizeof(ival);
        db_set_data(hDB, hKey, &ival, size, 1, TID_INT);
      } else {
        ival = 1;
        size = sizeof(ival);
        db_set_data(hDB, hKey, &ival, size, 1, TID_INT);
      }
    }
  }

  return SUCCESS;
}

/*-- EOR routine ---------------------------------------------------*/
INT sc_eor(INT run_number)
{
  INT nbin;
  INT odb_period;
  double period;
  char name[256], title[256];
  extern char runname[256];
  TH1F *hist;
  HNDLE hDB, hKey;
  INT status, size;

  // check if histo already exists, and if so remove them: (i) from the folder, and (ii) from the memory
  for (int i=0; i<N_SC_HIST; i++){
    sprintf(name, "%s", sc_ana_param.histotitles.titles[i]);
    if (!gHistoFolderStack->Last())
      hist = (TH1F *) gManaHistosFolder->FindObjectAny(name);
    else
      hist = (TH1F *) ((TFolder *) gHistoFolderStack->Last())->FindObjectAny(name);

    if (hist != NULL) {
      if (!gHistoFolderStack->Last())
        gManaHistosFolder->Remove(hist);
      else
        ((TFolder *) gHistoFolderStack->Last())->Remove(hist);

      // since object is now removed, it can be deleted
      delete hist;
    }
  }

  // get reading period
  cm_get_experiment_database(&hDB, NULL);
  status = db_find_key(hDB, 0, "/Equipment/SlowControl/Common/Period", &hKey);
  if (status != DB_SUCCESS) {
    cm_msg(MERROR, "sc_eor", "sc_eor: couldn't obtain the reading period, will assume 30 sec.");
    period = 30.0;
  } else {
    size = sizeof(INT);
    status = db_get_data(hDB, hKey, &odb_period, &size, TID_INT);
    if (status == DB_SUCCESS) {
      period = (double) odb_period / 1.0e3;
    } else {
      cm_msg(MERROR, "sc_eor", "sc_eor: couldn't obtain the reading period, will assume 30 sec.");
      period = 30.0;
    }
  }

  // create all the necessary slow control histograms
  for (int i=0; i<N_SC_HIST; i++) {
    sprintf(name, "%s", sc_ana_param.histotitles.titles[i]);
    sprintf(title, "%s Run %s_%04d", sc_ana_param.histotitles.titles[i], runname, runinfo.run_number);
    nbin  = (INT) sc_histo_vec[i].size();
    hSC_Hist[i] = new TH1F(name, title, nbin, -period*0.5, ((double)nbin+0.5)*period);
    if (hSC_Hist[i] == 0) {
      cm_msg(MERROR, "sc_eor", "sc_eor: hSC_Hist[%d]==0!! Cannot book the slow control histograms.", i);
      return SUCCESS;
    }

    // set the title
    hSC_Hist[i]->GetXaxis()->SetTitle("time (sec since SOR)");

    // fill the histos
    for (int j=0; j<sc_histo_vec[i].size(); j++) {
      hSC_Hist[i]->SetBinContent(j+1, sc_histo_vec[i][j]);
    }

    // add the slow control histograms to the SCAnaModule folder
    TFolder *scFolder = (TFolder *) gManaHistosFolder->FindObject("SCAnaModule");
    scFolder->Add(hSC_Hist[i]);
  }

  return SUCCESS;
}
/*-- event routine -------------------------------------------------*/
INT sc_ana(EVENT_HEADER *pheader, void *pevent)
{
  INT   n, ind;
  float *ptdc, field, current;
  HNDLE hDB, hKey;
  double meanModeratorHV, meanSampleHV, meanSampleT, meanSampleB;
  double errorModeratorHV, errorSampleHV, errorSampleT, errorSampleB;
  unsigned int i;

  // look for MMOD bank, ModCryo data
  n = bk_locate(pevent, "MMOD", &ptdc);
  if (n != 0 ){
    ind = sc_ana_param.mod_cf1_channel;
    sc_histo_vec[SC_MOD_CF].push_back(ptdc[ind]);
    ind = sc_ana_param.mod_xtc_channel;
    sc_histo_vec[SC_MOD_XTC].push_back(ptdc[ind]);
  }

  // look for MSAM bank, Sample Temperature data
  n = bk_locate(pevent, "MSAM", &ptdc);
  if (n != 0 ){
    ind = sc_ana_param.sample_cf_ctrl_channel;
    sc_histo_vec[SC_SAMPLE_CF].push_back(ptdc[ind]);
  }

  // look for MVAC bank, vacuum data
  n = bk_locate(pevent, "MVAC", &ptdc);
  if (n != 0 ){
    ind = sc_ana_param.sample_sc_gj_channel;
    sc_histo_vec[SC_PRESSURE_SC_GJ].push_back(ptdc[ind]);
  }

  // look for MHVT bank, transport HV data
  n = bk_locate(pevent, "MHVT", &ptdc);
  if (n != 0 ){
    ind = sc_ana_param.mod_hv_channel;
    sc_histo_vec[SC_MOD_HV].push_back(ptdc[ind]);
    ind = sc_ana_param.sample_hv_channel;
    sc_histo_vec[SC_SAMPLE_HV].push_back(ptdc[ind]);
  }

  // look for MVAC bank, vacuum data
  n = bk_locate(pevent, "M900", &ptdc);
  if (n != 0 ){
    ind = sc_ana_param.sample_zeroflux_channel;
    current = ptdc[ind]/10.*500.;
    field   = get_magnetic_field(current);
    sc_histo_vec[SC_SAMPLE_ZERO_FLUX].push_back(field);
  }

  // ----------------------------------------------------------------
  // Calculate mean values and errors of mean values for selected
  // slow control parameter.
  // Write mean and errors to /Equipment/SlowControl/Variables/MEAN.

  // calculate mean value for the moderator HV
  meanModeratorHV = 0.0;
  for (i=0; i<sc_histo_vec[SC_MOD_HV].size(); i++)
    meanModeratorHV += sc_histo_vec[SC_MOD_HV][i];
  meanModeratorHV /= sc_histo_vec[SC_MOD_HV].size();
  // calculate standard deviation for the moderator HV
  errorModeratorHV = 0.0;
  for (i=0; i<sc_histo_vec[SC_MOD_HV].size(); i++)
    errorModeratorHV += (meanModeratorHV - sc_histo_vec[SC_MOD_HV][i]) * (meanModeratorHV - sc_histo_vec[SC_MOD_HV][i]);
  errorModeratorHV = sqrt(errorModeratorHV/(float)(sc_histo_vec[SC_MOD_HV].size()-1));

  // calculate mean value for the sample HV
  meanSampleHV = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_HV].size(); i++)
    meanSampleHV += sc_histo_vec[SC_SAMPLE_HV][i];
  meanSampleHV /= sc_histo_vec[SC_SAMPLE_HV].size();
  // calculate standard deviation for the sample HV
  errorSampleHV = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_HV].size(); i++)
    errorSampleHV += (meanSampleHV - sc_histo_vec[SC_SAMPLE_HV][i]) * (meanSampleHV - sc_histo_vec[SC_SAMPLE_HV][i]);
  errorSampleHV = sqrt(errorSampleHV/(float)(sc_histo_vec[SC_SAMPLE_HV].size()-1));

  // calculate mean value for the sample temperature
  meanSampleT = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_CF].size(); i++)
    meanSampleT += sc_histo_vec[SC_SAMPLE_CF][i];
  meanSampleT /= sc_histo_vec[SC_SAMPLE_CF].size();
  // calculate standard deviation for the sample temperature
  errorSampleT = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_CF].size(); i++)
    errorSampleT += (meanSampleT - sc_histo_vec[SC_SAMPLE_CF][i]) * (meanSampleT - sc_histo_vec[SC_SAMPLE_CF][i]);
  errorSampleT = sqrt(errorSampleT/(float)(sc_histo_vec[SC_SAMPLE_CF].size()-1));

  // calculate mean value for the sample field
  meanSampleB = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_ZERO_FLUX].size(); i++)
    meanSampleB += sc_histo_vec[SC_SAMPLE_ZERO_FLUX][i];
  meanSampleB /= sc_histo_vec[SC_SAMPLE_ZERO_FLUX].size();
  // calculate standard deviation for the sample field
  errorSampleB = 0.0;
  for (i=0; i<sc_histo_vec[SC_SAMPLE_ZERO_FLUX].size(); i++)
    errorSampleB += (meanSampleB - sc_histo_vec[SC_SAMPLE_ZERO_FLUX][i]) * (meanSampleB - sc_histo_vec[SC_SAMPLE_ZERO_FLUX][i]);
  errorSampleB = sqrt(errorSampleB/(float)(sc_histo_vec[SC_SAMPLE_ZERO_FLUX].size()-1));

  mean.moderator_hv = (float) meanModeratorHV;
  mean.sample_hv    = (float) meanSampleHV;
  mean.sample_t     = (float) meanSampleT;
  mean.sample_b     = (float) meanSampleB;

  mean.var_moderator_hv = (float) errorModeratorHV;
  mean.var_sample_hv    = (float) errorSampleHV;
  mean.var_sample_t     = (float) errorSampleT;
  mean.var_sample_b     = (float) errorSampleB;

  cm_get_experiment_database(&hDB, NULL);
  db_find_key(hDB, 0, "/Equipment/SlowControl/Variables/MEAN", &hKey);
  db_set_record(hDB, hKey, &mean, sizeof(mean), 0);

  return SUCCESS;
}
/* ----------------------------------------------------------------------------*\
        EOF sc_ana_module.c
\* ----------------------------------------------------------------------------*/