tof_ana_module.cxx

Go to the documentation of this file.

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

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

  Contents:     Ana module to create time-of-flight and beam spot
                histograms for LEM experiment 
                
  December 2020: Fill only BCR3/BCR4 channels in TOF histograms, which are two Lyman alpha MCPs of the Mu-MASS/Lambshift experiment.
                 Channels BCL1-BCR2 are six BGO detectors.


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

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

/* standard includes */
#include <stdio.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 <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;
TOFANAMODULE_PARAM      tof_ana_param;
TOFANAMODULE_PARAM_STR(tof_ana_param_str);

INT masterChannel;

/*-- Module declaration --------------------------------------------*/
INT tof_init(void);
INT tof_ana(EVENT_HEADER*,void*);
INT tof_bor(INT run_number);
INT tof_eor(INT run_number);

ANA_MODULE tof_ana_module = {
  "TOFAnaModule",                /* module name           */
  "Thomas Prokscha",             /* author                */
  tof_ana,                       /* event routine         */
  tof_bor,                       /* BOR routine           */
  tof_eor,                       /* EOR routine           */
  tof_init,                      /* init routine          */
  NULL,                          /* exit routine          */
  &tof_ana_param,                /* parameter structure   */
  sizeof(tof_ana_param),         /* structure size        */
  (const char **)tof_ana_param_str, /* initial parameters    */
};

/*-- other global variables ----------------------------------------*/
#define N_TOFHIST  17
#define N_POSHIST  8
#ifndef HAVE_HBOOK
static TH1F* hTOFHist[2*N_TOFHIST]; // double number of histograms for Ext. stimulus ON/OFF
static TH1F* hProjHist[2*2];        // double number of histograms for Ext. stimulus ON/OFF
static TH2F* hPosHist[2*N_POSHIST]; // double number of histograms for Ext. stimulus ON/OFF
//the following objects are from mana.c
extern TFolder   *gManaHistosFolder;
extern TObjArray *gHistoFolderStack;
#endif

// rotation cosine and sine
static float rotCos, rotSin;

