decay_ana_module.cxx

Go to the documentation of this file.

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

  Name:         decay_ana_module.c
  Created by:   Thomas Prokscha
  Date:         01-Apr-2007

  Contents:     Analyzer module for LE-muSR decayhistograms

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

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

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

/* midas includes */
#include "midas.h"
#include "rmidas.h"
#include "v1190.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

/*-- ODB Parameters ------------------------------------------------*/
//extern RUNINFO        runinfo;
//extern EXP_PARAM      exp_param;
extern RUNINFO          runinfo;
extern TRIGGER_SETTINGS trigger_settings;
extern T0SHIFT_PARAM    t0shift_param;
DECAYANAMODULE_PARAM    decay_ana_param;
DECAYANAMODULE_PARAM_STR(decay_ana_param_str);

/*-- Module declaration --------------------------------------------*/
INT decay_init(void);
INT decay_ana(EVENT_HEADER*,void*);
INT decay_bor(INT run_number);
INT decay_eor(INT run_number);

ANA_MODULE decay_ana_module = {
  "DecayAnaModule",               /* module name           */
  "Thomas Prokscha",              /* author                */
  decay_ana,                      /* event routine         */
  decay_bor,                      /* BOR routine           */
  decay_eor,                      /* EOR routine           */
  decay_init,                     /* init routine          */
  NULL,                           /* exit routine          */
  &decay_ana_param,               /* parameter structure   */
  sizeof(decay_ana_param),        /* structure size        */
  (const char **)decay_ana_param_str, /* initial parameters    */
};

/*-- other global variables ----------------------------------------*/
#ifndef HAVE_HBOOK
static TH1F* hDecHist[4*N_DECAY_HISTS];
//the following objects are from mana.c
extern TFolder   *gManaHistosFolder;
extern TObjArray *gHistoFolderStack;
#endif

#define POS_FLAG_LD  (1<<0)    
#define POS_FLAG_TD  (1<<1)
#define POS_FLAG_RD  (1<<2)
#define POS_FLAG_BD  (1<<3)
#define POS_FLAG_LU  (1<<4)    
#define POS_FLAG_TU  (1<<5)
#define POS_FLAG_RU  (1<<6)
#define POS_FLAG_BU  (1<<7)

static DWORD dataWindowMuonDecay;  //needed for post-pileup rejection

/*-- INIT routine --------------------------------------------------*/
INT decay_init(void)
{
  extern char runname[256];
  char title[256], name[256];
  int   nbin;
  float xlow, xhigh;

  // decay histograms with external stimulus OFF, npp and ppc 
  for (int i=0; i<2*N_DECAY_HISTS; i++){
    if ( i > N_DECAY_HISTS-1)
     sprintf(name,  "hDecay%03d", i+1+N_OFFSET_PPC_HISTOGRAMS-N_DECAY_HISTS); //histograms with post-pileup correction: offset 20 in name
    else
     sprintf(name,  "hDecay%03d", i+1);
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", decay_ana_param.histotitles.decaytitles.titles[i],
            runname, runinfo.run_number);
    nbin  = decay_ana_param.histobinning.decaybin.histonbin[i];
    xlow  = decay_ana_param.histobinning.decaybin.histoxlow[i];
    xhigh = decay_ana_param.histobinning.decaybin.histoxup[i];
    hDecHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }
  
  // decay histograms with external stimulus ON, npp and ppc 
  for (int i=2*N_DECAY_HISTS; i<4*N_DECAY_HISTS; i++){
    if ( i > 3*N_DECAY_HISTS-1)
     sprintf(name,  "hDecay%03d", i+1+N_OFFSET_PPC_HISTOGRAMS+N_OFFSET_ONOFF_HISTOGRAMS-3*N_DECAY_HISTS); //histograms with post-pileup correction: offset 20 in name
    else
     sprintf(name,  "hDecay%03d", i+1+N_OFFSET_ONOFF_HISTOGRAMS-2*N_DECAY_HISTS);
    sprintf(title, "%s, Ext. ON, Run %s_%04d", decay_ana_param.histotitles.decaytitles.titles[i],
            runname, runinfo.run_number);
    nbin  = decay_ana_param.histobinning.decaybin.histonbin[i-2*N_DECAY_HISTS];
    xlow  = decay_ana_param.histobinning.decaybin.histoxlow[i-2*N_DECAY_HISTS];
    xhigh = decay_ana_param.histobinning.decaybin.histoxup[i-2*N_DECAY_HISTS];
    hDecHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }
  
  dataWindowMuonDecay = 0;
  
  return SUCCESS;
}

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

