Jpp  pmt_effective_area_update
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JMergeCalibrateK40.cc File Reference
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TKey.h"
#include "TString.h"
#include "TMath.h"
#include "JROOT/JRootFileReader.hh"
#include "JSupport/JMeta.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)
Author
mkarel Auxiliary program to merge K40 data.

Definition at line 28 of file JMergeCalibrateK40.cc.

29 {
30  using namespace std;
31  using namespace JPP;
32 
33  vector<string> inputFile;
34  string outputFile;
35  int debug;
36 
37  try {
38 
39  JParser<> zap("Auxiliary program to merge K40 data.");
40 
41  zap['f'] = make_field(inputFile, "input file (output from JCalibrateK40).");
42  zap['o'] = make_field(outputFile, "output file (input to JFitK40).") = "mergek40.root";
43  zap['d'] = make_field(debug, "debug.") = 1;
44 
45  zap(argc, argv);
46  }
47  catch(const exception &error) {
48  FATAL(error.what() << endl);
49  }
50 
51 
52  gErrorIgnoreLevel = kError;
53 
54  TH1* weights_hist = NULL;
55 
56  typedef map<TString, TH1*> map_type;
57 
58  map_type zmap;
59 
60 
61  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
62 
63  DEBUG("Processing " << *i << endl) ;
64 
65  TFile in(i->c_str(), "read");
66 
67  // check if file contains weights vector
68 
69  if (!in.GetListOfKeys()->Contains(weights_hist_t)) {
70  FATAL("Histogram does not contain weights histogram." << endl);
71  }
72 
73  TH1* weights_hist_i = (TH1*) in.Get(weights_hist_t);
74 
75  if (weights_hist == NULL) {
76 
77  weights_hist = (TH1*) weights_hist_i->Clone();
78 
79  } else {
80 
81  // check normalization
82 
83  TAxis* x_axis = weights_hist->GetXaxis();
84 
85  if (weights_hist_i->GetBinContent(x_axis->FindBin(WS_t)) !=
86  weights_hist ->GetBinContent(x_axis->FindBin(WS_t))) {
87  FATAL("Number of y-bins does not match for this file " << *i << endl);
88  }
89 
90  if (weights_hist_i->GetBinContent(x_axis->FindBin(WB_t)) !=
91  weights_hist ->GetBinContent(x_axis->FindBin(WB_t))) {
92  FATAL("TMax_ns is not the same for this file " << *i << endl);
93  }
94 
95  // add normalization
96 
97  weights_hist->SetBinContent(x_axis->FindBin(W1_overall_t),
98  weights_hist ->GetBinContent(x_axis->FindBin(W1_overall_t)) +
99  weights_hist_i->GetBinContent(x_axis->FindBin(W1_overall_t)));
100  }
101 
102 
103  TIter iter(in.GetListOfKeys());
104 
105  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
106 
107  if (TString(key->GetName()).EndsWith(_2S) ||
108  TString(key->GetName()).EndsWith(_1B) ||
109  TString(key->GetName()).EndsWith(_1L)) {
110 
111  TH1* h1 = dynamic_cast<TH1*>(key->ReadObj());
112 
113  map_type::iterator p = zmap.find(h1->GetName());
114 
115  if (p == zmap.end()) {
116 
117  DEBUG("Clone " << h1->GetName() << endl);
118 
119  p = zmap.insert(make_pair(h1->GetName(), (TH1*) h1->Clone())).first;
120 
121  } else {
122 
123  DEBUG("Add " << h1->GetName() << endl);
124 
125  p->second->Add(h1);
126  }
127  }
128  }
129 
130  weights_hist->SetDirectory(0);
131 
132  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
133  i->second->SetDirectory(0);
134  }
135 
136  in.Close();
137  }
138 
139 
140  if (weights_hist == NULL) {
141  FATAL("Missing histogram " << weights_hist_t << endl);
142  }
143 
144 
145  // write file
146 
147  TFile out(outputFile.c_str(), "recreate");
148 
149  const Double_t WS = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WS_t)) ; // [ns]
150  const Double_t WB = weights_hist->GetBinContent(weights_hist->GetXaxis()->FindBin(WB_t)) ; // [ns]
151 
152  for (map_type::iterator i = zmap.begin(); i != zmap.end(); ++i) {
153 
154  if (i->first.EndsWith(_2S)) {
155 
156  TH2D* h2s = (TH2D*) i->second;
157 
158  if (h2s->GetEntries() != 0) {
159 
160  TH1D* h1b = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1B))->second;
161  TH1D* h1L = (TH1D*) zmap.find(TString(i->first).ReplaceAll(_2S, _1L))->second;
162 
163  //----------------------------------------------------------
164  // normalisation and background estimation
165  //----------------------------------------------------------
166 
167  for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
168 
169  // normalised background estimation
170  // if live time is zero then set bin count to zero and error to some large value
171 
172  if (h1L->GetBinContent(ix) == 0.0) {
173 
174  h1b->SetBinContent(ix, 0.0);
175  h1b->SetBinError (ix, 0.0);
176 
177  for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
178  h2s->SetBinContent(ix, iy, 0.0);
179  h2s->SetBinError (ix, iy, 999.999);
180  }
181 
182  } else {
183 
184  const double Y = h1b->GetBinContent(ix) / h1L->GetBinContent(ix) / WB;
185  const double W = h1b->GetBinError (ix) / h1L->GetBinContent(ix) / WB;
186 
187  h1b->SetBinContent(ix, Y);
188  h1b->SetBinError (ix, W);
189 
190  for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
191 
192  double y = h2s->GetBinContent(ix, iy) / h1L->GetBinContent(ix) / WS;
193  double w = h2s->GetBinError (ix, iy) / h1L->GetBinContent(ix) / WS;
194 
195  if (w == 0.0) {
196 
197  // for chi2 minimisation; set accordingly error of empty bins
198 
199  w = 1.0 / h1L->GetBinContent(ix) / WS;
200  }
201 
202  h2s->SetBinContent(ix, iy, y - Y);
203  h2s->SetBinError (ix, iy, TMath::Sqrt(w*w + W*W));
204  }
205  }
206  }
207 
208  // reset under and overflows.
209 
210  for (int ix = 0; ix <= h2s->GetXaxis()->GetNbins() + 1; ++ix) {
211  h2s->SetBinContent(ix, 0, 0.0); h2s->SetBinContent(ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
212  h2s->SetBinError (ix, 0, 0.0); h2s->SetBinError (ix, h2s->GetYaxis()->GetNbins() + 1, 0.0);
213  }
214 
215  for (int iy = 0; iy <= h2s->GetYaxis()->GetNbins() + 1; ++iy) {
216  h2s->SetBinContent(0, iy, 0.0); h2s->SetBinContent(h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
217  h2s->SetBinError (0, iy, 0.0); h2s->SetBinError (h2s->GetXaxis()->GetNbins() + 1, iy, 0.0);
218  }
219 
220  NOTICE("Biased coincidence rate [Hz] " << h2s->GetName() << ' ' << h2s->GetSumOfWeights() * WS << endl);
221 
222  h2s->SetName(TString(i->first).ReplaceAll(_2S, _2R));
223 
224  h2s->Write();
225  }
226  }
227  }
228 
229  putObject(&out, JMeta(argc, argv));
230 
231  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
232  JMeta::copy(i->c_str(), out);
233  }
234 
235  out.Write();
236  out.Close();
237 }
Utility class to parse command line options.
Definition: JParser.hh:1500
data_type w[N+1][M+1]
Definition: JPolint.hh:741
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
bool putObject(TDirectory *dir, const TObject &object)
Write object to ROOT directory.
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
string outputFile
static const char *const WS_t
Named bin for time residual bin width.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
static const char *const _1B
Name extension for 1D background.
#define NOTICE(A)
Definition: JMessage.hh:64
static const char *const _2R
Name extension for 2D rate measured.
int debug
debug level
Definition: JSirene.cc:63
#define FATAL(A)
Definition: JMessage.hh:67
static const char *const weights_hist_t
Histogram naming.
static const char *const _1L
Name extension for 1D live time.
static const char *const W1_overall_t
Named bin for duration of the run.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:40
static const char *const _2S
Name extension for 2D counts.
do if[[!-f $ACOUSTICS_WORKDIR/${KEY}.txt]]
Definition: JAcoustics.sh:38