-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathModuleCalibration.cc
More file actions
6273 lines (5387 loc) · 331 KB
/
Copy pathModuleCalibration.cc
File metadata and controls
6273 lines (5387 loc) · 331 KB
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
//------------------------------------------------------------//
// //
// PROGRAM FOR matrix calibration USING DT5740 //
// //
//------------------------------------------------------------//
// compile with cmake and make
// this program is meant to read a "short" acquisition data and find the limits that will be used by the reconstruction algorithm to
// do the images. Limits will be searched on u,v,w and e, i.e
// u = weighted average in the x direction
// v = weighted average in the y direction
// w = ratio highest charge / total charge collected -> converted to z
// e = energy of the incident particle, to select photopeaks
// CAREFUL: this program is meant to run on a single module, so if the acquisition is done in parallel on the two modules, it is
// supposed to be run twice, with different configuration files that will select the proper input channels, and with, as an output,
// different text files
// CAREFUL: at the moment this program is used to characterize different modules. So instead of really looking for the limits, it takes
// limits (for the crystal positions in u,v) and computes some relevat paramenters, especially the peak position, energy resolutions and
// width (FWHM) of the w histograms for each crystal.
// The user can input the crystal limits in the config file, only crystals that are specified in the config file will be analized (see newBoardConfig.cfg)
// PROGRAM STRUCTURE
// This file performes the relevant analysis, while the module structure is defined by the Element, Module, Mppc and Crystal classes. Module, Mppc and Crystal
// inherit from Element, but also have some specific method
// Config files are read in various parts of the program, via the ConfigFile class, to input the user parameters
// So ModuleCalibration does:
// - some initial check
// - read the input, via the InputFile class (see those files). This will read the Tchain of input root files and creat a specific TTree that will be used for the analysis
// - read config parameters
// - Create the modules, mppcs, crystals, according to the user requests
// - Run on all the elements and produce the desided plots. Mainly 2D and 3D histograms, crystal and mppc spectra, w histograms and so on (see in the code below)
// - Creates some Canvases for nice plot storing
// - Creates an output file with a directory structure and saves the plots
// - Optionally, saves the analysis TTree in a root file
#include "TROOT.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3I.h"
#include "TString.h"
#include "TApplication.h"
#include "TLegend.h"
#include "TTree.h"
#include "TFile.h"
#include "TF2.h"
#include "TGraph2D.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "TSpectrum2.h"
#include "TTreeFormula.h"
#include "TMath.h"
#include "TChain.h"
#include "TCut.h"
#include "TLine.h"
#include "TError.h"
#include "TEllipse.h"
#include "TFormula.h"
#include "TGraphErrors.h"
#include "TGraph2DErrors.h"
#include "TMultiGraph.h"
#include "TCutG.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TProfile.h"
#include "TH1D.h"
#include "TPaveText.h"
#include "TGraphDelaunay.h"
#include "TFitResult.h"
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <algorithm>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <cmath>
#include <assert.h>
#include <cstddef>
// #include <thread>
#include <getopt.h>
#include "ConfigFile.h"
#include "InputFile.h"
#include "Element.h"
#include "Crystal.h"
#include "Module.h"
#include "Mppc.h"
#include "Utilities.h"
#include "./libraries/Extract.h"
// #include <omp.h>
//test
#define ENERGY_RESOLUTION 0.07
#define ENERGY_RESOLUTION_SATURATION_CORRECTION 0.12
int main (int argc, char** argv)
{
// gStyle->SetOptStat(0);
//----------------------------------------------------------//
// Check input //
//----------------------------------------------------------//
if(argc<2)
{
std::cout << " Usage: " << std::endl;
std::cout << " ModuleCalibration [-c config file ] <input-files> " << std::endl;
std::cout << " note: -c option is optional, but if you use it, it has to be the first argument" << std::endl;
std::cout << " without it, a configFile.cfg will be assumed" << std::endl;
return 0;
}
//----------------------------------------------------------//
// Say hello to the user //
//----------------------------------------------------------//
std::cout<<"\n"<<std::endl;
std::cout<<"###########################################################"<<std::endl;
std::cout<<"# #"<<std::endl;
std::cout<<"# Module calibration #"<<std::endl;
std::cout<<"# #"<<std::endl;
std::cout<<"###########################################################"<<std::endl;
std::cout<<"\n\n"<<std::endl;
std::cout<<"=====> C O N F I G U R A T I O N <====\n"<<std::endl;
//----------------------------------------------------------//
// Enable Implicit MultiThreading of ROOT //
//----------------------------------------------------------//
// unsigned int nthreads = std::thread::hardware_concurrency();
// std::cout << nthreads << " concurrent threads are supported " << std::endl;
// ROOT::EnableImplicitMT(nthreads-1);
// std::cout << nthreads-1 << " concurrent threads will be used " << std::endl;
//----------------------------------------------------------//
// Import input and parse the config file //
//----------------------------------------------------------//
// Set a default config file name
std::string ConfigFileName = "config.cfg";
// and assume as default that there is no config file name from command line
// then check
if(std::string(argv[1]) == std::string("-c")) // first argument is -c, then the config file name is passed by command line
{
ConfigFileName = argv[2];
std::cout << "Configuration file: '" << argv[2] << "'"<< std::endl;
}
else // the config file was indeed the default one
{
std::cout << "Configuration file set to default: config.cfg "<< std::endl;
}
//PUT THE ENTIRE CONFIG FILE IN A STRINGSTREAM, TO SAVE IT LATER IN THE OUTPUT FILE
std::stringstream streamConfigFile;
std::string strConfig;
std::ifstream inConfig;
inConfig.open (ConfigFileName.c_str(), std::ifstream::in);
while (std::getline(inConfig, strConfig))
{
streamConfigFile << strConfig << std::endl;
}
inConfig.close();
ConfigFile config(ConfigFileName); // create a ConfigFile object
InputFile input(config); // read the input chain of root files, passing the inputs and the config object
/////////////////////////////////////////////
// READ CONFIG FILE //
////////////////////////////////////////////7
// std::string chainName = config.read<std::string>("chainName","adc");
// int digitizerTotalCh = config.read<int>("digitizerTotalCh");
// int digitizerType = config.read<int>("digitizerType",0);
std::string fname = config.read<std::string>("chainName","adc");
int ncrystalsx = config.read<int>("ncrystalsx",2); // number of crystals in x direction per mppc - default to 2 if the key is not found in the config file
int ncrystalsy = config.read<int>("ncrystalsy",2); // number of crystals in y direction per mppc - default to 2 if the key is not found in the config file
int nmppcx = config.read<int>("nmppcx",2); // number of mppc in x direction per mppc - default to 2 if the key is not found in the config file
int nmppcy = config.read<int>("nmppcy",2); // number of mppc in y direction per mppc - default to 2 if the key is not found in the config file
int nmodulex = config.read<int>("nmodulex",1); // number of modules in x direction per mppc - default to 1 if the key is not found in the config file
int nmoduley = config.read<int>("nmoduley",1); // number of modules in y direction per mppc - default to 1 if the key is not found in the config file
int histo1Dmax = config.read<int>("histo1Dmax",25000); // max of the 1D charge histograms (in ADC channels)
int histo1Dbins = config.read<int>("histo1Dbins",250); // number of bins of the 1D charge histograms
int histoLYmax = config.read<int>("histoLYmax",20000); // max of the 1D charge histograms (in ADC channels)
int histoLYbins = config.read<int>("histoLYbins",500); // number of bins of the 1D charge histograms
float qe = config.read<float>("qe",0.35); // mppc QE
float gainMPPC = config.read<float>("gainMPPC",1.25e6); // mppc gain
float chargeBinningADC = config.read<float>("chargeBinningADC",156e-15); // adc charge binning
float sourceMeV = config.read<float>("sourceMeV",0.511); // gamma peak in MeV
int histo2DchannelBin = config.read<int>("histo2DchannelBin",250); // number of bins of the 2D flood histograms, for single channels
int histo2DglobalBins = config.read<int>("histo2DglobalBins",1000); // number of bins of the 2D flood histograms, for entire module
int histo3DglobalBins = config.read<int>("histo3DglobalBins",100); // number of bins of the 3D flood histograms, for entire module
int taggingPeakMin = config.read<int>("taggingPeakMin",0); // min range of tagging crystal photopeak, in ADC channels - to help TSpectrum
int taggingPeakMax = config.read<int>("taggingPeakMax",0); // max range of tagging crystal photopeak, in ADC channels - to help TSpectrum
// int clusterLevelPrecision = config.read<int>("clusterLevelPrecision",10); // precision of the level search when separating the cluster of 3D points
float taggingPosition = config.read<float>("taggingPosition",0); // position of the tagging bench in mm
bool usingTaggingBench = config.read<bool>("usingTaggingBench",0); // true if the input is using tagging bench, false if not
int taggingCrystalChannel = config.read<int>("taggingCrystalChannel",16); // input channel where the tagging crystal information is stored
bool calcDoiResWithDelta = config.read<bool>("calcDoiResWithDelta",0); // alternative calcolation of doi res, based on deltas. only if it's usingTaggingBench
std::string calcDoiFileName = config.read<std::string>("calcDoiFileName",""); // name (and path if not in this folder) of calibration_params.txt
int pointsFromDoi = config.read<int>("pointsFromDoi",0); // points measured in file calcDoiFileName. this is MANDATORY if calcDoiFileName is specified!
bool correctingSaturation = config.read<bool>("correctingSaturation"); // true if saturation correction is applied, false if it's not
bool correctingForDOI = config.read<bool>("correctingForDOI",0); // true if the energy correction using DOI info is computed
float energyResolution = config.read<float>("expectedEnergyResolution",0); // energy resolution input by the user, if any, otherwise 0
bool usingRealSimData = config.read<bool>("usingRealSimData",0);
float moduleLateralSideX = config.read<float>("moduleLateralSideX",7.0); //
float moduleLateralSideY = config.read<float>("moduleLateralSideY",7.0); //
bool backgroundRun = config.read<bool>("backgroundRun",0); // whether this is a background run or not
bool lateralRun = config.read<bool>("lateralRun",0); // whether this is a lateral irradiation run or not
float userBroadCut = config.read<float>("userBroadCut",0.0); // if in backgroundRun, cut to get rid of low energy events is not done on photopeak search but by user input (default 0ch)
float thresholdKev = config.read<float>("thresholdKev",0.0);
float wThreshold = config.read<float>("wThreshold",0.1); // Threshold for w plots limits
float crystalz = config.read<float>("crystalz",15);
// double DoiResolutionVsIJmax = config.read<double>("DoiResolutionVsIJmax",10); // max of the 2d DoiResolution values plot (starts from 0) - it's mm
float EnergyResolutionVsIJmax= config.read<float>("EnergyResolutionVsIJmax",0.3); // max of the 2d EnergyResolution values plot (starts from 0)
float LYvsIJmax = config.read<float>("LYvsIJmax",40000);
float PeakPositionVsIJmax = config.read<float>("PeakPositionVsIJmax",12000); // max of the 2d PeakPosition values plot (starts from 0) - it's ADC channels
int wHistogramsBins = config.read<int>("wHistogramsBins",250);
int doiColumnOffset = config.read<int>("doiColumnOffset",0); // for DOI output, fix the column i by adding this quantity. if not stated, 0 by default
float energyCorrectionMin = config.read<float>("energyCorrectionMin",0.25); // once the wmin and wmax are found for each w histo, choose at which point to start and to stop the linear fitting
float energyCorrectionMax = config.read<float>("energyCorrectionMax",0.75); // (as percentage from min to max)
float lambda511 = config.read<float>("lambda511",12.195); //everything in mm
// bool wAllChannels = config.read<bool>("wAllChannels",0); // whether we use the sum of all channels to compute w (true = 1) of just the neighbours (false = 0). Deafult to false.
bool comptonAnalysis = config.read<bool>("comptonAnalysis",0); //whether to perform or not the compton recovery analysis part. Default to false
bool lightYieldComputation = config.read<bool>("lightYieldComputation",0); //whether to perform or not the light yield calculation. Default to false
float peakSearchRangeMin = config.read<int>("peakSearchRangeMin",0); //lower limit for search of 511KeV peak - wide limitation to help peak search
float peakSearchRangeMax = config.read<int>("peakSearchRangeMax",histo1Dmax); //upper limit for search of 511KeV peak - wide limitation to help peak search
bool saturationRun = config.read<bool>("saturationRun",0); //whether this is a saturation run or not
float histoSingleChargeMax = config.read<float>("histoSingleChargeMax",0); // max in the histograms of charge for saturationRun
float histoSingleChargeBin = config.read<float>("histoSingleChargeBin",0); // max in the histograms of charge for saturationRun
float histoSumChargeMax = config.read<float>("histoSumChargeMax",0); // max in the histograms of charge for saturationRun
float histoSumChargeBin = config.read<float>("histoSumChargeBin",0); // max in the histograms of charge for saturationRun
float saturationPeakFractionLow = config.read<float>("saturationPeakFractionLow",0.06); // lower limit for saturation peak fitting, expressed in fraction of peak position. limit will be = (peak - peak*saturationPeakFractionLow)
float saturationPeakFractionHigh = config.read<float>("saturationPeakFractionHigh",0.06); // upper limit for saturation peak fitting, expressed in fraction of peak position. limit will be = (peak + peak*saturationPeakFractionHigh)
bool performSaturationPeakSearch = config.read<bool>("performSaturationPeakSearch",1); //perform of nor satuartion peak search
bool backgroundSaturationRun = config.read<bool>("backgroundSaturationRun",0);
// set output file name
std::string outputFileName = config.read<std::string>("output");
std::string digitizer_s = config.read<std::string>("digitizer");
std::string saturationPeakEnergy_s = config.read<std::string>("saturationPeakEnergy","");
std::string saturationPeakMin_s = config.read<std::string>("saturationPeakMin","");
std::string saturationPeakMax_s = config.read<std::string>("saturationPeakMax","");
std::string loadAnalysisTreeName = config.read<std::string>("loadAnalysisTreeName","0"); // look for input analysis ttree in config file. if no input, it will produced by this program
std::string saveAnalysisTreeName = config.read<std::string>("saveAnalysisTreeName","0"); // look for filename to save analysis ttree in config file. if no filename, analysis tree won't be saved
bool saveAnalysisTree = config.read<bool>("saveAnalysisTree",0); // choice to save or not the analysis TTree, in a file (name chosen above)
bool calcDoiResWithCalibration = config.read<bool>("calcDoiResWithCalibration",0);
std::string calibrationFileName = config.read<std::string>("calibrationFileName","0");
bool usingAllChannelsForEnergySpectra = config.read<bool>("usingAllChannelsForEnergySpectra",0); // whether to use all digitizer channels for energy spectra (1) or just neighbours (0)
bool usingAllChannels = config.read<bool>("usingAllChannels",0);
bool wAllChannels = config.read<bool>("wAllChannels",0); // whether we use the sum of all channels to compute w (true = 1) of just the neighbours (false = 0). Deafult to false.
int taggingCrystalBins = config.read<int>("taggingCrystalBins",1200);
float taggingSpectrumMin = config.read<float>("taggingSpectrumMin",0.0);
float taggingSpectrumMax = config.read<float>("taggingSpectrumMax",12000.0);
bool TagEdgeCalculation = config.read<bool>("tagEdgeCalculation",0);
std::string TagEdgeCalculationLabels_s = config.read<std::string>("tagEdgeCalculationLabels","");
int digitizerType = config.read<int>("digitizerType",0); // type of digitizer. 0 = desktop, 1 = vme
int amplitudeData = config.read<int>("amplitudeData",0);
bool cuttingOnTagPhotopeak = config.read<bool>("cuttingOnTagPhotopeak",1);
int CTRbins = config.read<int>("CTRbins",500);
float CTRmin = config.read<float>("CTRmin",-5e-9);
float CTRmax = config.read<float>("CTRmax",5e-9);
float histo3Dmin = config.read<float>("histo3Dmin",0);
float histo3Dmax = config.read<float>("histo3Dmax",1);
int DeltaTimeBins = config.read<int>("DeltaTimeBins",500);
float DeltaTimeMin = config.read<float>("DeltaTimeMin",-5e-9);
float DeltaTimeMax = config.read<float>("DeltaTimeMax",5e-9);
bool taggingForTiming = config.read<bool>("taggingForTiming",0);
float photopeakSigmasMin = config.read<float>("photopeakSigmasMin",2.0); // how many sigmas far from mean is the lower bound of cut on photopeak - default = 2.0
float photopeakSigmasMax = config.read<float>("photopeakSigmasMax",4.0); // how many sigmas far from mean is the upper bound of cut on photopeak - default = 4.0
float TaggingPhotopeakSigmasMin = config.read<float>("TaggingPhotopeakSigmasMin",1.5); // how many sigmas far from mean is the lower bound of cut on Tagging photopeak - default = 1.5
float TaggingPhotopeakSigmasMax = config.read<float>("TaggingPhotopeakSigmasMax",2.0); // how many sigmas far from mean is the upper bound of cut on Tagging photopeak - default = 2.0
int WrangeBinsForTiming = config.read<int>("WrangeBinsForTiming",10);
bool smearTaggingTime = config.read<bool>("smearTaggingTime",0);// whether to smear the time stamp of external tagging. Needed for simulations, where the tagging time stamp is always 0 (i.e. the gamma emission time) - default = 0
float tagFitLowerFraction = config.read<float>("tagFitLowerFraction",0.06); // enRes = 2.355*sigma/peak --> sigma = enRes*Peak/2.355 EnRes about 0.15 --> sigma = 0.06* peak -> limits -1sigma +2 sigma
float tagFitUpperFraction = config.read<float>("tagFitUpperFraction",0.12);
bool timingCorrection = config.read<bool>("timingCorrection",1);// perform or not the timing correction (for example, it could make no sense to perform it fully in polished arrays) - deafult = 1
bool timingCorrectionForPolished = config.read<bool>("timingCorrectionForPolished",1);// produce the plots for the simple combination of timestamps that can be done for polished crystals (i.e. when DOI info is not available)
int adcChannels = config.read<int>("digitizerTotalCh");// number of output channels in the adc in use (32 for the CAEN DT5740)
float noiseSigmas = config.read<float>("noiseSigmas",1.0); // how many sigmas of noise far from "0" to cut the dataset
int taggingCrystalTimingChannel = config.read<int>("taggingCrystalTimingChannel",16); // input timing channel where the tagging crystal information is stored - default = 16
float marginWZgraph = config.read<float>("marginWZgraph",0.1); // distance from crystal limit for calculation of beginW and endW
float tagCrystalPeakResolutionFWHM = config.read<float>("tagCrystalPeakResolutionFWHM",0.07);
float photopeakFitRangeMin = config.read<float>("photopeakFitRangeMin",1.2);
float photopeakFitRangeMax = config.read<float>("photopeakFitRangeMax",1.2);
int TimeCorrectionFitFunction = config.read<int>("TimeCorrectionFitFunction",0); // function to be used on time correction slice fitting. 0 = crystalball, 1 = gauss+exp
bool applyNoiseCut = config.read<bool>("applyNoiseCut",1); // apply or not the noise cut - default = 1 (true)
bool applyMaxAllowedCut = config.read<bool>("applyMaxAllowedCut",1); // apply or not the cut on max allowed charge, to avoid overflows - default = 1 (true)
int noiseCutLevel = config.read<int>("noiseCutLevel",0); // number of channels that can be equal to 0 for an event to be accepted. So noiseCutLevel = 0 means ALL channels have to be different to 0 for an event to be accepted, 1 means that 1 channel can be 0, etc...
bool apply3Dcut = config.read<bool>("apply3Dcut",1); // apply or not the 3D cut - default = 1 (true)
// bool likelihoodCorrection = config.read<bool>("likelihoodCorrection",0); // perform likelihood correction
// bool noTimeFitting = = config.read<bool>("noTimeFitting",0); // avoid fitting of time function, just take
// bool calcDOIresInternally = config.read<bool>("calcDOIresInternally",0); // calculate DOi res running on events and computing distance to calibrationGraph
bool essentialOnly = config.read<bool>("essentialOnly",1); // do only essential plots, to minimize computation time - default = 1 (true)
float spectrumSearchMin = config.read<float>("spectrumSearchMin",1);
float spectrumSearchMax = config.read<float>("spectrumSearchMax",histo1Dmax);
int neighCTRbins = config.read<int>("neighCTRbins",500); // number of bins in neighbour CTR plots - default = 500
float neighCTRmin = config.read<float>("neighCTRmin",-1e-9); // min of neighbour CTR plots - default = -1e-9
float neighCTRmax = config.read<float>("neighCTRmax",10e-9); // max of neighbour CTR plots - default = 10e-9
bool lowStat = config.read<bool>("lowStat",1); // if low statistics, apply fits to slices - default 1 (true)
// bool gexp = config.read<bool>("gexp",0); // use the EMG gaussian fit - default 0 (false)
bool dumpImages = config.read<bool>("dumpImages",0); // dump images of some plots for quick check - default 0 (false)
bool applyTriggerCut = config.read<bool>("applyTriggerCut",1); // apply trigger cut (channel in analysis has biggest charge in the event)
// FIXME hardcoded
bool manualCutG = config.read<bool>("manualCutG",0); // read the cutg from external file - default = 0 (false). FIXME this is HIGHLY hardcode fo the moment, will work just on 1-to-1 coupling crystal mppc and fo rjust 1 crystal per ModuleCalibration run!!!
std::string cutgsFileName = config.read<std::string> ("cutgsFileName","cutgs.root"); // name of root file with manual cutgs - default = cutgs.root
float minTimestamp = config.read<float>("minTimestamp",-70e-9);
float maxTimestamp = config.read<float>("maxTimestamp",-40e-9);
bool dumpSinglePeaks = config.read<bool>("dumpSinglePeaks",0); // dump a string with single peaks positions, to use to equalize 511 kev positons (hence gains) in following ModuleCalibration run
// it works like this.
// 1. run standard modulecalibration without manualCutG option
// 2. maybe no cluster found, or you don't like the cluster
// 3. take a look at the 2d projections of Flood Histogram 3D for the mppc under analysis (the _zx and _zy ones)
// 4. create a graphical cut for one of them
// a) open toolbar in gui
// b) click on scissors symbol (Graphical Cut)
// c) draw a cut by left clicking on points, then close it with double click (cut has to be a closed line)
// d) retrieve from ROOT "world" and change name
// TCutG *mycutg1;
// mycutg1 = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
// mycutg1->SetName("cutg_0");
// e) do the same as c) for the other plot
// f) same as d) for this plot
// TCutG *mycutg2;
// mycutg2 = (TCutG*)gROOT->GetListOfSpecials()->FindObject("CUTG");
// mycutg2->SetName("cutg_1");
// g) save both in a tfile (cutgs.root)
// TFile *fOut = new TFile("cutgs.root","RECREATE")
// fOut->cd()
// mycutg1->Write()
// mycutg2->Write()
// fOut->Close()
// h) re-run ModuleCalibration by adding the manualCutG = 1 key
//
std::ofstream singlePeaksFile;
std::map<std::string, float> mapGain;
std::string mppc_s; // input string of the mppc key from config file
std::vector <std::string> mppc_f; // tokenization of above strings
std::vector <std::string> mppc_label; // above tokenized string transformed in the proper variable types
mppc_s = config.read<std::string>("mppc");
config.split( mppc_f, mppc_s, "," );
for(unsigned int i = 0 ; i < mppc_f.size() ; i++)
{
config.trim(mppc_f[i]);
mppc_label.push_back(mppc_f[i]);
}
if(dumpSinglePeaks)
{
singlePeaksFile.open("singlePeaksFile.txt", std::ofstream::out);
}
// channels to exclude from time correction (will affect only polished correction)
std::string excludeChannels_s = config.read<std::string>("excludeChannels",""); //channels to exclude from time correction (will affect only polished correction, the others have to be specified in timeAnalysis)
std::vector<std::string> excludeChannels_f;
std::vector<int> excludeChannels;
config.split( excludeChannels_f, excludeChannels_s, "," );
for(unsigned int i = 0 ; i < excludeChannels_f.size() ; i++)
{
config.trim(excludeChannels_f[i]);
excludeChannels.push_back(atoi(excludeChannels_f[i].c_str()));
}
std::vector<std::string> tagEdgeCalculationLabels_f;
std::vector<std::string> tagEdgeCalculationLabels;
config.split( tagEdgeCalculationLabels_f, TagEdgeCalculationLabels_s, "," );
for(unsigned int i = 0 ; i < tagEdgeCalculationLabels_f.size() ; i++)
{
config.trim(tagEdgeCalculationLabels_f[i]);
tagEdgeCalculationLabels.push_back(tagEdgeCalculationLabels_f[i].c_str());
}
for(unsigned int i = 0 ; i <tagEdgeCalculationLabels.size(); i++ )
{
std::cout << tagEdgeCalculationLabels[i] << " " ;
}
std::cout << std::endl;
//get a more readable variable for essential plots
bool completeAnalysis = true; // when this is true, all plots are done and saved. When it's false, only essentials
if(essentialOnly)
{
completeAnalysis = false;
}
std::vector<tagEdge_t> tagEdge;
// GET INFO ON PC AND FILES
char hostname[1024];
hostname[1023] = '\0';
gethostname(hostname, 1023);
std::string HostNameString(hostname);
char cwd[1024];
getcwd(cwd, sizeof(cwd));
std::string PWDstring(cwd);
// std::cout << cwd << std::endl;
ROOT::v5::TFormula::SetMaxima(10000,10000,10000);
// list of all input files
std::stringstream InputFiles;
if(std::string(argv[1]) == std::string("-c")) // first argument is -c, then the config file name is passed by command line
{
for (int i = 3; i < argc ; i++) // run on the remaining arguments to add all the input files
{
// std::cout << "Adding file " << argv[i] << std::endl;
InputFiles << argv[i] << std::endl;
}
}
else // the config file was indeed the default one
{
for (int i = 1; i < argc ; i++) // run on the remaining arguments to add all the input files
{
// std::cout << "Adding file " << argv[i] << std::endl;
InputFiles << argv[i] << std::endl;
}
}
// std::cout << InputFiles.str() << std::endl;
//set outputFileName
// change name if parallel keys are given
std::string parallelOutput = config.read<std::string>("parallelOutput","0");
if (parallelOutput.compare("0") != 0)
{
outputFileName = parallelOutput;
}
int parallelCrystal = config.read<int>("parallelCrystal",-1);
if (parallelCrystal != -1)
{
outputFileName = parallelOutput;
}
outputFileName += ".root";
TFile *calibrationFile = NULL;
if(calcDoiResWithCalibration)
{
calibrationFile = new TFile(calibrationFileName.c_str());
}
//----------------------------------------------------------//
// Creating the module and elements //
//----------------------------------------------------------//
// create the elements of the module
Module*** module;
Mppc*** mppc;
Crystal*** crystal;
// make an array of module pointers
module = new Module**[nmodulex];
for(int j = 0; j < nmodulex ; j++) module[j] = new Module* [nmoduley];
// make an array of mppc pointers
mppc = new Mppc**[nmodulex*nmppcx];
for(int j = 0; j < nmodulex*nmppcx ; j++) mppc[j] = new Mppc* [nmoduley*nmppcy];
// make an array of crystal pointers
crystal = new Crystal**[nmodulex*nmppcx*ncrystalsx];
for(int j = 0; j < nmodulex*nmppcx*ncrystalsx ; j++) crystal[j] = new Crystal* [nmoduley*nmppcy*ncrystalsy];
// fill the elements
input.FillElements(module,mppc,crystal);
//----------------------------------------------------------//
//----------------------------------------------------------//
// Plots and spectra //
//----------------------------------------------------------//
//properly set timing correction flags
// the problem is back compatibility
// timingCorrection is used to enable the complete time correciton analysis
// timingCorrectionForPolished enables only the correction for polished crystals
// but the second is part of the first, so the second controls and if statement that
// would not allow the first to happen. Therefore, if the first is set on but the second is off
// the second needs to be enabled too
// if(timingCorrection == true && timingCorrectionForPolished == false)
// {
// timingCorrectionForPolished = true;
// } // -- FIXME
// doi bench specific part
std::ofstream doiFile;
std::ofstream AltDoiFile;
std::ifstream altDoiCalcFile;
// saturation dataset part
std::ofstream saturationFile;
if(usingTaggingBench && calcDoiResWithDelta)
{
if( (calcDoiFileName.compare("") == 0) )
{
std::cout << "No file with doi Calibration provided!" << std::endl;
calcDoiResWithDelta = false;
}
if( pointsFromDoi == 0)
{
std::cout << "Need to specify the number of points in doi calibration file!" << std::endl;
calcDoiResWithDelta = false;
}
}
//saturation run part
std::vector<SaturationPeak_t> saturationPeak;
if(saturationRun)
{
saturationFile.open("saturationData.txt", std::ofstream::out);
if(histoSingleChargeMax == 0)// calc histoSingleChargeMax from the histo1Dmax,
{
histoSingleChargeMax = histo1Dmax * chargeBinningADC / 2.0 ; // the /2.0 it's because the histo1D is sum spectrum
}
if(histoSingleChargeBin == 0)// set the histoSingleChargeBin to histo1Dbins
{
histoSingleChargeBin = histo1Dbins;
}
if(histoSumChargeMax == 0)// set it to the same as single
{
histoSumChargeMax = histoSingleChargeMax ; // the /2.0 it's because the histo1D is sum spectrum
}
if(histoSumChargeBin == 0)// set it to the same as single
{
histoSumChargeBin = histoSingleChargeBin;
}
//read the config file to set the search areas for peaks
// std::string saturationPeakEnergy_s,saturationPeakMin_s,saturationPeakMax_s;
std::vector <std::string> saturationPeakEnergy_f,saturationPeakMin_f,saturationPeakMax_f;
std::vector <float> saturationPeakEnergy,saturationPeakMin,saturationPeakMax;
config.split( saturationPeakEnergy_f, saturationPeakEnergy_s, "," );
config.split( saturationPeakMin_f, saturationPeakMin_s, "," );
config.split( saturationPeakMax_f, saturationPeakMax_s, "," );
for(unsigned int i = 0 ; i < saturationPeakEnergy_f.size() ; i++)
{
config.trim(saturationPeakEnergy_f[i]);
saturationPeakEnergy.push_back(atof(saturationPeakEnergy_f[i].c_str()));
}
for(unsigned int i = 0 ; i < saturationPeakMin_f.size() ; i++)
{
config.trim(saturationPeakMin_f[i]);
saturationPeakMin.push_back(atof(saturationPeakMin_f[i].c_str()));
}
for(unsigned int i = 0 ; i < saturationPeakMax_f.size() ; i++)
{
config.trim(saturationPeakMax_f[i]);
saturationPeakMax.push_back(atof(saturationPeakMax_f[i].c_str()));
}
// store data in saturationPeak vector of struct
for(unsigned int i = 0 ; i < saturationPeakEnergy.size() ; i++)
{
SaturationPeak_t tempPeak;
tempPeak.energy = saturationPeakEnergy[i];
tempPeak.peakMin = saturationPeakMin[i];
tempPeak.peakMax = saturationPeakMax[i];
saturationPeak.push_back(tempPeak);
}
}
//points form DOI input file
std::vector<inputDoi_t> inputDoi;
inputDoi_t tempInputDoi(pointsFromDoi);
if(usingTaggingBench)
{
doiFile.open("doiData.txt", std::ofstream::out);
if(calcDoiResWithDelta)
{
altDoiCalcFile.open(calcDoiFileName.c_str(), std::ios::in);
if(!altDoiCalcFile)
{
std::cout << "Doi calibration file not found!" << std::endl;
calcDoiResWithDelta = false;
}
else //finally, everything is fine...
{
AltDoiFile.open("altDoiRes.txt", std::ofstream::out);
while(altDoiCalcFile >> tempInputDoi)
{
inputDoi.push_back(tempInputDoi);
tempInputDoi.clear();
}
}
}
}
// simulation dataset spectific part
TCut SingleCrystalInteraction;
TCut SingleEnergyDeposition;
if(usingRealSimData) // only if this is a sim dataset
{
SingleCrystalInteraction = "CrystalsHit == 1";
SingleEnergyDeposition = "NumbOfInteractions == 1";
}
//variable that will be used only if the user requires the TagEdgeCalculation
float tagPeakHgEntries = 0.0;
// MAIN LOOP v2
// get the TTree, to plot the spectra
// TTree* tree = input.GetTree();
// fill the tchain with input files
//----------------------------------------------------------//
// Get TChain of input files //
//----------------------------------------------------------//
TChain* tree = new TChain(fname.c_str()); // create the input tchain and the analysis ttree
// first, create the adc channels variables and branches
// ChainAdcChannel = new Int_t [adcChannels];
ChainDesktopAdcChannel = new Short_t [adcChannels]; // input from ADC desktop
ChainVMEadcChannel = new UShort_t [adcChannels]; // input from VME
// ChainVMEadcChannel = new Float_t [adcChannels]; // input from VME
ChainVMEAmplitudeChannel = new UShort_t [adcChannels]; // input amplitude from VME
ChainTimeStamp = new Float_t[adcChannels];
// ChainTimeStamp = new Long_t[adcChannels];
// TDCBinning = new Float_t[adcChannels];
// DigitizerChannelOn = new bool[adcChannels];
bChainAdcChannel = new TBranch* [adcChannels];
bChainTimeStamp = new TBranch* [adcChannels];
bChainAmplChannel = new TBranch* [adcChannels];
if(std::string(argv[1]) == std::string("-c")) // first argument is -c, then the config file name is passed by command line
{
for (int i = 3; i < argc ; i++) // run on the remaining arguments to add all the input files
{
std::cout << "Adding file " << argv[i] << std::endl;
tree->Add(argv[i]);
}
}
else // the config file was indeed the default one
{
for (int i = 1; i < argc ; i++) // run on the remaining arguments to add all the input files
{
std::cout << "Adding file " << argv[i] << std::endl;
tree->Add(argv[i]);
}
}
// set branches for reading the input files
tree->SetBranchAddress("ExtendedTimeTag", &ChainExtendedTimeTag, &bChainExtendedTimeTag);
tree->SetBranchAddress("DeltaTimeTag", &ChainDeltaTimeTag, &bChainDeltaTimeTag);
if(usingRealSimData)
{
tree->SetBranchAddress("RealX", &RealX, &bRealX);
tree->SetBranchAddress("RealY", &RealY, &bRealY);
tree->SetBranchAddress("RealZ", &RealZ, &bRealZ);
// fchain->SetBranchAddress("Tagging", &simTaggingCharge, &bsimTaggingCharge);
// fchain->SetBranchAddress("TaggingTimeStamp", &simTaggingTime, &bsimTaggingTime);
tree->SetBranchAddress("CrystalsHit",&CrystalsHit, &bCrystalsHit);
tree->SetBranchAddress("NumbOfInteractions",&NumbOfInteractions, &bNumbOfInteractions);
// fchain->SetBranchAddress("TotalCryEnergy",&TotalCryEnergy, &bTotalCryEnergy);
}
for(int i=0; i<adcChannels; i++)
{
if(digitizerType == 0)
{
std::stringstream sname;
sname << "ch" << i;
tree->SetBranchAddress(sname.str().c_str(), &ChainDesktopAdcChannel[i], &bChainAdcChannel[i]);
sname.str("");
}
if(digitizerType == 1)
{
std::stringstream sname;
sname << "ch" << i;
tree->SetBranchAddress(sname.str().c_str(), &ChainVMEadcChannel[i], &bChainAdcChannel[i]);
sname.str("");
sname << "t" << i;
tree->SetBranchAddress(sname.str().c_str(), &ChainTimeStamp[i],&bChainTimeStamp[i]);
sname.str("");
if(amplitudeData == 1)
{
sname << "ampl" << i;
tree->SetBranchAddress(sname.str().c_str(), &ChainVMEAmplitudeChannel[i], &bChainAmplChannel[i]);
sname.str("");
}
}
}
TH1F* TaggingCrystalSpectrum = NULL;
TH1F *TriggerSpectrumHighlight = NULL;
TList formulas;
std::vector<detector_t> detector;
// Loop on modules, mppcs and crystal
for(int iModule = 0; iModule < nmodulex ; iModule++)
{
for(int jModule = 0; jModule < nmoduley ; jModule++)
{
// different variables logic, trying to use directly the original files:
// almost everything need to be produced channel by channel (otherwise the
// definition of neighbours is impossible)
// so the variables in the original files are just:
//
// - BoardTimeTags ... ch0 ... chN , t0 ... tN
//
// and with the pre-processing we were genetating some additional variables:
//
// - TriggerChannel = the channel with highest charge among all channels in the event
// - FloodX,Y = weighted average on charge of detector positions, with multiple coices (i.e. summing on all detectors or only on trigger channel + 8 neighbours)
// - FloodZ = ratio (charge in trigger channel)/(sum charge) where sum chage can be the sum of all channels or the sum of trigger channel + 8 neighbours
// - BadEvent = counter of events with charge higher of saturation parameter for that specific channel. these events are "bad" because given the logarithmic formula of saturation correction, it would be mathematically impossible to correct them. It's just a counter that provides a feedback to user (typically if all makes sense, badEvents need to be just a small fraction of all events)
// - Tagging = charge measured in the external tagging crystal, if any
// - TaggingTimeStamp = time stamp of the external tagging crystal, if any
// - ZPosition = z position of the tagging bench, for DOI measurements, for this event
//
// These variables can be computed on the fly directly on the original data (with some limitations, see later)
//
// write down the channel collection
// std::stringstream TriggerChannel;
// TriggerChannel
//useful strings etc
// First cut is a general, global cut
// it will be composed of all cuts that will apply generally to all the events
// For example
// 1. A global XYZ cut that allows us to consider only the events with "reasonable" uvw coordinates
// 2. A cut on photopeak of the tagging crystal, if there is any such crystal involved
std::cout << "Generating global spectra..." << std::endl;
std::stringstream var;
std::stringstream sname;
std::stringstream cut;
TCut BasicCut;
//FIXME what do we do with the "old" CutXYZ? <--------------------
// std::stringstream CutXYZstream;
// CutXYZstream << "FloodX > " << -moduleLateralSideX << " && FloodX < " << moduleLateralSideX << "&& FloodY > " << -moduleLateralSideY << " && FloodY < " << moduleLateralSideY << "&& FloodZ > " << histo3Dmin << " && FloodZ < " << histo3Dmax;
// TCut CutXYZ = CutXYZstream.str().c_str();
// get a std::vector of all module channels
// std::vector<int> moduleChannels = module[iModule][jModule]->GetChannels();
//get the detector data
detector = module[iModule][jModule]->GetDetector();
// // //DEBUG OUTPUT
// // output the detector[] data for check
// // struct detector_t
// // // {
// // // int digitizerChannel;
// // // int timingChannel;
// // // std::string label;
// // // int i;
// // // int j;
// // // float saturation;
// // // int plotPosition;
// // // float xPosition;
// // // float yPosition;
// // // float pedestal;
// // // int OnForDOI;
// // // std::vector<int> neighbourChannels;
// // // bool OnForModular;
// // // // bool operator<(const masks_t& rhs) const { meanx < rhs.meanx; }
// // // };
// std::cout << "////////////////////////////////////////////////" << std::endl;
// std::cout << "// Detector Structure //" << std::endl;
// std::cout << "////////////////////////////////////////////////" << std::endl;
// for(unsigned int iDet = 0; iDet < detector.size() ; iDet++)
// {
// std::cout << "------------------------------------------------" << std::endl;
// std::cout << "Index = " << iDet << std::endl;
// std::cout << "digitizerChannel = " << detector[iDet].digitizerChannel << std::endl;
// std::cout << "timingChannel = " << detector[iDet].timingChannel << std::endl;
// std::cout << "label = " << detector[iDet].label << std::endl;
// std::cout << "i = " << detector[iDet].i << std::endl;
// std::cout << "j = " << detector[iDet].j << std::endl;
// std::cout << "saturation = " << detector[iDet].saturation << std::endl;
// std::cout << "plotPosition = " << detector[iDet].plotPosition << std::endl;
// std::cout << "xPosition = " << detector[iDet].xPosition << std::endl;
// std::cout << "yPosition = " << detector[iDet].yPosition << std::endl;
// std::cout << "pedestal = " << detector[iDet].pedestal << std::endl;
// std::cout << "OnForDOI = " << detector[iDet].OnForDOI << std::endl;
// std::cout << "Neighbours = ";
// for(unsigned int iNeigh = 0; iNeigh < detector[iDet].neighbourChannels.size(); iNeigh++)
// {
// std::cout << detector[iDet].neighbourChannels[iNeigh] << " " ;
// }
// std::cout << std::endl;
// std::cout << "OnForModular = " << detector[iDet].OnForModular << std::endl;
// std::cout << std::endl;
// }
// std::cout << "////////////////////////////////////////////////" << std::endl;
// std::cout << "// End of Detector Structure //" << std::endl;
// std::cout << "////////////////////////////////////////////////" << std::endl;
// // //END of DEBUG OUTPUT
std::vector<int> allModuleChID;
std::vector<float> saturationData;
std::vector<float> pedestalData;
std::vector<float> gainData;
std::vector<int> detChData;
for(unsigned int iDet = 0; iDet < detector.size() ; iDet++)
{
allModuleChID.push_back(iDet); //dummy... we need to get all the indexes of detector[] that identify all the channels in the array, and of course these are all the indexes of the detector[] vector.
detChData.push_back(detector[iDet].digitizerChannel);
saturationData.push_back(detector[iDet].saturation);
pedestalData.push_back(detector[iDet].pedestal);
gainData.push_back(detector[iDet].gain);
}
module[iModule][jModule]->SetDetChData(detChData);
module[iModule][jModule]->SetSaturationData(saturationData);
module[iModule][jModule]->SetPedestalData(pedestalData);
module[iModule][jModule]->SetGainData(gainData);
// //DEBUG OUTPUT
// std::cout << std::endl;
// for(unsigned int iDet = 0; iDet < moduleChannels.size() ; iDet++)
// {
// std::cout << moduleChannels[iDet] << " ";
// }
// std::cout << std::endl;
// BasicCut += CutXYZ; // add the module physical constrains the to global cut of accepted events
//
TCut taggingPhotopeakCut = "" ;
// int Tagging =
// prepare GLOBAL SPECTRA
// Flood histogram
// TString nameModule;
// std::stringstream varModule;
sname << "Flood Histogram 2D - " << module[iModule][jModule]->GetName();
// varModule << "FloodY:FloodX >> " << nameModule;
TH2F *spectrum2dModule = new TH2F(sname.str().c_str(),sname.str().c_str(),histo2DglobalBins,-moduleLateralSideX,moduleLateralSideX,histo2DglobalBins,-moduleLateralSideY,moduleLateralSideY);
// tree->Draw(varModule.str().c_str(),"","COLZ");
spectrum2dModule->SetName(sname.str().c_str());
spectrum2dModule->SetTitle(sname.str().c_str());
spectrum2dModule->GetXaxis()->SetTitle("U");
spectrum2dModule->GetYaxis()->SetTitle("V");
// module[iModule][jModule]->SetFloodMap2D(spectrum2dModule);
sname.str("");
if(usingTaggingBench || calcDoiResWithCalibration || taggingForTiming)//trigger spectrum
{
if(taggingPeakMax == 0)
{
taggingPeakMax = taggingSpectrumMax;
}
if(taggingPeakMin == 0)
{
taggingPeakMin = taggingSpectrumMin;
}
std::stringstream tagStream;
tagStream << "ch" << taggingCrystalChannel; // if using original files
TaggingCrystalSpectrum = new TH1F("TaggingCrystalSpectrum","TaggingCrystalSpectrum",taggingCrystalBins,taggingSpectrumMin,taggingSpectrumMax);
var << tagStream.str().c_str() << ">> TaggingCrystalSpectrum";
// std::cout << nameModule << " ... ";
TaggingCrystalSpectrum->SetName("TaggingCrystalSpectrum");
TaggingCrystalSpectrum->GetXaxis()->SetTitle("ADC");
TaggingCrystalSpectrum->GetYaxis()->SetTitle("Counts");
TaggingCrystalSpectrum->SetTitle("Spectrum of Tagging Crystal");
tree->Draw(var.str().c_str(),"");
var.str("");
//restrict the region where to look for peaks. Fix for tspectrum...
TaggingCrystalSpectrum->GetXaxis()->SetRangeUser(taggingPeakMin,taggingPeakMax);
//find peak in the tagging crystal
TSpectrum *sTagCrystal;
sTagCrystal = new TSpectrum(5);
Int_t TagCrystalPeaksN = sTagCrystal->Search(TaggingCrystalSpectrum,1,"",0.5);
Double_t *TagCrystalPeaks = sTagCrystal->GetPositionX();
Double_t *TagCrystalPeaksY = sTagCrystal->GetPositionY();
// float saturationPeakFraction
//delete s;
// float distPeak = INFINITY;
float maxPeak = 0.0;
int peakID = 0;
for (int peakCounter = 0 ; peakCounter < TagCrystalPeaksN ; peakCounter++ )
{
// //DEBUG
// std::cout << peakCounter << "\t" << TagCrystalPeaks[peakCounter] << "\t" << TagCrystalPeaksY[peakCounter] << std::endl;
// if( fabs(CrystalPeaks[peakCounter] - 0.5*(saturationPeak[iSaturation].peakMin+saturationPeak[iSaturation].peakMax)) < distPeak)//take closest peak to the center of search range selected by the user
if(TagCrystalPeaksY[peakCounter] > maxPeak)
{
// distPeak = CrystalPeaks[peakCounter];
maxPeak = TagCrystalPeaksY[peakCounter];
peakID = peakCounter;
}
}
// //DEBUG
// std::cout << "chosen " << peakID << "\t" << TagCrystalPeaks[peakID] << "\t" << TagCrystalPeaksY[peakID] << std::endl;
TF1 *gaussTag = new TF1("gaussTag",
"gaus",
TagCrystalPeaks[peakID] - tagFitLowerFraction*TagCrystalPeaks[peakID],
TagCrystalPeaks[peakID] + tagFitUpperFraction*TagCrystalPeaks[peakID]);
// TF1 *gaussTag = new TF1("gaussTag",
// "gaus",
// taggingPeakMin,
// taggingPeakMax);
gaussTag->SetParameter(0,TagCrystalPeaksY[peakID]); // heigth
gaussTag->SetParameter(1,TagCrystalPeaks[peakID]); //mean
gaussTag->SetParameter(2,(tagCrystalPeakResolutionFWHM*TagCrystalPeaks[peakID])/2.355); // sigma
//
// TaggingCrystalSpectrum->Fit("gaussTag",
// "Q",
// "",
// taggingPeakMin,
// taggingPeakMax);
TaggingCrystalSpectrum->Fit("gaussTag",
"Q",
"",
TagCrystalPeaks[peakID] - tagFitLowerFraction*TagCrystalPeaks[peakID],
TagCrystalPeaks[peakID] + tagFitUpperFraction*TagCrystalPeaks[peakID]);
// // DEBUG
// std::cout << "fit pars " << gaussTag->GetParameter(1) << "\t" << gaussTag->GetParameter(2) << std::endl;
TaggingCrystalSpectrum->GetXaxis()->SetRangeUser(taggingSpectrumMin,taggingSpectrumMax);
//define a TCut for this peak
float tagPhotopeakMin = gaussTag->GetParameter(1) - TaggingPhotopeakSigmasMin*gaussTag->GetParameter(2);
float tagPhotopeakMax = gaussTag->GetParameter(1) + TaggingPhotopeakSigmasMax*gaussTag->GetParameter(2);
std::stringstream tagString;
tagString << tagStream.str().c_str() << " > " << tagPhotopeakMin << "&& " << tagStream.str().c_str() <<" < " << tagPhotopeakMax;
taggingPhotopeakCut = tagString.str().c_str();
//highlighted spectrum
TriggerSpectrumHighlight = new TH1F("TriggerSpectrumHighlight","",taggingCrystalBins,taggingSpectrumMin,taggingSpectrumMax);
var.str("");
var << tagStream.str().c_str() << " >> TriggerSpectrumHighlight";
TriggerSpectrumHighlight->SetLineColor(3);
TriggerSpectrumHighlight->SetFillColor(3);
TriggerSpectrumHighlight->SetFillStyle(3001);
tree->Draw(var.str().c_str(),taggingPhotopeakCut);
var.str("");
if(TagEdgeCalculation)
{
tagPeakHgEntries = TriggerSpectrumHighlight->Integral();
}
if(cuttingOnTagPhotopeak) // only if tagging on photopeak is chosen
{
BasicCut += taggingPhotopeakCut; // add the cut on photopeak of the tagging crystal.
}
//but always save the tagging cut in the output file
taggingPhotopeakCut.SetName("taggingPhotopeakCut");
module[iModule][jModule]->SetTaggingPhotopeakCut(taggingPhotopeakCut); // energy cut, events in the photopeak
TTreeFormula* FormulaTag = new TTreeFormula("FormulaTag",taggingPhotopeakCut,tree);
formulas.Add(FormulaTag);
module[iModule][jModule]->SetFormulaTaggingPhotopeakCut(FormulaTag); // energy cut, events in the photopeak
// delete sTagCrystal;
delete gaussTag;
// std::cout << " done" << std::endl;