hit_analyse_v2/Scripts_20180425/hit_analyse.c

#define hit_analyse_cxx
#include "hit_analyse.h"

int main(int argc, char **argv){

  opendatafiles(argc, argv);
  histograms(argc, argv);
  analyse(argc, argv);
  closedatafiles();
  return 0;
}

int opendatafiles(int argc, char ** argv){
  if (argc>2){
    //open bpm data file
    filename = Form("%s%s.dat",argv[1],argv[2]);
    file.open(filename, ifstream::in | ifstream::binary);
    fileframesize = getFileSize(filename) / ( 4*sizeof(BufferData) );
    if (fileframesize>0){     std::cout << "Number of frames in data file: " << fileframesize << std::endl;}
    else { std:cout << "BPM .dat dile not found." << endl;  return -1;}

    //bpm data timestamps
    timestampfilename =  Form("%s%s_timestamp.csv",argv[1],argv[2]); 
    timestampfile.open(timestampfilename, ifstream::in);
    if (!timestampfile.is_open()){
      printf("timestamp.csv did not open.\n");
      ethercat = false;
      //return -1;	//file could not be opened
    }


  //open ethercat file
  if (argc>3){
    ethercatfile = argv[3];
    tree2 = new TTree("t2", "t2");
    std::cout << " Loading Ethercat data." << std::endl;
    tree2->ReadFile(ethercatfile, "RELTIME2/D:IC1/D:MW1_POSX/D:MW1_POSY/D:ANALOG_IN1/D:ENERGY_INDEX/D:INTENSITY_INDEX/D:ION-SORT/D:TIME2/D", '\t');
    std::cout << "Ethercat data loaded." << std::endl;
    tree2->Print();
  }

  //open amplitude offset correction file
  if (argc>4){
    offsetfilename = Form("%s",argv[4]);
    offsetfile.open(offsetfilename, ifstream::in);
    if (!offsetfile.is_open()){
      printf("no offset.txt file found\n");
      //   return -1;	//file could not be opened
    }
  }
    
     
    
    dataptr = new BufferData();

    
    return 1;
}

int closedatafiles(int argc, char ** argv){

  if (file.is_open()) file.close();
  if (timestampfile.is_open()) timestampfile.close();
  if (offsetfile.is_open())  offsetfile.close();

  rootFile->Write();
  rootFile->Close();
}


int analyse(int argc, char **argv)
{
  int bkg_frames = 1000;
  set_background_v1(bkg_frames);

  for (int frame = 0; frame< fileframesize - bkg_frames; frame++ ){
    if (frame%10000==0) std::cout << "Frame: " <<  frame << " (" <<double(frame)/double(fileframesize)*100.0 << "%)" << std::endl;

    //must read all boards to keep read correct position in data file
    for (int boardnumber = 0; boardnumber<4;boardnumber++){
      board_b[boardnumber] = readboard(frame, boardnumber); //read in the frame
      BPMbeamrecon[boardnumber] = beamreconstruction(board_b[boardnumber], 50.); // do the linear regression fit of the beam;
      
    }

    
  }
  return 1;
}

bpm_frame_v1 readboard(int frame, int boardnumber){

  bpm_frame_v1 board;
  board.integratedsignalamp = 0.;
  
  file.seekg(boardnumber*sizeof(BufferData)+4*frame*sizeof(BufferData));
  file.read ((char*)dataptr ,sizeof(BufferData));
  if (dataptr->sync_frame.device_nr==boardnumber){
    for (int j = 1; j<128;j++){
      //subtract the background from the data
      board.channel_amp[j] = dataptr->sensor_data[j] -  board_b_bkg[boardnumber].channel_amp[j];
      //     std::cout <<  j << " " << board.channel_amp[j] << "    "  << dataptr->sensor_data[j] << std::endl;

      //sum the signal across channels
      board.integratedsignalamp +=  board.channel_amp[j]; 

      //find the peak channel
      if (board.channel_amp[j]> board.maxchannel_amp) {
	board.maxchannel = j;
	board.maxchannel_amp = board.channel_amp[j];
	//	cout << maxchannel_b0 << " " <<maxchannelamp_b0 << endl;
      }
    }
  } 
  else std::cerr << "Error reading board data." << std::endl;
  
  

  
  return board;
}

void histograms(){
   
  //open output root file
  rootfilename = Form("%s/root/%s.root",argv[1],argv[2]); 
  rootFile = new TFile(rootfilename,"RECREATE");
  if ( rootFile->IsOpen() ) printf("ROOT file opened successfully\n");
}
 else return -1;
 TTree *rootTree = new TTree("t","HIT Data Root Tree");
 
}

