-
Notifications
You must be signed in to change notification settings - Fork 0
/
sfaplot.C
executable file
·3822 lines (3274 loc) · 136 KB
/
sfaplot.C
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
//sfa7
// group of functions for loading and analysing the data produced by vmelist2root
#define NDETS 24
#include "TROOT.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TCutG.h"
#include "TFormula.h"
#include "TRandom.h"
#include "TFunction.h"
#include "TMethodCall.h"
#include "TObjString.h"
#include "TError.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TColor.h"
#include "TDiamond.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TGaxis.h"
#include "TVector.h"
#include "TMatrix.h"
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <cstring>
#include <vector>
using namespace std;
struct event{
Int_t mod; // module number - bit useless
Int_t trig_det; // the detector number that had the highest energy of the last adc module... bit useless
Double_t hi_tim; // high presision time. Multiply by 8E-6 to convert to seconds. Cycles every few hours?
Double_t time; // time of event in unix time
Double_t energy[NDETS]; // energy of detectors
Int_t status[4]; // temperature1, temperature2, humidity, flowmeter
Int_t flag; // event flag. Reject if flag==1;
Int_t veto; // used to store TTL flag.
Int_t runno; // the squential run number that the event belongs to
Int_t groupno; // number to which the calibration group the event belongs. starts at 0
};
#define EVTREE2buffer sprintf(buffer,\
"mod/I:trig_det/I:hi_tim/D:time/D:energy[%d]/D:status[4]/I:flag/I:veto/I:runno/I:groupno/I",\
NDETS);
struct run_stats{
Int_t firstevno; // event number of the first event in the run
Int_t starttime; // time of the start of the run in unix time
Int_t groupno; // the group that the run belongs to
Int_t flag[NDETS]; // flags used to reject runs. reject if flag[det-1]==1
Int_t flaggedtime[NDETS]; // used for small cuts... not yet used.
Int_t nevents; // number of saved events per run
Int_t nvetoevents; // the number of HV veto events
Int_t neventsall; // actual number of events in daq file. Does not including nvetoedevents. Deadtime = (neventsall+nevetoedevents)* deadTimePerEvent
Int_t det_on[NDETS]; // whether the detector is on for this run. This number is set by the daq. Use data if det_on==1
Int_t hv[NDETS]; // the high voltage setting reported by the daq
Int_t amp[NDETS]; // the amplifacation setting (0-127) reported by the daq
Int_t thr[NDETS]; // the thresholds setting (0-127) reported by the daq
Int_t runlen; // the run length in seconds (usually 3600s unless short runs are accepted when compiling data)
Double_t lo_thresh[NDETS]; // the threshold above which the data can be analysied (keV)
};
#define RSTREE2buffer sprintf(buffer,\
"firstevno/I:starttime/I:groupno/I:flag[%d]/I:flaggedtime[%d]/I:nevents/I:nvetoedevents/I:neventsall/I:det_on[%d]/I:hv[%d]/I:amp[%d]/I:thr[%d]/I:runlen/I:lo_thresh[%d]/D",\
NDETS,NDETS,NDETS,NDETS,NDETS,NDETS,NDETS);
struct group_stats{
Int_t firstrunno; // run number of first run in group
Int_t starttime[NDETS]; // start time of group (for each detector) // not much use
Int_t nruns; // number of runs in a group
Int_t runs_toobig; // number of runs rejected due to file size cut
Int_t det_on[NDETS]; // specifies if a detector is on or off. This is set by the user in group_stats
Double_t hi_thresh[NDETS]; // upper energy threshold of data
Double_t lo_thresh[NDETS]; // lower threshold can be set high by user in groups_stats. This is then moved over to rs.lo_thresh in sfaplot.
Double_t cal_eqn[4][NDETS]; // ={slope,error,constant,error}
Double_t res_eqn[4][NDETS]; // ={slope,error,constant,error}
Int_t group_name[50]; // name of the group. Replace each number with its corespondic ASCI character to reviel.
Int_t flag; // to flag whole group and groups sellections. Useful for comparing results from one group against another.
Int_t nevents; // number of events in the group
};
#define GSTREE2buffer sprintf(buffer,\
"firstrunno/I:starttime[%d]/I:nruns/I:runs_toobig/I:det_on[%d]/I:hi_thresh[%d]/D:lo_thresh[%d]/D:cal_eqn[4][%d]/D:res_eqn[4][%d]/D:group_name[50]/I:flag/I:nevents/I",\
NDETS,NDETS,NDETS,NDETS,NDETS,NDETS);
struct eppstruct{
vector <Double_t> epp; // number of events per period
vector <Double_t> time; // time of start of period
vector <Int_t> runno; // sequential run number of first run in period
vector <Int_t> used_run_cnt; // number of runs used in period
} epp;
struct CutStats{
Double_t mean;
Double_t error;
Double_t probX2;
Int_t availruns;
Int_t extracut;
} *cs;
// #ifdef __CINT__
// #pragma link C++ class std::vector<Int_t>+;
// #pragma link C++ class std::vector<Double_t>+;
// #endif
#define SKIPFLAGS() \
while((gs[ev.groupno].flag || !gs[ev.groupno].det_on[detn-1]) && i<nentries){ \
i += gs[ev.groupno].nevents; \
evtree->GetEvent(i); \
} \
\
while( (!rs[ev.runno].det_on || (use_flags[1]!='_' && rs[ev.runno].flag[detn-1]==(use_flags[1]=='0'))) && i<nentries){ \
i += rs[ev.runno].nevents; \
evtree->GetEvent(i); \
} \
// Global variables
event ev;
run_stats *rs;
group_stats *gs;
Double_t gbl_deadtimeperevent=1.2E-4; // this is the approx deadtime per event Uncertainty? What impact will this have? Check!!!
Double_t samescalef=1;
Long64_t nentries;
Int_t ngroups, nruns;
Int_t ndetectors = NDETS;
Int_t *flag;
Int_t *BFlag;
TCanvas *c1;
TGraph *g_waterfall;
TGraph *g_handle, *g_handle2;
TH1D *h_handle, *h_handle2;
TH2D *h2_handle, *h2_handle2;
Char_t *fname;
TFile *ssfile;
TTree *evtree;
Int_t fit_locut =0, fit_hicut = 0; // used for output of fit_poisson
Int_t gblperiod =1; // global user period, ie the number of runs per period
Int_t gbl_verbose=1;
Int_t gbl_draw=1;
TF1 *gblfit; // used for fit_poisson. needs to be global if I want to be able to get information of fit externally.
int gblpause=0;
Int_t gbldetx, gbldety;
vector <Double_t> gblresults;
// vector <Double_t> *gblcuts = new vector <Double_t> [NDETS];
// vector <Double_t> *newcuts = new vector <Double_t> [NDETS];
// decleration of some functions
Double_t SSgetdeadtime(Int_t det);
Double_t SSgetusedtime(Int_t detn, Int_t incl_flaggedtime=1, Int_t incl_deadtime=1);
//Int_t getentries(Int_t group=-1, Int_t return_output=0); outof date
void llfn_getgroupname(int groupno);
void SShelp();
void llfn_makewhite(TCanvas *c1);
void SSgroups_set(Char_t use_groups[20]=" ");
Int_t SSgroups_get();
void SSfit_poisson(Double_t mean =-1, TH1D *hist1 = h_handle, Int_t min=0, Int_t max=100);
void SSfit_gaus(Double_t mean =-1., Double_t sigma = -0.10, Double_t scale = -1., TH1D *hist1 = h_handle);
void SSfit_poly();
void llfn_poly1(Double_t *a, Double_t x1, Double_t y1, Double_t x2, Double_t y2);
//subfunctions
char *llfn_getname(FILE *fd);
Int_t llfn_time2runno(Double_t time);
Double_t llfn_time2unixtime(Double_t time);
void llfn_runno2runname(Int_t starttime, char *runlist);
Double_t SSgetuptime(Int_t detn);
Double_t SSgetflaggedtime(Int_t detn);
Double_t SSgetusedtime(Int_t detn, Int_t incl_flaggedtime, Int_t incl_deadtime);
Double_t SSgetusedtimeR(Int_t detn, Double_t Emin, Double_t Emax, Double_t deadtimeperevent=1.2E-4);
Double_t SSgetdeadtime(Int_t det);
Int_t SSgetdatetimedif(int i=0, int f=nentries-1);
void SSgetResEqn(Int_t group, Int_t det);
void SSaddleggend(char title1[],char title2[],TH1D *h_handle1, TH1D *h_handle2);
void SSprint(Int_t opt = 0, char name[64]="c1", TCanvas *c1=c1);
void llfn_setstyles(char OptStat[16]="e");
void SScheckorder();
// Double_t SSgetcuttime(Int_t detn);
// void llfn_addcuts();
int llfn_NOT(int val, int sign);
Int_t llfn_Decimal2Bit(Int_t decimal, Int_t bitN, Int_t Nbits);
/// ################################################ FUNDERMENTAL FUNCTIONS ##############################################
/// ---------------------------------------------- low level functions -----------------------------------------------
/// ############################################ END OF FUNDERMENTAL FUNCTIONS ###########################################
/// #################################################### MAIN FUNCTIONS ##################################################
/// #################################################### OLD FUNCTIONS ##################################################
// calibration functuion threshold number to channel number
// TF1 SSth2ch("th2ch","26.6*x+24+5",0,3000);
// TF1 SSch2th("ch2th","1./26.6*x-24/26.6",0,127);
TF1 SSth2ch("th2ch","26.6*x+24+26",0,3000);
TF1 SSch2th("ch2th","1./26.6*x-50/26.6",0,127);
Double_t SSe2ch(Double_t E, Int_t det, Int_t group=1) {
return E/gs[group-1].cal_eqn[0][det-1] - gs[group-1].cal_eqn[2][det-1]/gs[group-1].cal_eqn[0][det-1];
}
Double_t SSch2e(Double_t Ch, Int_t det, Int_t group=1) {
return gs[group-1].cal_eqn[0][det-1]*Ch + gs[group-1].cal_eqn[2][det-1];
}
Double_t SSget_coincident_usedtime(Int_t *dets) {
// this function gets the used time for periods when all of the "dets" are not flagged.
// The function only takes into account runs so far
int i,j;
TH1D *f1 = new TH1D("f1","f1",nruns,0,nruns);
TH1D *f2 = new TH1D("f2","f2",nruns,0,nruns);
// put the used runs of dets[0] into a histogram
for(j=0;j<nruns;++j) if(rs[j].flag[dets[0]] && gs[rs[j].groupno].flag) f1->Fill(j,1);
for(i=1;dets[i];++i) {
// put the used runs of the next dets[] into a histogram
for(j=0;j<nruns;++j) if(rs[j].flag[dets[i]] && !gs[rs[j].groupno].flag) f2->Fill(j,1);
// f2 to f1.
f1->Add(f2);
}
Int_t time=0;
for(j=0;j<nruns;++j) if(f1->GetBinContent(j+1)==0) ++time;
return Double_t(time);
}
void SScut_plot_coincedencesHelp() {
cout << "This function serves two porposes. \n"
<< "It can plot coinidences and or can flag coincidences of event across detectors.\n\n"
<< "Specify the detectors to the analysed by creating an integer array of the required\n"
<< " detector numbers setting the last one in the list to 0 which marks the end of the\n"
<< " list similar to the \\0 character in a string\n"
<< "eg: $ Int_t dets[]={2,3,5,7,8,0}; // selects the current LHS detectors only\n\n";
cout << "flagevents = 0 - Do not flag coincident events\n"
<< " = 1 - Do flag coincident events\n"
<< " = 2, 3, etc - Flag only events were there 2, 3 etc coincidences\n\n";
cout << "detn = 0, - Coincidences between all detectors are used\n"
<< " = det# - Coincidences between det# (replace with number) and all other\n"
<< " detectors is used\n\n";
cout << "plot_layer = 0 - Search for coincidences between all layers of detectors\n"
<< " (2D graph is not drawn) \n"
<< " = layer# - Search for coincidences only within specified layer \n"
<< " A 2D colour plot of all the coincidence counts is also plotted\n"
<< " where the detectors are numbered as they are in Gran Sasso\n"
<< " and the detector between coincidences are search, is marked\n\n";
cout << "Emin = 0 - Do not restrict lower energy range\n"
<< " = # (>0) - Restrict search above specified energy (keV)\n"
<< " NB The lower thresholds are always restricted by rs.lo_thresh\n\n";
cout << "Emax = 0 - Do not restrict uppper energy rangea\n"
<< " = # (>0) - restrict search below specfied energy (keV)\n\n";
cout << "Do not use the remaining options for now\n\n";
cout << "So what do the plots mean? I'll update that section soon!\n";
}
void SScut_plot_coincedences(Int_t *dets, Int_t flagevents=0, Char_t use_flags[3] = "000", Int_t detn=0, Int_t plot_layer=0, Double_t Emin=0, Double_t Emax=0, Int_t useDetnDefaultThresh = 0, Double_t cutsecsbefore=0., Double_t cutsecsafter=0.){
//This function flags events that are coincident with another event in the list of searched detectors "dets"
// NB if flag events > 1 then only coincidences of greater than flagevents are flagged.
// if detn==0 then coincidences between all 'dets' are used. Else coincidences between detn is only used.
// plots are only made when 'flagevents' is not zero.
// plot_layer: set to 0 to plot all layers (2D graph is not drawn) otherwise set the the layer that should be plotted.
// if dets=0 then using defaults:{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,0};, to use 'dets' create an array of detectors to be used
// set useDetnDefaultThresh: specify if lo_thresh (0) of the run or the trigger threshold (1) is used
int i, j, k, cnt;
char buffer[50];
TDiamond *diamond=0;
Double_t cur_time;
vector <Int_t> coincidence_counts(NDETS,0);
vector <Int_t> coincidence_counts_total(NDETS,0);
// for(j=0;j<NDETS;j++) newcuts[j].clear();
Int_t default_dets[]={ 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,*/ 0};
if(!dets) dets = default_dets;
int plot_energies=0;//!flagevents;
int count=3;
if(plot_energies) gblresults.clear();
Double_t depossitionE[NDETS];
for(i=0;i<NDETS;++i) depossitionE[i]=0;
// cout << "hi1\n";
//Int_t *Nhits = new Int_t[NDETS];
//for(i=0;i<NDETS;i++) Nhits[i] = 0;
h_handle = new TH1D("Coincidences","Coincidences",NDETS-1,1.,NDETS);
if(Emin>0) sprintf(buffer,"Number depositions between %1.0lf and %1.0lfkeV",Emin,Emax);
else sprintf(buffer,"Number depositions between");
h_handle->GetXaxis()->SetTitle(buffer);
h2_handle = new TH2D("CoinArray","Coincidence Array",4,0.5,4.5,4,0.5,4.5);
Double_t lower_thresh = 0;
int cnt2=0, cnt3=0, cnt4=0;
for(cnt=0,i=0,evtree->GetEvent(0);i<nentries;++i,evtree->GetEvent(i)){
for(j=0;dets[j];++j) {
if(useDetnDefaultThresh) lower_thresh = SSch2e(SSth2ch.Eval(rs[ev.runno].thr[dets[j]-1]),dets[j],ev.groupno+1);
else lower_thresh = rs[ev.runno].lo_thresh[dets[j]-1];
if(
((ev.energy[dets[j]-1] > Emin && Emax ==0) || (ev.energy[dets[j]-1] > Emin && ev.energy[dets[j]-1] < Emax)) &&
// check if within allowed thresholds
(ev.energy[dets[j]-1] > lower_thresh && (Emax == 0 || (Emax < gs[ev.groupno].hi_thresh[dets[j]-1]))) &&
rs[ev.runno].det_on[dets[j]-1] &&
(use_flags[0]=='_' || !rs[ev.runno].flag[dets[j]-1]) &&
gs[ev.groupno].det_on[dets[j]-1] &&
!gs[ev.groupno].flag &&
// old 12/03/2008 &&(use_flags==0 || (use_flags==1 && !flag[i] && !rs[ev.runno].flag[dets[j]-1])) ) {
(use_flags[0]=='_' || ((use_flags[0]=='0') && !flag[i]) || ((use_flags[0]=='1') && flag[i]) || (use_flags[0]==flag[i])) &&
(use_flags[2]=='_' || !llfn_Decimal2Bit(BFlag[i],detn,NDETS)==(use_flags[2]=='0'))
) {
// cout << cnt << endl;
++cnt; //check for coincidences
++coincidence_counts[dets[j]-1];
if(plot_energies) depossitionE[dets[j]-1]=ev.energy[dets[j]-1];
}
}
if(cnt>1) ++cnt2; // is cnt2 not used?
if(detn!=0 && coincidence_counts[detn-1]==0){ // check that detn fired, reject if it did not.
cnt=0;
for(k=0;k<NDETS;k++) coincidence_counts[k]=0;
}
else for(k=0;k<NDETS;k++) coincidence_counts_total[k]+=coincidence_counts[k];
h_handle->Fill(cnt);
if(plot_energies) {
if(cnt==count) {
gblresults.push_back(i);
for(j=0;j<NDETS;++j) {
if(depossitionE[j]>0) {
gblresults.push_back(j+1);
gblresults.push_back(depossitionE[j]);
}
}
gblresults.push_back(-1);
}
for(j=0;j<NDETS;++j) depossitionE[j] = 0;
}
if(flagevents > 0 && cnt>flagevents) {
cur_time = ev.time;
if(cur_time-cutsecsbefore <ev.time ) { // go back "cutsecsbefore" seconds
for(;(cur_time-cutsecsbefore)<ev.time && i>=0;i--,evtree->GetEvent(i));
i++;
evtree->GetEvent(i);
}
++cnt3;
for(;(cur_time+cutsecsafter)>=ev.time && i<nentries;i++,evtree->GetEvent(i)) {flag[i]=1, ++cnt4;};
--i;
}
//delete counters
cnt=0;
for(k=0;k<NDETS;k++) coincidence_counts[k]=0;
}
// cout << "hi5\n";
if(flagevents==0) {
c1 = new TCanvas("c1","c1",10,10,600,400);
c1->SetLogy();
c1->cd();
h_handle->Draw();
int k_start=(plot_layer-1)*16;
if(plot_layer>0) {
TCanvas *c2 = new TCanvas("c2","c2",10,10,600,400);
c2->cd();
for(i=4,k=k_start;i>0 && k<(k_start+16);--i) {
for(j=1;j<5 && k<(k_start+16);++j) {
if(k==detn-1) {
h2_handle->Fill(i,j,0);
if(int((detn-1)/16)+1==plot_layer) diamond = new TDiamond(i-0.25,j-0.25,i+0.25,j+0.25);
}
else h2_handle->Fill(i,j,coincidence_counts_total[k]);
k++;
}
}
h2_handle->Draw("colz");
if(int((detn-1)/16)+1==plot_layer) {
diamond->SetFillColor(1);
diamond->SetTextAlign(22);
diamond->SetTextSize(0);
diamond->SetBorderSize(0);
diamond->Draw();
}
}
}
// cout << "hi6\n";
if(plot_energies) {
TCanvas *c3 = new TCanvas("c3","c3",10,10,600,400);c3->cd();
TH2D *h2_handle3= new TH2D("h3"," ",4,0.5,4.5,4,0.5,4.5);
for(i=1;i<int(gblresults.size());++i) {
while(i< int(gblresults.size()) && gblresults[i]!=-1) {
h2_handle3->Fill( -1*(int(gblresults[i])-16)/4+1, (int(gblresults[i])%4!=0)*int(gblresults[i])%4+(int(gblresults[i])%4==0)*4, gblresults[i+1]);
cout << "det: " << gblresults[i] << " filling " << -1*(int(gblresults[i])-16)/4+1 << " " << (int(gblresults[i])%4!=0)*int(gblresults[i])%4+(int(gblresults[i])%4==0)*4 << " " << gblresults[i+1] << endl;
i+=2;
}
cout << endl;
++i;
h2_handle3->Draw("colz");
// c3->Update();
// gPad->WaitPrimitive();
sprintf(buffer,"coin%d.gif",i);
c3->Print(buffer);
h2_handle3->Reset();
}
}
// llfn_addcuts();
cout << cnt2 << " " << cnt3 << " " << cnt4 << endl;
cout << "NB in the event of a call for a flag, all consequative events with the same system time are flagged.\n";
cout << "To remove this feature the hi_time must be compared. This is essential for shorter time cuts anyway\n";
// cout << ".........\nNB this function does not work well across groups. (does not do range safety checks)\n";
cout << "\n\nUsed time = " << SSget_coincident_usedtime(dets) << " hours\n\n";
}
// because I don't trust SScut_plot_coincidences I'm writing a new simple function...
// the conclusion was that it does not work either. Basically this remove all the events from about 2 groups
void SScut_coincidences() {
// this function flags all data that is coincident with other events above the trigger threshold
// the only problem in doing this is that it is possible that a noise event that has been cut out with the cleaning is included in this
// search and could be coincident with a real event, which is then rejected.
int i, j, cnt;
for(i=0,evtree->GetEvent(0);i<nentries;++i,evtree->GetEvent(i)){
cnt = 0;
for(j=0;j<NDETS;++j) {
if( rs[ev.runno].det_on[j] &&
// gs[ev.groupno].det_on[j] &&
// !gs[ev.groupno].flag &&
(ev.energy[j] > SSch2e(SSth2ch.Eval(rs[ev.runno].thr[j]),j+1,ev.groupno+1))
) ++cnt;
}
if(cnt>1) flag[i]=1;
}
}
// NB itoa can be used to convert an integer to binay string
//#include <string.h>
int llfn_Binary2Decimal(char *bin){
int dec = 0, fac = 1;
char *s = bin + strlen(bin) -1;
while (s >= bin) {
if (*s == '1') {
dec += fac;
}
fac *= 2;
s--;
}
return dec;
}
void SSflag_greatest(Int_t plot_triggered=0){
int i, j;
char *expandedBFlag = new char[NDETS+3];
for(i=0;i<NDETS;i++) expandedBFlag[i] = '1'; expandedBFlag[NDETS] = '\0';
Double_t Egreatest;
Int_t DetNo = 0;
h_handle = new TH1D("Triggered","Triggered",NDETS,0.,NDETS);
for(i=0,evtree->GetEvent(0);i<nentries;i++,evtree->GetEvent(i)){
Egreatest = 0.;
for(j=0;j<NDETS;j++){
if(ev.energy[j] > Egreatest){
Egreatest = ev.energy[j];
DetNo = j;
}
}
if(plot_triggered) h_handle->Fill(DetNo);
expandedBFlag[DetNo]='0';
BFlag[i] = llfn_Binary2Decimal(expandedBFlag);
expandedBFlag[DetNo]='1';
// cout << DetNo << " " << BFlag[i] << " " << llfn_Decimal2Bit(BFlag[i],detn,NDETS)<< endl;
}
if(plot_triggered) h_handle->Draw();
delete expandedBFlag;
}
// TH1F *th1_thresholds;
void SSlocate_Thresholds(){
int i, j, k, largest_thresh=0, first_run=0, numberOfRuns=0;
for(i=0;i<ngroups;i++){
numberOfRuns += gs[rs[first_run].groupno].nruns;
// cout << "group " << i << endl;
for(j=0;j<NDETS;j++){
//locate largest threshold used in group
// cout << "locating largest threshold of det " << j << endl;
for(k=first_run,largest_thresh=0; k<numberOfRuns; ++k) if(rs[k].thr[j]>rs[largest_thresh].thr[j]) largest_thresh=k;
//set thresholds
// cout << "setting threshold\n";
// cout << "det = " << j+1 << " largest = " << largest_thresh << " ch = " << SSth2ch.Eval(rs[largest_thresh].thr[j]) << endl;
gs[i].lo_thresh[j] = gs[i].cal_eqn[0][j]*(SSth2ch.Eval(rs[largest_thresh].thr[j])) + gs[i].cal_eqn[2][j];
}
// cout << "setting first_run to " << first_run << endl;
first_run += gs[i].nruns;
}
}
void SSplot_settings_thr(char help='H'){
help = help;
cout << "SSplot_settings_thr plots the threshold settings as a function of time\n";
cout << "To select plotting formate enter mode number:\n";
cout << "0 ..... threshold numbers are plotted\n";
cout << "1 ..... thresholds in channels are plotted\n";
cout << "2 ..... thresholds in energy are plotted\n";
}
void SSplot_settings_thr(Int_t detn, int mode =0, int drawPlot=1) {
Double_t *run_time = new Double_t[nruns];
Double_t *run_thr = new Double_t[nruns];
for(int run_counter=0;run_counter<nruns;++run_counter) {
run_thr[run_counter] = rs[run_counter].thr[detn-1];
if(mode==1) run_thr[run_counter] = SSth2ch.Eval(rs[run_counter].thr[detn-1]);
if(mode==2) run_thr[run_counter] = gs[rs[run_counter].groupno].cal_eqn[0][detn-1]*(SSth2ch.Eval(rs[run_counter].thr[detn-1])) + gs[rs[run_counter].groupno].cal_eqn[2][detn-1];
run_time[run_counter] = rs[run_counter].starttime;/*cout << run_thr[run_counter] << " " << run_time[run_counter] << endl;*/
}
g_handle = new TGraph(nruns,run_time,run_thr);
// g_handle->SetStats(0);
if(drawPlot) {
g_handle->GetXaxis()->SetTimeDisplay(1);
g_handle->GetXaxis()->SetTimeFormat("");
g_handle->GetXaxis()->SetNdivisions(-510);
g_handle->GetXaxis()->SetTitle("Date");
g_handle->GetYaxis()->SetTitle("Threshold Number");
if(mode==1) g_handle->GetYaxis()->SetTitle("Threshold (Channels)");
if(mode==2) g_handle->GetYaxis()->SetTitle("Threshold (keV)");
g_handle->Draw("APl");
}
}
void SSplot_settings_userthr(Int_t detn, Char_t plot_opt[10]="same") {
Double_t time_scale=gblperiod;
if(time_scale==0) h_handle = new TH1D("thresholds","run thresholds",(rs[nruns-1].starttime-rs[0].starttime)/3600,rs[0].starttime,rs[nruns-1].starttime);
else h_handle = new TH1D("thresholds","run thresholds",(rs[nruns-1].starttime-rs[0].starttime)/3600,0,(rs[nruns-1].starttime-rs[0].starttime)/3600*3600./time_scale);
for(int run_counter=0;run_counter<nruns;++run_counter) {
if(time_scale==0) h_handle->Fill(rs[run_counter].starttime,rs[run_counter].lo_thresh[detn-1]);
else h_handle->Fill(Double_t(rs[run_counter].starttime-rs[0].starttime)/time_scale,rs[run_counter].lo_thresh[detn-1]);
}
if(time_scale==0) {
h_handle->GetXaxis()->SetTimeDisplay(1);
h_handle->GetXaxis()->SetTimeFormat("");
h_handle->GetXaxis()->SetNdivisions(-510);
h_handle->GetXaxis()->SetTitle("Date");
}
else h_handle->GetXaxis()->SetTitle("Time");
h_handle->GetYaxis()->SetTitle("Threshold (keV)");
h_handle->SetLineColor(kBlue);
h_handle->Draw(plot_opt);
}
void SSplot_settings_hv(Int_t detn, int drawPlot=1) {
Double_t *run_time = new Double_t[nruns];
Double_t *setting = new Double_t[nruns];
for(int run_counter=0;run_counter<nruns;++run_counter) {
setting[run_counter] = rs[run_counter].hv[detn-1];
run_time[run_counter] = rs[run_counter].starttime;/*cout << run_thr[run_counter] << " " << run_time[run_counter] << endl;*/
}
g_handle = new TGraph(nruns,run_time,setting);
// g_handle->SetStats(0);
c1 = new TCanvas("settings", "hv settings",14,30,700,500);
if(drawPlot) {
g_handle->GetXaxis()->SetTimeDisplay(1);
g_handle->GetXaxis()->SetTimeFormat("");
g_handle->GetXaxis()->SetNdivisions(-510);
g_handle->GetXaxis()->SetTitle("Date");
g_handle->GetYaxis()->SetTitle("Voltage (V)");
g_handle->Draw("APl");
}
}
void SSlocate_NoiseThresholds(int detn, double tolerance = 0.1){
int i, j;
Int_t cur_group=-1, cur_run=-1, lo_thresh=0;
int nbins=0, last_nruns=0;
int *fraction_ok_regions=NULL;
// function of typical approximate spectrum Normalised to Counts/10keV/hr
TF1 SSGausFunc("gausfunc","1.11*TMath::Gaus(x,110,98)",0,500);
TH1D * spec=NULL;
for(i=0,evtree->GetEvent(i); i<nentries ;i++,evtree->GetEvent(i)){
// SKIPFLAGS(); //skip flagged groups and runs (defined above)
// cout << "h1\n";
if(ev.groupno!=cur_group){
// test and set thresholds
if(fraction_ok_regions) {
for(j=nbins-1;j>=0;j--) {
if(double(fraction_ok_regions[j])/last_nruns < tolerance && gs[cur_group].lo_thresh[detn-1] < spec->GetBinCenter(j+1) + 5) {
gs[cur_group].lo_thresh[detn-1] = spec->GetBinCenter(j+1) + 5;
break;
}
}
}
for(j=0;j<nbins;j++) printf("%3.2lf ",double(fraction_ok_regions[j])/last_nruns); cout << endl;
for(j=0;j<nbins;j++) printf("%3.0lf ",spec->GetBinCenter(j+1)); cout << endl;
// cout << "h2\n";
cur_group=ev.groupno;
if(int(gs[ev.groupno].lo_thresh[detn-1])%10) lo_thresh=Int_t(gs[ev.groupno].lo_thresh[detn-1])+10;
else lo_thresh=Int_t(gs[ev.groupno].lo_thresh[detn-1]);
if(lo_thresh<=400) nbins = (400-lo_thresh)/10 ;
else nbins = 1;
last_nruns = gs[ev.groupno].nruns;
if(fraction_ok_regions) delete fraction_ok_regions;
// cout << nbins << endl;
fraction_ok_regions = new int[nbins];
// cout << "h6\n";
for(j=0;j<nbins;j++) fraction_ok_regions[j] = 0;
}
if(ev.runno!=cur_run) {
// cout << "h3\n";
cur_run = ev.runno;
// count ok regions if counts are less than 4 sigmas away from standard spectrum
if(spec) for(j=0;j<nbins;j++) if(spec->GetBinContent(j+1) < SSGausFunc.Eval(spec->GetBinCenter(j+1)) + 4.*sqrt(SSGausFunc.Eval(spec->GetBinCenter(j+1)))) fraction_ok_regions[j]++;
// if(spec) cout << spec->GetBinContent(3) << " " << GausFunc.Eval(spec->GetBinCenter(3)) + 4.*sqrt(GausFunc.Eval(spec->GetBinCenter(3))) << endl;
// cout << "h4\n";
// create new test spectrum
if(spec) spec->Delete();
// cout << "h5\n";
spec = new TH1D("spec","spec",nbins,lo_thresh,400);
}
spec->Fill(ev.energy[detn-1]);
}
if(fraction_ok_regions)
for(j=nbins-1;j>=0;j--)
if(double(fraction_ok_regions[j])/last_nruns < tolerance)
if(gs[cur_group].lo_thresh[detn-1] > spec->GetBinCenter(j+1) + 10) gs[cur_group].lo_thresh[detn-1] = spec->GetBinCenter(j+1) + 5;
for(j=0;j<nbins;j++) printf("%1.2lf ",double(fraction_ok_regions[j])/last_nruns); cout << endl;
for(j=0;j<nbins;j++) printf("%3.0lf ",spec->GetBinCenter(j+1)); cout << endl;
}
void SSplot_spec2(char help='H'){
help = help;
cout << "Use plot_spec2 to plot the spectrum of detector \"detn\" using the lower threshold rs.lo_thresh\n";
cout << " where the life time is calculated for each keV bin. The spectrum is plotted with h_handle. \n";
cout << " eg you can use h_handle to redraw (h_handle->Draw()). The lifetime is saved to h_handle2.\n\n";
cout << "rebin: define binning. NB histogramed is rescaled so that it is still in bins/keV\n";
cout << "use_flags: flag options: _ = ignore, 0 = do not use flagged data, 1 = only use flagged data\n";
cout << " 1st number is ev flags. That is the flag fag for each event (blind to detector number)\n";
cout << " 2nd number is run flags. An hour of data is rejected if flagged.\n";
cout << " 3rd number is Binary ev flag. With this flag the specific detector numbers can be flagged\n";
cout << "sumw2: set this to draw spectrum with errors\n";
cout << "rawdata: set this if you want the scale in Counts/keV. NB rebinned data is also in /keV.\n";
cout << "reset_used_groups: set this to reset groups sellection with \"groups_set();\"\n";
}
void SSplot_spec2(Int_t detn, Int_t rebin = 1, Char_t use_flags[3] = "000", Int_t use_sumw2=0, Int_t rawdata=0, Double_t undercut=10.){
// This function plots a spectrum like plot_spec() but uses thresholds from each run and life times per keV bin to
// generate the spectrum. This way the full data can be used.
// undercut is the amount bellow the upper threshold that the data is used and avoids inclusion of saturation events.
Char_t buffer[512];
if(gbl_draw) {
c1 = new TCanvas("cspec", "cspec",14,30,700,500);
c1->SetLogy();
}
int i, j, k, j_start;
int specUpperRange = 12000;
// initialise plot ranges
sprintf(buffer,"D%d",detn);
h_handle = new TH1D("hspec",buffer,specUpperRange,0.,specUpperRange);
sprintf(buffer,"D%d time",detn);
h_handle2 = new TH1D("spec_time",buffer,specUpperRange,0.,specUpperRange);
i=0;
evtree->GetEvent(i);
for(i=0,j=0,j_start=0,evtree->GetEvent(j);i<nruns;++i) { //loop over runs
// if run or group is flagged, skip to next run (set j to i+=event in run
// skip run if...
if( !gs[rs[i].groupno].det_on[detn-1] || // skip off detectors for whole group
gs[rs[i].groupno].flag || // skip flagged groups
!rs[i].det_on[detn-1] || // skip off detectors
(use_flags[1]!='_' && rs[i].flag[detn-1]==(use_flags[1]=='0')) // skip flagged runs
) {
j = j_start+rs[i].nevents;
j_start = j;
// maybe check that it is still in sync (ie that the previous event if of the previous run
// idealy it would go like this: j = rs[i+1].evnumber which would take it to the first event of the i+1'th run
continue;
}
for(evtree->GetEvent(j);j<j_start+rs[i].nevents;++j,evtree->GetEvent(j)) { // loop over events in run (could have been j<rs[i+1].evnumber)
// check if event is flagged
if( int(rs[i].lo_thresh[detn-1]) < ev.energy[detn-1] && // Fill if above lower run threshold
(gs[rs[i].groupno].hi_thresh[detn-1]-undercut > ev.energy[detn-1]) && // Fill if below upper group threshold (so far there is no upper run threshold as it is unlikely to be used)
(use_flags[0]=='_' || !flag[j]==(use_flags[0]=='0')) && // Fill if using evflags and is not evflagged
(use_flags[2]=='_' || !llfn_Decimal2Bit(BFlag[j],detn,NDETS)==(use_flags[2]=='0')) ) { // Fill if using evBitFlagged and is not evBitFlagged (evBitFlag is a binary number which is decomposed into a flag for each detector)
// Add event to spectrum
h_handle->Fill(ev.energy[detn-1]+(gRandom->Uniform()-0.5));
}
}
// add lifetime in seconds to runningtime histogram
for(k=int(rs[i].lo_thresh[detn-1]);Double_t(k)<gs[rs[i].groupno].hi_thresh[detn-1]-undercut && k<specUpperRange;++k)
h_handle2->Fill(k,3600.-(rs[i].neventsall+rs[i].nvetoevents)*gbl_deadtimeperevent); //
j = j_start+rs[i].nevents;
j_start = j;
evtree->GetEvent(j);
//if(i+1!=ev.runno) cout << "Warning, out of sync!\n";
}
//export?
// Add error bars are required
if(use_sumw2) {
h_handle->Sumw2();
h_handle2->Sumw2();
}
// normalise / rebin
// before 28/02/2008 a detector mass of 5.78g was assumed so any plots must be rescaled accordingly.
// an average mass of 6.5g is a better approximation. although idealy the precise mass should be used.
if(!rawdata) h_handle->Divide(h_handle,h_handle2,3600.*24./6.5e-3);
//rebin
h_handle->Rebin(rebin); // I check, these functions work as predicted. Only point is that the error on 1 and and 0 is sqrt(1).
h_handle->Scale(1./rebin);
// label and plot if required
if(gs[0].cal_eqn[0][detn-1]!=1.) h_handle->GetXaxis()->SetTitle("Energy (keV)");
else h_handle->GetXaxis()->SetTitle("Channels");
if(rawdata && gs[0].cal_eqn[0][detn-1]!=1.) h_handle->GetYaxis()->SetTitle("Counts/keV");
else if(rawdata && gs[0].cal_eqn[0][detn-1]==1.) h_handle->GetYaxis()->SetTitle("Counts/Channel");
else if(!rawdata && gs[0].cal_eqn[0][detn-1]!=1.) h_handle->GetYaxis()->SetTitle("Counts/keV/kg/day");
else if(!rawdata && gs[0].cal_eqn[0][detn-1]==1.) h_handle->GetYaxis()->SetTitle("Counts/Channel/kg/day");
if(gbl_draw) h_handle->Draw();
}
void SSplot_spec(char help='H'){
help = help;
cout << "Use plot_spec to plot the spectra of detector number \"detn\"\n\n";
cout << "rebin: define binning. NB histogramed is rescaled so that it is still in bins/keV\n";
cout << "Emin: Lower energy on x-axis for spectrum. Group used only if between gs[i].hi_thresh[detn-1] and gs[i].lo_thresh[detn-1]\n";
cout << "Emax: Higher energy boundery. Set to zero to use maximum but note that different groups mays have different end points.\n";
cout << "use_flags: flag options: _ = ignore, 0 = do not use flagged data, 1 = only use flagged data\n";
cout << " 1st number is ev flags. That is the flag fag for each event (blind to detector number)\n";
cout << " 2nd number is run flags. An hour of data is rejected if flagged.\n";
cout << " 3rd number is Binary ev flag. With this flag the specific detector numbers can be flagged\n";
cout << "rawdata: set this if you want the scale in Counts/keV. NB rebin data is also in /keV.\n";
cout << "sumw2: set this to draw spectrum with errors\n";
cout << "close_file: set this to save spectrum to spec.C, instead of plotting\n";
cout << "reset_used_groups: set this to reset groups sellection with \"groups_set();\"\n";
}
void SSplot_spec(Int_t detn, Int_t rebin = 1, Double_t Emin = 120, Double_t Emax = 10000, Char_t use_flags[3] = "000", Int_t rawdata=0, Int_t use_sumw2=0, Int_t closefile=0, Int_t reset_used_groups =1){
Char_t buffer[512];
//TROOT simple("simple","nTuple");
gROOT->SetStyle("Plain");
//sprintf(fname,"%s.root", argv[1]);
// TFile *hfile = new TFile("hist.root","RECREATE","hist",9);
TFile *hfile = new TFile("hist.root","UPDATE","hist",9);
if(gbl_draw) {
c1 = new TCanvas("spec", "spec",14,30,700,500);
c1->SetLogy();
}
// llfn_makewhite(spec);
int i, j;
int *used_group = new int[ngroups];
for(i=0;i<ngroups;i++) used_group[i] = 0;
Double_t used_time =0;
sprintf(buffer,"D%d",detn);
TH1D *h_spec = new TH1D("spec",buffer,Int_t(Emax-Emin),Emin,Emax);
if(reset_used_groups) SSgroups_set();
for(i=0,evtree->GetEvent(i); i<nentries ;i++,evtree->GetEvent(i)){
//for(j=1;j<=24;j++) cout << llfn_Decimal2Bit(BFlag[i],j,24); cout << endl;
SKIPFLAGS(); //skip flagged groups and runs (defined above)
if((/*gs[ev.groupno].lo_thresh[detn-1] <= Emin && */gs[ev.groupno].hi_thresh[detn-1] >= Emax) // don't include groups that have smaller ranged (this is now redundent now that rs.lo_thresh has been introduced)
&& (use_flags[0]=='_' || !flag[i]==(use_flags[0]=='0'))
&& (use_flags[2]=='_' || !llfn_Decimal2Bit(BFlag[i],detn,NDETS)==(use_flags[2]=='0'))
&& (Emin < ev.energy[detn-1])
/// NB I just added the following line and that the life time is calcalated with the new function...
&& (Emin > rs[ev.runno].lo_thresh[detn-1])){
h_spec->Fill(ev.energy[detn-1]+(gRandom->Uniform()-0.5));
used_group[ev.groupno] = 1; //could use a more efficient way then this as
}
}
//flag skipped groups
for(i=0;i<ngroups;i++) if(!used_group[i]) gs[i].flag = 1;
if(use_sumw2) h_spec->Sumw2();
h_spec->Rebin(rebin);
h_spec->Scale(1./rebin);
//Get Running Time.....
used_time = SSgetusedtimeR(detn,Emin,Emax);
if(use_flags[1]!='0') cout << "WARNING: don't trust the running time\n";
// Scale to g/keV/day
if(!rawdata) {
// h_spec->Scale(3600.*24./5.78e-3/used_time);
// h_spec->SetBinContent(0,5.78e-3*used_time/3600./24.); // changed 28/02/2008
h_spec->Scale(3600.*24./6.5e-3/used_time);
h_spec->SetBinContent(0,6.5e-3*used_time/3600./24.);
}
else h_spec->SetBinContent(0,used_time);
//h_spec->SetEntries(used_time);
// sprintf(buffer,"D%d %1.2f hours. (Deadtime of %1.2f hours subtracted)",detn,used_time/3600.,getdeadtime(detn)/3600);
sprintf(buffer,"D%d %1.2f hours.",detn,used_time/3600.);
h_spec->SetTitle(buffer);
if(gs[0].cal_eqn[0][detn-1]!=1.) h_spec->GetXaxis()->SetTitle("Energy (keV)");
else h_spec->GetXaxis()->SetTitle("Channels");
if(rawdata && gs[0].cal_eqn[0][detn-1]!=1.) h_spec->GetYaxis()->SetTitle("Counts/keV");
else if(rawdata && gs[0].cal_eqn[0][detn-1]==1.) h_spec->GetYaxis()->SetTitle("Counts/Channel");
else if(!rawdata && gs[0].cal_eqn[0][detn-1]!=1.) h_spec->GetYaxis()->SetTitle("Counts/keV/kg/day");
else if(!rawdata && gs[0].cal_eqn[0][detn-1]==1.) h_spec->GetYaxis()->SetTitle("Counts/Channel/kg/day");
if(gbl_draw) h_spec->Draw();
h_handle = h_spec;
if(gblpause) gPad->WaitPrimitive();
hfile->Write();
if(closefile) hfile->Close();
if(gbl_verbose){
// print out which groups were used and total timed used.
cout << "The following groups were used for this plot:\n";
j = SSgroups_get();
cout << j << "/" << ngroups << " groups used. Total " << used_time/3600 << " hours of data.\n\n";
cout << "includes deadtime = " << SSgetdeadtime(detn)/3600 << " hours. (" << SSgetdeadtime(detn)/used_time*100 << "%)\n";
}
}
// use detn = 0 to analyse all detectors
void SSplot_epp(Int_t detn, Double_t period = 1E-3, Int_t flag_finds = 0, Double_t Emin=150, Double_t Emax=0, Double_t range_Xmax = 100){
// if(time_period>=3600) {
// cout << "This function is only for period searches
// plot events per given time period. If trig_filter = 0 then applies to all detectors.
// Energy range only works for detn != 0 currantly. This could be changed to "if any detector detectors energy > 0" with a fors loop and break statements.
//!!!! WARNING: if search periods of greater than 1 hour then the periods will not be complete before a group (or set) is complete! ie, every
c1 = new TCanvas("c1", "c1",14,30,700,500);
c1->SetLogy();
int i, j, k, l;
int eventCounter;
Double_t period_start, period_start_hitime, skip_steps;//, hi_time_start;
Char_t buffer[100];
// save the event number of multiple events so that they can be flagged
vector <int> flagged_event_log;
// ensure that 3600 is devisable by period
if(period>3600.) period = 3600.;
period = 3600./int(3600./period);
sprintf(buffer,"Events per %f seconds. D%d",period,detn);
h_handle = new TH1D("Events per hour",buffer,Int_t(range_Xmax),0,range_Xmax);
h_handle->GetXaxis()->SetTitle("number of events");
h_handle->GetYaxis()->SetTitle("frequecy");
for(i=0;i<ngroups;++i) { // loop over groups
// skip group if flagged
if( !gs[i].det_on[detn-1] || gs[i].flag /*|| Emax > gs[i].hi_thresh[detn-1] */) continue;
for(j=gs[i].firstrunno;j<nruns && j<gs[i+1].firstrunno;++j) { // loop over runs
// skip runs if flagged
if( !rs[j].det_on[detn-1] || rs[j].flag[detn-1] /*|| Emin < rs[j].lo_thresh[detn-1]*//*SSch2e(SSth2ch.Eval(rs[j].thr[detn-1]),detn)*/ ) continue;
// initialise everything... counting is not done accross runs with this tool
eventCounter = 0;
flagged_event_log.clear();
period_start = rs[j].starttime;
evtree->GetEvent(rs[j].firstevno);
period_start_hitime = ev.hi_tim + (rs[j].starttime - ev.time)/8e-6; // set to the value the hi_time would have had at start of run
for(k=rs[j].firstevno;k<nentries && k<rs[j+1].firstevno;++k,evtree->GetEvent(k)) { // loop over event in run
if( flag[k] ||
llfn_Decimal2Bit(BFlag[k],detn,NDETS) ||
ev.energy[detn-1] < Emin ||
ev.energy[detn-1] > Emax ||
ev.energy[detn-1] < rs[j].lo_thresh[detn-1] ||
ev.energy[detn-1] > gs[i].hi_thresh[detn-1] ||
ev.energy[detn-1] < SSth2ch.Eval(rs[j].thr[detn-1])
) continue;
if(period >= 1) { // use unix time
if(ev.time - period_start > period) {
h_handle->Fill(eventCounter);
if(flag_finds && eventCounter>flag_finds) for(l=0;l<int(flagged_event_log.size());++l) flag[flagged_event_log[l]]='p';
eventCounter = 0;
flagged_event_log.clear();
period_start += period;
--k; // count this event in next period
}
else {
++eventCounter;
flagged_event_log.push_back(k);
}
}
else { // time stamp
if((ev.hi_tim - period_start_hitime)*8e-6 > period) {
if(eventCounter>0) {
h_handle->Fill(eventCounter);
if(flag_finds && int(flagged_event_log.size())>flag_finds) for(l=0;l<int(flagged_event_log.size());++l) flag[flagged_event_log[l]]='p';
flagged_event_log.clear();
eventCounter = 0;
period_start_hitime += period/8e-6;
}
else { // it is important not to recheck un filled regions to save time...
skip_steps = (ev.hi_tim - period_start_hitime)*8e-6/period;
for(l=0;l<skip_steps;++l) h_handle->Fill(0);
period_start_hitime += skip_steps*period/8e-6;
}
--k; // count this event in next period
}
else if((ev.hi_tim - period_start_hitime)*8e-6 < 0) { // check if counter has restarted
period_start_hitime = period_start_hitime - pow(2.,32) -1; // hopefully this sets the timer to the new time
}
else {
++eventCounter;
flagged_event_log.push_back(k);