-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathaw_lukas_backup_2019_06_30.cxx
5398 lines (5123 loc) · 216 KB
/
aw_lukas_backup_2019_06_30.cxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/**************************************
Initial code by [email protected] and [email protected]
Written by [email protected]
modified for single board / self triggered readout
**************************************/
#include <stdlib.h>
#include <stdio.h>
#include <cmath>
#include <math.h>
#include <string.h>
#include "aw_lib.h"
#include "tagger_lib.h"
#include <vector>
#include <fstream>
#include <cstddef>
#include <chrono> // measuring high-res execution time
#include "Math/Interpolator.h"
#include "Math/Polynomial.h"
#include "TROOT.h"
#include "TMath.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TRandom.h"
#include "TTree.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TDirectory.h"
#include "TStyle.h"
#include "TList.h"
#include "TExec.h"
#include "TText.h"
#include "TGraphPainter.h"
#include "TSpectrum.h" // For spectrum and peak analysis
#include "TVirtualFitter.h" // fitting
#include "TPaveStats.h"
// RooFit Framework
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "TLatex.h"
#include "RooNovosibirsk.h"
using namespace RooFit ;
using namespace std;
using namespace chrono;
using namespace RooFit;
int doof = 0;
// global constants
int TRACELEN = 250;
int MATRIX_J = 1;
int MATRIX_K = 2;
char VERBOSE[200];
int CHANNELS = 8;
int TAG_CHANNELS = 16;
int TAGGER_WINDOW = 100;
int CENTRAL = 4;
int BOARDS = 1;
int CHANNELS_EFF;
int SAMPLE_t = 10; // Sample frequency in ns
char MODE[100];
double SAMPLE_t_eff; // Effective sample frequency
Int_t N_E_WINDOW = 12;
Int_t E_WINDOW_LEFT = 0;
Int_t E_WINDOW_RIGHT = 6000;
Int_t E_WINDOW_LENGTH = (int) (E_WINDOW_RIGHT-E_WINDOW_LEFT)/N_E_WINDOW;
int COINC_LEVEL = 1; // number of coincidences needed per event
int BASELINE_CUT = 70;
int ENERGY_WINDOW_MAX = 100;
Int_t N_INTPOL_SAMPLES = 10; // must be even
Int_t NB = 2;
Int_t ENERGY_NORM = 2500; // Cosmis peak will be normed to this channel
int NB_ACT_CHANNELS = 0; // Number of active channels
double GENERAL_SCALING = 1.; // General histogram scaling for energy calibration
int EXTRACT_PROTO=1; // Extract proto trace and save it to file
double LOWER_RATIO = 0.;
double UPPER_RATIO = 100000.;
// Multi sampling calibration
int MULTIS = 1;
bool MULTIS_CALIB_MODE = false;
// global settings
double THRESHOLD_MULTIPLICY= 1;
bool GLITCH_FILTER = false;
bool SATURATION_FILTER = false;
int GLITCH_FILTER_RANGE = 0;
double CFD_fraction = 0.5;
int L=5; // Length of moving average intervals, centered around current value (mus be odd)
int DELAY=5; // DELAY for the CFD
int M=5; // Window for MWD
double TAU=1.; // Impact value for MA part of MWD
int nTabs = 20; // Number of tabs for the FIR filter
vector<double> FIR_COEF; // Array for the FIR filter coefficients
double LIN_COMP = 0.; // Fit Parameter for linearizing the detector energy sum
// global counters
unsigned int NOE=0;
ReadSystem_class DETECTOR;
// Taggerfile_class tagger;
struct mapping_struct
{
int multis = 1; // multiplication for the multisampling (multisampling)
int polarity = +1; // Polarity of the input signal
// if multisampling is true, save start and end h_channel
int h_channel_start = 0;
int h_channel_stop = 0;
int s_channel = 0; // software channel address
int board_nb = 0; // board number
};
struct hist_struct_TH1D
{
TF1 *fit;
TH1D *hist;
vector<Double_t> params;
};
struct hist_struct_TH2D
{
TF1 *fit;
TH2D *hist;
vector<Double_t> params;
};
struct calib_struct
{
vector<double> multis; // calibration for inter sampling mode
vector<double> RAW_energy; // calibration between different xtals
vector<double> MA_energy; // calibration between different xtals
vector<double> MWD_energy; // calibration between different xtals
vector<double> TMAX_energy; // calibration between different xtals
vector<double> NMO_energy; // calibration between different xtals
hist_struct_TH1D h_RAW_energy; // hist the calibration parameters
hist_struct_TH1D h_MA_energy;
hist_struct_TH1D h_MWD_energy;
hist_struct_TH1D h_TMAX_energy;
hist_struct_TH1D h_NMO_energy;
};
struct multis_norm_struct
{
hist_struct_TH1D h_hist;
Double_t ratio = 0.;
Double_t ratio_err = 0.;
};
struct time_struct
{
vector<hist_struct_TH1D> h_timing;
vector<double> timing; // Is N_E_WINDOWS long when built
};
struct tagger_energy
{
hist_struct_TH1D h_energy; // General energy tagged histogram
hist_struct_TH1D h_energy_m; // Histogram w/o multiples
hist_struct_TH1D h_energy_mt; // Histogram w/o multiples and w/ timing cut
double energy = 0.0;
double energy_m = 0.0;
double energy_mt = 0.0;
//
hist_struct_TH1D h_integral; // General integral tagged histogram
hist_struct_TH1D h_integral_m; // Histogram w/o multiples
hist_struct_TH1D h_integral_mt; // Histogram w/o multiples and w/ timing cut
double integral = 0.0;
double integral_m = 0.0;
double integral_mt = 0.0;
};
struct baseline_struct
{
// Container for storing the baseline traces
vector<Double_t> trace;
// Variables for storing the baseline statistics
double mean = 0.;
double std = 0.;
// Store the thresholds
double TH = 0.;
double TH_multiplicity = 0.;
hist_struct_TH1D h_mean;
hist_struct_TH1D h_std;
hist_struct_TH1D h_samples;
};
struct CFD_struct
{
vector<Double_t> trace; // CFD trace
vector<double> x_interpol; //x- interpolated section of the signal zero crossing
vector<double> y_interpol; //y- interpolated section of the signal zero crossing
double max = 0.;
int Xzero = 0;
int Xzero_int = 0;
double int_x0 = 0.; // the important time information
double int_b = 0.;
double int_m = 0.;
};
struct signal_struct
{
// Trace is saved
vector<Double_t> trace;
// weighted sum of all traces, not to be reset
hist_struct_TH2D TH2D_proto_trace;
vector<Double_t> proto_trace;
vector<Double_t> proto_trace_fit;
vector<Double_t> proto_trace_maxbin;
vector<hist_struct_TH1D> TH1D_proto_trace;
// Baseline is copied and saved separatly here
baseline_struct base;
// CFD conversion of trace
CFD_struct CFD;
// Energy of the signal by pulsehight
double energy = 0.;
int energy_n = 0; // sample number of peak
// Energy of the signal
double integral = 0.;
double ratio = 0.; // ratio of integral/pulseheight
// Marker for signal or non-signal
int is_signal = 0;
// General energy histogram
hist_struct_TH1D h_energy;
// General integral histogram
hist_struct_TH1D h_integral;
hist_struct_TH1D h_ratio; // ratio of integral ober pulseheight
// Tagged energy histograms and energies
vector<tagger_energy> tagged;
// Timing content
vector<time_struct> time;
// Some information about the harware characteritics of the channel
bool is_valid = false; // if this channel is not used for feature extraction, then config file will set this to false
int hardware_addr = 0;
int clock_speed = 100; // in MHz
double sample_t = 10.; // Sampling time in ns
int multis = 1; // Multisampling, number of hardware channels per software channel
int tracelen = 250; // number of samples in trace
int polarity = 1; // Polarity of input signal
bool is_raw = false; // Flag for the RAW container
};
struct tagger_struct
{
double time[16]; // Time information of tagger (arbitrary information)
int counts[16]; // Counts per tagger energy
int multiples_per_count[16]; // Multiples per counts
int multiples_per_channel[16]; // Multiples per tagger channel
vector<hist_struct_TH1D> t_hist; // Histogram for tagger timing distribution
hist_struct_TH2D h_tagger_vs_energy; // Tagger time vs energy in ECAL
int cut[16]; // Mean tagger time for cutting times off the mean time
double energy[16] = {41.3,44.8,51.6,69.14,79.3,99.8,149.37,200.,249.8,350.0,450.2,550.37,650.3,699.7,724.4,725.5}; // Energies corresponding to each tagger channel
};
// Construct containers for storing all wave forms + histograms + various informations
vector<signal_struct> RAW;
vector<signal_struct> RAW_CALIB;
vector<signal_struct> MA;
vector<signal_struct> MWD;
vector<signal_struct> TMAX;
vector<signal_struct> NMO; // NMO: Nelder-Mead Optimization Algorithm
//
vector<signal_struct> ECAL25; // Sum of all channels, only one element in vector
//
// Proto trace
vector<signal_struct> PROTO;
// Struct to save all multi sampling renormalisation histograms and parameters
// Dimensions: [eff. channels][MULTIS][MULTIS]
vector<vector<vector<multis_norm_struct> > > MULTIS_NORM;
// Build root file
TFile *hfile;
// output stream for writing the proto trace
ofstream *proto_out;
// Build the calibration struct
calib_struct CALIB;
// Build the Mapping struct
vector<vector<mapping_struct> > MAPPING;
// Initialize a tagger
tagger_struct TAGGER;
// Container for costum multiplicities
vector<double> TH_MULTIPLICITY;
// Definition of functions
void extraction();
void multis_calib();
void plot_waves(vector<signal_struct> &array, char const *name, char const *modus);
void plot_waves_compare(char const *name, char const *modus);
void plot_interpol(vector<double> &x, vector<double> &y);
void plot_time_energy(time_struct &array);
void plot_energy_hist(vector<signal_struct> &array, char const *path, const char *mode);
void plot_tagger_hist(vector<signal_struct> &array, char const *path, const char *mode);
void plot_timing_hist(vector<signal_struct> &signal, char const *path);
void plot_multis_hist();
void plot_TH2D_hist(TH2D *hist, char const *path, const char *name);
void plot_TH2D_hist_graph(TH2D *hist, vector<double> trace, char const *path, const char *name);
void plot_energy_vs_tagged(signal_struct &signal, char const *path, const char *name);
void plot_energy_resolution(signal_struct &signal, char const *path, const char *name);
void plot_trace(vector<double> &trace, char const *name, char const *modus);
double randit(int ini=0);
double array_mean(const vector<double> &array, int start, int end);
double array_std(const vector<double> &array, int start, int end, double mean);
int array_largest(vector<double> &array, int lower, int upper);
int array_zero_xing(vector<double> &array, int lower, int upper);
double array_compare(vector<double> &array1, vector<double> &array2, vector<double> &weigths, vector<double> &par, int start, int end, int debug);
vector<double> array_adjust(vector<double> &array, vector<double> &x, int debug);
vector<double> array_sum(vector<double> &array1, vector<double> &array2, double factor);
vector<double> array_simulate_proto();
vector<double> array_smooth(vector<double> &array, int s, int L);
void interpolate(vector<signal_struct> &signal);
void time_compare(vector<signal_struct> &signal);
vector<Double_t> fit_hist(TH1D *hist, TF1 *fit, char const *func, Double_t lower = 0, Double_t upper = 1, int verbose = 0);
vector<Double_t> fit_graph_err(TGraphErrors *graph, char const *func, Double_t lower = 0, Double_t upper = 1, int verbose = 0);
void print_usage();
void build_structure();
void print_final_statistics();
void print_energy_statistics(vector<signal_struct> &array, const char *name);
void print_energy_calib();
void print_timing_statistics(time_struct &array, Int_t total_coincidents, const char *name);
void print_stat_multis_calib();
void init_intersamp_hist(int channels);
void init_signal(vector<signal_struct> &signal, int channels, bool is_raw = false);
void init_times(vector<vector<time_struct> > &array, int channels);
void reset_signal(vector<signal_struct> &signal);
void reset_times(vector<vector<time_struct> > &array);
void init_multis_norm(vector<vector<vector<multis_norm_struct> > > &array, int channels);
void fill_hists();
void init_hists(int channels);
bool read_file(string file);
bool read_config(char const *file);
bool linreg(vector<double> &x, vector<double> &y, double *m, double *b);
// Double_t fpeaks(Double_t *x, Double_t *par);
Double_t langaufun(Double_t *x, Double_t *par);
TF1 *langaufit(TH1 *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF, bool silent);
Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM);
Int_t largest_1Dbin(TH1D *hist, Int_t lower, Int_t upper);
vector<double> FIR_filter(vector<double> &trace, double calib);
vector<double> MA_filter(vector<double> &trace, double calib);
vector<double> MWD_filter(vector<double> &trace, double calib);
vector<double> CFD(vector<double> &trace, double multiplier);
void print_detector_config();
bool is_in_string(char const *character, char const *letter);
bool is_glitch(vector<double> &trace, double TH, int n);
bool is_saturation(signal_struct &signal, int n);
bool baseline_weird(signal_struct &signal);
int is_valid_max(signal_struct &signal, int n);
double signal_integral(signal_struct &signal, int debug);
double polnx(double x, vector<double> &par);
double log3x(double x, vector<double> &par);
double ExpDecay1(double x, vector<double> &par);
double ExpGro1(double x, vector<double> &par);
Double_t SIPMpixel(Double_t *x, Double_t *par);
Double_t SIPMlinearize(Double_t x, Double_t A);
Double_t resolution(Double_t *x, Double_t *par);
#include "simplex.h"
int main(int argc, char *argv[])
{
// Timer for program execution time measurement
high_resolution_clock::time_point t_total_begin = high_resolution_clock::now();
///////////////////////////////////////////////////////////////
// Read all existing command line argunments and file names
///////////////////////////////////////////////////////////////
char inputfile[200];
char outputfile[200];
char outputfile_proto_trace[200];
char configfile[200];
unsigned int no_of_events=0, do_no_of_events=0;
int realign_first_NOE=0;
// Parse the input file
if(argc<=1){
printf("No datafile set!!!\n");
print_usage();
return(-1);
}
strcpy(inputfile, argv[1]);
// Parse the output file
if(argc<=2){
printf("No outputfile declared.\n");
sprintf(outputfile,"%s.root",argv[1]);
printf("outputfile set to: %s\n",outputfile);
}
else strcpy(outputfile, argv[2]);
// Parse the config file
bool conf_exists = false;
if(argc<=3){
printf("No config file set, using standard or command line parameters.\n");
}
else {
conf_exists = true;
strcpy(configfile, argv[3]);
printf("Config file read in: %s\n",configfile);
}
// Now parse all command line options
for(int n=0; n<argc; n++){
if(strstr(argv[n],"-n")!=NULL){ // set stop after # of counts
n++;
printf("Max count enabled!\n");
if(n<argc){
do_no_of_events=1;
no_of_events=atoi(argv[n]);
}
else{
printf("Missing max. number of events!\n");
return(-1);
}
}
// Reading of Digital threshold
if(strstr(argv[n],"-t")!=NULL){
n++;
if(n<argc){
THRESHOLD_MULTIPLICY=atoi(argv[n]);
printf("Software threshold multiplicy set to: %3.2f !\n", THRESHOLD_MULTIPLICY);
}
else{
printf("Missing treshold multiplicy!\n");
return(-1);
}
}
// Reading of glitch filter settings
if(strstr(argv[n],"-g")!=NULL){ // set stop after # of counts
n++;
if(n<argc){
GLITCH_FILTER_RANGE=atoi(argv[n]);
GLITCH_FILTER = true;
printf("Glitch filter range set to: %i !\n", GLITCH_FILTER_RANGE);
}
else{
printf("Missing glitch filter parameter!\n");
return(-1);
}
}
if(strstr(argv[n],"-L")!=NULL){ // set stop after # of counts
n++;
if(n<argc){
L=atoi(argv[n]);
printf("Moving average range set to: %i !\n", L);
}
// Check if L is odd
if (L%2 == 0){
printf("SETTINGS ERROR: MA window must be odd!\n");
return(-1);
}
}
if(strstr(argv[n],"-d")!=NULL){ // set stop after # of counts
n++;
if(n<argc){
DELAY=atoi(argv[n]);
printf("DELAY for CFD set to: %i !\n", DELAY);
}
}
if(strstr(argv[n],"-f")!=NULL){ // set stop after # of counts
n++;
if(n<argc){
CFD_fraction=atof(argv[n]);
if ( CFD_fraction < 1. && CFD_fraction > 0. ){
printf("Fraction for CFD set to: %f !\n", CFD_fraction);
}
else{
printf("SETTINGS ERROR: CFD fraction must be set between 0 and 1!\n");
return(-1);
}
}
}
if(strstr(argv[n],"-M")!=NULL){ // set stop after # of counts
n++;
if(n<argc){
M=atoi(argv[n]);
printf("Window length for MWD set to: %i !\n", M);
}
// Check if M is odd
if (M%2 == 0){
printf("SETTINGS ERROR: MWD window must be odd!\n");
return(-1);
}
}
if(strstr(argv[n],"-v")!=NULL){ // turn verbose output on
n++;
if(n<argc){
strcpy(VERBOSE, argv[n]);
printf("Verbose options set to: %s !\n", VERBOSE);
}
}
if(strstr(argv[n],"-I")!=NULL){ // set the number of inter samples
// Check if MULTIS is 1, 2, or 4
n++;
if (atoi(argv[n]) == 1 || atoi(argv[n]) == 2 || atoi(argv[n]) == 4){
if(n<argc){
MULTIS=atoi(argv[n]);
printf("Window length for MWD set to: %i !\n", M);
}
}
else{
printf("SETTINGS ERROR: Intersampling number not 1,2, or 4!\n");
return(-1);
}
}
if(strstr(argv[n],"-r")!=NULL){ // realign_first_NOE
realign_first_NOE=1;
}
if(strstr(argv[n],"-m")!=NULL){ // multiple files with startfile
realign_first_NOE=42;
}
if(strstr(argv[n],"-e")!=NULL){ // turn verbose output on
MULTIS_CALIB_MODE = true;
printf("\n\n++++++++++++++++++++++++++++++++\n");
printf("+++ ENERGY CALIBRATION MODE +++\n");
printf("++++++++++++++++++++++++++++++++\n");
}
if(strstr(argv[n],"-h")!=NULL){ // print help page
print_usage();
return(-1);
}
}
// Check if MA interval length is larger than sample of baseline cut
if (L > BASELINE_CUT){
printf("SETTINGS ERROR: Moving average interval larger than baseline! Not allowed!\n");
return(-1);
}
// Check if MWD interval length is larger than sample of baseline cut
if (M > BASELINE_CUT){
printf("SETTINGS ERROR: MWD moving average interval larger than baseline! Not allowed!\n");
return(-1);
}
//////////////////////////////////////////////////////////
// Initialize detector readout
//////////////////////////////////////////////////////////
// Standard detector analysis mode is cosmics mode, is changed by config file
sprintf(MODE, "COSMICS"); // standard setting
// Read the config file if file is given
if(conf_exists){
bool conf_healthy = read_config(configfile);
if (conf_healthy == false){
printf("ERROR: Config file error!\n");
return(-1);
}
}
// Print mapping if verbose flag is set
if (is_in_string(VERBOSE, "c")) print_detector_config();
// Set up detector read out
DETECTOR.init_readout(inputfile, realign_first_NOE);
//TRACELEN=DETECTOR.get_maxevents(1);
randit(1);
//////////////////////////////////////////////////////////
// Construct TFile for storing all plot/histograms and define folder structure
//////////////////////////////////////////////////////////
//TFile hfile(outputfile,"RECREATE","NTEC analysis");
hfile = new TFile(outputfile,"RECREATE","NTEC analysis");
hfile->SetCompressionLevel(1);
// Create folder structure
build_structure();
// Proto signal is to be extracted
if (EXTRACT_PROTO == 1){
sprintf(outputfile_proto_trace, "%s_proto_trace.dat", outputfile);
proto_out = new ofstream(outputfile_proto_trace);
if (proto_out->is_open()){
printf("NOTICE(main): Opened proto trace output file %s.\n", outputfile_proto_trace);
}
else printf("WARNING(statistics): Unable to write proto trace file.\n");
}
// array_simulate_proto();
////////////////////////////////////////////////////////
// BUILD HISTOGRAMS / LOAD PROTO TRACES FROM FILE
////////////////////////////////////////////////////////
// Reset/initialize signal container
// If in normal extraction mode, use CHANNELS_EFF
if (MULTIS_CALIB_MODE == false){
// Initialize the signal containers
init_signal(RAW, CHANNELS, "true"); // Has to be CHANNELS
init_signal(RAW_CALIB, CHANNELS_EFF);
init_signal(MA, CHANNELS_EFF);
init_signal(MWD, CHANNELS_EFF);
init_signal(TMAX, CHANNELS_EFF);
init_signal(NMO, CHANNELS_EFF);
// Init the Calorimeter sum (only one element per filter type in the vector)
init_signal(ECAL25, 5);
// Initialize histograms, done in extra functionS
init_hists(CHANNELS_EFF);
}
else{
init_signal(RAW, CHANNELS);
// Initialize histograms, done in extra function
// Initialize the MULTIS container
init_multis_norm(MULTIS_NORM, CHANNELS);
init_hists(CHANNELS);
}
// Save the Proto Trace
if (EXTRACT_PROTO == 0){
for (int i = 0; i<(int)RAW_CALIB.size(); i++){
RAW_CALIB[i].proto_trace_fit = PROTO[i].proto_trace_fit;
}
}
// Save costum threshold multiplicies
if (TH_MULTIPLICITY.size() != RAW_CALIB.size() && THRESHOLD_MULTIPLICY == -1){
printf("ERROR: Check your costum threshold multiplicies from the ini file!\n");
return(-1);
}
if (THRESHOLD_MULTIPLICY == -1){
for (int i = 0; i<(int)RAW_CALIB.size(); i++){
RAW_CALIB[i].base.TH_multiplicity = TH_MULTIPLICITY[i];
}
}
else{
for (int i = 0; i<(int)RAW_CALIB.size(); i++){
RAW_CALIB[i].base.TH_multiplicity = THRESHOLD_MULTIPLICY;
}
}
printf("INITIALIZATION COMPLETED! BEGIN READOUT!\n");
/////////////////////////////////////////////////////
// BEGIN SIGNAL READOUT
/////////////////////////////////////////////////////
// Initialize the loop condition
int m=1;
// Reset NOE number of events counter
// loop timer
high_resolution_clock::time_point t_loop_begin = high_resolution_clock::now();
NOE=0;
do{
// Keep reading as long there are unread events
m=DETECTOR.read_one_event(NOE);
// Print some verbose information
if (is_in_string(VERBOSE, "p")){ // "p" for progress report
if (no_of_events != 0) {
if(NOE%((int)no_of_events/10)==0){
high_resolution_clock::time_point t_loop_end = high_resolution_clock::now();
// duration<double, milli> dur = ( t_loop_end - t_loop_begin );
auto duration = duration_cast<milliseconds>( t_loop_end - t_loop_begin ).count();
cout << "Analysing event: " << NOE << " (" << duration/1000 << "s per cycle)" << endl;
}
}
else {
if(NOE%1000 ==0) {
high_resolution_clock::time_point t_loop_end = high_resolution_clock::now();
auto duration = duration_cast<milliseconds>( t_loop_end - t_loop_begin ).count();
cout << "Analysing event: " << NOE << " (" << duration/1000 << "s per cycle)" << endl;
}
}
// reset timer
t_loop_begin = high_resolution_clock::now();
}
// Increase Global number of events counter
NOE++;
// If there is an event, do the extraction
if(m==1){ // && t!=-13){ // not eof for either of these files
// MULTISAMPLING CALIBRATION MODE
// Depending on the program mode, do either calibration or final extraction
//
if(strcmp(MODE, "MULTIS") == 0 || MULTIS_CALIB_MODE == true) {
// Extract calibration information
multis_calib();
// Fill histograms with result, respecting the correct program mode
fill_hists();
}
// FEATURE EXTRACTION MODE
// Depending on the program mode, do either calibration or final extraction
//
else{
// Extract energy and timing information
extraction();
// Fill histograms with result, respecting the correct program mode
fill_hists();
}
// Every 1000 events plot a random event
if (NOE%1000==0){
// Plot every 100th raw signal
hfile->cd("WAVE_FORMS/RAW");
plot_waves(RAW, "Signal_RAW", "TRACE");
}
}
if(do_no_of_events==1){
if(NOE>=no_of_events) m=0;
}
}while(m==1);
/////////////////////////////////////////////////////
// END SIGNAL READOUT
/////////////////////////////////////////////////////
// Print final statistics
if(strcmp(MODE, "MULTIS") == 0 || MULTIS_CALIB_MODE == true) {
print_stat_multis_calib();
plot_multis_hist();
}
else{
print_final_statistics();
}
/////////////////////////////////////////////////////
// END PHYSICS PROGRAM
/////////////////////////////////////////////////////
// End programm timer
high_resolution_clock::time_point t_total_end = high_resolution_clock::now();
auto duration = duration_cast<milliseconds>( t_total_end - t_total_begin ).count();
if (NB_ACT_CHANNELS == 0) NB_ACT_CHANNELS++; // Avoid deviding by zero
cout << endl << "Program exectuion time: " << (duration/1000)
<< "s (" << duration/NB_ACT_CHANNELS/1000 << "s per active channel)" << endl << endl;
// Write and save and delete the root file element in program
printf("Writing root file...\n");
hfile->Write();
// If put stuff after deleting hfile, software might crash. Beware!
delete hfile;
printf("\nRoot file written! Program ends here.\n");
// If put stuff after deleting hfile, software might crash. Beware!
/////////////////////////////////////////////////////
// END OF PROGRAM
/////////////////////////////////////////////////////
}
/////////////////////////////////////////////////////
// ENERGY AND TIMING EXTRACITON
/////////////////////////////////////////////////////
void extraction(){
/////////////////////////////////////////////////////
// INITIALIZE SIGNAL CONTAINERS
/////////////////////////////////////////////////////
reset_signal(RAW); // raw has to be of lentgh CHANNELS
reset_signal(RAW_CALIB);
reset_signal(MA);
reset_signal(MWD);
reset_signal(TMAX);
reset_signal(NMO);
reset_signal(ECAL25);
/////////////////////////////////////////////////////
// FETCH RAW CONTENT FROM DATA
/////////////////////////////////////////////////////
// Create dummy container for fetching the signal
int entry = 0;
unsigned int dumm_cont[CHANNELS][TRACELEN];
// Initialization of the boards (has to start with 1, since 0 is usually the IOL board)
// and fetching the ADC content
for(int board=1; board<=BOARDS; board++){
for(int channel=0; channel<8; channel++){
// Convert the board channel to overall channel
entry = (board-1)*8 + channel;
// Fetch the channels
DETECTOR.get_trace(board, channel, dumm_cont[entry]);
// Save the entries in the RAW container
for(int n=0; n<TRACELEN; n++){
double dumm = (double) dumm_cont[entry][n]; // Has no function but is needed to work (dunno why)
RAW[entry].trace.push_back((double) dumm_cont[entry][n]);
dumm++; // Has no function but is needed to work (dunno why)
}
// Calculate the baseline information for the RAW and RAW_Calib
RAW[entry].base.mean = array_mean(RAW[entry].trace, 0,BASELINE_CUT);
// Already subtract the baseline from the signals
for (int n = 0; n<TRACELEN; n++){
RAW[entry].trace[n] -= RAW[entry].base.mean;
}
// Recalculate the baseline information for the RAW and RAW_Calib
RAW[entry].base.mean = array_mean(RAW[entry].trace, 0,BASELINE_CUT);
// Now calculate the std and TH for the baselines
RAW[entry].base.std = array_std(RAW[entry].trace,0, BASELINE_CUT,RAW[entry].base.mean);
RAW[entry].base.TH = THRESHOLD_MULTIPLICY * RAW[entry].base.std;
}
}
// According to polarity and baseline, transform negative signals into positive signals
for (int a = 0; a<(int)MAPPING.size(); a++){
for (int b = 0; b<(int)MAPPING[a].size(); b++){
// Mapped hardware channel in RAW container
int h = MAPPING[a][b].board_nb * 8 + MAPPING[a][b].h_channel_start;
// If it's a solo channel
if (MAPPING[a][b].multis == 1 && MAPPING[a][b].polarity == -1){
for (int n = 0; n < (int)RAW[h].trace.size(); n++){
RAW[h].trace[n] = 2 * RAW[h].base.mean + (RAW[h].trace[n]*(-1));
}
}
// If not:
if (MAPPING[a][b].multis == 2 && MAPPING[a][b].polarity == -1){
for (int n = 0; n < (int)RAW[h].trace.size(); n++){
RAW[h].trace[n] = 2 * RAW[h].base.mean + (RAW[h].trace[n]*(-1));
RAW[h+1].trace[n] = 2 * RAW[h+1].base.mean + (RAW[h+1].trace[n]*(-1));
}
}
if (MAPPING[a][b].multis == 4 && MAPPING[a][b].polarity == -1){
for (int n = 0; n < (int)RAW[h].trace.size(); n++){
RAW[h].trace[n] = 2 * RAW[h].base.mean + (RAW[h].trace[n]*(-1));
RAW[h+1].trace[n] = 2 * RAW[h+1].base.mean + (RAW[h+1].trace[n]*(-1));
RAW[h+2].trace[n] = 2 * RAW[h+2].base.mean + (RAW[h+2].trace[n]*(-1));
RAW[h+3].trace[n] = 2 * RAW[h+3].base.mean + (RAW[h+3].trace[n]*(-1));
}
}
}
}
// i is the software channel, including all multisampling channels. mapping from hardware channels has to be done
for (int a = 0; a<(int)MAPPING.size(); a++){
for (int b = 0; b<(int)MAPPING[a].size(); b++){
// Chrystal channel (Channel of matrix)
int ch = b*(int)MAPPING.size()+a;
// Mapped hardware channel in RAW container
int h = MAPPING[a][b].board_nb * 8 + MAPPING[a][b].h_channel_start;
// Check if hardware channel is valid channel, otherwise leave that channel out and set the flag
if (MAPPING[a][b].board_nb != 99) {
RAW_CALIB[ch].is_valid = true;
MA[ch].is_valid = true;
MWD[ch].is_valid = true;
TMAX[ch].is_valid = true;
NMO[ch].is_valid = true;
// Aditionally for the RAW channels
RAW[h].is_valid = true;
if (MAPPING[a][b].multis == 2){
RAW[h+1].is_valid = true;
}
if (MAPPING[a][b].multis == 4){
RAW[h+1].is_valid = true;
RAW[h+2].is_valid = true;
RAW[h+3].is_valid = true;
}
}
else continue;
// Now fill the software channels according to the mapping and multisampling
for(int n=0; n<TRACELEN; n++){
if (MAPPING[a][b].multis == 1){
RAW_CALIB[ch].trace.push_back( ((double) RAW[h].trace[n] )* CALIB.multis[h] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
// Update the channel info
RAW_CALIB[ch].hardware_addr = h;
RAW_CALIB[ch].polarity = MAPPING[a][b].polarity;
}
else if (MAPPING[a][b].multis == 2){
// Save the samples correctly
RAW_CALIB[ch].trace.push_back((double) RAW[h+1].trace[n] * CALIB.multis[h+1] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
RAW_CALIB[ch].trace.push_back((double) RAW[h].trace[n] * CALIB.multis[h] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
// Update the channel information
RAW_CALIB[ch].multis = 2;
RAW_CALIB[ch].clock_speed = 200; // in MHz
RAW_CALIB[ch].sample_t = 5.; // in MHz
RAW_CALIB[ch].tracelen = TRACELEN * 2;
RAW_CALIB[ch].hardware_addr = h;
RAW_CALIB[ch].polarity = MAPPING[a][b].polarity;
}
else if (MAPPING[a][b].multis == 4){
// Save the samples correctly
RAW_CALIB[ch].trace.push_back((double) RAW[h+3].trace[n] * CALIB.multis[h+3] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
RAW_CALIB[ch].trace.push_back((double) RAW[h+2].trace[n] * CALIB.multis[h+2] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
RAW_CALIB[ch].trace.push_back((double) RAW[h+1].trace[n] * CALIB.multis[h+1] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
RAW_CALIB[ch].trace.push_back((double) RAW[h].trace[n] * CALIB.multis[h] * CALIB.RAW_energy[ch] * GENERAL_SCALING);
// Update the channel information
RAW_CALIB[ch].multis = 4;
RAW_CALIB[ch].clock_speed = 400; // in MHz
RAW_CALIB[ch].sample_t = 2.5; // in MHz
RAW_CALIB[ch].tracelen = TRACELEN * 4;
RAW_CALIB[ch].hardware_addr = h;
RAW_CALIB[ch].polarity = MAPPING[a][b].polarity;
}
}
}
}
for(int board=1; board<=BOARDS; board++){
for(int channel=0; channel<8; channel++){
// Convert the board channel to overall channel
entry = (board-1)*8 + channel;
RAW[entry].CFD.trace.clear();
RAW[entry].CFD.trace = CFD(RAW[entry].trace, 1);
}
}
for (int i = 0; i < (int)RAW_CALIB.size(); i++){
/////////////////////////////////////////////////////
// APPLICATION OF THE FILTER TO RAW_CALIB
/////////////////////////////////////////////////////
MA[i].trace.clear();
MWD[i].trace.clear();
TMAX[i].trace.clear();
MA[i].trace = MA_filter(RAW_CALIB[i].trace, CALIB.MA_energy[i]);
MWD[i].trace = MWD_filter(RAW_CALIB[i].trace, CALIB.MWD_energy[i]);
TMAX[i].trace = FIR_filter(RAW_CALIB[i].trace, CALIB.TMAX_energy[i]);
for (int n = 0; n < (int)TMAX[i].trace.size(); n++){
// printf("%3.1f %3.1f %3.1f\n", RAW_CALIB[i].trace[n], TMAX[i].trace[n], CALIB.TMAX_energy[i] );
}
/////////////////////////////////////////////////////
// APPLICATION OF THE CONSTANT FRACTION DISCRIMINATOR
/////////////////////////////////////////////////////
RAW_CALIB[i].CFD.trace.clear();
MA[i].CFD.trace.clear();
MWD[i].CFD.trace.clear();
TMAX[i].CFD.trace.clear();
RAW_CALIB[i].CFD.trace = CFD(RAW_CALIB[i].trace, 1);
MA[i].CFD.trace = CFD(MA[i].trace, 1);
MWD[i].CFD.trace = CFD(MWD[i].trace, 1);
TMAX[i].CFD.trace = CFD(TMAX[i].trace, 1);
}
/////////////////////////////////////////////////////
// CALCULATE BASELINE STATISTICS
/////////////////////////////////////////////////////
// Calculate the baseline statistics
for (int i = 0; i < (int)RAW_CALIB.size(); i++){
// Check if channel is valid
if (RAW_CALIB[i].is_valid == false) continue;
// otherwise do the calculation
RAW_CALIB[i].base.mean = array_mean(RAW_CALIB[i].trace, 0,BASELINE_CUT*RAW_CALIB[i].multis);
RAW_CALIB[i].base.std = array_std(RAW_CALIB[i].trace, 0, BASELINE_CUT*RAW_CALIB[i].multis,RAW_CALIB[i].base.mean);
RAW_CALIB[i].base.TH = RAW_CALIB[i].base.TH_multiplicity * RAW_CALIB[i].base.std;
//
MA[i].base.mean = array_mean(MA[i].trace, 0,BASELINE_CUT*MA[i].multis);
MA[i].base.std = array_std(MA[i].trace, 0, BASELINE_CUT*MA[i].multis,MA[i].base.mean);
MA[i].base.TH = RAW_CALIB[i].base.TH_multiplicity * MA[i].base.std;
//
MWD[i].base.mean = array_mean(MWD[i].trace, 0,BASELINE_CUT*MWD[i].multis);
MWD[i].base.std = array_std(MWD[i].trace, 0, BASELINE_CUT*MWD[i].multis,MWD[i].base.mean);
MWD[i].base.TH = RAW_CALIB[i].base.TH_multiplicity * MWD[i].base.std;
//
TMAX[i].base.mean = array_mean(TMAX[i].trace, 0,BASELINE_CUT*TMAX[i].multis);
TMAX[i].base.std = array_std(TMAX[i].trace, 0, BASELINE_CUT*TMAX[i].multis,TMAX[i].base.mean);
TMAX[i].base.TH = RAW_CALIB[i].base.TH_multiplicity * TMAX[i].base.std;
}
/////////////////////////////////////////////////////
// EXTRACT FEATURES FROM TRACES
/////////////////////////////////////////////////////
// First the RAW traces
for(int i=0; i<(int)RAW.size(); i++){
// Extract the maximum
RAW[i].energy = RAW[i].trace[array_largest(RAW[i].trace, BASELINE_CUT, ENERGY_WINDOW_MAX)];
// Look for CFD Zero Crossing
RAW[i].CFD.Xzero = array_zero_xing(RAW[i].CFD.trace, BASELINE_CUT, ENERGY_WINDOW_MAX);
// Already fill the samples into a histogram for baseline noise estimation
for(int n=0; n < BASELINE_CUT; n++){
RAW[i].base.h_samples.hist->Fill(RAW[i].trace[n]);
}
}
int plot = 0;
//
// Then the Filtered traces
//
// In beam mode first check if the central crystal is healthy
int central_healty = 0;
// Left and Right borders for maximum detection
int left = BASELINE_CUT * RAW_CALIB[0].multis;
int right = ENERGY_WINDOW_MAX * RAW_CALIB[0].multis;
// Already extract features of RAW_CALIB
if (strcmp(MODE, "BEAM")==0){
RAW_CALIB[CENTRAL].energy_n = array_largest(RAW_CALIB[CENTRAL].trace, left, right); // sample number of largest sample
RAW_CALIB[CENTRAL].energy = RAW_CALIB[CENTRAL].trace[RAW_CALIB[CENTRAL].energy_n];
RAW_CALIB[CENTRAL].integral = signal_integral(RAW_CALIB[CENTRAL], 0);
RAW_CALIB[CENTRAL].ratio = RAW_CALIB[CENTRAL].integral / RAW_CALIB[CENTRAL].energy;
if ( is_valid_max( RAW_CALIB[CENTRAL], RAW_CALIB[CENTRAL].energy_n ) == 0 ){
if ( RAW_CALIB[CENTRAL].ratio < UPPER_RATIO && RAW_CALIB[CENTRAL].ratio > LOWER_RATIO){
RAW_CALIB[CENTRAL].is_signal = 1;
central_healty = 1;
}
else central_healty = 0;
}
else central_healty = 0;
}
else central_healty = 1;
// Now loop through all the channels
for(int i=0; i<(int)RAW_CALIB.size(); i++){
// plot = 1;
// Only extract features from valid channels
if (RAW_CALIB[i].is_valid == false) continue;
// Left and Right borders for maximum detection
int left = BASELINE_CUT * RAW_CALIB[i].multis;
int right = ENERGY_WINDOW_MAX * RAW_CALIB[i].multis;