Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JFitSinglesRates.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <map>
5
6#include "TROOT.h"
7#include "TFile.h"
8#include "TH1D.h"
9#include "TH2D.h"
10#include "TF1.h"
11#include "TF2.h"
12#include "TMath.h"
13
14#include "JROOT/JMinimizer.hh"
15
21
22#include "JTools/JRange.hh"
23
24#include "Jeep/JParser.hh"
25#include "Jeep/JMessage.hh"
26
27/**
28 * \file
29 *
30 * Auxiliary program to fit singles rate distributions. The fitted parameters are the mean singles rate and the spread.
31 * \author mkarel
32 */
33int main(int argc, char **argv)
34{
35 using namespace std;
36 using namespace JPP;
37
38 string inputFile;
39 string outputFile;
40 string detectorFile;
41 bool setDefaultLimits;
42 double peakFraction;
43 string option;
44 bool writeFits;
45 int debug;
46
47 //Parser module handles input parameters
48 try {
49
50 JParser<> zap("Auxiliary program to fit singles rate distributions.");
51
52 zap['f'] = make_field(inputFile);
53 zap['o'] = make_field(outputFile) = "fit.root";
54 zap['a'] = make_field(detectorFile);
55 zap['L'] = make_field(setDefaultLimits);
56 zap['p'] = make_field(peakFraction) = 0.1;
57 zap['O'] = make_field(option) = "R W M";
58 zap['w'] = make_field(writeFits);
59 zap['d'] = make_field(debug) = 1;
60
61 zap(argc, argv);
62 }
63 catch(const exception &error) {
64 FATAL(error.what() << endl);
65 }
66
67
68 using namespace KM3NETDAQ;
69
71
72 try {
73 load(detectorFile, detector);
74 }
75 catch(const JException& error) {
76 FATAL(error);
77 }
78
79 TFile* in = TFile::Open(inputFile.c_str(), "exist");
80
81 if (in == NULL || !in->IsOpen()) {
82 FATAL("File: " << inputFile << " not opened." << endl);
83 }
84
85 TFile out(outputFile.c_str(), "recreate");
86
87 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
88
89 DEBUG("Module " << module->getID() << endl);
90 const JString title = JString(module->getID());
91
92 TH2D* h2s = (TH2D*) in->Get((title+".2S").c_str());
93
94 if (h2s == NULL) {
95
96 WARNING("No data in histogram " << title << endl);
97
98 continue;
99 }
100
101 const double factor = 1.0/1000 ;
102
103 TH1D mean((title+".1mean").c_str(), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5) ;
104 TH1D sigma((title+".1sigma").c_str(), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5) ;
105
106 for (int i = 1; i <= h2s->GetXaxis()->GetNbins(); ++i) {
107
108 TH1D slice((title+Form(".%i.1S", i-1)).c_str(), NULL, h2s->GetYaxis()->GetNbins(), JDAQRate::getData(factor)) ;
109
110 slice.Sumw2() ;
111
112 for (int j = 1 ; j <= h2s->GetNbinsY() ; ++j) {
113
114 slice.SetBinContent(j, h2s->GetBinContent(i, j)) ;
115 slice.SetBinError (j, sqrt(h2s->GetBinContent(i, j))) ;
116
117 }
118
119 if (slice.Integral() <= 0.0) {
120
121 WARNING("No data in PMT " << i-1 << " of module " << title << ", skipping fit" << endl) ;
122
123 continue ;
124 }
125
126 slice.Scale(1.0/slice.Integral()) ;
127
128 // determine fit range
129 int binmax = slice.GetMaximumBin() ;
130 double ratemax = slice.GetBinContent(binmax) ;
131
132 int bin_l = binmax ;
133 int bin_r = binmax ;
134 while ((slice.GetBinContent(bin_l)==0 || slice.GetBinContent(bin_l) >= peakFraction*ratemax) && bin_l > 0) { bin_l-- ; } ; bin_l++ ;
135 while ((slice.GetBinContent(bin_r)==0 || slice.GetBinContent(bin_r) >= peakFraction*ratemax) && bin_r < slice.GetNbinsX()) { bin_r++ ; } ; bin_r-- ;
136
137 JRange<double> fitrange(slice.GetXaxis()->GetBinCenter(bin_l), slice.GetXaxis()->GetBinCenter(bin_r)) ;
138
139 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(sqrt(2*pi)*[2])");
140
141 // default parameters for start of fit
142 f1.SetParameter(0, 0.2) ;
143 f1.SetParameter(1, 6.0) ; // [kHz]
144 f1.SetParameter(2, 0.3) ; // [kHz]
145
146 if (setDefaultLimits) {
147 f1.SetParLimits(0, 0.0, 10.0);
148 f1.SetParLimits(1, 0.0, 20.0); // [kHz]
149 f1.SetParLimits(2, 0.0, 1.0); // [kHz]
150 }
151
152 if (debug < JEEP::debug_t) {
153 option += " Q";
154 }
155
156 DEBUG("Fit histogram " << slice.GetName() << " in range " << fitrange << " kHz " << endl ) ;
157
158 slice.Fit(&f1, option.c_str(), "same", max(0.0, fitrange.first), fitrange.second) ;
159
160 DEBUG( f1.GetParameter(0)<<" "<<f1.GetParameter(1)<<" "<<f1.GetParameter(2)<<endl ) ;
161
162 mean.SetBinContent(i, f1.GetParameter(1)) ;
163 sigma.SetBinContent(i, f1.GetParameter(2)) ;
164
165 if (writeFits) {
166 slice.Write() ;
167 }
168
169 }
170
171 mean.Write() ;
172 sigma.Write() ;
173
174 }
175
176 out.Close() ;
177
178}
179
180
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define WARNING(A)
Definition JMessage.hh:65
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
Auxiliary class to define a range between two values.
Detector data structure.
Definition JDetector.hh:96
General exception.
Definition JException.hh:24
Wrapper class around STL string class.
Definition JString.hh:29
Utility class to parse command line options.
Definition JParser.hh:1698
JKey_t first
Definition JPair.hh:128
JValue_t second
Definition JPair.hh:129
Range of values.
Definition JRange.hh:42
static const double * getData(const double factor=1.0)
Get abscissa values.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
@ debug_t
debug
Definition JMessage.hh:29
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
int j
Definition JPolint.hh:801
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Detector file.
Definition JHead.hh:227