leverington
/
hit_analyse_v2
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
data analysis scripts
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
20 Commits
1 Branch
0 Tags
15 MiB
 
 
 
 
Tree: 3c2c696fa6
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '3c2c696fa6'
${ noResults }
hit_analyse_v2/Scripts_20180425/convert_clean_fix.c

1306 lines
48 KiB
Raw Blame History

/// to do:
/// 1. correct for gap between arrays
#define convert_cxx
#include "convert_clean.h"
int init(int argc, char **argv);
int analyse(int argc, char **argv);
//using namespace std;
/// compile with:
//// $ make clean; make
/// run with:
//// $ ./convert <path to data> <run> <ethercatfile> <offsetfile>
///// $convert_clean /work/leverington/beamprofilemonitor/hitdata/HIT_2018-04-25/ run41 /work/leverington/beamprofilemonitor/hitdata/HIT_2018-04-25/ethercat/with_timestamp/run41.csv ./run41_offsetcor.txt
int main(int argc, char **argv){
init(argc, argv); // initialise some variables, create histograms and root tree
analyse(argc, argv);
}
int init(int argc, char **argv){
smooth_on = true; // boxcar smoothing of data
doFit = false; // gaussian fitting (slow)
doLinInt = true; // Linear Interpolation (fast)
ethercat = true; // align with ethercat data (IC and MWPC)
f_sp_air->SetParameters(303.746, -0.873697,1.01335); //protons between 50 and 300 MeV
f_sp_ps->SetParameters(361.936, -0.892255, 1.19133); //protons between 50 and 220 MeV
f_h2sp_air->SetParameters(1236.51, -0.880225, 4.17041); //helium between 50 and 300 MeV
f_h2sp_ps->SetParameters(1387.08, -0.882686,4.60833); //heelium between 50 and 300 MeV
f_c12sp_air->SetParameters(13291.2, -0.925464,45.3876); //carbon between 80 and 480 MeV
f_c12sp_ps->SetParameters(14538.9, -0.922423,49.2859); //carbon between 80 and 480 MeV
f_o16sp_air->SetParameters(24624.9, -0.925425,80.6828); //oxygen between 80 and 480 MeV
f_o16sp_ps->SetParameters(26465.6, -0.928309,88.6728); //oxygen between 80 and 480 MeV
// TVirtualFitter::SetDefaultFitter("Minuit2");
// TVirtualFitter::SetPrecision(0.1);
// std::cout << "Default Fitter:" << TVirtualFitter::GetDefaultFitter() << std::endl;
//initialise some values;
for (int i = 0;i<128;i++){
board_b0_ch_bkg[i] = 0.;
board_b1_ch_bkg[i] = 0.;
board_b2_ch_bkg[i] = 0.;
board_b3_ch_bkg[i] = 0.;
channellist_gsl_b0[i] = 0;
channellist_gsl_b1[i] = 0;
channellist_gsl_b2[i] = 0;
channellist_gsl_b3[i] = 0;
calibration_b0[i] = 1.0;
calibration_b1[i] = 1.0;
calibration_b2[i] = 1.0;
calibration_b3[i] = 1.0;
errorx[i]=17.;
errory[i]=17.;
if (i<64) {pos[i] = double(i);}
else {pos[i] = double(i) + 0.2; }
}
// if (argc > 3){
//open the file names specified in the command line
ethercatfile = argv[3];
// argv[1]= "/work/leverington/beamprofilemonitor/hitdata/HIT_17_12_2017/";
// argv[2]= "Run11"
filename = Form("%s%s.dat",argv[1],argv[2]);
offsetfilename = Form("%s",argv[4]);
timestampfilename = Form("%s%s_timestamp.csv",argv[1],argv[2]);
// modify location of rootfilename to somewhere writeable
rootfilename = Form("%s/root/%s.root",argv[1],argv[2]);
file.open(filename, ifstream::in | ifstream::binary);
timestampfile.open(timestampfilename, ifstream::in);
offsetfile.open(offsetfilename, ifstream::in);
if (!file.is_open()){
printf(".dat did not open.\n");
return -1; //file could not be opened
}
else {
printf("Data file opened.\n");
}
if (!timestampfile.is_open()){
printf("timestamp.csv did not open.\n");
ethercat = false;
//return -1; //file could not be opened
}
if (!offsetfile.is_open()){
printf("no offset.txt file found\n");
// return -1; //file could not be opened
for (int ii = 0; ii<30;ii++) {
signaloffset_b0[ii] = 0.; signaloffset_b1[ii] = 0.; signaloffset_b2[ii] = 0.; signaloffset_b3[ii] = 0.;
}
}
else {
for (int ii = 0; ii<30;ii++) offsetfile >> signaloffset_b0[ii] >> signaloffset_b1[ii] >>signaloffset_b2[ii] >>signaloffset_b3[ii] ;
}
///read in offsets produced by offsetcor.C script (needs to be run after convert_clean, then repeat convert_clean.)
///some variables to get the data
// int number_of_records_to_read = 4*10;
// BufferData* buffer = new BufferData[number_of_records_to_read];
//if(argc > 4) {
// framestart = atoi(argv[3]);
// frameend = atoi(argv[4]);
// }
board_b = -1;
fileframesize = getFileSize(filename) / ( 4*sizeof(BufferData) );
dataptr = new BufferData();
if (fileframesize>0){
std::cout << "Number of frames:" << fileframesize << std::endl;
//open ROOT file for saving histos and TTree.
rootFile = new TFile(rootfilename,"RECREATE");
if ( rootFile->IsOpen() ) printf("ROOT file opened successfully\n");
th2_signal_vs_channel_b0 = new TH2D("th2_signal_vs_channel_b0","th2_signal_vs_channel_b0",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_b1 = new TH2D("th2_signal_vs_channel_b1","th2_signal_vs_channel_b1",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_b2 = new TH2D("th2_signal_vs_channel_b2","th2_signal_vs_channel_b2",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_b3 = new TH2D("th2_signal_vs_channel_b3","th2_signal_vs_channel_b3",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_sub_b0 = new TH2D("th2_signal_vs_channel_sub_b0","th2_signal_vs_channel_sub_b0",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_sub_b1 = new TH2D("th2_signal_vs_channel_sub_b1","th2_signal_vs_channel_sub_b1",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_sub_b2 = new TH2D("th2_signal_vs_channel_sub_b2","th2_signal_vs_channel_sub_b2",128,0,128,2200,-1000,10000);
th2_signal_vs_channel_sub_b3 = new TH2D("th2_signal_vs_channel_sub_b3","th2_signal_vs_channel_sub_b3",128,0,128,2200,-1000,10000);
th1_signal_b0 = new TH1D("th1_signal_b0","th1_signal_b0",2200,-10000,50000);
th1_signal_b1 = new TH1D("th1_signal_b1","th1_signal_b1",2200,-10000,50000);
th1_signal_b2 = new TH1D("th1_signal_b2","th1_signal_b2",2200,-10000,50000);
th1_signal_b3 = new TH1D("th1_signal_b3","th1_signal_b3",2200,-10000,50000);
h_regresidual_b0 = new TH1D("h_regresidual_b0","h_regresidual_b0", 200,-.20,.20);
h_regresidual_b1 = new TH1D("h_regresidual_b1","h_regresidual_b1", 200,-.20,.20);
h_regresidual_b2 = new TH1D("h_regresidual_b2","h_regresidual_b2", 200,-.20,.20);
h_regresidual_b3 = new TH1D("h_regresidual_b3","h_regresidual_b3", 200,-.20,.20);
graph_bkg_b0 = new TGraph();
graph_bkg_b1 = new TGraph();
graph_bkg_b2 = new TGraph();
graph_bkg_b3 = new TGraph();
th2_signal_vs_channel_bkg_b0 = new TH2D("th2_signal_vs_channel_bkg_b0","th2_signal_vs_channel_bkg_b0",128,0,128,1000,-500,500);
th2_signal_vs_channel_bkg_b1 = new TH2D("th2_signal_vs_channel_bkg_b1","th2_signal_vs_channel_bkg_b1",128,0,128,1000,-500,500);
th2_signal_vs_channel_bkg_b2 = new TH2D("th2_signal_vs_channel_bkg_b2","th2_signal_vs_channel_bkg_b2",128,0,128,1000,-500,500);
th2_signal_vs_channel_bkg_b3 = new TH2D("th2_signal_vs_channel_bkg_b3","th2_signal_vs_channel_bkg_b3",128,0,128,1000,-500,500);
th2_resid_vs_channel_bkg_b0 = new TH2D("th2_resid_vs_channel_bkg_b0","th2_resid_vs_channel_bkg_b0",128,0,128,300,-.15,.15);
th2_resid_vs_channel_bkg_b1 = new TH2D("th2_resid_vs_channel_bkg_b1","th2_resid_vs_channel_bkg_b1",128,0,128,300,-.15,.15);
th2_resid_vs_channel_bkg_b2 = new TH2D("th2_resid_vs_channel_bkg_b2","th2_resid_vs_channel_bkg_b2",128,0,128,300,-.15,.15);
th2_resid_vs_channel_bkg_b3 = new TH2D("th2_resid_vs_channel_bkg_b3","th2_resid_vs_channel_bkg_b3",128,0,128,300,-.15,.15);
for (int iii=0; iii<=29; iii++){
sprintf(histname,"h_beamSignal_b0_%d",iii);
h_beamSignal_b0[iii] = new TH1D(histname,histname,200,0,100000);
sprintf(histname,"h_beamSignal_b1_%d",iii);
h_beamSignal_b1[iii] = new TH1D(histname,histname,200,0,100000);
sprintf(histname,"h_beamSignal_b2_%d",iii);
h_beamSignal_b2[iii] = new TH1D(histname,histname,200,0,100000);
sprintf(histname,"h_beamSignal_b3_%d",iii);
h_beamSignal_b3[iii] = new TH1D(histname,histname,200,0,100000);
}
for (int iii=0; iii<=29; iii++){
sprintf(histname,"h_beamSignal_cor_b0_%d",iii);
h_beamSignal_cor_b0[iii] = new TH1D(histname,histname,200,0,10);
sprintf(histname,"h_beamSignal_cor_b1_%d",iii);
h_beamSignal_cor_b1[iii] = new TH1D(histname,histname,200,0,10.);
sprintf(histname,"h_beamSignal_cor_b2_%d",iii);
h_beamSignal_cor_b2[iii] = new TH1D(histname,histname,200,0,10.);
sprintf(histname,"h_beamSignal_cor_b3_%d",iii);
h_beamSignal_cor_b3[iii] = new TH1D(histname,histname,200,0,10.);
sprintf(histname,"h_beamSignal_ic1_%d",iii);
h_beamSignal_ic1[iii] = new TH1D(histname,histname,200,0,1000.);
}
rootTree = new TTree("t","HIT Data Root Tree");
rootTree->Branch("beamPosX_reg_b0",&beamPosX_reg_b0,"beamPosX_reg_b0/D");
rootTree->Branch("beamPosX_reg_b1",&beamPosX_reg_b1,"beamPosX_reg_b1/D");
rootTree->Branch("beamPosX_reg_b2",&beamPosX_reg_b2,"beamPosX_reg_b2/D");
rootTree->Branch("beamPosX_reg_b3",&beamPosX_reg_b3,"beamPosX_reg_b3/D");
rootTree->Branch("beamFocusX_reg_b0",&beamFocusX_reg_b0,"beamFocusX_reg_b0/D");
rootTree->Branch("beamFocusX_reg_b1",&beamFocusX_reg_b1,"beamFocusX_reg_b1/D");
rootTree->Branch("beamFocusX_reg_b2",&beamFocusX_reg_b2,"beamFocusX_reg_b2/D");
rootTree->Branch("beamFocusX_reg_b3",&beamFocusX_reg_b3,"beamFocusX_reg_b3/D");
rootTree->Branch("beamPeakX_reg_b0",&beamPeakX_reg_b0,"beamPeakX_reg_b0/D");
rootTree->Branch("beamPeakX_reg_b1",&beamPeakX_reg_b1,"beamPeakX_reg_b1/D");
rootTree->Branch("beamPeakX_reg_b2",&beamPeakX_reg_b2,"beamPeakX_reg_b2/D");
rootTree->Branch("beamPeakX_reg_b3",&beamPeakX_reg_b3,"beamPeakX_reg_b3/D");
rootTree->Branch("beamRsqr_b0",&beamRsqr_b0,"beamRsqr_b0/D");
rootTree->Branch("beamRsqr_b1",&beamRsqr_b1,"beamRsqr_b1/D");
rootTree->Branch("beamRsqr_b2",&beamRsqr_b2,"beamRsqr_b2/D");
rootTree->Branch("beamRsqr_b3",&beamRsqr_b3,"beamRsqr_b3/D");
rootTree->Branch("sp2_ps",&sp2_ps,"sp2_ps/D");
rootTree->Branch("sp2_air",&sp2_air,"sp2_air/D");
rootTree->Branch("mybeta",&mybeta,"mybeta/D");
rootTree->Branch("beamSkewX_b0",&beamSkewX_b0,"beamSkewX_b0/D");
rootTree->Branch("beamSkewX_b1",&beamSkewX_b1,"beamSkewX_b1/D");
rootTree->Branch("beamSkewX_b2",&beamSkewX_b2,"beamSkewX_b2/D");
rootTree->Branch("beamSkewX_b3",&beamSkewX_b3,"beamSkewX_b3/D");
rootTree->Branch("beamKurtX_b0",&beamKurtX_b0,"beamKurtX_b0/D");
rootTree->Branch("beamKurtX_b1",&beamKurtX_b1,"beamKurtX_b1/D");
rootTree->Branch("beamKurtX_b2",&beamKurtX_b2,"beamKurtX_b2/D");
rootTree->Branch("beamKurtX_b3",&beamKurtX_b3,"beamKurtX_b3/D");
rootTree->Branch("beamNumX_b0",&numtocalc_b0,"beamNumX_b0/I");
rootTree->Branch("beamNumX_b1",&numtocalc_b1,"beamNumX_b1/I");
rootTree->Branch("beamNumX_b2",&numtocalc_b2,"beamNumX_b2/I");
rootTree->Branch("beamNumX_b3",&numtocalc_b3,"beamNumX_b3/I");
rootTree->Branch("beamSignal_b0",&signal_b0,"beamSignal_b0/D");
rootTree->Branch("beamSignal_b1",&signal_b1,"beamSignal_b1/D");
rootTree->Branch("beamSignal_b2",&signal_b2,"beamSignal_b2/D");
rootTree->Branch("beamSignal_b3",&signal_b3,"beamSignal_b3/D");
rootTree->Branch("beamSignal_cor_b0",&signal_cor_b0,"beamSignal_cor_b0/D");
rootTree->Branch("beamSignal_cor_b1",&signal_cor_b1,"beamSignal_cor_b1/D");
rootTree->Branch("beamSignal_cor_b2",&signal_cor_b2,"beamSignal_cor_b2/D");
rootTree->Branch("beamSignal_cor_b3",&signal_cor_b3,"beamSignal_cor_b3/D");
rootTree->Branch("eventID",&eventID,"eventID/I");
//rootTree->Branch("board_b0_ch",&board_b0_ch,"board_b0_ch[128]/D");
//rootTree->Branch("board_b1_ch",&board_b1_ch,"board_b1_ch[128]/D");
//rootTree->Branch("board_b2_ch",&board_b2_ch,"board_b2_ch[128]/D");
//rootTree->Branch("board_b3_ch",&board_b3_ch,"board_b3_ch[128]/D");
rootTree->Branch("arrayavg_b0",&arrayavg_b0,"arrayavg_b0/D");
rootTree->Branch("arrayavg_b1",&arrayavg_b1,"arrayavg_b1/D");
rootTree->Branch("arrayavg_b2",&arrayavg_b2,"arrayavg_b2/D");
rootTree->Branch("arrayavg_b3",&arrayavg_b3,"arrayavg_b3/D");
rootTree->Branch("M1elements_b0",&M1elements_b0,"M1elements_b0[9]/D");
rootTree->Branch("M1elements_b1",&M1elements_b1,"M1elements_b1[9]/D");
rootTree->Branch("M1elements_b2",&M1elements_b2,"M1elements_b2[9]/D");
rootTree->Branch("M1elements_b3",&M1elements_b3,"M1elements_b3[9]/D");
rootTree->Branch("goodevents_b0",&goodevents_b0,"goodevents_b0/I");
rootTree->Branch("goodevents_b1",&goodevents_b1,"goodevents_b1/I");
rootTree->Branch("goodevents_b2",&goodevents_b2,"goodevents_b2/I");
rootTree->Branch("goodevents_b3",&goodevents_b3,"goodevents_b3/I");
double mytime;
// import and align matching ethercat data
TTree *tree2 = new TTree("t2", "t2");
if(ethercat){
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();
}
if(ethercat){
rootTree->Branch("ic1_avg", &ic1_avg, "ic1_avg/D");
rootTree->Branch("mw1_posx", &mw1_posx_avg, "mw1_posx/D");
rootTree->Branch("mw1_posy", &mw1_posy_avg, "mw1_posy/D");
rootTree->Branch("energy_index", &energy_index, "energy_index/D");
rootTree->Branch("intensity_index", &intensity_index, "intensity_index/D");
rootTree->Branch("ionsort", &ionsort, "ionsort/D");
timeoffset = -0.050; //offset between ic and bpm readouts
mwoffset = 2; // offset for timestamped event. MW chamber position correlation seems to be better in other windows
timewindow = -0.0999; //should be a negative number. window size in time to average over.
timeoffset2 = -0.00; //offset between ic and bpm readouts
timewindow2 = -0.0999; //should be a negative number. window size in time to average over.
//////////////////////////////////////////////////////////////////
///// ETHERCAT DATA
//////////////////////////////////////////////////////////////////
// tree->SetBranchAddress("time", &time);
tree2->SetBranchAddress("RELTIME2", &rel_time2);
tree2->SetBranchAddress("TIME2", &time2);
tree2->SetBranchAddress("IC1", &ic1);
// tree2->SetBranchAddress("IC2", &ic2);
//tree2->SetBranchAddress("MW1_FOCUSX", &mw1_focusx);
//tree2->SetBranchAddress("MW1_FOCUSY", &mw1_focusy);
///tree2->SetBranchAddress("MW2_FOCUSX", &mw2_focusx);
///tree2->SetBranchAddress("MW2_FOCUSY", &mw2_focusy);
tree2->SetBranchAddress("MW1_POSX", &mw1_posx);
tree2->SetBranchAddress("MW1_POSY", &mw1_posy);
//tree2->SetBranchAddress("MW2_POSX", &mw2_posx);
//tree2->SetBranchAddress("MW2_POSY", &mw2_posy);
tree2->SetBranchAddress("ENERGY_INDEX", &energy_index);
tree2->SetBranchAddress("INTENSITY_INDEX", &intensity_index);
tree2->SetBranchAddress("ION-SORT", &ionsort);
tree2->SetBranchAddress("ANALOG_IN1", &analog_in1);
}
int currentEntryTree2 = 1;
// int nevents = tree->GetEntries();
int nevents2 = tree2->GetEntries();
int icCounter = 0;
int count = 0;
int count2 = 0;
}
}
int analyse(int argc, char **argv)
{
//loop through recorded data frames
for (int frame = 0; frame<fileframesize; frame++){
// if (icCounter>10000) break;
// if(frame>2000&&frame<120000) continue;
eventID = frame;
if (ethercat){
if (timestampfile) timestampfile >> time ;
//printf("%i %f\n", eventID, time);
count= 0;
count2 = 0;
ic1_avg = 0.;
mw1_posx_avg = 0.;
mw1_posy_avg = 0.;
tree2->GetEntry(currentEntryTree2);
/* if (frame % 100 == 0)
{
printf("merging event %d ,", frame);
printf("Time %f \n", time);
printf("Entry hit %d ,", currentEntryTree2);
printf("Time hit %f \n", time2);
}
*/
while (time2 < mytime + timeoffset && currentEntryTree2 < nevents2 )
{
if (time2 - mytime - timeoffset > timewindow)
{
tree2->GetEntry(currentEntryTree2);
if (ic1>0.0) ic1_avg += ic1;
if (ic1>0.0) count++;
// if (frame % 200 == 0) printf("%i %2.3f %2.3f %2.3f %2.3f %2.3f %i \n", count, ic1, time2, mytime, timeoffset, time2 - mytime - timeoffset, mwoffset);
}
tree2->GetEntry(currentEntryTree2);
if ( time2 - mytime - timeoffset2 > timewindow2)
{
tree2->GetEntry(currentEntryTree2 + mwoffset); //why currentEtryTree2-4?
mw1_posx_avg += mw1_posx;
mw1_posy_avg += mw1_posy;
count2++;
// if (i % 2000 == 0) printf("%i %2.3f %2.3f %2.3f %2.3f\n", count2, mytime, time2, ic1, mw1_posx);
//if (i % 2001 == 0) printf("%i %2.3f %2.3f %2.3f %2.3f\n", count2, mytime, time2, ic1, mw1_posx);
// if (i % 2002 == 0) printf("%i %2.3f %2.3f %2.3f %2.3f\n", count2, mytime, time2, ic1, mw1_posx);
// printf("%i %2.3f %2.3f %2.3f %2.3f\n", count2, mytime, time2, ic1, mw1_posx);
}
// currentEntryTree2++;
tree2->GetEntry(currentEntryTree2);
currentEntryTree2++;
if (count2>0){
mw1_posx_avg /= double(count2); //the positions are weighted by the charge
mw1_posy_avg /= double(count2);
}
if(count>0){
ic1_avg /= double(count);
if (ic1_avg>1.) icCounter++;
// if (frame % 2000 == 0) printf("%i %f.2 %i \n", count, ic1_avg, icCounter);
}
// std::cout << ic1_avg << " " << icCounter << std::endl;
}
/////////end of ethercat matching
if (floor(energy_index/10.)>=0&&floor(energy_index/10)<=26){energy_indexbin = floor(energy_index/10.);}
else {energy_indexbin=27;}
if (ionsort==3){
//protons
energy = energy_list_p[ energy_indexbin];
intensity = intensity_list_p[ int(intensity_index) ];
sp_air = f_sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_sp_air->Eval(energy);
sp2_ps = f_sp_ps->Eval(energy);
lorentz_gamma = energy/938.272 +1.; // lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 1.;
}
else if (ionsort==1){
//helium
energy = energy_list_he[ energy_indexbin ];
intensity = intensity_list_he[ int(intensity_index) ];
sp_air = f_h2sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_h2sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_h2sp_air->Eval(energy);
sp2_ps = f_h2sp_ps->Eval(energy);
lorentz_gamma = energy/932.1055 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 4.;
}
else if (ionsort==4){
//carbon
energy = energy_list_c[ energy_indexbin];
intensity = intensity_list_c[ int(intensity_index) ];
sp_air = f_c12sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_c12sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_c12sp_air->Eval(energy);
sp2_ps = f_c12sp_ps->Eval(energy);
lorentz_gamma = energy/932.3539 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 36.;
}
else if (ionsort==2){
//oxygen
energy = energy_list_o[ energy_indexbin];
intensity = intensity_list_o[ int(intensity_index) ];
sp_air = f_o16sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_o16sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_o16sp_air->Eval(energy);
sp2_ps = f_o16sp_ps->Eval(energy);
lorentz_gamma = energy/931.4418 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 64.;
}
else {
sp_air = -1.;
sp_ps = -1.;
sp2_air = -1.;
sp2_ps = -1.;
lorentz_gamma = 0.;
mybeta = 0.;
zsqr = -1.;
}
// intcorr = (Intensity /19.8E6) / (totalcurrent/sp2_air/totaltime);
intcorr = (intensity / 19.8E6) / (ic1_avg/sp2_air)/ 3.11;
// if (frame>2000&&ic1_avg<0.5) continue;
}
///////////////////////////////////
//// HIT DATA
///////////////////////////////////
if (frame>2000 && beam_wason == true && beam_off == true) { beamoffcounter++;
// cout << beamoffcounter<< endl;
}
if (samplenoise==false && frame> 2000 && beam_wason == true && beam_off == true && beamoffcounter==5000) {
samplenoise = true;
}
//board_b 0
board_b=0;
signal_b0 = 0.;
maxchannelamp_b0 = 0.;
file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
file.read ((char*)dataptr ,sizeof(BufferData));
if (dataptr->sync_frame.device_nr==board_b){
// std::cout << frame << " " << dataptr->sync_frame.device_nr << " " << dataptr->sync_frame.local_ctr << " " << dataptr->sync_frame.global_ctr << " " << dataptr->sync_frame.sma_state << std::endl;
if (frame%10000==0) std::cout << "Frame: " << frame << " (" <<double(frame)/double(fileframesize)*100.0 << "%)" << std::endl;
arrayavg_b0 = 0;
//shift the events of the rolling average of the signal one left
for (int j = 1; j<length;j++){
array_b0[j-1] = array_b0[j];
arrayavg_b0 += array_b0[j-1];
// cout << array_b0[j-1] << " ";
}
array_b0[length-1] = 0.;
// cout << endl;
///loop over the channels
for (int i = 0;i<128;i++){
board_b0_ch[i] = dataptr->sensor_data[i];
if (frame<1000 || samplenoise==true){
if (samplenoise==true&& beamoffcounter==5000) {
board_b0_ch_bkg[i]=0.;
if(i==0) cout << frame << " Resampling noise " << beamoffcounter << endl;
}
board_b0_ch_bkg[i] += board_b0_ch[i]/1000.; //find baseline from the average of first 1000 events
if (samplenoise==true && beamoffcounter==6000) {
beam_wason = false;
samplenoise =false;
cout << frame << " End sampling." << endl;
}
}
else if (frame>=1000&& frame<2000){
board_b0_ch[i] -= board_b0_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
th2_signal_vs_channel_bkg_b0->Fill(i,board_b0_ch[i]);
if (i==1) graph_bkg_b0->SetPoint(graph_bkg_b0->GetN(), frame, board_b0_ch[i]);
}
else if (frame>=2000 && samplenoise==false) {
board_b0_ch[i]-=board_b0_ch_bkg[i]; // the rest background subtracted
// board_b0_ch[i]*=calibration_b0[i]; //calibration factor
th2_signal_vs_channel_b0->Fill(i,board_b0_ch[i]);
//find the peak channel
if (board_b0_ch[i]> maxchannelamp_b0) {
maxchannel_b0 = i;
maxchannelamp_b0 = board_b0_ch[i];
// cout << maxchannel_b0 << " " <<maxchannelamp_b0 << endl;
}
///add new event to the rolling average.
if( board_b0_ch[i]>50.0){
array_b0[length-1] += board_b0_ch[i];
}
////////////////////////////////////////////
}
}//end of 128
// th1_signal_b0->Fill(signal_b0);
if (frame>2000) {
arrayavg_b0 += array_b0[length-1]/length;
if (arrayavg_b0>= 20000) {
beamoncounter++;
}
if (beamoncounter>100 && beam_wason == false && arrayavg_b0>= 20000) {
beam_off = false;
beam_wason = true;
cout << frame << " Beam on " << arrayavg_b0 << endl;
}
if (arrayavg_b0< 10000) {
if(beam_wason== true && beam_off==false ) {
cout << frame << "Beam went OFF!" << arrayavg_b0 << endl;
beamoffcounter = 0;
}
beam_off = true;
goodevents_b0=0;
goodevents_b1=0;
goodevents_b2=0;
goodevents_b3=0;
}
}
// cout << frame << " Rolling average: " << arrayavg_b0 << endl;// " ic: " << ic1_avg << endl;
}
else {
std::cout << "Error." << std::endl;
}
//board_b 1
board_b=1;
signal_b1 = 0.;
maxchannelamp_b1 = 0.;
file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
file.read ((char*)dataptr ,sizeof(BufferData));
if (dataptr->sync_frame.device_nr==board_b){
// std::cout << frame << " " << dataptr->sync_frame.device_nr << " " << dataptr->sync_frame.local_ctr << " " << dataptr->sync_frame.global_ctr << " " << dataptr->sync_frame.sma_state << std::endl;
arrayavg_b1 = 0;
//shift the events of the rolling average of the signal one left
for (int j = 1; j<length;j++){
array_b1[j-1] = array_b1[j];
arrayavg_b1 += array_b1[j-1];
// cout << array_b1[j-1] << " ";
}
array_b1[length-1] = 0.;
// cout << endl;
///loop over the channels
for (int i = 0;i<128;i++){
board_b1_ch[i] = dataptr->sensor_data[i];
if (frame<1000){
board_b1_ch_bkg[i] += board_b1_ch[i]/1000.; //find baseline from the average of first 1000 events
}
else if (frame>=1000&& frame<2000){
board_b1_ch[i] -= board_b1_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
th2_signal_vs_channel_bkg_b1->Fill(i,board_b1_ch[i]);
if (i==1) graph_bkg_b1->SetPoint(graph_bkg_b1->GetN(), frame, board_b1_ch[i]);
}
else if (frame>=2000) {
board_b1_ch[i]-=board_b1_ch_bkg[i]; // the rest background subtracted
// board_b1_ch[i]*=calibration_b1[i]; //calibration factor
th2_signal_vs_channel_b1->Fill(i,board_b1_ch[i]);
//find the peak channel
if (board_b1_ch[i]> maxchannelamp_b1) {
maxchannel_b1 = i;
maxchannelamp_b1 = board_b1_ch[i];
// cout << maxchannel_b1 << " " <<maxchannelamp_b1 << endl;
}
///add new event to the rolling average.
if( board_b1_ch[i]>50.0){
array_b1[length-1] += board_b1_ch[i];
}
////////////////////////////////////////////
}
}//end of 128
// th1_signal_b1->Fill(signal_b1);
arrayavg_b1 += array_b1[length-1]/length;
// cout << frame << " Rolling average: " << arrayavg_b1 << endl;// " ic: " << ic1_avg << endl;
}
else {
std::cout << "Error." << std::endl;
}
//board_b 2
board_b=2;
signal_b2 = 0.;
maxchannelamp_b2 = 0.;
file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
file.read ((char*)dataptr ,sizeof(BufferData));
if (dataptr->sync_frame.device_nr==board_b){
// std::cout << frame << " " << dataptr->sync_frame.device_nr << " " << dataptr->sync_frame.local_ctr << " " << dataptr->sync_frame.global_ctr << " " << dataptr->sync_frame.sma_state << std::endl;
arrayavg_b2 = 0;
//shift the events of the rolling average of the signal one left
for (int j = 1; j<length;j++){
array_b2[j-1] = array_b2[j];
arrayavg_b2 += array_b2[j-1];
// cout << array_b2[j-1] << " ";
}
array_b2[length-1] = 0.;
// cout << endl;
///loop over the channels
for (int i = 0;i<128;i++){
board_b2_ch[i] = dataptr->sensor_data[i];
if (frame<1000){
board_b2_ch_bkg[i] += board_b2_ch[i]/1000.; //find baseline from the average of first 1000 events
}
else if (frame>=1000&& frame<2000){
board_b2_ch[i] -= board_b2_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
th2_signal_vs_channel_bkg_b2->Fill(i,board_b2_ch[i]);
if (i==1) graph_bkg_b2->SetPoint(graph_bkg_b2->GetN(), frame, board_b2_ch[i]);
}
else if (frame>=2000) {
board_b2_ch[i]-=board_b2_ch_bkg[i]; // the rest background subtracted
// board_b2_ch[i]*=calibration_b2[i]; //calibration factor
th2_signal_vs_channel_b2->Fill(i,board_b2_ch[i]);
//find the peak channel
if (board_b2_ch[i]> maxchannelamp_b2) {
maxchannel_b2 = i;
maxchannelamp_b2 = board_b2_ch[i];
// cout << maxchannel_b2 << " " <<maxchannelamp_b2 << endl;
}
///add new event to the rolling average.
if( board_b2_ch[i]>50.0){
array_b2[length-1] += board_b2_ch[i];
}
////////////////////////////////////////////
}
}//end of 128
// th1_signal_b2->Fill(signal_b2);
arrayavg_b2 += array_b2[length-1]/length;
// cout << frame << " Rolling average: " << arrayavg_b2 << endl;// " ic: " << ic1_avg << endl;
}
else {
std::cout << "Error." << std::endl;
}
//board_b 3
board_b=3;
signal_b3 = 0.;
maxchannelamp_b3 = 0.;
file.seekg(framestart*sizeof(BufferData)+board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
file.read ((char*)dataptr ,sizeof(BufferData));
if (dataptr->sync_frame.device_nr==board_b){
// std::cout << frame << " " << dataptr->sync_frame.device_nr << " " << dataptr->sync_frame.local_ctr << " " << dataptr->sync_frame.global_ctr << " " << dataptr->sync_frame.sma_state << std::endl;
arrayavg_b3 = 0;
//shift the events of the rolling average of the signal one left
for (int j = 1; j<length;j++){
array_b3[j-1] = array_b3[j];
arrayavg_b3 += array_b3[j-1];
// cout << array_b3[j-1] << " ";
}
array_b3[length-1] = 0.;
// cout << endl;
///loop over the channels
for (int i = 0;i<128;i++){
board_b3_ch[i] = dataptr->sensor_data[i];
if (frame<1000){
board_b3_ch_bkg[i] += board_b3_ch[i]/1000.; //find baseline from the average of first 1000 events
}
else if (frame>=1000&& frame<2000){
board_b3_ch[i] -= board_b3_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
th2_signal_vs_channel_bkg_b3->Fill(i,board_b3_ch[i]);
if (i==1) graph_bkg_b3->SetPoint(graph_bkg_b3->GetN(), frame, board_b3_ch[i]);
}
else if (frame>=2000) {
board_b3_ch[i]-=board_b3_ch_bkg[i]; // the rest background subtracted
// board_b3_ch[i]*=calibration_b3[i]; //calibration factor
th2_signal_vs_channel_b3->Fill(i,board_b3_ch[i]);
//find the peak channel
if (board_b3_ch[i]> maxchannelamp_b3) {
maxchannel_b3 = i;
maxchannelamp_b3 = board_b3_ch[i];
// cout << maxchannel_b3 << " " <<maxchannelamp_b3 << endl;
}
///add new event to the rolling average.
if( board_b3_ch[i]>50.0){
array_b3[length-1] += board_b3_ch[i];
}
////////////////////////////////////////////
}
}//end of 128
// th1_signal_b3->Fill(signal_b3);
arrayavg_b3 += array_b3[length-1]/length;
// cout << frame << " Rolling average: " << arrayavg_b3 << endl;// " ic: " << ic1_avg << endl;
}
else {
std::cout << "Error." << std::endl;
}
////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
//start the signal analysis
////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
if (frame>=2000&&beam_off == false&&beam_wason == true){
////////////////////////////////////////////////////////////
//find the approx FWHM
////////////////////////////////////////////////////////////
nfwhm_b0 =0;
nfwhm_b1 =0;
nfwhm_b2 =0;
nfwhm_b3 =0;
for (int i = 0;i<128;i++){
if (board_b0_ch[i] > maxchannelamp_b0/3.) nfwhm_b0++; //switched to a lower threshold 1/3 instead of 1/2
if (board_b1_ch[i] > maxchannelamp_b1/3.) nfwhm_b1++;
if (board_b2_ch[i] > maxchannelamp_b2/3.) nfwhm_b2++;
if (board_b3_ch[i] > maxchannelamp_b3/3.) nfwhm_b3++;
signal_gsl_b0[i] = 0.;
signal_gsl_b1[i] = 0.;
signal_gsl_b2[i] = 0.;
signal_gsl_b3[i] = 0.;
}
////////////////////////////////////////////////////////////
//integrate under the approximate peak
//build the channel list for statistical analysis
////////////////////////////////////////////////////////////
numtocalc_b0 = 0;
numtocalc_b1 = 0;
numtocalc_b2 = 0;
numtocalc_b3 = 0;
signal_b0 = 0. - signaloffset_b0[energy_indexbin];
// cout << energy_indexbin << endl;
for (int i = maxchannel_b0-nfwhm_b0 ; i <= maxchannel_b0 + nfwhm_b0; i++){
if (i>=0 && i<=127 && board_b0_ch[i]>0){
signal_b0 +=board_b0_ch[i];
signal_gsl_b0[numtocalc_b0]=board_b0_ch[i];
// channellist_gsl_b0[numtocalc_b0] = i;
channellist_gsl_b0[numtocalc_b0] = pos[i];
numtocalc_b0++;
}
}
th1_signal_b0->Fill(signal_b0); //cout << signal_b0 << endl;
if (ic1_avg>0.1){
h_beamSignal_b0[energy_indexbin]->Fill(signal_b0);
goodevents_b0++;
}
signal_b1 = 0. - signaloffset_b1[energy_indexbin];;
for (int i = maxchannel_b1-nfwhm_b1 ; i <= maxchannel_b1 + nfwhm_b1; i++){
if (i>=0 && i<=127&&board_b1_ch[i]>0){
signal_b1 +=board_b1_ch[i] ;
signal_gsl_b1[numtocalc_b1]=board_b1_ch[i];
// channellist_gsl_b1[numtocalc_b1] = i;
channellist_gsl_b1[numtocalc_b1] = pos[i];
numtocalc_b1++;
}
}
th1_signal_b1->Fill(signal_b1);
if (ic1_avg>0.1){
h_beamSignal_b1[energy_indexbin]->Fill(signal_b1);
goodevents_b1++;
}
signal_b2 = 0. - signaloffset_b2[energy_indexbin];;
for (int i = maxchannel_b2-nfwhm_b2 ; i <= maxchannel_b2 + nfwhm_b2; i++){
if (i>=0 && i<=127&&board_b2_ch[i]>0){
signal_b2 +=board_b2_ch[i];
signal_gsl_b2[numtocalc_b2]=board_b2_ch[i];
// channellist_gsl_b2[numtocalc_b2] = i;
channellist_gsl_b2[numtocalc_b2] = pos[i];
numtocalc_b2++;
}
}
th1_signal_b2->Fill(signal_b2);
if (ic1_avg>0.1) {
h_beamSignal_b2[energy_indexbin]->Fill(signal_b2);
goodevents_b2++;
}
signal_b3 = 0. - signaloffset_b3[energy_indexbin];
for (int i = maxchannel_b3-nfwhm_b3 ; i <= maxchannel_b3 + nfwhm_b3; i++){
if (i>=0 && i<=127&&board_b3_ch[i]>0){
signal_b3 +=board_b3_ch[i];
signal_gsl_b3[numtocalc_b3]=board_b3_ch[i];
//channellist_gsl_b3[numtocalc_b3] = i;
channellist_gsl_b3[numtocalc_b3] = pos[i];
numtocalc_b3++;
}
}
th1_signal_b3->Fill(signal_b3);
if (ic1_avg>0.1) {
h_beamSignal_b3[energy_indexbin]->Fill(signal_b3);
h_beamSignal_ic1[energy_indexbin]->Fill(ic1_avg);
goodevents_b3++;
}
signal_cor_b0 = signal_b0*intcorr/(intensity/1.E6)/sp2_ps;
signal_cor_b1 = signal_b1*intcorr/(intensity/1.E6)/sp2_ps;
signal_cor_b2 = signal_b2*intcorr/(intensity/1.E6)/sp2_ps;
signal_cor_b3 = signal_b3*intcorr/(intensity/1.E6)/sp2_ps;
if (signal_cor_b0 >0.05 && signal_cor_b0 < 10.) {h_beamSignal_cor_b0[energy_indexbin]->Fill(signal_cor_b0);}
if (signal_cor_b1 >0.05 && signal_cor_b1 < 10.) {h_beamSignal_cor_b1[energy_indexbin]->Fill(signal_cor_b1);}
if (signal_cor_b2 >0.05 && signal_cor_b2 < 10.) {h_beamSignal_cor_b2[energy_indexbin]->Fill(signal_cor_b2);}
if (signal_cor_b3 >0.05 && signal_cor_b3 < 10.) {h_beamSignal_cor_b3[energy_indexbin]->Fill(signal_cor_b3);}
///////////////// linear regression using Integration by parts of gaussian function.
if (doLinInt&&numtocalc_b0>=3&&numtocalc_b1>=3&&numtocalc_b2>=3&&numtocalc_b3>=3){
for (int j = 0; j<=4; j++){
if (j== 0) {
copy(signal_gsl_b0,signal_gsl_b0+128, signal_reg);
copy(channellist_gsl_b0,channellist_gsl_b0+128, channellist_reg);
numtocalc_reg = numtocalc_b0;
}
else if (j== 1) {
copy(signal_gsl_b1,signal_gsl_b1+128, signal_reg);
copy(channellist_gsl_b1,channellist_gsl_b1+128, channellist_reg);
numtocalc_reg = numtocalc_b1;
}
else if (j== 2) {
copy(signal_gsl_b2,signal_gsl_b2+128, signal_reg);
copy(channellist_gsl_b2,channellist_gsl_b2+128, channellist_reg);
numtocalc_reg = numtocalc_b2;
}
else if (j== 3) {
copy(signal_gsl_b3,signal_gsl_b3+128, signal_reg);
copy(channellist_gsl_b3,channellist_gsl_b3+128, channellist_reg);
numtocalc_reg = numtocalc_b3;
}
SumY = 0.;
SumS = 0.;
SumT = 0.;
SumS2 = 0.;
SumST = 0.;
SumT2 = 0.;
SumYS = 0.;
SumYT = 0.;
b_den = 0.;
b_num = 0.;
b = 0.;
SumYYM = 0.;
SumYYP = 0.;
MeanY = 0.;
for(int k=0; k<numtocalc_b0;k++){
if (k==0){
S[k]=0.; T[k]=0.;
}
else{
S[k] = S[k-1]+0.5*( signal_reg[k] + signal_reg[k-1] ) * ( channellist_reg[k] - channellist_reg[k-1] );
T[k] = T[k-1]+0.5*( channellist_reg[k] * signal_reg[k] + channellist_reg[k-1] * signal_reg[k-1] ) * ( channellist_reg[k] - channellist_reg[k-1] );
}
// cout << S[k] << " " << T[k] << endl;
SumS += S[k]; SumT += T[k];
SumY += signal_reg[k];
SumS2 += S[k]*S[k]; SumST += S[k]*T[k]; SumT2 += T[k]*T[k];
SumYS += signal_reg[k]*S[k];
SumYT += signal_reg[k]*T[k];
MeanY+=signal_reg[k];
}
MeanY/=numtocalc_reg;
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) = numtocalc_reg;
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<numtocalc_reg;k++){
b_num += exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) * signal_reg[k];
b_den += exp(-2*p*(channellist_reg[k]-c)*(channellist_reg[k]-c));
}
b = b_num/b_den;
//calculate R-squared
///residual plots (should be normally distributed). They are normalised to the standard error of the baseline (17 ADC units)
for(int k=0; k<numtocalc_reg;k++){
SumYYM+= (signal_reg[k]-MeanY)*(signal_reg[k]-MeanY);
SumYYP+= (b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - MeanY )*(b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - MeanY );
if (j == 0) {
h_regresidual_b0->Fill((b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
th2_resid_vs_channel_bkg_b0->Fill(channellist_reg[k],(b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
}
if (j == 1) {
h_regresidual_b1->Fill((b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
th2_resid_vs_channel_bkg_b1->Fill(channellist_reg[k],(b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
}
if (j == 2) {
h_regresidual_b2->Fill((b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
th2_resid_vs_channel_bkg_b2->Fill(channellist_reg[k],(b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
}
if (j == 3) {
h_regresidual_b3->Fill((b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
th2_resid_vs_channel_bkg_b3->Fill(channellist_reg[k],(b*exp(-p*(channellist_reg[k]-c)*(channellist_reg[k]-c)) - signal_reg[k])/signal_reg[k]);
}
}
// cout << "R-squared = " << SumYYP/SumYYM << endl;
//final results
if (j == 0) {
beamFocusX_reg_b0 = 2.3548/sqrt(2*p);
beamPosX_reg_b0 = -ABC(1)/ ABC(0);
beamPeakX_reg_b0 = b;
beamRsqr_b0 = SumYYP/SumYYM;
beamSkewX_b0 = gsl_stats_wskew_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b0,beamFocusX_reg_b0/2.3548); //skewness (symmetry)
beamKurtX_b0 = gsl_stats_wkurtosis_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b0,beamFocusX_reg_b0/2.3548); //excess kurtosis (well behaved tails)
M1elements_b0[0] = SumT2;
M1elements_b0[1] = SumST;
M1elements_b0[2] = SumT;
M1elements_b0[3] = SumST;
M1elements_b0[4] = SumS2;
M1elements_b0[5] = SumS;
M1elements_b0[6] = SumT;
M1elements_b0[7] = SumS;
M1elements_b0[8] = numtocalc_reg;
}
else if (j == 1) {
beamFocusX_reg_b1 = 2.3548/sqrt(2*p);
beamPosX_reg_b1 = -ABC(1)/ ABC(0);
beamPeakX_reg_b1 = b;
beamRsqr_b1 = SumYYP/SumYYM;
beamSkewX_b1 = gsl_stats_wskew_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b1,beamFocusX_reg_b1/2.3548); //skewness (symmetry)
beamKurtX_b1 = gsl_stats_wkurtosis_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b1,beamFocusX_reg_b1/2.3548); //excess kurtosis (well behaved tails)
M1elements_b0[0] = SumT2;
M1elements_b1[1] = SumST;
M1elements_b1[2] = SumT;
M1elements_b1[3] = SumST;
M1elements_b1[4] = SumS2;
M1elements_b1[5] = SumS;
M1elements_b1[6] = SumT;
M1elements_b1[7] = SumS;
M1elements_b1[8] = numtocalc_reg;
}
else if (j == 2) {
beamFocusX_reg_b2 = 2.3548/sqrt(2*p);
beamPosX_reg_b2 = -ABC(1)/ ABC(0);
beamPeakX_reg_b2 = b;
beamRsqr_b2 = SumYYP/SumYYM;
beamSkewX_b2 = gsl_stats_wskew_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b2,beamFocusX_reg_b2/2.3548); //skewness (symmetry)
beamKurtX_b2 = gsl_stats_wkurtosis_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b2,beamFocusX_reg_b2/2.3548); //excess kurtosis (well behaved tails)
M1elements_b2[0] = SumT2;
M1elements_b2[1] = SumST;
M1elements_b2[2] = SumT;
M1elements_b2[3] = SumST;
M1elements_b2[4] = SumS2;
M1elements_b2[5] = SumS;
M1elements_b2[6] = SumT;
M1elements_b2[7] = SumS;
M1elements_b2[8] = numtocalc_reg;
}
else if (j == 3) {
beamFocusX_reg_b3 = 2.3548/sqrt(2*p);
beamPosX_reg_b3 = -ABC(1)/ ABC(0);
beamPeakX_reg_b3 = b;
beamRsqr_b3 = SumYYP/SumYYM;
beamSkewX_b3 = gsl_stats_wskew_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b3,beamFocusX_reg_b3/2.3548); //skewness (symmetry)
beamKurtX_b3 = gsl_stats_wkurtosis_m_sd(signal_reg,1,channellist_reg,1,numtocalc_reg,beamPosX_reg_b3,beamFocusX_reg_b3/2.3548); //excess kurtosis (well behaved tails)
M1elements_b0[0] = SumT2;
M1elements_b3[1] = SumST;
M1elements_b3[2] = SumT;
M1elements_b3[3] = SumST;
M1elements_b3[4] = SumS2;
M1elements_b3[5] = SumS;
M1elements_b3[6] = SumT;
M1elements_b3[7] = SumS;
M1elements_b3[8] = numtocalc_reg;
}
}
}
else if (numtocalc_b0>0&&numtocalc_b1>0) {
beamPosX_reg_b0 = gsl_stats_wmean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0); //calculate the weighted mean
beamFocusX_reg_b0 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_reg_b0); //Standard Deviation
beamFocusX_reg_b0 *=2.3548;//SD-->FWHM
beamPosX_reg_b1 = gsl_stats_wmean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1); //calculate the weighted mean
beamFocusX_reg_b1 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_reg_b1); //Standard Deviation
beamNumX_b1 = numtocalc_b1;
beamFocusX_reg_b1 *=2.3548;//SD-->FWHM
beamPosX_reg_b2 = gsl_stats_wmean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2); //calculate the weighted mean
beamFocusX_reg_b2 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_reg_b2); //Standard Deviation
beamNumX_b2 = numtocalc_b2;
beamFocusX_reg_b2 *=2.3548;//SD-->FWHM
beamPosX_reg_b3 = gsl_stats_wmean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3); //calculate the weighted mean
beamFocusX_reg_b3 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_reg_b3); //Standard Deviation
beamNumX_b3 = numtocalc_b3;
beamFocusX_reg_b3 *=2.3548;//SD-->FWHM
beamRsqr_b0 = .0;
beamRsqr_b1 = .0;
beamRsqr_b2 = .0;
beamRsqr_b3 = .0;
}
else
{
beamFocusX_reg_b0 = 0.;
beamFocusX_reg_b1 = 0.;
beamFocusX_reg_b1 = 0.;
beamFocusX_reg_b1 = 0.;
beamPosX_reg_b0 = 0.;
beamPosX_reg_b1 = 0.;
beamPosX_reg_b2 = 0.;
beamPosX_reg_b3 = 0.;
beamRsqr_b0 = .0;
beamRsqr_b1 = .0;
beamRsqr_b2 = .0;
beamRsqr_b3 = .0;
}
rootTree->Fill();
}
if (frameend>0 && frame+framestart>=frameend) break;
}//end of loop over frames
graph_bkg_b0->SetName("graph_bkg_b0"); // graph_bkg_b0->Write();
graph_bkg_b1->SetName("graph_bkg_b1"); // graph_bkg_b1->Write();
graph_bkg_b2->SetName("graph_bkg_b2"); // graph_bkg_b2->Write();
graph_bkg_b3->SetName("graph_bkg_b3"); // graph_bkg_b3->Write();
// bic1->Fill();
// bmw1_posx->Fill();
//bmw1_posy->Fill();
///make reduced signal graph.
for (int jjj = 0; jjj<23;jjj++) {
// beamSignal_cor_b0_mean[jjj] = h_beamSignal_cor_b0[jjj]->GetMean();
// beamSignal_cor_b1_mean[jjj] = h_beamSignal_cor_b1[jjj]->GetMean();
// beamSignal_cor_b2_mean[jjj] = h_beamSignal_cor_b2[jjj]->GetMean();
// beamSignal_cor_b3_mean[jjj] = h_beamSignal_cor_b3[jjj]->GetMean();
//// get the median (0.5 quantile)
q = 0.5; // 0.5 for "median"
if (ionsort==3){
//protons
energy = energy_list_p[ jjj];
intensity = intensity_list_p[ int(intensity_index) ];
sp_air = f_sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_sp_air->Eval(energy);
sp2_ps = f_sp_ps->Eval(energy);
lorentz_gamma = energy/938.272 +1.; // lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 1.;
}
else if (ionsort==1){
//helium
energy = energy_list_he[ jjj ];
intensity = intensity_list_he[ int(intensity_index) ];
sp_air = f_h2sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_h2sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_h2sp_air->Eval(energy);
sp2_ps = f_h2sp_ps->Eval(energy);
lorentz_gamma = energy/932.1055 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 4.;
}
else if (ionsort==4){
//carbon
energy = energy_list_c[ jjj];
intensity = intensity_list_c[ int(intensity_index) ];
sp_air = f_c12sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_c12sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_c12sp_air->Eval(energy);
sp2_ps = f_c12sp_ps->Eval(energy);
lorentz_gamma = energy/932.3539 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 36.;
}
else if (ionsort==2){
//oxygen
energy = energy_list_o[ jjj];
intensity = intensity_list_o[ int(intensity_index) ];
sp_air = f_o16sp_air->Eval(energy)*(intensity/1.0E6);
sp_ps = f_o16sp_ps->Eval(energy)*(intensity/1.0E6);
sp2_air = f_o16sp_air->Eval(energy);
sp2_ps = f_o16sp_ps->Eval(energy);
lorentz_gamma = energy/931.4418 +1.;// lorentz_gamma = E_u/M_u +1
mybeta = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
zsqr = 64.;
}
else {
sp_air = -1.;
sp_ps = -1.;
sp2_air = -1.;
sp2_ps = -1.;
lorentz_gamma = 0.;
mybeta = 0.;
zsqr = -1.;
}
h_beamSignal_ic1[jjj]->ComputeIntegral(); //just a precaution
h_beamSignal_ic1[jjj]->GetQuantiles(1,&median, &q);
//median = h_beamSignal_ic1[jjj]->GetMean();
// intcorr = (Intensity /19.8E6) / (totalcurrent/sp2_air/totaltime);
if (median>0&&energy>0) {intcorr = (intensity / 19.8E6) / (median/sp2_air)/ 3.11;}
// else if (median>0&&energy>0) {intcorr = 0.; median = 0.;}
else {intcorr = 0.; median = 0.;}
cout << "energy: " << energy << endl;
cout << "intensity: " << intensity << endl;
cout << "sp2_air: " << sp2_air << endl;
cout << "median: " << median << endl;
cout << "intcorr: " << intcorr << endl;
mybeta_graph[jjj] = sqrt(1. - 1./(lorentz_gamma*lorentz_gamma));
h_beamSignal_b0[jjj]->ComputeIntegral(); //just a precaution
h_beamSignal_b0[jjj]->GetQuantiles(1,&median, &q);
beamSignal_cor_b0_mean[jjj] = median*intcorr /(intensity/1.E6)/sp2_ps;
//beamSignal_cor_b0_mean[jjj] = h_beamSignal_b0[jjj]->GetMean()*intcorr /(intensity/1.E6)/sp2_ps;
cout << beamSignal_cor_b0_mean[jjj] << endl;
h_beamSignal_b1[jjj]->ComputeIntegral(); //just a precaution
h_beamSignal_b1[jjj]->GetQuantiles(1,&median, &q);
beamSignal_cor_b1_mean[jjj] = median*intcorr /(intensity/1.E6)/sp2_ps;
//beamSignal_cor_b1_mean[jjj] = h_beamSignal_b1[jjj]->GetMean()*intcorr /(intensity/1.E6)/sp2_ps;
h_beamSignal_b2[jjj]->ComputeIntegral(); //just a precaution
h_beamSignal_b2[jjj]->GetQuantiles(1,&median, &q);
beamSignal_cor_b2_mean[jjj] = median*intcorr /(intensity/1.E6)/sp2_ps;
// beamSignal_cor_b2_mean[jjj] = h_beamSignal_b2[jjj]->GetMean()*intcorr /(intensity/1.E6)/sp2_ps;
h_beamSignal_b3[jjj]->ComputeIntegral(); //just a precaution
h_beamSignal_b3[jjj]->GetQuantiles(1,&median, &q);
beamSignal_cor_b3_mean[jjj] = median*intcorr /(intensity/1.E6)/sp2_ps;
//beamSignal_cor_b3_mean[jjj] = h_beamSignal_b3[jjj]->GetMean()*intcorr /(intensity/1.E6)/sp2_ps;
}
TGraph * graph_b0 = new TGraph(23,mybeta_graph, beamSignal_cor_b0_mean);
TGraph * graph_b1 = new TGraph(23,mybeta_graph, beamSignal_cor_b1_mean);
TGraph * graph_b2 = new TGraph(23,mybeta_graph, beamSignal_cor_b2_mean);
TGraph * graph_b3 = new TGraph(23,mybeta_graph, beamSignal_cor_b3_mean);
graph_b0->SetName("graph_b0");graph_b0->Write();
graph_b1->SetName("graph_b1");graph_b1->Write();
graph_b2->SetName("graph_b2");graph_b2->Write();
graph_b3->SetName("graph_b3");graph_b3->Write();
rootFile->Write();
rootFile->Close();
std::cout << eventID << " frames analysed." << std::endl;
return 0;
// }
// else {
// std::cerr << "Error 1" << std::endl;
// return 1;
// }
}