Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
Functions | Variables
JFitMultiplicity.cc File Reference

Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF2.h"
#include "TMath.h"
#include "TFitResult.h"
#include "JROOT/JMinimizer.hh"
#include "JPhysics/JConstants.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMeta.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JFitK40.hh"
#include "JROOT/JManager.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Variables

const char * livetime_t = "LiveTime"
 

Detailed Description

Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.

Author
mlincett

Definition in file JFitMultiplicity.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 44 of file JFitMultiplicity.cc.

44 {
45 using namespace std;
46 using namespace JPP;
47
48 string inputFile;
49 string outputFile;
50 string summaryFile;
51 string option;
52 int debug;
53 bool fitBackground;
54 bool nominalLiveTime;
55
56 try {
57 JParser<> zap("Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.");
58
59 zap['f'] = make_field(inputFile, "input file");
60 zap['o'] = make_field(outputFile, "output file") = "fitmultiplicity.root";
61 zap['O'] = make_field(option, "fitting option") = "Q";
62 zap['B'] = make_field(fitBackground, "get signal by fit and subtraction from background from the integral");
63 zap['L'] = make_field(nominalLiveTime, "use nominal live time from LiveTime histogram, instead of individual DOMs");
64 zap['d'] = make_field(debug) = 1;
65
66
67 zap(argc, argv);
68 }
69 catch(const exception &error) {
70 FATAL(error.what() << endl);
71 }
72
73 TFile* in = TFile::Open(inputFile.c_str(), "exist");
74
75 if (in == NULL || !in->IsOpen()) {
76 FATAL("File: " << inputFile << " not opened." << endl);
77 }
78
79
80 //----------------------------------------------------------
81 // general histograms
82 //----------------------------------------------------------
83
84 const int nx = 1 + NUMBER_OF_PMTS;
85 const double xmin = -0.5;
86 const double xmax = nx - 0.5;
87
88 typedef JManager<TString, TH1D> JManager_t;
89
90 // inclusive
91 JManager_t HTF(new TH1D("%F", "", nx, xmin, xmax));
92 // exclusive
93 JManager_t HTFE(new TH1D("%FE", "", nx, xmin, xmax));
94
95 HTF->Sumw2(kFALSE);
96
97 //----------------------------------------------------------
98 // get overall normalization constants
99 //----------------------------------------------------------
100
101 TH1D* hLiveTime = (TH1D*) in->Get(livetime_t);
102
103 if (hLiveTime == NULL) {
104 FATAL("File " << inputFile << " does not contain livetime histogram." << endl);
105 }
106
107 //----------------------------------------------------------
108 // initialise output file
109 //----------------------------------------------------------
110
111 TFile out(outputFile.c_str(), "recreate");
112
113 //----------------------------------------------------------
114 // loop over bins of livetime histogram
115 //----------------------------------------------------------
116
117 TDirectory *fdir = out.mkdir("FIT");
118 TDirectory *pdir = out.mkdir("MUL");
119
120 // nominal live time
121 double nT = hLiveTime->GetBinContent(1);
122
123 for (int bin = 2; bin <= hLiveTime->GetNbinsX(); bin++) {
124
125 TString prefix = hLiveTime->GetXaxis()->GetBinLabel(bin);
126
127 DEBUG("Processing data from " << prefix << endl);
128
129 double T = nominalLiveTime ? nT : hLiveTime->GetBinContent(bin);
130
131 TString h1D_label = prefix + TString("_P" );
132 TString h2D_label = prefix + TString("_HT");
133
134 if (T > 0) {
135
136 TH2D* h1D = (TH2D*) in->Get(h1D_label);
137 TH2D* h2D = (TH2D*) in->Get(h2D_label);
138
139 if (h1D != NULL) {
140
141 h1D->Scale(1.0 / T);
142
143 pdir->cd();
144
145 h1D->Write();
146
147 } else {
148 DEBUG(h1D_label << " not found" << endl);
149 }
150
151 if (h2D != NULL) {
152
153 TDirectory* cdir = fdir->mkdir(h2D_label); cdir->cd();
154
155 for (int M = 2; M <= NUMBER_OF_PMTS; M++) {
156
157 int Mbin = h2D->GetXaxis()->FindBin(M);
158
159 TH1D* hI = h2D->ProjectionY(h2D_label + Form("_%d", M), Mbin, NUMBER_OF_PMTS);
160
161 hI->Sumw2();
162
163 if (hI->GetEntries() > 0) {
164
165 const double W = (hI->GetXaxis()->GetXmax() - hI->GetXaxis()->GetXmin());
166 const double N = hI->GetNbinsX();
167 const double V = W / N;
168 const double minRMS = 0.1;
169
170 double inclusiveCountVal = 0, inclusiveCountErr = 0;
171
172 if (hI->GetRMS() > minRMS) {
173
174 hI->Scale(1.0 / V);
175
176 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
177
178 f1.SetParameter(0, hI->GetMaximum());
179 f1.SetParameter(1, hI->GetMean());
180 f1.SetParameter(2, hI->GetRMS() * 0.25);
181 f1.SetParameter(3, hI->GetMinimum());
182
183 hI->Fit(&f1, option.c_str(), "same");
184
185 double sg_v = f1.GetParameter(0);
186 double sg_e = f1.GetParError( 0);
187
188 double bg_v = f1.GetParameter(3);
189 double bg_e = f1.GetParError( 3);
190
191 // 2-sigma acceptance criteria for fitted values
192 if (bg_v / bg_e <= 2.0) { bg_v = bg_e = 0; }
193 if (sg_v / sg_e <= 2.0) { sg_v = sg_e = 0; }
194
195 double integral = hI->GetSumOfWeights();
196
197 if (fitBackground) {
198
199 inclusiveCountVal = V * (integral - (bg_v * N));
200 inclusiveCountErr = V * bg_e * N;
201
202 } else {
203
204 inclusiveCountErr = sg_e;
205 inclusiveCountVal = sg_v;
206
207 }
208
209 hI->Scale(1.0 / T);
210
211 // re-fit after scaling
212 hI->Fit(&f1, option.c_str(), "same");
213
214 }
215
216 if (inclusiveCountVal <= 0) { inclusiveCountVal = 0; }
217
218 // --------------------------------
219 // write fitted histogram
220 // --------------------------------
221
222 hI->Write();
223
224 // --------------------------------
225 // build inclusive rates histograms
226 // --------------------------------
227
228 int bin = HTF[h2D_label]->FindBin(M);
229
230 HTF[h2D_label]->SetBinContent(bin, inclusiveCountVal);
231
232 double binError = sqrt(inclusiveCountVal + pow(inclusiveCountErr, 2));
233
234 HTF[h2D_label]->SetBinError(bin, binError);
235
236 // --------------------------------
237 // build exclusive rates histograms
238 // --------------------------------
239
240 double exclusiveCountVal = 0;
241
242 if (M == NUMBER_OF_PMTS) {
243 exclusiveCountVal = inclusiveCountVal;
244 } else {
245 int prevBin = HTF[h2D_label]->GetXaxis()->FindBin(M+1);
246 exclusiveCountVal = inclusiveCountVal - HTF[h2D_label]->GetBinContent(prevBin);
247 }
248
249 HTFE[h2D_label]->SetBinContent(HTFE[h2D_label]->FindBin(M), exclusiveCountVal);
250
251
252 }
253
254 }
255
256 } else {
257 DEBUG(h1D_label << " not found" << endl);
258 }
259
260 // HTF[hLabel]->Sumw2(kTRUE);
261
262 HTFE[h2D_label]->Sumw2();
263
264 HTF[h2D_label ]->Scale(1.0 / T);
265 HTFE[h2D_label]->Scale(1.0 / T);
266
267 }
268 }
269
270 fdir = out.mkdir("HTF" ); fdir->cd();
271
272 HTF.Write(*fdir);
273
274 fdir = out.mkdir("HTFE"); fdir->cd();
275
276 HTFE.Write(*fdir);
277
278 out.Close();
279
280}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
livetime_t
const char * livetime_t
Definition JFitMultiplicity.cc:35
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
GAUSS_LEGENDRE::xmax
const double xmax
Definition JQuadrature.cc:24
GAUSS_LEGENDRE::xmin
const double xmin
Definition JQuadrature.cc:23
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43

Variable Documentation

◆ livetime_t

const char* livetime_t = "LiveTime"

Definition at line 35 of file JFitMultiplicity.cc.

Generated by doxygen 1.12.0