void set_background_v1(int max_frames){

  for (int j = 0; j<128; j++){
    for (int k = 0; k<4; k++){
      board_b_bkg[k].channel_amp[j] = 0.;
    }
  }  
  for (int i = 0;i<max_frames;i++){
    //must read all boards to keep read correct position in data file
    for (int boardnumber = 0; boardnumber<4; boardnumber++){

      file.seekg(boardnumber*sizeof(BufferData)+4*i*sizeof(BufferData));
      file.read ((char*)dataptr ,sizeof(BufferData));
      if (dataptr->sync_frame.device_nr==boardnumber){
	for (int j = 1; j<128;j++){
	  board_b_bkg[boardnumber].channel_amp[j] += double(dataptr->sensor_data[j]) / double(max_frames);
	  //  std::cout <<  j << " " << board.channel_amp[j] << "    "  << dataptr->sensor_data[j] << std::endl;
	}
      } 
      else std::cerr << "Error reading board data." << std::endl;
    }
  }
}




beamRecon beamreconstruction(bpm_frame_v1 frametoanalyse, double threshold = 50.){

  ///////////////// linear regression using Integration by parts of gaussian function.

  beamRecon beam;
  double SumT, SumS, SumS2, SumST, SumT2, SumY, SumYS, SumYT, sigmaABC, muABC,p,c, b, b_den, b_num, SumYYP, SumYYM, MeanY;
  TMatrixD M1(3,3);
  TMatrixD M1inv(3,3);
  TVectorD ABC(3);
  TVectorD M2(3);
  vector<double> signal_list;
  vector<double> channel_list;

  SumY = 0.;
  SumS = 0.;
  SumT = 0.;
  SumS2 = 0.;
  SumST = 0.;
  SumT2 = 0.;
  SumYS = 0.;
  SumYT = 0.;
  b_den = 0.;
  b_num = 0.;
  b = 0.;
  p = 0.;
  c = 0.;
  SumYYM = 0.;
  SumYYP = 0.;
  MeanY = 0.;

  
  //  const int array_length = sizeof(frametoanalyse.channel_amp)/sizeof(double);
  const int array_length = 128;
  for (int i = 0; i< array_length; i++){
    if (frametoanalyse.channel_amp[i]>=threshold) {
	  signal_list.push_back(frametoanalyse.channel_amp[i]);
	  channel_list.push_back(i);//correct for actual detector position 
    }
  }
  const int vector_length = channel_list.size();
  if (vector_length<=3) return beam;
  
  double S[vector_length];
  double T[vector_length];
  
  

  
  for(int k=0; k<vector_length;k++){
    if (k==0){
      S[k]=0.;		T[k]=0.;
    }
    else{
      S[k] = S[k-1]+0.5*( signal_list[k] + signal_list[k-1] ) *  ( channel_list[k] - channel_list[k-1] );
      T[k] = T[k-1]+0.5*( channel_list[k] * signal_list[k] + channel_list[k-1] * signal_list[k-1] ) * ( channel_list[k] - channel_list[k-1] );
    }
    //  cout << S[k] << " " << T[k] << endl;
    SumS += S[k];	      SumT += T[k];
    SumY += signal_list[k];
    SumS2 += S[k]*S[k]; 	      SumST += S[k]*T[k];	      SumT2 += T[k]*T[k];
    SumYS += signal_list[k]*S[k];
    SumYT += signal_list[k]*T[k];
    MeanY+=signal_list[k];
  }
  MeanY/=vector_length;

	      
  M1(0,0) = SumT2;	      M1(0,1) = SumST;	      M1(0,2) = SumT;	      M1(1,0) = SumST;	      M1(1,1) = SumS2;
  M1(1,2) = SumS;	      M1(2,0) = SumT;	      M1(2,1) = SumS;
  M1(2,2) = vector_length;
  
  M2(0) = SumYT;	      M2(1) = SumYS;	      M2(2) = SumY;
  M1inv = M1.Invert();	      ABC = M1inv * M2;

	     

	      
  //calculate b,p,c ---> y = b*exp(-p*(x-c)*(x-c))
  p = -ABC(0)/2.;	     c = -ABC(1)/ABC(0);

  for(int k=0; k<vector_length;k++){
    b_num += exp(-p*(channel_list[k]-c)*(channel_list[k]-c)) * signal_list[k];
    b_den +=  exp(-2*p*(channel_list[k]-c)*(channel_list[k]-c));
		
  }
  b = b_num/b_den;

  beam.Position = -ABC(1)/ ABC(0);
  beam.Focus = 2.3548/sqrt(2*p);
  beam.Peak = b;
  beam.Rsqr = SumYYP/SumYYM;
  beam.Skew = gsl_stats_wskew_m_sd(&signal_list[0],1,&channel_list[0],1,vector_length,beam.Position,beam.Focus/2.3548); //skewness (symmetry)
  beam.Kurtosis = gsl_stats_wkurtosis_m_sd(&signal_list[0],1,&channel_list[0],1,vector_length,beam.Position,beam.Focus/2.3548); //excess kurtosis (well behaved tails)
  
  return beam;
  
}