Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFitMultiplicity.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <map>
6 #include <cmath>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TH2D.h"
12 #include "TF2.h"
13 #include "TMath.h"
14 #include "TFitResult.h"
15 
16 #include "JROOT/JMinimizer.hh"
17 
18 #include "JPhysics/JConstants.hh"
21 
22 #include "JROOT/JRootToolkit.hh"
23 #include "JTools/JRange.hh"
24 #include "JSupport/JMeta.hh"
25 
27 #include "JCalibrate/JFitK40.hh"
28 
29 #include "JROOT/JManager.hh"
30 
31 #include "Jeep/JPrint.hh"
32 #include "Jeep/JParser.hh"
33 #include "Jeep/JMessage.hh"
34 
35 const char* livetime_t = "LiveTime";
36 
37 /**
38  * \file
39  *
40  * Auxiliary program to fit multiplicity rates from JMonitorMultiplicity output.
41  * \author mlincett
42  */
43 
44 int main(int argc, char **argv) {
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 }
281 
282 
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1711
int main(int argc, char *argv[])
Definition: Main.cc:15
Dynamic ROOT object management.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
const JPolynome f1(1.0, 2.0, 3.0)
Function.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
do set_variable OUTPUT_DIRECTORY $WORKDIR T
ROOT I/O of application specific meta data.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
Physics constants.
const char * livetime_t
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
const double xmin
Definition: JQuadrature.cc:23
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
Auxiliary class to define a range between two values.
then fatal The output file must have the wildcard in the e g root fi 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
Definition: JCanberra.sh:48
Utility class to parse command line options.
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62