31 int main(
int argc,
char **argv)
39 bool setDefaultLimits;
48 JParser<> zap(
"Auxiliary program to fit singles rate distributions.");
61 catch(
const exception &error) {
62 FATAL(error.what() << endl);
66 using namespace KM3NETDAQ;
77 TFile*
in = TFile::Open(inputFile.c_str(),
"exist");
79 if (
in == NULL || !
in->IsOpen()) {
80 FATAL(
"File: " << inputFile <<
" not opened." << endl);
85 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
87 DEBUG(
"Module " << module->getID() << endl);
90 TH2D* h2s = (TH2D*)
in->Get((title+
".2S").c_str());
94 WARNING(
"No data in histogram " << title << endl);
99 const double factor = 1.0/1000 ;
104 for (
int i = 1; i <= h2s->GetXaxis()->GetNbins(); ++i) {
106 TH1D slice((title+Form(
".%i.1S", i-1)).c_str(), NULL, h2s->GetYaxis()->GetNbins(), JDAQRate::getData(factor)) ;
110 for (
int j = 1 ;
j <= h2s->GetNbinsY() ; ++
j) {
112 slice.SetBinContent(
j, h2s->GetBinContent(i,
j)) ;
113 slice.SetBinError (
j, sqrt(h2s->GetBinContent(i,
j))) ;
117 if (slice.Integral() <= 0.0) {
119 WARNING(
"No data in PMT " << i-1 <<
" of module " << title <<
", skipping fit" << endl) ;
124 slice.Scale(1.0/slice.Integral()) ;
127 int binmax = slice.GetMaximumBin() ;
128 double ratemax = slice.GetBinContent(binmax) ;
132 while ((slice.GetBinContent(bin_l)==0 || slice.GetBinContent(bin_l) >= peakFraction*ratemax) && bin_l > 0) { bin_l-- ; } ; bin_l++ ;
133 while ((slice.GetBinContent(bin_r)==0 || slice.GetBinContent(bin_r) >= peakFraction*ratemax) && bin_r < slice.GetNbinsX()) { bin_r++ ; } ; bin_r-- ;
135 JRange<double> fitrange(slice.GetXaxis()->GetBinCenter(bin_l), slice.GetXaxis()->GetBinCenter(bin_r)) ;
137 TF1
f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(sqrt(2*pi)*[2])");
140 f1.SetParameter(0, 0.2) ;
141 f1.SetParameter(1, 6.0) ;
142 f1.SetParameter(2, 0.3) ;
144 if (setDefaultLimits) {
145 f1.SetParLimits(0, 0.0, 10.0);
146 f1.SetParLimits(1, 0.0, 20.0);
147 f1.SetParLimits(2, 0.0, 1.0);
154 DEBUG(
"Fit histogram " << slice.GetName() <<
" in range " << fitrange <<
" kHz " << endl ) ;
156 slice.Fit(&f1, option.c_str(),
"same", max(0.0, fitrange.first), fitrange.
second) ;
158 DEBUG( f1.GetParameter(0)<<
" "<<f1.GetParameter(1)<<
" "<<f1.GetParameter(2)<<endl ) ;
160 mean.SetBinContent(i, f1.GetParameter(1)) ;
161 sigma.SetBinContent(i, f1.GetParameter(2)) ;
Utility class to parse command line options.
int main(int argc, char *argv[])
Wrapper class around STL string class.
Data structure for detector geometry and calibration.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
General purpose messaging.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Utility class to parse command line options.
do set_variable DETECTOR_TXT $WORKDIR detector
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
#define DEBUG(A)
Message macros.