Jpp  19.1.0-rc.1
the software that should make you happy
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 
16 #include "JPhysics/JConstants.hh"
19 #include "JDetector/JDetector.hh"
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  */
33 int 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:69
#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:2158
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:1714
JKey_t first
Definition: JPair.hh:128
JValue_t second
Definition: JPair.hh:129
const double sigma[]
Definition: JQuadrature.cc:74
const JPolynome f1(1.0, 2.0, 3.0)
Function.
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:792
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Definition: JSTDTypes.hh:14
Detector file.
Definition: JHead.hh:227