/*-- INIT routine --------------------------------------------------*/
INT tof_init(void)
{
  HNDLE  hDB;
  POSI_BANK_STR(posi_bank_str);
  extern char runname[256];
  char title[256], name[256];
  int   nbin, i;
  float xlow, xhigh;

  cm_get_experiment_database(&hDB, NULL);
  db_create_record(hDB, 0, "/Equipment/Trigger/Variables/POSI", strcomb(posi_bank_str));

  // histograms with external stimulus OFF 
  for (i=0; i<N_TOFHIST; i++) {
    sprintf(name,  "hTOF%02d", i);
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", tof_ana_param.histotitles.toftitles.titles[i], 
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.tofbin.histonbin[i];
    xlow  = tof_ana_param.histobinning.tofbin.histoxlow[i];
    xhigh = tof_ana_param.histobinning.tofbin.histoxup[i];
    hTOFHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }
  
  // histograms with external stimulus ON
  for (i=N_TOFHIST; i<2*N_TOFHIST; i++) {
    sprintf(name,  "hTOF%02d", i+N_OFFSET_ONOFF_HISTOGRAMS-N_TOFHIST);
    sprintf(title, "%s, Ext. ON, Run %s_%04d", tof_ana_param.histotitles.toftitles.titles[i-N_TOFHIST], 
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.tofbin.histonbin[i-N_TOFHIST];
    xlow  = tof_ana_param.histobinning.tofbin.histoxlow[i-N_TOFHIST];
    xhigh = tof_ana_param.histobinning.tofbin.histoxup[i-N_TOFHIST];
    hTOFHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }

  // histograms with external stimulus OFF 
  for (i=0; i<2; i++) {
    sprintf(name,  "hProj%02d", i);
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", tof_ana_param.histotitles.projtitles.titles[i],
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.projbin.histonbin[i];
    xlow  = tof_ana_param.histobinning.projbin.histoxlow[i];
    xhigh = tof_ana_param.histobinning.projbin.histoxup[i];
    hProjHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }
  
  // histograms with external stimulus ON
  for (i=2; i<4; i++) {
    sprintf(name,  "hProj%02d", i+N_OFFSET_ONOFF_HISTOGRAMS-2);
    sprintf(title, "%s, Ext. ON, Run %s_%04d", tof_ana_param.histotitles.projtitles.titles[i-2],
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.projbin.histonbin[i-2];
    xlow  = tof_ana_param.histobinning.projbin.histoxlow[i-2];
    xhigh = tof_ana_param.histobinning.projbin.histoxup[i-2];
    hProjHist[i] = H1_BOOK(name, title, nbin, xlow, xhigh);
  }

  // histograms with external stimulus OFF 
  for (i=0; i<N_POSHIST; i++) {
    sprintf(name,  "hPos%02d", i);
    sprintf(title, "%s, EXT. OFF, Run %s_%04d", tof_ana_param.histotitles.positiontitles.titles[i],
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.positionbin.histonbin[i];
    xlow  = tof_ana_param.histobinning.positionbin.histoxlow[i];
    xhigh = tof_ana_param.histobinning.positionbin.histoxup[i];
    hPosHist[i] = H2_BOOK(name, title, nbin, xlow, xhigh, nbin, xlow, xhigh);
  }

  // histograms with external stimulus ON
  for (i=N_POSHIST; i<2*N_POSHIST; i++) {
    sprintf(name,  "hPos%02d", i+N_OFFSET_ONOFF_HISTOGRAMS-N_POSHIST);
    sprintf(title, "%s, EXT. ON, Run %s_%04d", tof_ana_param.histotitles.positiontitles.titles[i-N_POSHIST],
            runname, runinfo.run_number);
    nbin  = tof_ana_param.histobinning.positionbin.histonbin[i-N_POSHIST];
    xlow  = tof_ana_param.histobinning.positionbin.histoxlow[i-N_POSHIST];
    xhigh = tof_ana_param.histobinning.positionbin.histoxup[i-N_POSHIST];
    hPosHist[i] = H2_BOOK(name, title, nbin, xlow, xhigh, nbin, xlow, xhigh);
  }

  rotCos = rotSin = 0.;
   
  return SUCCESS;
}

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

INT tof_bor(INT run_number)
{
  int i;
  char title[256];
  extern char runname[256];

  // histograms with external stimulus OFF 
  for (i=0; i<N_TOFHIST; i++) {
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", tof_ana_param.histotitles.toftitles.titles[i],
            runname, runinfo.run_number);
    hTOFHist[i]->SetTitle(title);
  }
  
  // histograms with external stimulus ON 
  for (i=N_TOFHIST; i<2*N_TOFHIST; i++) {
    sprintf(title, "%s, Ext. ON, Run %s_%04d", tof_ana_param.histotitles.toftitles.titles[i-N_TOFHIST],
            runname, runinfo.run_number);
    hTOFHist[i]->SetTitle(title);
  }
  
  // histograms with external stimulus OFF 
  for (i=0; i<2; i++) {
    sprintf(title, "%s, Ext. OFF, Run %s_%04d", tof_ana_param.histotitles.projtitles.titles[i],
            runname, runinfo.run_number);
    hProjHist[i]->SetTitle(title);
  }
  
  // histograms with external stimulus ON 
  for (i=2; i<4; i++) {
    sprintf(title, "%s, Ext. ON, Run %s_%04d", tof_ana_param.histotitles.projtitles.titles[i-2],
            runname, runinfo.run_number);
    hProjHist[i]->SetTitle(title);
  }
  
  // histograms with external stimulus OFF 
  for (i=0; i<N_POSHIST; i++) {
    if ( (i == 2) || (i == 3) )
      sprintf(title, "%s, Ext. OFF, %d - %d Run %s_%04d",
        tof_ana_param.histotitles.positiontitles.titles[i],
        (int) tof_ana_param.tofm2fmin,
        (int) tof_ana_param.tofm2fmax,
        runname, runinfo.run_number);
    else
      sprintf(title, "%s, Ext. OFF, Run %s_%04d", tof_ana_param.histotitles.positiontitles.titles[i],
              runname, runinfo.run_number);

    hPosHist[i]->SetTitle(title);
  }

  // histograms with external stimulus ON 
  for (i=N_POSHIST; i<2*N_POSHIST; i++) {
    if ( (i == 2+N_POSHIST) || (i == 3+N_POSHIST) )
      sprintf(title, "%s, Ext. ON, %d - %d Run %s_%04d",
        tof_ana_param.histotitles.positiontitles.titles[i-N_POSHIST],
        (int) tof_ana_param.tofm2fmin,
        (int) tof_ana_param.tofm2fmax,
        runname, runinfo.run_number);
    else
      sprintf(title, "%s, Ext. ON, Run %s_%04d", tof_ana_param.histotitles.positiontitles.titles[i-N_POSHIST],
              runname, runinfo.run_number);

    hPosHist[i]->SetTitle(title);
  }

  if ( trigger_settings.event_1_settings.td_master )
    masterChannel = TD;
  else
    masterChannel = MCP2F;

  rotCos = cos(tof_ana_param.rotationangle/180.*3.14);
  rotSin = sin(tof_ana_param.rotationangle/180.*3.14);
  
  return SUCCESS;
}

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

INT tof_eor(INT run_number)
{

  return SUCCESS;
}

/*-- event routine -------------------------------------------------*/
INT tof_ana(EVENT_HEADER *pheader, void *pevent)
{
  INT       n, i, j, channel, counter[N_TDC_CHANNELS];
  DWORD     *ptdc, tdc_data[N_TDC_CHANNELS][DATA_N_HITS];
  float     tof_TD_M2F, tof_TD_M2R;
  float   x, y, x1, x2, y1, y2, timesum;
  float     rotX, rotY; //for beam spot rotation
  POSI_BANK *pposi;
  BOOL    TIMEcut, SUMcut, TOFcut, foundTD, foundM2, foundBC; /* flag for time cuts on x,y position */
  LEM_DATA  *lem_data;
  INT        difftime_last, difftime_now;

  /* 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  EVENT_0_TYPE (Muons on MCP1)
  if ( ptdc[0] & (EVENT_0_TYPE<<16) ) return 1;

  TIMEcut = SUMcut = TOFcut = foundM2 = foundTD = foundBC = FALSE;
  memset(&counter, 0x00, sizeof(counter));
  memset(&tdc_data, 0x00, sizeof(tdc_data));

  for (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:
        foundM2 = TRUE;
        break;

      case BCL1: case BCL2: case BCL3: case BCL4: case BCR1: case BCR2: case BCR3: case BCR4:
        foundBC = TRUE;
        break;

      default:
        break;
    }
  }

  if ( foundM2 ) {
   // don't consider pileup counts here
   if (foundTD) {
     tof_TD_M2F = (float) (tdc_data[MCP2F][0]  - tdc_data[TD][0]);
     tof_TD_M2R = (float) (tdc_data[MCP2R][0]  - tdc_data[TD][0]);
   }
   x1         = (float) (tdc_data[MCP2X1][0] - tdc_data[MCP2F][0]);
   x2         = (float) (tdc_data[MCP2X2][0] - tdc_data[MCP2F][0]);
   y1         = (float) (tdc_data[MCP2Y1][0] - tdc_data[MCP2F][0]);
   y2         = (float) (tdc_data[MCP2Y2][0] - tdc_data[MCP2F][0]);
   timesum    = x1 + x2 + y1 + y2;
   if ( ptdc[0] & EXTOFF1_FLAG) {
     if (foundTD) {
       hTOFHist[0]->Fill(tof_TD_M2F, 1.);
       hTOFHist[1]->Fill(tof_TD_M2R, 1.);
     }
     hTOFHist[2]->Fill(x1, 1.);
     hTOFHist[3]->Fill(x2, 1.);
     hTOFHist[4]->Fill(y1, 1.);
     hTOFHist[5]->Fill(y2, 1.);
     hTOFHist[6]->Fill(timesum, 1.);
   }
   if ( ptdc[0] & EXTON1_FLAG) {
     if (foundTD) {
       hTOFHist[0+N_TOFHIST]->Fill(tof_TD_M2F, 1.);
       hTOFHist[1+N_TOFHIST]->Fill(tof_TD_M2R, 1.);
     }
     hTOFHist[2+N_TOFHIST]->Fill(x1, 1.);
     hTOFHist[3+N_TOFHIST]->Fill(x2, 1.);
     hTOFHist[4+N_TOFHIST]->Fill(y1, 1.);
     hTOFHist[5+N_TOFHIST]->Fill(y2, 1.);
     hTOFHist[6+N_TOFHIST]->Fill(timesum, 1.);
   }
   
   x    = (x1 - x2) * tof_ana_param.scalefactor + tof_ana_param.xoffset_mm_; // scale to mm and add offset in mm
   y    = (y1 - y2) * tof_ana_param.scalefactor + tof_ana_param.yoffset_mm_; // scale to mm and add offset in mm
   rotX = rotCos*x - rotSin*y;
   rotY = rotSin*x + rotCos*y;

   bk_create(pevent, "POSI", TID_STRUCT, (void**)&pposi);
   pposi->x = rotX;
   pposi->y = rotY;
   bk_close(pevent, pposi+1);
   
   if ((tof_ana_param.timemin <= x1 && x1 <= tof_ana_param.timemax) &&
      (tof_ana_param.timemin <= x2 && x2 <= tof_ana_param.timemax) &&
      (tof_ana_param.timemin <= y1 && y1 <= tof_ana_param.timemax) &&
      (tof_ana_param.timemin <= y2 && y2 <= tof_ana_param.timemax) )
     TIMEcut = TRUE;
   else
     TIMEcut = FALSE;

   if ( tof_ana_param.timesummin <= timesum && timesum <= tof_ana_param.timesummax)
     SUMcut = TRUE;
   else
     SUMcut = FALSE;

   if ( tof_ana_param.tofm2fmin <= tof_TD_M2F && tof_TD_M2F <= tof_ana_param.tofm2fmax)
     TOFcut = TRUE;
   else
     TOFcut = FALSE;

   if ( TIMEcut && SUMcut ) {
     if ( ptdc[0] & EXTOFF1_FLAG) {
       hPosHist[0]->Fill(pposi->x, pposi->y, 1.); // "fine" binning
       hPosHist[1]->Fill(pposi->x, pposi->y, 1.); // "course" binning
       hPosHist[4]->Fill(x1+x2, x1-x2, 1.);          // control histogram
       hPosHist[5]->Fill(y1+y2, y1-y2, 1.);          // control histogram
       hPosHist[6]->Fill(x1, x2, 1.);        // control histogram
       hPosHist[7]->Fill(y1, y2, 1.);        // control histogram
     }
     if ( ptdc[0] & EXTON1_FLAG) {
       hPosHist[0+N_POSHIST]->Fill(pposi->x, pposi->y, 1.); // "fine" binning
       hPosHist[1+N_POSHIST]->Fill(pposi->x, pposi->y, 1.); // "course" binning
       hPosHist[4+N_POSHIST]->Fill(x1+x2, x1-x2, 1.);        // control histogram
       hPosHist[5+N_POSHIST]->Fill(y1+y2, y1-y2, 1.);        // control histogram
       hPosHist[6+N_POSHIST]->Fill(x1, x2, 1.);              // control histogram
       hPosHist[7+N_POSHIST]->Fill(y1, y2, 1.);              // control histogram
     }    
   }
   if ( TIMEcut && SUMcut && TOFcut ) {
     if ( ptdc[0] & EXTOFF1_FLAG ) {
       hPosHist[2]->Fill(pposi->x, pposi->y, 1.); // "course" binning
       hPosHist[3]->Fill(pposi->x, pposi->y, 1.);
       hProjHist[0]->Fill(pposi->x, 1.);
       hProjHist[1]->Fill(pposi->y, 1.);
     }
     if ( ptdc[0] & EXTON1_FLAG ) {
       hPosHist[2+N_POSHIST]->Fill(pposi->x, pposi->y, 1.); // "course" binning
       hPosHist[3+N_POSHIST]->Fill(pposi->x, pposi->y, 1.);
       hProjHist[0+2]->Fill(pposi->x, 1.);
       hProjHist[1+2]->Fill(pposi->y, 1.);
     }
    }
  }

  if ( foundBC ) {
  // find the next most distant BC hit with respect to TD and fill it in the TOF histogram
/*
    difftime_now  = difftime_last = 0;
    if ( foundTD ){
     for ( i=BCL1; i<=BCR4; i++){
      if ( counter[i] == 0 ) continue;
      for (j = 0; j<counter[i]; j++){
          difftime_now = tdc_data[TD][0] - tdc_data[i][j];
          if ( (difftime_now > 0) &&
             (difftime_now < trigger_settings.tof_bc_window) &&
             (difftime_now > difftime_last))
             difftime_last = difftime_now;
      }
     }
     hTOFHist[7]->Fill( (float) difftime_last, 1.);
*/    
     // this is for forward timing, when "BC" signal occurs after TD
     difftime_now  = difftime_last = 0.;
     if ( foundTD ) {
       for ( i=BCR3; i<=BCR4; i++ ) {
        if ( counter[i] == 0 ) continue;
        if ( counter[i] > 0) {
          for (j = 0; j<counter[i]; j++) {
            difftime_last = tdc_data[i][j] - tdc_data[TD][0];

            if ( ptdc[0] & EXTOFF1_FLAG ) {
              hTOFHist[7]->Fill( (float) difftime_last, 1.);
              switch (i) {
                case BCR3:
                  hTOFHist[9]->Fill( (float) difftime_last, 1.);
                  break;

                case BCR4:
                  hTOFHist[10]->Fill( (float) difftime_last, 1.);
                  break;

                default:
                  break;
              }
            }

            hTOFHist[7+N_TOFHIST]->Fill( (float) difftime_last, 1.);
          }
        }
      }
    }
    
    if ( foundM2 ) {
      for ( i=BCR3; i<=BCR4; i++) {
        if ( counter[i] == 0 ) continue;
        if (counter[i] > 0) {
          for (j = 0; j<counter[i]; j++) {
            //difftime_last = tdc_data[MCP2F][0] - tdc_data[i][j];
            difftime_last = tdc_data[i][j]-tdc_data[MCP2F][0]; //changed to check for MCP delay

            if ( ptdc[0] & EXTOFF1_FLAG ) {
              hTOFHist[8]->Fill( (float) difftime_last, 1.);
              switch (i) {
                case BCR3:
                  hTOFHist[13]->Fill( (float) difftime_last, 1.);
                  break;
         
                case BCR4:
                  hTOFHist[14]->Fill( (float) difftime_last, 1.);
                  break;
         
                default:
                  break;
              }
            }
            if ( ptdc[0] & EXTON1_FLAG )
              hTOFHist[8+N_TOFHIST]->Fill( (float) difftime_last, 1.);
          }
        }
      }
    }    
  }
  
  //added by Gianluca
  if (foundM2) {
    float difftime_TDBack = 0;
    
    for(int l=0; l<counter[MCP2F]; l++) {
      difftime_TDBack = (float)(tdc_data[MCP2F][l] - tdc_data[TD][0]);
      
      if (difftime_TDBack < 1800 && difftime_TDBack > 800) {
        hTOFHist[11]->Fill(difftime_TDBack, 1.);
              
        for (i=BCR3; i<=BCR4; i++) {
          float difftime_TDLya = 0;
          float difftime_MCP2Lya = 0;
        
          for (j=0; j<counter[i]; j++) {
            difftime_TDLya = (float)(tdc_data[i][j]-tdc_data[TD][0])-66560;
            difftime_MCP2Lya = 66560-(float)(tdc_data[i][j]-tdc_data[MCP2F][l]);
          
            if (difftime_TDLya > 700 && difftime_TDLya < 900) {
              if (difftime_MCP2Lya > 400 && difftime_MCP2Lya < 800) {

                if ( ptdc[0] & EXTOFF1_FLAG ) {
                  hTOFHist[12]->Fill(difftime_TDBack, 1.);
                }
              }
            }
          }
        }
      }
    }
  }

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