/*
* $Id$
*
* Author: Sarah-Jane Lonsdale
* Date: 02-Dec-2015
* Version: 5.1 (8-Aug-2017)
* Update: Claudia Lederer-Woods
* Date: 31-May-2018
* Version: 6.0 (31-May-2018)
* Update: Nikolay Sosnin
* Date: 24-Nov-2022
* Version: 6.1 (24-Nov-2022)
*/
#include <math.h>
#include <string.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <vector>
#include <TFile.h>
#include <TH1F.h>
#include <TString.h>
using namespace std;
#include "MWDdetector.h"
void FormatHist(TH1D* h, TString name, TString title, TString xtitle, TString ytitle, int color, int width, int marker_color, int marker_style);
bool MWDDetector::parseConfigLine(char* line, const char* settings_file)
{
cout << "Config: " << line << endl;
if (!Detector::parseConfigLine(line)) return false;
// threshold
char* pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid Threshold" << endl;
return false;
}
threshold = atof(pch);
// polarity
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid negative polarity" << endl;
return false;
}
polarity = atoi(pch);
if (polarity>=0)
polarity = 1;
else
polarity = -1;
// pole zero correction
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid pole zero correction" << endl;
return false;
}
pz = atof(pch);
// deconvolution window
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid deconvolution window" << endl;
return false;
}
m = atoi(pch);
// average window
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid average window" << endl;
return false;
}
l = atoi(pch);
// presample
pch = strtok(NULL," "); //NS: 27.07.2023
if (!pch) {
cerr << "UserInput: Invalid presample" << endl;
return false;
}
presample = atoi(pch);
//averaging windows 2
pch = strtok(NULL," "); //NS: 27.07.2023
if (!pch) {
cerr << "UserInput: Invalid presample" << endl;
return false;
}
window = atoi(pch);
// presample/averaging window 2
/*if (!parse2real(&presample, &window)) {
cerr << "UserInput: Invalid averager 1/2" << endl;
return false;
}*/
// gain/offset
if (!parse2real(&gain, &offset)) {
cerr << "UserInput: Invalid gain/offset" << endl;
return false;
}
// gamma flash search start
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid gamma_threshold" << endl;
return false;
}
g_threshold = atoi(pch);
// minimum amp_threshold
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid amp_threshold" << endl;
return false;
}
amp_threshold = atoi(pch);
// gamma flash primary (additional) offset
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: Invalid fixed Dead time(ns)" << endl;
return false;
}
gamma_time_primary = atoi(pch);//SL 08/07/17
// time different for baseline determination
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: time_diff_baselne" << endl;
return false;
}
time_diff_baseline = atoi(pch);//SL 08/07/17
// time after g flash which has undershoot
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: tailtime" << endl;
return false;
}
tailtime = atoi(pch);//SL 08/07/17
// Time constant
pch = strtok(NULL," ");
if (!pch) {
cerr << "UserInput: time_const" << endl;
return false;
}
time_const = atoi(pch); // 28/08/23
return true;
} // parseConfigLine
int MWDDetector::analysis(
ntof::lib::ReaderStructEVEH& eveh, // EVEH event information
ntof::lib::ReaderStructMODH& modh, // MODH header information
ntof::lib::ReaderStructACQC& acqc, // ACQC pulse record
PulseVector* pulsevec, // vector of pulses
int movie_number,
bool html) // I: movie number
{
Detector::analysis(eveh, modh, acqc, pulsevec, movie_number, html);
double rate = modh.getSampleRate();
int NofPeaks = 0;
TString name = modh.getDetectorType();
//double g_threshold = 2000; //move to .h file //CHANGED
double tdiffsig = gamma_time_primary;
double tdiffbase = time_diff_baseline; // consider moving to h file
int aver = presample; //
//double tailtime = 100000; // time up to which there is an undershoot after gflash
double* x = new double[length_of_movie];
double* xsmooth = new double[length_of_movie]; //13.3.18
double* xsmooth_ma = new double[length_of_movie]; //NS 22.05.2023
double* y = new double[length_of_movie];
double* z = new double[length_of_movie];
double* zdiff = new double[length_of_movie]; //13.3.18
double* mwd_m = new double[length_of_movie];
double* ma_l = new double[length_of_movie]; // moving average array
double* mwd_deriv = new double[length_of_movie]; // mwd derivative
double timeScale = 1000.0 / rate; // in ns/Sample
for(int i = 0; i < length_of_movie; i++){
x[i] = polarity * (acqc[i] * gain + offset);
y[i] = z[i] = zdiff[i] = xsmooth[i] = 0.0; //13.3.18
}
//CHANGED
int startofevent = aver / 2;
//Averaging preamplifier output
for(int i = startofevent; i < length_of_movie; i++){
xsmooth[i] = 0.;
for(int j = -1 * (aver - 1) / 2; j <= (aver - 1) / 2; j++){
xsmooth[i] += x[i + j];
}
xsmooth[i] /= aver;
}
//Extra moving average test
//const double window = 39.;
int startofevent2 = window / 2; //NS 22.05.2023
for(int i = startofevent2; i < length_of_movie; i++){
xsmooth_ma[i] = 0.;
for(int j = -1 * (window - 1) / 2; j <= (window - 1) / 2; j++){
xsmooth_ma[i] += xsmooth[i + j];
}
xsmooth_ma[i] /= window;
}
// Locating maximum and minimum of derivative
if(verbose){cout << "Begin amplitude extraction." << endl;}
// Parameters for semi-gauss discriminator
//const double time_const = 0; // SL 12/08/16 Emmanuel
const double pole_zero = 5.e5;
double a0, a1, b1;
b1 = exp(-1. / int(time_const));
a0 = (1. + b1) / 2.0;
a1 = -1. * (1. + b1) / 2.0;
// Single pole high pass with pz correction
for(int i = 1; i < length_of_movie; i++){
//y[i] = b1 * y[i - 1] + a0 * xsmooth[i] + a1 * xsmooth[i - 1] + xsmooth[i - 1] / pole_zero;
y[i] = b1 * y[i - 1] + a0 * xsmooth_ma[i] + a1 * xsmooth_ma[i - 1] + xsmooth_ma[i - 1] / pz;
}
// Single pole low pass filter
for(int i = 1; i < length_of_movie; i++){
z[i] = b1 * z[i - 1] + a0 * y[i];
}
// devirative of filters
for(int i = 1; i < length_of_movie; i++){
zdiff[i] = -1. * z[i - 1] + z[i];
}
//bool beamType = false;
//if(eveh.getBeamType() != 1){beamType = true;}
double twait = 16000.; // introduce different time window for baseline depending on ded or par CLW 31/05/18
/*if(eveh.getBeamType() == 2){twait = 13000.;} // dedicated
if(eveh.getBeamType() == 3){twait = 19000.;} // parasitic
else{twait = 16000.;}*/
////////////////////////////////////////
// MWD and filtering //
////////////////////////////////////////
// moving window deconvolution
for(int i = startofevent + m; i < length_of_movie; i++){ // SL 10/08/16 // 13/03/18 MWD on smoothed
//double d_m = xsmooth[i] - xsmooth[i - m];
double d_m = xsmooth_ma[i] - xsmooth_ma[i - m];
double ma_m = 0.;
for(int j = (i - m); j < (i - 1); j++) {
//ma_m += xsmooth[j];
ma_m += xsmooth_ma[j];
}
mwd_m[i] = d_m + ma_m / pz;
}
// moving average
//for(int i = l + m + gamma_flash; i < length_of_movie; i++){ // SL 10/08/16
for(int i = l + m; i < length_of_movie; i++){ // NS 31.07.2023
ma_l[i] = 0.;
... 285 more lines ...
|