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.

914 lines
38 KiB

  1. /// to do:
  2. /// 1. correct for gap between arrays
  3. #define convert_cxx
  4. #include "hitutils.h"
  5. #include <string>
  6. #include <stdio.h>
  7. #include <iostream>
  8. #include <vector>
  9. #include <utility>
  10. #include <TH2.h>
  11. #include <TStyle.h>
  12. #include <TCanvas.h>
  13. #include <TFile.h>
  14. #include <TTree.h>
  15. #include <TSystemDirectory.h>
  16. #include <gsl/gsl_statistics.h>
  17. #include <math.h>
  18. #include <gsl/gsl_errno.h>
  19. #include <gsl/gsl_fft_complex.h>
  20. #include <TF1.h>
  21. #include <TGraphErrors.h>
  22. #include <gsl/gsl_sort.h>
  23. #include "GaussSmoothen.h"
  24. #include <TVector.h>
  25. using namespace std;
  26. bool smooth_on = true;
  27. int eventID = 0;
  28. double board_b0_ch[128];
  29. double board_b1_ch[128];
  30. double board_b2_ch[128];
  31. double board_b3_ch[128];
  32. double board_b0_ch_bkg[128];
  33. double board_b1_ch_bkg[128];
  34. double board_b2_ch_bkg[128];
  35. double board_b3_ch_bkg[128];
  36. double signal_b0 = 0.;
  37. double signal_b1 = 0.;
  38. double signal_b2 = 0.;
  39. double signal_b3 = 0.;
  40. double signal_b0_left = 0.;
  41. double signal_b1_left = 0.;
  42. double signal_b2_left = 0.;
  43. double signal_b3_left = 0.;
  44. double signal_b0_right = 0.;
  45. double signal_b1_right = 0.;
  46. double signal_b2_right = 0.;
  47. double signal_b3_right = 0.;
  48. double signal_gsl_b0[128];
  49. double signal_gsl_b1[128];
  50. double signal_gsl_b2[128];
  51. double signal_gsl_b3[128];
  52. double channellist_gsl_b0[128];
  53. double channellist_gsl_b1[128];
  54. double channellist_gsl_b2[128];
  55. double channellist_gsl_b3[128];
  56. double channellist[128];
  57. double pos[128];
  58. double maxchannelamp_b0 =0.;
  59. double maxchannelamp_b1 =0.;
  60. double maxchannelamp_b2 =0.;
  61. double maxchannelamp_b3 =0.;
  62. int maxchannel_b0 =0;
  63. int maxchannel_b1 =0;
  64. int maxchannel_b2 =0;
  65. int maxchannel_b3 =0;
  66. int numtocalc_b0 =0;
  67. int numtocalc_b1 =0;
  68. int numtocalc_b2 =0;
  69. int numtocalc_b3 =0;
  70. int nfwhm_b0 = 0;
  71. int nfwhm_b1 = 0;
  72. int nfwhm_b2 = 0;
  73. int nfwhm_b3 = 0;
  74. double beamPosX_b0,beamPosX_b1,beamPosX_b2,beamPosX_b3;
  75. double beamFocusX_b0,beamFocusX_b1,beamFocusX_b2,beamFocusX_b3;
  76. double beamSkewX_b0,beamSkewX_b1,beamSkewX_b2,beamSkewX_b3;
  77. double beamKurtX_b0,beamKurtX_b1,beamKurtX_b2,beamKurtX_b3;
  78. double beamNumX_b0,beamNumX_b1,beamNumX_b2,beamNumX_b3;
  79. double beamSidebandNoise_b0, beamSidebandNoise_b1, beamSidebandNoise_b2, beamSidebandNoise_b3 ;
  80. int sidenumtocalc_b0,sidenumtocalc_b1,sidenumtocalc_b2,sidenumtocalc_b3;
  81. size_t size = 5;
  82. vector<double> boxcar;
  83. vector<double> boxcarsort;
  84. vector<double> boxcarweight{1,3,5,3,1};
  85. vector<double> data(128);
  86. TVector sumvector_b0(128);
  87. TVector sumvector_b1(128);
  88. TVector sumvector_b2(128);
  89. TVector sumvector_b3(128);
  90. TVector * sumvector_b0_ptr = &sumvector_b0;
  91. TVector * sumvector_b1_ptr = &sumvector_b1;
  92. TVector * sumvector_b2_ptr = &sumvector_b2;
  93. TVector * sumvector_b3_ptr = &sumvector_b3;
  94. const int length = 100; //length of the rolling average
  95. double array_b0[length][128] = {{0.}};
  96. double array_b1[length][128] = {{0.}};
  97. double array_b2[length][128] = {{0.}};
  98. double array_b3[length][128] = {{0.}};
  99. double arrayavg_b0[128] = {0.};
  100. double arrayavg_b1[128] = {0.};
  101. double arrayavg_b2[128] = {0.};
  102. double arrayavg_b3[128] = {0.};
  103. double board_b0_smooth_ch[128];
  104. double board_b1_smooth_ch[128];
  105. double board_b2_smooth_ch[128];
  106. double board_b3_smooth_ch[128];
  107. bool graphsaved_b0 = false;
  108. bool graphsaved_b1 = false;
  109. bool graphsaved_b2 = false;
  110. bool graphsaved_b3 = false;
  111. TVector beamontime(0,3,0.,0.,0.,0.,"END");
  112. TVector * beamontime_ptr = &beamontime;
  113. double calibration_b0[128] = {0.680457, 0.70036, 0.729569, 0.768141, 0.776347, 0.785101, 0.799768, 0.819888, 0.817106, 0.830481, 0.872981, 0.866716, 0.898705, 0.8508, 0.851149, 0.913431, 0.908141, 0.920829, 0.954078, 1.0219, 1.05272, 1.01594, 1.09155, 1.07093, 1.06138, 1.04821, 1.04332, 1.07227, 1.09676, 1.13826, 1.16939, 1.16456, 1.11643, 1.16456, 1.08821, 1.11285, 1.0864, 1.05949, 1.10727, 1.13911, 1.13338, 1.1812, 1.14498, 1.15755, 1.20076, 1.1486, 1.15815, 1.13221, 1.18893, 1.20203, 1.21786, 1.28368, 1.08994, 1.09368, 1.12954, 1.11233, 1.11923, 1.13531, 1.13921, 1.12821, 1.10786, 1.14957, 1.29386, 1.3473, 1.01231, 0.992629, 0.983174, 0.987423, 0.990555, 0.987429, 1.01338, 1.02899, 1.03465, 1.0364, 1.03649, 1.04556, 1.0285, 1.0183, 1.0018, 1.00728, 1.01542, 1.00067, 0.995893, 0.998084, 0.969798, 0.982129, 0.957553, 0.975154, 0.971954, 0.975809, 0.945857, 0.947371, 0.930662, 0.946869, 0.950247, 0.939276, 0.925114, 0.944206, 0.934665, 0.921425, 0.938742, 0.931046, 0.924208, 0.927826, 0.906987, 0.89058, 0.884395, 0.871027, 0.878559, 0.862506, 0.8569, 0.829157, 0.820562, 0.833909, 0.815283, 0.806924, 0.823096, 0.811991, 0.808052, 0.79241, 0.812052, 0.817675, 0.814948, 0.829207, 0.839157, 0.846599, 0.829502, 0.787898};
  114. double calibration_b1[128] = {0.688689, 0.775108, 0.887893, 0.993394, 0.898088, 0.878858, 1.20321, 0.945858, 0.898906, 0.872124, 0.878262, 0.860341, 0.860317, 0.895689, 0.906266, 0.879183, 0.871434, 0.868314, 1.07037, 0.9266, 0.90419, 0.918598, 0.911545, 0.940526, 0.947977, 0.967936, 0.936727, 0.932638, 0.960146, 0.93711, 0.936098, 0.943812, 0.9285, 0.921318, 0.9103, 0.921049, 0.914352, 0.91389, 0.92357, 0.927336, 0.929008, 0.924783, 0.93331, 0.967505, 0.958619, 1.06722, 1.14702, 1.14996, 1.09687, 0.987663, 1.01112, 1.02958, 1.00921, 1.08554, 1.09195, 1.01458, 1.02371, 1.01004, 1.00361, 1.0078, 1.0087, 1.00997, 1.03622, 1.04314, 0.894689, 0.842492, 0.912018, 0.926507, 0.936238, 0.939629, 0.959568, 0.982508, 1.31579, 0.966401, 0.977061, 0.971446, 0.989141, 1.03233, 1.05314, 1.06946, 1.04801, 1.09925, 1.54443, 1.37196, 1.39673, 1.31932, 1.11508, 1.00852, 1.01871, 1.153605, 1.06298, 0.997849, 1.00636, 1.03707, 1.03591, 1.02363, 1.01388, 1.06183, 1.10148, 1.21118, 1.26586, 1.07663, 1.07481, 1.07749, 1.0698, 1.04575, 1.06815, 1.06138, 1.06487, 1.04332, 1.01416, 1.0149, 0.989772, 0.989825, 0.985936, 0.985193, 0.959173, 0.943177, 0.993867, 0.936754, 0.934994, 0.948134, 0.922357, 0.911016, 0.805087, 0.737094, 1.03925, 0.921076};
  115. double calibration_b2[128] = {0.792171, 0.859978, 0.879552, 0.89414, 0.898794, 0.921472, 0.867588, 0.880401, 0.902657, 0.89956, 0.913535, 0.903335, 0.904012, 0.87603, 0.892562, 0.920291, 0.93814, 0.912946, 0.917767, 0.936739, 0.953897, 0.953534, 0.982874, 0.974959, 0.987859, 1.01291, 1.06464, 1.02794, 1.04071, 1.00659, 1.00576, 0.982194, 0.972739, 1.00489, 1.01822, 1.00853, 1.0207, 1.03832, 1.04594, 1.02243, 1.03193, 1.03358, 1.03278, 1.04483, 1.04611, 1.09238, 1.06066, 1.06584, 1.06113, 1.06598, 1.07107, 1.07641, 1.0793, 1.07425, 1.06559, 1.06963, 1.07979, 1.0577, 1.06135, 1.08137, 1.09443, 1.06233, 1.07256, 1.05033, 1.03563, 1.04431, 1.04057, 1.05344, 1.05814, 1.03095, 1.04926, 1.05164, 1.03729, 1.02816, 1.04004, 1.02974, 1.04938, 1.04055, 1.0446, 1.01919, 1.03717, 0.987029, 1.00179, 1.01179, 0.997439, 0.991328, 0.997722, 0.987031, 0.992538, 0.983989, 0.984387, 0.992612, 1.01848, 0.995258, 1.003, 0.993957, 0.992327, 0.990547, 1.00074, 0.996346, 1.01041, 1.01222, 1.00646, 1.03651, 1.00847, 1.01119, 0.999438, 0.979914, 1.02522, 0.973713, 0.987926, 0.957933, 0.958882, 0.950941, 0.941317, 0.933169, 0.92743, 0.936562, 0.932922, 0.913018, 0.920427, 0.909679, 0.918784, 0.933422, 0.937952, 0.950281, 0.908778, 0.859974};
  116. double calibration_b3[128] = {0.870684, 0.722589, 0.801319, 0.899014, 0.943942, 0.951251, 0.962286, 0.943771, 0.955227, 0.950627, 0.971831, 0.974907, 0.977956, 0.986517, 0.954527, 0.96526, 0.939636, 1.11263, 0.957583, 0.977304, 0.99436, 0.995887, 0.993831, 0.989441, 1.02161, 1.01742, 1.07226, 1.05602, 1.03457, 1.04393, 1.0073, 0.992433, 0.990369, 1.018, 0.983336, 1.00309, 0.996907, 0.987621, 1.00206, 0.984202, 0.990234, 0.978864, 0.982685, 0.988023, 1.00838, 0.991885, 1.01678, 0.984361, 1.01042, 0.979532, 1.00935, 0.995128, 1.00761, 1.00347, 0.999018, 0.992036, 0.994325, 0.999203, 0.991379, 0.966761, 0.923005, 0.98845, 0.944313, 0.917033, 0.966649, 0.934711, 0.934409, 0.97089, 1.00189, 1.01703, 1.01652, 1.02944, 1.03673, 1.01944, 1.04403, 1.02387, 1.04477, 1.0339, 1.03241, 1.02265, 0.987512, 1.0377, 1.01763, 1.01628, 1.02939, 1.00738, 1.02215, 1.01216, 1.00725, 1.01186, 1.02182, 1.0041, 1.01231, 1.03678, 1.02598, 1.10276, 1.03272, 1.04435, 1.03799, 1.03239, 1.05349, 1.06976, 1.0274, 1.04304, 1.03917, 1.04833, 1.04097, 1.04549, 1.06443, 1.02653, 1.03372, 1.01791, 1.02118, 1.01925, 1.01674, 1.02208, 0.990651, 1.0075, 1.00273, 1.00169, 0.999678, 1.00178, 0.981252, 0.996679, 0.960323, 0.841463, 0.724784, 0.873291};
  117. /// compile with:
  118. //// $ make clean; make
  119. /// run with:
  120. //// $ ./convert <path to data> <run>
  121. //// $ ./convert /work/leverington/beamprofilemonitor/hitdata/HIT_17_12_2017/ run2
  122. int main(int argc, char **argv)
  123. {
  124. //initialise some values;
  125. for (int i = 0;i<128;i++){
  126. board_b0_ch_bkg[i] = 0.;
  127. board_b1_ch_bkg[i] = 0.;
  128. board_b2_ch_bkg[i] = 0.;
  129. board_b3_ch_bkg[i] = 0.;
  130. channellist_gsl_b0[i] = 0;
  131. channellist_gsl_b1[i] = 0;
  132. channellist_gsl_b2[i] = 0;
  133. channellist_gsl_b3[i] = 0;
  134. sumvector_b0[i] = 0;
  135. sumvector_b1[i] = 0;
  136. sumvector_b2[i] = 0;
  137. sumvector_b3[i] = 0;
  138. if (i<64) {pos[i] = double(i);}
  139. else {pos[i] = double(i) + 0.2; }
  140. }
  141. if (argc > 1){
  142. //open the file names specified in the command line
  143. string filename;
  144. string rootfilename;
  145. // argv[1]= "/work/leverington/beamprofilemonitor/hitdata/HIT_17_12_2017/";
  146. // argv[2]= "Run11"
  147. filename = Form("%s%s.dat",argv[1],argv[2]);
  148. rootfilename = Form("%s/root/%s.root",argv[1],argv[2]); // Giulia needs to change the path %s/root/ of the written root file.
  149. ifstream file(filename, ifstream::in | ifstream::binary);
  150. if (!file.is_open()){
  151. printf(".dat did not open.\n");
  152. return -1; //file could not be opened
  153. }
  154. ///some variables to get the data
  155. // int number_of_records_to_read = 4*10;
  156. // BufferData* buffer = new BufferData[number_of_records_to_read];
  157. if(argc > 3) int block = atoi(argv[3]);
  158. int board_b = -1;
  159. long int fileframesize = getFileSize(filename.c_str()) / ( 4*sizeof(BufferData) );
  160. BufferData* dataptr = new BufferData();
  161. if (fileframesize>0){
  162. //open ROOT file for saving histos and TTree.
  163. TFile *rootFile = new TFile(rootfilename.c_str(),"RECREATE");
  164. if ( rootFile->IsOpen() ) printf("ROOT file opened successfully\n");
  165. TH2D * th2_signal_vs_channel_b0 = new TH2D("th2_signal_vs_channel_b0","th2_signal_vs_channel_b0",128,0,128,2200,-1000,10000);
  166. TH2D * th2_signal_vs_channel_b1 = new TH2D("th2_signal_vs_channel_b1","th2_signal_vs_channel_b1",128,0,128,2200,-1000,10000);
  167. TH2D * th2_signal_vs_channel_b2 = new TH2D("th2_signal_vs_channel_b2","th2_signal_vs_channel_b2",128,0,128,2200,-1000,10000);
  168. TH2D * th2_signal_vs_channel_b3 = new TH2D("th2_signal_vs_channel_b3","th2_signal_vs_channel_b3",128,0,128,2200,-1000,10000);
  169. TH2D * 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);
  170. TH2D * 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);
  171. TH2D * 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);
  172. TH2D * 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);
  173. TH1D * th1_signal_b0 = new TH1D("th1_signal_b0","th1_signal_b0",2200,-10000,50000);
  174. TH1D * th1_signal_b1 = new TH1D("th1_signal_b1","th1_signal_b1",2200,-10000,50000);
  175. TH1D * th1_signal_b2 = new TH1D("th1_signal_b2","th1_signal_b2",2200,-10000,50000);
  176. TH1D * th1_signal_b3 = new TH1D("th1_signal_b3","th1_signal_b3",2200,-10000,50000);
  177. TGraph * graph_bkg_b0 = new TGraph();
  178. TGraph * graph_bkg_b1 = new TGraph();
  179. TGraph * graph_bkg_b2 = new TGraph();
  180. TGraph * graph_bkg_b3 = new TGraph();
  181. TGraph * gr_b0;
  182. TGraph * gr_sm_b0;
  183. TGraph * gr_b1;
  184. TGraph * gr_sm_b1;
  185. TGraph * gr_b2;
  186. TGraph * gr_sm_b2;
  187. TGraph * gr_b3;
  188. TGraph * gr_sm_b3;
  189. TH2D * 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);
  190. TH2D * 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);
  191. TH2D * 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);
  192. TH2D * 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);
  193. TTree *rootTree = new TTree("t","HIT Data Root Tree");
  194. rootTree->Branch("beamPosX_b0",&beamPosX_b0,"beamPosX_b0/D");
  195. rootTree->Branch("beamPosX_b1",&beamPosX_b1,"beamPosX_b1/D");
  196. rootTree->Branch("beamPosX_b2",&beamPosX_b2,"beamPosX_b2/D");
  197. rootTree->Branch("beamPosX_b3",&beamPosX_b3,"beamPosX_b3/D");
  198. rootTree->Branch("beamFocusX_b0",&beamFocusX_b0,"beamFocusX_b0/D");
  199. rootTree->Branch("beamFocusX_b1",&beamFocusX_b1,"beamFocusX_b1/D");
  200. rootTree->Branch("beamFocusX_b2",&beamFocusX_b2,"beamFocusX_b2/D");
  201. rootTree->Branch("beamFocusX_b3",&beamFocusX_b3,"beamFocusX_b3/D");
  202. rootTree->Branch("beamSkewX_b0",&beamSkewX_b0,"beamSkewX_b0/D");
  203. rootTree->Branch("beamSkewX_b1",&beamSkewX_b1,"beamSkewX_b1/D");
  204. rootTree->Branch("beamSkewX_b2",&beamSkewX_b2,"beamSkewX_b2/D");
  205. rootTree->Branch("beamSkewX_b3",&beamSkewX_b3,"beamSkewX_b3/D");
  206. rootTree->Branch("beamKurtX_b0",&beamKurtX_b0,"beamKurtX_b0/D");
  207. rootTree->Branch("beamKurtX_b1",&beamKurtX_b1,"beamKurtX_b1/D");
  208. rootTree->Branch("beamKurtX_b2",&beamKurtX_b2,"beamKurtX_b2/D");
  209. rootTree->Branch("beamKurtX_b3",&beamKurtX_b3,"beamKurtX_b3/D");
  210. rootTree->Branch("beamNumX_b0",&beamNumX_b0,"beamNumX_b0/D");
  211. rootTree->Branch("beamNumX_b1",&beamNumX_b1,"beamNumX_b1/D");
  212. rootTree->Branch("beamNumX_b2",&beamNumX_b2,"beamNumX_b2/D");
  213. rootTree->Branch("beamNumX_b3",&beamNumX_b3,"beamNumX_b3/D");
  214. rootTree->Branch("beamSignal_b0",&signal_b0,"beamSignal_b0/D");
  215. rootTree->Branch("beamSignal_b1",&signal_b1,"beamSignal_b1/D");
  216. rootTree->Branch("beamSignal_b2",&signal_b2,"beamSignal_b2/D");
  217. rootTree->Branch("beamSignal_b3",&signal_b3,"beamSignal_b3/D");
  218. rootTree->Branch("beamSignal_b0_left",&signal_b0_left,"beamSignal_b0_left/D");
  219. rootTree->Branch("beamSignal_b1_left",&signal_b1_left,"beamSignal_b1_left/D");
  220. rootTree->Branch("beamSignal_b2_left",&signal_b2_left,"beamSignal_b2_left/D");
  221. rootTree->Branch("beamSignal_b3_left",&signal_b3_left,"beamSignal_b3_left/D");
  222. rootTree->Branch("beamSignal_b0_right",&signal_b0_right,"beamSignal_b0_right/D");
  223. rootTree->Branch("beamSignal_b1_right",&signal_b1_right,"beamSignal_b1_right/D");
  224. rootTree->Branch("beamSignal_b2_right",&signal_b2_right,"beamSignal_b2_right/D");
  225. rootTree->Branch("beamSignal_b3_right",&signal_b3_right,"beamSignal_b3_right/D");
  226. rootTree->Branch("eventID",&eventID,"eventID/I");
  227. rootTree->Branch("board_b0_ch",&board_b0_ch,"board_b0_ch[128]/D");
  228. rootTree->Branch("board_b1_ch",&board_b1_ch,"board_b1_ch[128]/D");
  229. rootTree->Branch("board_b2_ch",&board_b2_ch,"board_b2_ch[128]/D");
  230. rootTree->Branch("board_b3_ch",&board_b3_ch,"board_b3_ch[128]/D");
  231. rootTree->Branch("beamSidebandNoise_b0",&beamSidebandNoise_b0,"beamSidebandNoise_b0/D");
  232. rootTree->Branch("beamSidebandNoise_b1",&beamSidebandNoise_b1,"beamSidebandNoise_b1/D");
  233. rootTree->Branch("beamSidebandNoise_b2",&beamSidebandNoise_b2,"beamSidebandNoise_b2/D");
  234. rootTree->Branch("beamSidebandNoise_b3",&beamSidebandNoise_b3,"beamSidebandNoise_b3/D");
  235. rootTree->Branch("arrayavg_b0",&arrayavg_b0,"arrayavg_b0[128]/D");
  236. rootTree->Branch("arrayavg_b1",&arrayavg_b1,"arrayavg_b1[128]/D");
  237. rootTree->Branch("arrayavg_b2",&arrayavg_b2,"arrayavg_b2[128]/D");
  238. rootTree->Branch("arrayavg_b3",&arrayavg_b3,"arrayavg_b3[128]/D");
  239. //loop through recorded data frames
  240. for (int frame = 0; frame<fileframesize; frame++){
  241. eventID = frame;
  242. // cout << eventID << endl;
  243. //board_b 0
  244. board_b=0;
  245. signal_b0 = 0.;
  246. maxchannelamp_b0 = 0.;
  247. file.seekg(board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
  248. file.read ((char*)dataptr ,sizeof(BufferData));
  249. if (dataptr->sync_frame.device_nr==board_b){
  250. // cout << frame << " " << dataptr->sync_frame.device_nr << endl;
  251. if (frame%1000==0) cout << "Frame: " << frame << endl;
  252. for (int i = 0;i<128;i++){
  253. board_b0_ch[i] = dataptr->sensor_data[i];
  254. if (frame<1000){
  255. board_b0_ch_bkg[i] += board_b0_ch[i]/1000.; //find baseline from the average of first 1000 events
  256. }
  257. else if (frame>=1000&& frame<2000){
  258. board_b0_ch[i] -= board_b0_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
  259. th2_signal_vs_channel_bkg_b0->Fill(i,board_b0_ch[i]);
  260. if (i==1) graph_bkg_b0->SetPoint(graph_bkg_b0->GetN(), frame, board_b0_ch[i]);
  261. }
  262. else if (frame>=2000) {
  263. board_b0_ch[i]-=board_b0_ch_bkg[i]; // the rest background subtracted
  264. board_b0_ch[i]*=calibration_b0[i]; //calibration factor
  265. th2_signal_vs_channel_b0->Fill(i,board_b0_ch[i]);
  266. // signal_b0 +=board_b0_ch[i] ;
  267. if (board_b0_ch[i]> maxchannelamp_b0) {
  268. maxchannel_b0 = i;
  269. maxchannelamp_b0 = board_b0_ch[i];
  270. }
  271. //calculate a rolling average of the signal
  272. arrayavg_b0[i] = 0;
  273. for (int j = 1; j<length;j++){
  274. array_b0[j-1][i] = array_b0[j][i];
  275. arrayavg_b0[i] += array_b0[j-1][i]/double(length);
  276. }
  277. if( board_b0_ch[i]>-1000){
  278. array_b0[length-1][i] = board_b0_ch[i];
  279. arrayavg_b0[i] += array_b0[length-1][i]/double(length);
  280. }
  281. ////////////////////////////////////////////
  282. }
  283. }
  284. // th1_signal_b0->Fill(signal_b0);
  285. }
  286. else {
  287. cout << "Error." << endl;
  288. }
  289. //board_b 1
  290. board_b=1;
  291. signal_b1=0.;
  292. maxchannelamp_b1 = 0.;
  293. file.seekg(board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
  294. file.read ((char*)dataptr ,sizeof(BufferData));
  295. if (dataptr->sync_frame.device_nr==board_b){
  296. // cout << frame << " " << dataptr->sync_frame.device_nr << endl;
  297. for (int i = 0;i<128;i++){
  298. board_b1_ch[i] = dataptr->sensor_data[i];
  299. if (frame<1000){
  300. board_b1_ch_bkg[i] += board_b1_ch[i]/1000.; //find baseline from the average of first 1000 events
  301. }
  302. else if (frame>=1000&& frame<2000){
  303. board_b1_ch[i] -= board_b1_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
  304. th2_signal_vs_channel_bkg_b1->Fill(i,board_b1_ch[i]);
  305. if (i==1) graph_bkg_b1->SetPoint(graph_bkg_b1->GetN(), frame, board_b1_ch[i]);
  306. }
  307. else if (frame>=2000) {
  308. board_b1_ch[i]-=board_b1_ch_bkg[i]; // the rest are background subtracted
  309. board_b1_ch[i]*=calibration_b1[i]; //calibration factor
  310. th2_signal_vs_channel_b1->Fill(i,board_b1_ch[i]);
  311. // signal_b1 +=board_b1_ch[i] ;
  312. if (board_b1_ch[i]> maxchannelamp_b1) {
  313. maxchannel_b1 = i;
  314. maxchannelamp_b1 = board_b1_ch[i];
  315. }
  316. //calculate a rolling average of the signal
  317. arrayavg_b1[i] = 0;
  318. for (int j = 1; j<length;j++){
  319. array_b1[j-1][i] = array_b1[j][i];
  320. arrayavg_b1[i] += array_b1[j-1][i]/double(length);
  321. }
  322. if( board_b1_ch[i]>-1000){
  323. array_b1[length-1][i] = board_b1_ch[i];
  324. arrayavg_b1[i] += array_b1[length-1][i]/double(length);
  325. }
  326. ////////////////////////////////////////////
  327. }
  328. }
  329. // th1_signal_b1->Fill(signal_b1);
  330. }
  331. else {
  332. cout << "Error." << endl;
  333. }
  334. //board_b 2
  335. board_b=2;
  336. signal_b2=0.;
  337. maxchannelamp_b2 = 0.;
  338. file.seekg(board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
  339. file.read ((char*)dataptr ,sizeof(BufferData));
  340. if (dataptr->sync_frame.device_nr==board_b){
  341. // cout << frame << " " << dataptr->sync_frame.device_nr << endl;
  342. for (int i = 0;i<128;i++){
  343. board_b2_ch[i] = dataptr->sensor_data[i];
  344. if (frame<1000){
  345. board_b2_ch_bkg[i] += board_b2_ch[i]/1000.; //find baseline from the average of first 1000 events
  346. }
  347. else if (frame>=1000&& frame<2000){
  348. board_b2_ch[i] -= board_b2_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
  349. th2_signal_vs_channel_bkg_b2->Fill(i,board_b2_ch[i]);
  350. if (i==1) graph_bkg_b2->SetPoint(graph_bkg_b2->GetN(), frame, board_b2_ch[i]);
  351. }
  352. else if (frame>=2000) {
  353. board_b2_ch[i]-=board_b2_ch_bkg[i]; // the rest background subtracted
  354. board_b2_ch[i]*=calibration_b2[i]; //calibration factor
  355. th2_signal_vs_channel_b2->Fill(i,board_b2_ch[i]);
  356. // signal_b2 +=board_b2_ch[i] ;
  357. if (board_b2_ch[i]> maxchannelamp_b2) {
  358. maxchannel_b2 = i;
  359. maxchannelamp_b2 = board_b2_ch[i];
  360. }
  361. //calculate a rolling average of the signal
  362. arrayavg_b2[i] = 0;
  363. for (int j = 1; j<length;j++){
  364. array_b2[j-1][i] = array_b2[j][i];
  365. arrayavg_b2[i] += array_b2[j-1][i]/double(length);
  366. }
  367. if( board_b2_ch[i]>-1000){
  368. array_b2[length-1][i] = board_b2_ch[i];
  369. arrayavg_b2[i] += array_b2[length-1][i]/double(length);
  370. }
  371. ////////////////////////////////////////////
  372. }
  373. }
  374. // th1_signal_b2->Fill(signal_b2);
  375. }
  376. else {
  377. cout << "Error." << endl;
  378. }
  379. //board_b 3
  380. board_b=3;
  381. signal_b3=0.;
  382. maxchannelamp_b3 = 0.;
  383. file.seekg(board_b*sizeof(BufferData)+4*frame*sizeof(BufferData));
  384. file.read ((char*)dataptr ,sizeof(BufferData));
  385. if (dataptr->sync_frame.device_nr==board_b){
  386. // cout << frame << " " << dataptr->sync_frame.device_nr << endl;
  387. for (int i = 0;i<128;i++){
  388. board_b3_ch[i] = dataptr->sensor_data[i];
  389. if (frame<1000){
  390. board_b3_ch_bkg[i] += board_b3_ch[i]/1000.; //find baseline from the average of first 1000 events
  391. }
  392. else if (frame>=1000&& frame<2000){
  393. board_b3_ch[i] -= board_b3_ch_bkg[i]; //histogram the subtracted baseline in the next 1000 events
  394. th2_signal_vs_channel_bkg_b3->Fill(i,board_b3_ch[i]);
  395. if (i==1) graph_bkg_b3->SetPoint(graph_bkg_b3->GetN(), frame, board_b3_ch[i]);
  396. }
  397. else if (frame>=2000) {
  398. board_b3_ch[i]-=board_b3_ch_bkg[i]; // the rest of the events are background subtracted
  399. board_b3_ch[i]*=calibration_b3[i]; //with a calibration factor
  400. th2_signal_vs_channel_b3->Fill(i,board_b3_ch[i]);
  401. // signal_b3 +=board_b3_ch[i] ;
  402. if (board_b3_ch[i]> maxchannelamp_b3) {
  403. maxchannel_b3 = i;
  404. maxchannelamp_b3 = board_b3_ch[i];
  405. }
  406. //calculate a rolling average of the signal
  407. arrayavg_b3[i] = 0;
  408. for (int j = 1; j<length;j++){
  409. array_b3[j-1][i] = array_b3[j][i];
  410. arrayavg_b3[i] += array_b3[j-1][i]/double(length);
  411. }
  412. if( board_b3_ch[i]>-1000){
  413. array_b3[length-1][i] = board_b3_ch[i];
  414. arrayavg_b3[i] += array_b3[length-1][i]/double(length);
  415. }
  416. ////////////////////////////////////////////
  417. }
  418. }
  419. // th1_signal_b3->Fill(signal_b3);
  420. }
  421. else {
  422. cout << "Error." << endl;
  423. }
  424. ////////////////////////////////////////////////////////////
  425. ////////////////////////////////////////////////////////////
  426. //start the signal analysis
  427. ////////////////////////////////////////////////////////////
  428. ////////////////////////////////////////////////////////////
  429. if (frame>=2000){
  430. ////////////////////////////////////////////////////////////
  431. //boxcar smoothing filter
  432. ////////////////////////////////////////////////////////////
  433. if (smooth_on){
  434. maxchannelamp_b0 = 0.0;
  435. maxchannelamp_b1 = 0.0;
  436. maxchannelamp_b2 = 0.0;
  437. maxchannelamp_b3 = 0.0;
  438. }
  439. boxcar.clear();
  440. boxcar.push_back(board_b0_ch[0]);
  441. boxcar.push_back(board_b0_ch[1]);
  442. // cout << frame << endl;
  443. for (int i = 0;i<128;i++){
  444. channellist[i] = i;
  445. if (i<126) {
  446. boxcar.push_back(board_b0_ch[i+2]);
  447. }
  448. if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
  449. board_b0_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
  450. if (smooth_on && board_b0_smooth_ch[i]> maxchannelamp_b0) {
  451. maxchannel_b0 = i;
  452. maxchannelamp_b0 = board_b0_smooth_ch[i];
  453. }
  454. }
  455. if (eventID>2001&&maxchannelamp_b0>400&&!graphsaved_b0) {
  456. gr_b0 = new TGraph(128, channellist, board_b0_ch);
  457. gr_b0->SetTitle("RawData_b0");
  458. gr_b0->SetName("RawData_b0");
  459. gr_b0->Write();
  460. gr_sm_b0 = new TGraph(128, channellist, board_b0_smooth_ch);
  461. gr_sm_b0->SetTitle("SmoothedData_b0");
  462. gr_sm_b0->SetName("SmoothedData_b0");
  463. gr_sm_b0->Write();
  464. graphsaved_b0=true;
  465. }
  466. boxcar.clear();
  467. boxcar.push_back(board_b1_ch[0]);
  468. boxcar.push_back(board_b1_ch[1]);
  469. // cout << frame << endl;
  470. for (int i = 0;i<128;i++){
  471. channellist[i] = i;
  472. if (i<126) {
  473. boxcar.push_back(board_b1_ch[i+2]);
  474. }
  475. if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
  476. board_b1_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
  477. if (smooth_on && board_b1_smooth_ch[i]> maxchannelamp_b1) {
  478. maxchannel_b1 = i;
  479. maxchannelamp_b1 = board_b1_smooth_ch[i];
  480. }
  481. }
  482. if (eventID>2001&&maxchannelamp_b1>400&&!graphsaved_b1) {
  483. gr_b1 = new TGraph(128, channellist, board_b1_ch);
  484. gr_b1->SetTitle("RawData_b1");
  485. gr_b1->SetName("RawData_b1");
  486. gr_b1->Write();
  487. gr_sm_b1 = new TGraph(128, channellist, board_b1_smooth_ch);
  488. gr_sm_b1->SetTitle("SmoothedData_b1");
  489. gr_sm_b1->SetName("SmoothedData_b1");
  490. gr_sm_b1->Write();
  491. graphsaved_b1=true;
  492. }
  493. boxcar.clear();
  494. boxcar.push_back(board_b2_ch[0]);
  495. boxcar.push_back(board_b2_ch[1]);
  496. // cout << frame << endl;
  497. for (int i = 0;i<128;i++){
  498. channellist[i] = i;
  499. if (i<126) {
  500. boxcar.push_back(board_b2_ch[i+2]);
  501. }
  502. if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
  503. board_b2_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
  504. if (smooth_on && board_b2_smooth_ch[i]> maxchannelamp_b2) {
  505. maxchannel_b2 = i;
  506. maxchannelamp_b2 = board_b2_smooth_ch[i];
  507. }
  508. }
  509. if (eventID>2001&&maxchannelamp_b2>400&&!graphsaved_b2) {
  510. gr_b2 = new TGraph(128, channellist, board_b2_ch);
  511. gr_b2->SetTitle("RawData_b2");
  512. gr_b2->SetName("RawData_b2");
  513. gr_b2->Write();
  514. gr_sm_b2 = new TGraph(128, channellist, board_b2_smooth_ch);
  515. gr_sm_b2->SetTitle("SmoothedData_b2");
  516. gr_sm_b2->SetName("SmoothedData_b2");
  517. gr_sm_b2->Write();
  518. graphsaved_b2=true;
  519. }
  520. boxcar.clear();
  521. boxcar.push_back(board_b3_ch[0]);
  522. boxcar.push_back(board_b3_ch[1]);
  523. // cout << frame << endl;
  524. for (int i = 0;i<128;i++){
  525. channellist[i] = i;
  526. if (i<126) {
  527. boxcar.push_back(board_b3_ch[i+2]);
  528. }
  529. if ( (i<126&&boxcar.size()>5) || (i>=126&&boxcar.size()>3) ) boxcar.erase(boxcar.cbegin());
  530. board_b3_smooth_ch[i] = gsl_stats_mean(boxcar.data(), 1, 5 );
  531. if (smooth_on && board_b3_smooth_ch[i]> maxchannelamp_b3) {
  532. maxchannel_b3 = i;
  533. maxchannelamp_b3 = board_b3_smooth_ch[i];
  534. }
  535. }
  536. if (eventID>2001&&maxchannelamp_b3>400&&!graphsaved_b3) {
  537. gr_b3 = new TGraph(128, channellist, board_b3_ch);
  538. gr_b3->SetTitle("RawData_b3");
  539. gr_b3->SetName("RawData_b3");
  540. gr_b3->Write();
  541. gr_sm_b3 = new TGraph(128, channellist, board_b3_smooth_ch);
  542. gr_sm_b3->SetTitle("SmoothedData_b3");
  543. gr_sm_b3->SetName("SmoothedData_b3");
  544. gr_sm_b3->Write();
  545. graphsaved_b3=true;
  546. }
  547. ////////////////////////////////////////////////////////////
  548. //find the approx FWHM
  549. ////////////////////////////////////////////////////////////
  550. nfwhm_b0 =0;
  551. nfwhm_b1 =0;
  552. nfwhm_b2 =0;
  553. nfwhm_b3 =0;
  554. if (smooth_on){
  555. for (int i = 0;i<128;i++){
  556. if (board_b0_smooth_ch[i] > maxchannelamp_b0/2.) nfwhm_b0++;
  557. if (board_b1_smooth_ch[i] > maxchannelamp_b1/2.) nfwhm_b1++;
  558. if (board_b2_smooth_ch[i] > maxchannelamp_b2/2.) nfwhm_b2++;
  559. if (board_b3_smooth_ch[i] > maxchannelamp_b3/2.) nfwhm_b3++;
  560. signal_gsl_b0[i] = 0.;
  561. signal_gsl_b1[i] = 0.;
  562. signal_gsl_b2[i] = 0.;
  563. signal_gsl_b3[i] = 0.;
  564. }
  565. }
  566. else {
  567. for (int i = 0;i<128;i++){
  568. if (board_b0_ch[i] > maxchannelamp_b0/2.) nfwhm_b0++;
  569. if (board_b1_ch[i] > maxchannelamp_b1/2.) nfwhm_b1++;
  570. if (board_b2_ch[i] > maxchannelamp_b2/2.) nfwhm_b2++;
  571. if (board_b3_ch[i] > maxchannelamp_b3/2.) nfwhm_b3++;
  572. signal_gsl_b0[i] = 0.;
  573. signal_gsl_b1[i] = 0.;
  574. signal_gsl_b2[i] = 0.;
  575. signal_gsl_b3[i] = 0.;
  576. }
  577. }
  578. ////////////////////////////////////////////////////////////
  579. //integrate the sidebands first to check for baseline shift
  580. ////////////////////////////////////////////////////////////
  581. /*
  582. beamSidebandNoise_b0 = 0.;
  583. sidenumtocalc_b0 = 0;
  584. for (int i =0; i<maxchannel_b0-nfwhm_b0; i++){
  585. if (i>=0 && i<=127){
  586. beamSidebandNoise_b0 += board_b0_ch[i]; //integrate the noise outside the peak
  587. sidenumtocalc_b0++;
  588. }
  589. }
  590. for (int i = maxchannel_b0+nfwhm_b0; i < 128; i++){
  591. if (i>=0 && i<=127){
  592. beamSidebandNoise_b0 += board_b0_ch[i]; //integrate the noise outside the peak
  593. sidenumtocalc_b0++;
  594. }
  595. }
  596. if (sidenumtocalc_b0>0) beamSidebandNoise_b0 = beamSidebandNoise_b0 /double(sidenumtocalc_b0); // channel baseline shift
  597. beamSidebandNoise_b1 = 0.;
  598. sidenumtocalc_b1 = 0;
  599. for (int i =0; i<maxchannel_b1-nfwhm_b1; i++){
  600. if (i>=0 && i<=127){
  601. beamSidebandNoise_b1 += board_b1_ch[i]; //integrate the noise outside the peak
  602. sidenumtocalc_b1++;
  603. }
  604. }
  605. for (int i = maxchannel_b1+nfwhm_b1; i < 128; i++){
  606. if (i>=0 && i<=127){
  607. beamSidebandNoise_b1 += board_b1_ch[i]; //integrate the noise outside the peak
  608. sidenumtocalc_b1++;
  609. }
  610. }
  611. if (sidenumtocalc_b1>0) beamSidebandNoise_b1 = beamSidebandNoise_b1 /double(sidenumtocalc_b1); // channel baseline shift
  612. beamSidebandNoise_b2 = 0.;
  613. sidenumtocalc_b2 = 0;
  614. for (int i =0; i<maxchannel_b2-nfwhm_b2; i++){
  615. if (i>=0 && i<=127){
  616. beamSidebandNoise_b2 += board_b2_ch[i]; //integrate the noise outside the peak
  617. sidenumtocalc_b2++;
  618. }
  619. }
  620. for (int i = maxchannel_b2+nfwhm_b2; i < 128; i++){
  621. if (i>=0 && i<=127){
  622. beamSidebandNoise_b2 += board_b2_ch[i]; //integrate the noise outside the peak
  623. sidenumtocalc_b2++;
  624. }
  625. }
  626. if (sidenumtocalc_b2>0) beamSidebandNoise_b2 = beamSidebandNoise_b2 /double(sidenumtocalc_b2); // channel baseline shift
  627. beamSidebandNoise_b3 = 0.;
  628. sidenumtocalc_b3 = 0;
  629. for (int i =0; i<maxchannel_b3-nfwhm_b3; i++){
  630. if (i>=0 && i<=127){
  631. beamSidebandNoise_b3 += board_b3_ch[i]; //integrate the noise outside the peak
  632. sidenumtocalc_b3++;
  633. }
  634. }
  635. for (int i = maxchannel_b3+nfwhm_b3; i < 128; i++){
  636. if (i>=0 && i<=127){
  637. beamSidebandNoise_b3 += board_b3_ch[i]; //integrate the noise outside the peak
  638. sidenumtocalc_b3++;
  639. }
  640. }
  641. if (sidenumtocalc_b3>0) beamSidebandNoise_b3 = beamSidebandNoise_b3 /double(sidenumtocalc_b3); // channel baseline shift
  642. */
  643. ////////////////////////////////////////////////////////////
  644. //integrate under the approximate peak
  645. //build the channel list for statistical analysis
  646. ////////////////////////////////////////////////////////////
  647. numtocalc_b0 = 0;
  648. numtocalc_b1 = 0;
  649. numtocalc_b2 = 0;
  650. numtocalc_b3 = 0;
  651. signal_b0_left = 0.;
  652. signal_b0_right = 0.;
  653. for (int i = maxchannel_b0-nfwhm_b0 ; i <= maxchannel_b0 + nfwhm_b0; i++){
  654. if (i>=0 && i<=127 && board_b0_ch[i]>50){
  655. signal_b0 +=board_b0_ch[i]-beamSidebandNoise_b0 ;
  656. if (i<64) {signal_b0_left +=board_b0_ch[i]-beamSidebandNoise_b0 ;}
  657. else {signal_b0_right +=board_b0_ch[i]-beamSidebandNoise_b0 ;}
  658. signal_gsl_b0[numtocalc_b0]=board_b0_ch[i]-beamSidebandNoise_b0 ;
  659. // channellist_gsl_b0[numtocalc_b0] = i;
  660. channellist_gsl_b0[numtocalc_b0] = pos[i];
  661. numtocalc_b0++;
  662. }
  663. }
  664. th1_signal_b0->Fill(signal_b0);
  665. signal_b1_left = 0.;
  666. signal_b1_right = 0.;
  667. for (int i = maxchannel_b1-nfwhm_b1 ; i <= maxchannel_b1 + nfwhm_b1; i++){
  668. if (i>=0 && i<=127&&board_b1_ch[i]>50){
  669. signal_b1 +=board_b1_ch[i]-beamSidebandNoise_b1 ;
  670. if (i<64) {signal_b1_left +=board_b1_ch[i]-beamSidebandNoise_b1 ;}
  671. else {signal_b1_right +=board_b1_ch[i]-beamSidebandNoise_b1 ;}
  672. signal_gsl_b1[numtocalc_b1]=board_b1_ch[i]-beamSidebandNoise_b1 ;
  673. // channellist_gsl_b1[numtocalc_b1] = i;
  674. channellist_gsl_b1[numtocalc_b1] = pos[i];
  675. numtocalc_b1++;
  676. }
  677. }
  678. th1_signal_b0->Fill(signal_b1);
  679. signal_b2_left = 0.;
  680. signal_b2_right = 0.;
  681. for (int i = maxchannel_b2-nfwhm_b2 ; i <= maxchannel_b2 + nfwhm_b2; i++){
  682. if (i>=0 && i<=127&&board_b2_ch[i]>50){
  683. signal_b2 +=board_b2_ch[i]-beamSidebandNoise_b2 ;
  684. if (i<64) {signal_b2_left +=board_b2_ch[i]-beamSidebandNoise_b2 ;}
  685. else {signal_b2_right +=board_b2_ch[i]-beamSidebandNoise_b2 ;}
  686. signal_gsl_b2[numtocalc_b2]=board_b2_ch[i]-beamSidebandNoise_b2 ;
  687. // channellist_gsl_b2[numtocalc_b2] = i;
  688. channellist_gsl_b2[numtocalc_b2] = pos[i];
  689. numtocalc_b2++;
  690. }
  691. }
  692. th1_signal_b0->Fill(signal_b2);
  693. signal_b3_left = 0.;
  694. signal_b3_right = 0.;
  695. for (int i = maxchannel_b3-nfwhm_b3 ; i <= maxchannel_b3 + nfwhm_b3; i++){
  696. if (i>=0 && i<=127&&board_b3_ch[i]>50){
  697. signal_b3 +=board_b3_ch[i]-beamSidebandNoise_b3 ;
  698. if (i<64) {signal_b3_left +=board_b3_ch[i]-beamSidebandNoise_b3 ;}
  699. else {signal_b3_right +=board_b3_ch[i]-beamSidebandNoise_b3 ;}
  700. signal_gsl_b3[numtocalc_b3]=board_b3_ch[i]-beamSidebandNoise_b3 ;
  701. //channellist_gsl_b3[numtocalc_b3] = i;
  702. channellist_gsl_b3[numtocalc_b3] = pos[i];
  703. numtocalc_b3++;
  704. }
  705. }
  706. th1_signal_b0->Fill(signal_b3);
  707. for (int i=0;i<128;i++){
  708. if(signal_b0>7000) {
  709. th2_signal_vs_channel_sub_b0->Fill(i,board_b0_ch[i]-beamSidebandNoise_b0);
  710. sumvector_b0[i] += board_b0_ch[i];//-beamSidebandNoise_b0;
  711. }
  712. if(signal_b1>7000) {
  713. th2_signal_vs_channel_sub_b1->Fill(i,board_b1_ch[i]-beamSidebandNoise_b1);
  714. sumvector_b1[i] += board_b1_ch[i];//-beamSidebandNoise_b1;
  715. }
  716. if(signal_b2>7000){
  717. th2_signal_vs_channel_sub_b2->Fill(i,board_b2_ch[i]-beamSidebandNoise_b2);
  718. sumvector_b2[i] += board_b2_ch[i];//-beamSidebandNoise_b2;
  719. }
  720. if(signal_b3>7000) {
  721. th2_signal_vs_channel_sub_b3->Fill(i,board_b3_ch[i]-beamSidebandNoise_b3);
  722. sumvector_b3[i] += board_b3_ch[i];//-beamSidebandNoise_b3;
  723. }
  724. }
  725. if(signal_b0>7000) beamontime[0]+=1.0;
  726. if(signal_b1>7000) beamontime[1]+=1.0;
  727. if(signal_b2>7000) beamontime[2]+=1.0;
  728. if(signal_b3>7000) beamontime[3]+=1.0;
  729. ///add gsl stuff here.
  730. /* cout << maxchannel_b0 << " " << maxchannel_b1 << " "<< maxchannel_b2 << " "<< maxchannel_b3 << " " << endl;
  731. cout << maxchannelamp_b0 << " " << maxchannelamp_b1 << " "<< maxchannelamp_b2 << " "<< maxchannelamp_b3 << " " << endl;
  732. cout << nfwhm_b0 << " " << nfwhm_b1 << " " << nfwhm_b2 << " " << nfwhm_b3 << " " << endl;
  733. cout << endl;*/
  734. beamPosX_b0 = gsl_stats_wmean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0); //calculate the weighted mean
  735. beamFocusX_b0 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0); //Standard Deviation
  736. beamSkewX_b0 = gsl_stats_wskew_m_sd(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0,beamFocusX_b0); //skewness (symmetry)
  737. beamKurtX_b0 = gsl_stats_wkurtosis_m_sd(signal_gsl_b0,1,channellist_gsl_b0,1,numtocalc_b0,beamPosX_b0,beamFocusX_b0); //excess kurtosis (well behaved tails)
  738. beamNumX_b0 = numtocalc_b0;
  739. beamFocusX_b0 *=2.3548;//SD-->FWHM
  740. beamPosX_b1 = gsl_stats_wmean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1); //calculate the weighted mean
  741. beamFocusX_b1 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1); //Standard Deviation
  742. beamSkewX_b1 = gsl_stats_wskew_m_sd(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1,beamFocusX_b1); //skewness (symmetry)
  743. beamKurtX_b1 = gsl_stats_wkurtosis_m_sd(signal_gsl_b1,1,channellist_gsl_b1,1,numtocalc_b1,beamPosX_b1,beamFocusX_b1); //excess kurtosis (well behaved tails)
  744. beamNumX_b1 = numtocalc_b1;
  745. beamFocusX_b1 *=2.3548;//SD-->FWHM
  746. beamPosX_b2 = gsl_stats_wmean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2); //calculate the weighted mean
  747. beamFocusX_b2 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2); //Standard Deviation
  748. beamSkewX_b2 = gsl_stats_wskew_m_sd(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2,beamFocusX_b2); //skewness (symmetry)
  749. beamKurtX_b2 = gsl_stats_wkurtosis_m_sd(signal_gsl_b2,1,channellist_gsl_b2,1,numtocalc_b2,beamPosX_b2,beamFocusX_b2); //excess kurtosis (well behaved tails)
  750. beamNumX_b2 = numtocalc_b2;
  751. beamFocusX_b2 *=2.3548;//SD-->FWHM
  752. beamPosX_b3 = gsl_stats_wmean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3); //calculate the weighted mean
  753. beamFocusX_b3 = gsl_stats_wsd_with_fixed_mean(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3); //Standard Deviation
  754. beamSkewX_b3 = gsl_stats_wskew_m_sd(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3,beamFocusX_b3); //skewness (symmetry)
  755. beamKurtX_b3 = gsl_stats_wkurtosis_m_sd(signal_gsl_b3,1,channellist_gsl_b3,1,numtocalc_b3,beamPosX_b3,beamFocusX_b3); //excess kurtosis (well behaved tails)
  756. beamNumX_b3 = numtocalc_b3;
  757. beamFocusX_b3 *=2.3548;//SD-->FWHM
  758. /////////////////////////////////////////////////
  759. /////Add fitting algorithm here
  760. /////////////////////////////////////////////////
  761. /*
  762. use board_b0_ch[i] to fill a TGraph.
  763. fit with a gaussian
  764. subtract the difference
  765. fill a THF2 with the difference, delta_amp vs channel
  766. */
  767. rootTree->Fill();
  768. }
  769. }//end of loop over frames
  770. /* th2_signal_vs_channel_b0->Write();
  771. th2_signal_vs_channel_b1->Write();
  772. th2_signal_vs_channel_b2->Write();
  773. th2_signal_vs_channel_b3->Write();
  774. th2_signal_vs_channel_bkg_b0->Write();
  775. th2_signal_vs_channel_bkg_b1->Write();
  776. th2_signal_vs_channel_bkg_b2->Write();
  777. th2_signal_vs_channel_bkg_b3->Write();
  778. th1_signal_b0->Write();
  779. th1_signal_b1->Write();
  780. th1_signal_b2->Write();
  781. th1_signal_b3->Write();*/
  782. graph_bkg_b0->SetName("graph_bkg_b0"); // graph_bkg_b0->Write();
  783. graph_bkg_b1->SetName("graph_bkg_b1"); // graph_bkg_b1->Write();
  784. graph_bkg_b2->SetName("graph_bkg_b2"); // graph_bkg_b2->Write();
  785. graph_bkg_b3->SetName("graph_bkg_b3"); // graph_bkg_b3->Write();
  786. sumvector_b0_ptr->Write("sumvector_b0");
  787. sumvector_b1_ptr->Write("sumvector_b1");
  788. sumvector_b2_ptr->Write("sumvector_b2");
  789. sumvector_b3_ptr->Write("sumvector_b3");
  790. beamontime_ptr->Write("beamontime");
  791. rootFile->Write();
  792. rootFile->Close();
  793. }
  794. cout << eventID << " frames analysed." << endl;
  795. return 0;
  796. }
  797. else {
  798. cerr << "Error 1" << endl;
  799. return 1;
  800. }
  801. }