INT decay_bor(INT run_number)
{
  char title[256];
  extern char runname[256];

  // decay histograms with external stimulus OFF
  for (int i=0; i<2*N_DECAY_HISTS; i++){
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", decay_ana_param.histotitles.decaytitles.titles[i], 
            runname, runinfo.run_number);
    hDecHist[i]->SetTitle(title);
  }
  
  // decay histograms with external stimulus ON
  for (int i=2*N_DECAY_HISTS; i<4*N_DECAY_HISTS; i++){
    sprintf(title, "%s, Ext. ON, Run %s_%04d", decay_ana_param.histotitles.decaytitles.titles[i-2*N_DECAY_HISTS], 
            runname, runinfo.run_number);
    hDecHist[i]->SetTitle(title);
  }

  // get data window for post-pileup rejection from ODB
  if ( trigger_settings.events.event_1_active)
   dataWindowMuonDecay = trigger_settings.event_1_settings.data_window;
  else
   dataWindowMuonDecay = trigger_settings.event_2_settings.data_window; 
  
  return SUCCESS;
}

/*-- EOR routine ---------------------------------------------------*/

INT decay_eor(INT run_number)
{
  return SUCCESS;
}

/*-- event routine -------------------------------------------------*/
INT decay_ana(EVENT_HEADER *pheader, void *pevent)
{
  INT       n, channel, counter[N_TDC_CHANNELS];
  INT     masterCH;
  DWORD     *ptdc, tdc_data[N_TDC_CHANNELS][DATA_N_HITS];
  BOOL    foundTD;      /* flag indicating that we got a start signal */
  BOOL      foundMCP2;
  DWORD     posHit_flag; /* flag indicating which e+ detector were hit*/
  LEM_DATA  *lem_data;
  INT        difftime_last;

  /* look for TDC0 bank */
  n = bk_locate(pevent, "TDC0", &ptdc);
  if (n == 0 )
    return 1;

  // return if we don't have the NEW_EVENT_FLAG present
  if ( !(ptdc[0] & NEW_EVENT_MASK)) return 1;

  // return if not EVENT_0_TYPE (Muons on MCP1)
  if ( ptdc[0] & (EVENT_0_TYPE<<16) ) return 1;

  masterCH    = TD;
  foundTD     = FALSE;
  foundMCP2   = FALSE;
  posHit_flag = 0;
  memset(&counter, 0x00, sizeof(counter));
  memset(&tdc_data, 0x00, sizeof(tdc_data));

  for (int i=1; i<n; i++){ // i=0 is a control word, including the event length
    lem_data = (LEM_DATA *) &ptdc[i];
    channel    = lem_data->channel;
    if ( channel > N_TDC_CHANNELS-1 ) continue; //not understood how this can happen, but it can happen
                                  //after a HV break at the moderator (all channel get "many" hits, like
                                  //a "flash" in the detectors
    tdc_data[channel][counter[channel]] = lem_data->data;
    counter[channel]++;

    switch ( channel ){
      case TD:
        foundTD = TRUE;
        break;

      case MCP2F:
        foundMCP2 = TRUE;
        break;

      case POSLDI: case POSLDO:
        posHit_flag = posHit_flag | POS_FLAG_LD;
        break;
      case POSLUI: case POSLUO:
        posHit_flag = posHit_flag | POS_FLAG_LU;
        break;

      case POSTDI: case POSTDO:
        posHit_flag = posHit_flag | POS_FLAG_TD;
        break;
      case POSTUI: case POSTUO:
        posHit_flag = posHit_flag | POS_FLAG_TU;
        break;

      case POSRDI: case POSRDO:
        posHit_flag = posHit_flag | POS_FLAG_RD;
        break;
      case POSRUI: case POSRUO:
        posHit_flag = posHit_flag | POS_FLAG_RU;
        break;

      case POSBDI: case POSBDO:
        posHit_flag = posHit_flag | POS_FLAG_BD;
        break;
      case POSBUI: case POSBUO:
        posHit_flag = posHit_flag | POS_FLAG_BU;
        break;

      default:
        break;
    }
  }

  if ( !foundTD && !foundMCP2 ) return SUCCESS;
  if ( !foundTD ) masterCH = MCP2F;

  // fill decay histograms without further cut; the "inner" detector
  // determines the timing, "downstream/forward" detectors first
  if ( tdc_data[POSLDI][0] > 0 && tdc_data[POSLDO][0] > 0 ) {
    difftime_last = tdc_data[POSLDI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[0]->Fill( (float) difftime_last + t0shift_param.delta_t0[0], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[16]->Fill( (float) difftime_last + t0shift_param.delta_t0[0], 1.);
  }
  if ( tdc_data[POSTDI][0] > 0 && tdc_data[POSTDO][0] > 0 ) {
    difftime_last = tdc_data[POSTDI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[1]->Fill( (float) difftime_last + t0shift_param.delta_t0[1], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[17]->Fill( (float) difftime_last + t0shift_param.delta_t0[1], 1.);
  }
  if ( tdc_data[POSRDI][0] > 0 && tdc_data[POSRDO][0] > 0 ) {
    difftime_last = tdc_data[POSRDI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[2]->Fill( (float) difftime_last + t0shift_param.delta_t0[2], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[18]->Fill( (float) difftime_last + t0shift_param.delta_t0[2], 1.);
  }
  if ( tdc_data[POSBDI][0] > 0 && tdc_data[POSBDO][0] > 0 ) {
    difftime_last = tdc_data[POSBDI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[3]->Fill( (float) difftime_last + t0shift_param.delta_t0[3], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[19]->Fill( (float) difftime_last + t0shift_param.delta_t0[3], 1.);
  }
  // fill decay histograms without further cut; the "inner" detector
  // determines the timing, "upstream/backward" detectors
  if ( tdc_data[POSLUI][0] > 0 && tdc_data[POSLUO][0] > 0 ) {
    difftime_last = tdc_data[POSLUI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[4]->Fill( (float) difftime_last + t0shift_param.delta_t0[4], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[20]->Fill( (float) difftime_last + t0shift_param.delta_t0[4], 1.);
  }
  if ( tdc_data[POSTUI][0] > 0 && tdc_data[POSTUO][0] > 0 ) {
    difftime_last = tdc_data[POSTUI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[5]->Fill( (float) difftime_last + t0shift_param.delta_t0[5], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[21]->Fill( (float) difftime_last + t0shift_param.delta_t0[5], 1.);
  }
  if ( tdc_data[POSRUI][0] > 0 && tdc_data[POSRUO][0] > 0 ) {
    difftime_last = tdc_data[POSRUI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[6]->Fill( (float) difftime_last + t0shift_param.delta_t0[6], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[22]->Fill( (float) difftime_last + t0shift_param.delta_t0[6], 1.);
  }
  if ( tdc_data[POSBUI][0] > 0 && tdc_data[POSBUO][0] > 0 ) {
    difftime_last = tdc_data[POSBUI][0] - tdc_data[masterCH][0];
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[7]->Fill( (float) difftime_last + t0shift_param.delta_t0[7], 1.);
    if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[23]->Fill( (float) difftime_last + t0shift_param.delta_t0[7], 1.);
  }

  // check for post-pileup in masterCH: if post-pileup, return
  if ( tdc_data[masterCH][1] != 0 && 
      (tdc_data[masterCH][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay) return SUCCESS;
  
  // check, if we had only 1 e+ detector fired, and that for this detector there was
  // no 2nd e+ during data window
  
  // "Left" detectors 
  if ( posHit_flag == POS_FLAG_LD ){
   if ( tdc_data[POSLDI][1] != 0 ){
    if ( (tdc_data[POSLDI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&
          abs(tdc_data[POSLDI][1]-tdc_data[POSLDO][1]) < 50 ) return SUCCESS; //2nd e+
   }
   if ( ptdc[0] & EXTOFF1_FLAG)
     hDecHist[8]->Fill( (float)(tdc_data[POSLDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[0], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[24]->Fill( (float)(tdc_data[POSLDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[0], 1.);
   return SUCCESS;
  }
  if ( posHit_flag == POS_FLAG_LU ){
   if ( tdc_data[POSLUI][1] != 0 ){
    if ( (tdc_data[POSLUI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&
          abs(tdc_data[POSLUI][1]-tdc_data[POSLUO][1]) < 50 ) return SUCCESS; //2nd e+
   }
   if ( ptdc[0] & EXTOFF1_FLAG)
     hDecHist[12]->Fill( (float)(tdc_data[POSLUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[4], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[28]->Fill( (float)(tdc_data[POSLUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[4], 1.);
   return SUCCESS;
  }
 
  // "Top" detectors
  if ( posHit_flag == POS_FLAG_TD ){   
   if ( tdc_data[POSTDI][1] != 0 ){
    if ( (tdc_data[POSTDI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&  
          abs(tdc_data[POSTDI][1]-tdc_data[POSTDO][1]) < 50) return SUCCESS; //2nd e+
   }
    if ( ptdc[0] & EXTOFF1_FLAG)
      hDecHist[9]->Fill( (float)(tdc_data[POSTDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[1], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
      hDecHist[25]->Fill( (float)(tdc_data[POSTDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[1], 1.);
   return SUCCESS;
  }
  if ( posHit_flag == POS_FLAG_TU ){   
   if ( tdc_data[POSTUI][1] != 0 ){
    if ( (tdc_data[POSTUI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&  
          abs(tdc_data[POSTUI][1]-tdc_data[POSTUO][1]) < 50) return SUCCESS; //2nd e+
   }
   if ( ptdc[0] & EXTOFF1_FLAG)
     hDecHist[13]->Fill( (float)(tdc_data[POSTUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[5], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[29]->Fill( (float)(tdc_data[POSTUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[5], 1.);
   return SUCCESS;
  }
  
  // "Right" detectors
  if ( posHit_flag == POS_FLAG_RD ){   
   if ( tdc_data[POSRDI][1] != 0 ){
    if ( (tdc_data[POSRDI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&  
          abs(tdc_data[POSRDI][1]-tdc_data[POSRDO][1]) < 50) return SUCCESS; //2nd e+
   }
   if ( ptdc[0] & EXTOFF1_FLAG)   
     hDecHist[10]->Fill( (float)(tdc_data[POSRDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[2], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[26]->Fill( (float)(tdc_data[POSRDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[2], 1.);
   return SUCCESS;
  }
  if ( posHit_flag == POS_FLAG_RU ){   
   if ( tdc_data[POSRUI][1] != 0 ){
    if ( (tdc_data[POSRUI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&  
          abs(tdc_data[POSRUI][1]-tdc_data[POSRUO][1]) < 50) return SUCCESS; //2nd e+
   }
   if ( ptdc[0] & EXTOFF1_FLAG)   
     hDecHist[14]->Fill( (float)(tdc_data[POSRUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[6], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[30]->Fill( (float)(tdc_data[POSRUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[6], 1.);
   return SUCCESS;
  }
  
  // "Bottom" detectors
  if ( posHit_flag == POS_FLAG_BD ){
   if ( tdc_data[POSBDI][1] != 0 ){ 
    if ( (tdc_data[POSBDI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&
        abs(tdc_data[POSBDI][1]-tdc_data[POSBDO][1]) < 50) return SUCCESS;
   }
   if ( ptdc[0] & EXTOFF1_FLAG)
     hDecHist[11]->Fill( (float)(tdc_data[POSBDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[3], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[27]->Fill( (float)(tdc_data[POSBDI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[3], 1.);
   return SUCCESS;
  }
  if ( posHit_flag == POS_FLAG_BU ){
   if ( tdc_data[POSBUI][1] != 0 ){ 
    if ( (tdc_data[POSBUI][1] - tdc_data[masterCH][0]) < dataWindowMuonDecay &&
        abs(tdc_data[POSBUI][1]-tdc_data[POSBUO][1]) < 50) return SUCCESS;
   }
   if ( ptdc[0] & EXTOFF1_FLAG)   
     hDecHist[15]->Fill( (float)(tdc_data[POSBUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[7], 1.);
   if ( ptdc[0] & EXTON1_FLAG)
     hDecHist[31]->Fill( (float)(tdc_data[POSBUI][0] - tdc_data[masterCH][0]) + t0shift_param.delta_t0[7], 1.);
   return SUCCESS;
  }
  
  return SUCCESS;
}
/* ----------------------------------------------------------------------------*\
        EOF decay_ana_module.c
\* ----------------------------------------------------------------------------*/