44int main(
int argc,
char **argv) {
57 JParser<> zap(
"Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.");
59 zap[
'f'] =
make_field(inputFile,
"input file");
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");
69 catch(
const exception &error) {
70 FATAL(error.what() << endl);
73 TFile* in = TFile::Open(inputFile.c_str(),
"exist");
75 if (in == NULL || !in->IsOpen()) {
76 FATAL(
"File: " << inputFile <<
" not opened." << endl);
84 const int nx = 1 + NUMBER_OF_PMTS;
85 const double xmin = -0.5;
86 const double xmax = nx - 0.5;
88 typedef JManager<TString, TH1D> JManager_t;
91 JManager_t HTF(
new TH1D(
"%F",
"", nx, xmin, xmax));
93 JManager_t HTFE(
new TH1D(
"%FE",
"", nx, xmin, xmax));
101 TH1D* hLiveTime = (TH1D*) in->Get(
livetime_t);
103 if (hLiveTime == NULL) {
104 FATAL(
"File " << inputFile <<
" does not contain livetime histogram." << endl);
117 TDirectory *fdir = out.mkdir(
"FIT");
118 TDirectory *pdir = out.mkdir(
"MUL");
121 double nT = hLiveTime->GetBinContent(1);
123 for (
int bin = 2; bin <= hLiveTime->GetNbinsX(); bin++) {
125 TString prefix = hLiveTime->GetXaxis()->GetBinLabel(bin);
127 DEBUG(
"Processing data from " << prefix << endl);
129 double T = nominalLiveTime ? nT : hLiveTime->GetBinContent(bin);
131 TString h1D_label = prefix + TString(
"_P" );
132 TString h2D_label = prefix + TString(
"_HT");
136 TH2D* h1D = (TH2D*) in->Get(h1D_label);
137 TH2D* h2D = (TH2D*) in->Get(h2D_label);
148 DEBUG(h1D_label <<
" not found" << endl);
153 TDirectory* cdir = fdir->mkdir(h2D_label); cdir->cd();
155 for (
int M = 2; M <= NUMBER_OF_PMTS; M++) {
157 int Mbin = h2D->GetXaxis()->FindBin(M);
159 TH1D* hI = h2D->ProjectionY(h2D_label + Form(
"_%d", M), Mbin, NUMBER_OF_PMTS);
163 if (hI->GetEntries() > 0) {
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;
170 double inclusiveCountVal = 0, inclusiveCountErr = 0;
172 if (hI->GetRMS() > minRMS) {
176 TF1 f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
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());
183 hI->Fit(&f1, option.c_str(),
"same");
185 double sg_v = f1.GetParameter(0);
186 double sg_e = f1.GetParError( 0);
188 double bg_v = f1.GetParameter(3);
189 double bg_e = f1.GetParError( 3);
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; }
195 double integral = hI->GetSumOfWeights();
199 inclusiveCountVal = V * (integral - (bg_v * N));
200 inclusiveCountErr = V * bg_e * N;
204 inclusiveCountErr = sg_e;
205 inclusiveCountVal = sg_v;
212 hI->Fit(&f1, option.c_str(),
"same");
216 if (inclusiveCountVal <= 0) { inclusiveCountVal = 0; }
228 int bin = HTF[h2D_label]->FindBin(M);
230 HTF[h2D_label]->SetBinContent(bin, inclusiveCountVal);
232 double binError = sqrt(inclusiveCountVal + pow(inclusiveCountErr, 2));
234 HTF[h2D_label]->SetBinError(bin, binError);
240 double exclusiveCountVal = 0;
242 if (M == NUMBER_OF_PMTS) {
243 exclusiveCountVal = inclusiveCountVal;
245 int prevBin = HTF[h2D_label]->GetXaxis()->FindBin(M+1);
246 exclusiveCountVal = inclusiveCountVal - HTF[h2D_label]->GetBinContent(prevBin);
249 HTFE[h2D_label]->SetBinContent(HTFE[h2D_label]->FindBin(M), exclusiveCountVal);
257 DEBUG(h1D_label <<
" not found" << endl);
262 HTFE[h2D_label]->Sumw2();
264 HTF[h2D_label ]->Scale(1.0 / T);
265 HTFE[h2D_label]->Scale(1.0 / T);
270 fdir = out.mkdir(
"HTF" ); fdir->cd();
274 fdir = out.mkdir(
"HTFE"); fdir->cd();