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