common_convert.C
Go to the documentation of this file.
1 /// common code for use in the ntuple converter code
2 ///
3 /// \author Robert Hatcher <rhatcher \at fnal.gov>
4 /// Fermi National Accelerator Laboratory
5 ///
6 /// \created 2012-11-07
7 /// \version $Id: common_convert.C,v 1.3 2012-11-15 09:09:26 rhatcher Exp $
8 ///==========================================================================
9 
10 #include <iostream>
11 #include <iomanip>
12 #include <cassert>
13 #include <map>
14 #include <cstring>
15 #include <fstream>
16 #include <stdlib.h>
17 
18 #include <float.h> // FLT_EPSILON and DBL_EPSILON definitions used for
19  // floating point comparisons
20 
21 #include "TFile.h"
22 #include "TTree.h"
23 #include "TRandom3.h"
24 #include "TH1D.h"
25 #include "TString.h"
26 
27 // dk2nu headers
28 #include "tree/dk2nu.h"
29 #include "tree/dkmeta.h"
30 /// functions for reading location text file and calculating
31 /// energy and weight vectors for locations
32 #include "tree/readWeightLocations.h"
33 #include "tree/calcLocationWeights.h"
34 
35 // some globals
36 TRandom3* rndm = 0;
39 TFile* treeFile = 0;
40 TTree* dk2nuTree = 0;
41 TTree* dkmetaTree = 0;
42 int myjob = 0;
43 int pots = 0;
44 
45 double frac_diff_tolerance = 5.0e-4;
46 double obs_frac_diff_max = 0;
47 
48 bool histCompare(double a, double b, string s); // forward declaration
49 
50 //____________________________________________________________________________
51 
52 void ConvertInit(std::string locfilename = "$(DK2NU)/etc/locations.txt")
53 {
54  ///-----------------------------------------------------------------------
55  ///
56  /// equivalent to NumiAnalysis::NumiAnalysis() in g4numi
57  ///
58  ///-----------------------------------------------------------------------
59 
60  // create the objects used for writing entries
61  dk2nu = new bsim::Dk2Nu;
62  dkmeta = new bsim::DkMeta;
63 
64  rndm = new TRandom3();
65 
66  // read the text file for locations, fill the dkmeta object
67  /*bsim::*/readWeightLocations(locfilename,dkmeta);
68 
69  // print out what we have for locations
70  /*bsim::*/printWeightLocations(dkmeta);
71 
72 }
73 
74 void ConvertBookNtuple(std::string ofname = "test_dk2nu.root")
75 {
76  ///-----------------------------------------------------------------------
77  ///
78  /// equivalent to NumiAnalysis::book() in g4numi
79  ///
80  ///-----------------------------------------------------------------------
81 
82  // create file, book tree, set branch address to created object
83  treeFile = new TFile(ofname.c_str(),"RECREATE");
84 
85  dk2nuTree = new TTree("dk2nuTree","neutrino ntuple");
86  dk2nuTree->Branch("dk2nu","bsim::Dk2Nu",&dk2nu,32000,99);
87 
88  dkmetaTree = new TTree("dkmetaTree","neutrino ntuple metadata");
89  dkmetaTree->Branch("dkmeta","bsim::DkMeta",&dkmeta,32000,99);
90 }
91 
93 {
94  ///-----------------------------------------------------------------------
95  ///
96  /// equivalent to NumiAnalysis::finish() in g4numi
97  ///
98  ///-----------------------------------------------------------------------
99 
100  /// fill the rest of the metadata (locations filled above)
101  //no! would clear location info // dkmeta->clear();
102  dkmeta->job = myjob; // needs to match the value in each dk2nu entry
103  dkmeta->pots = pots; // ntuple represents this many protons-on-target
104  //leave to user// dkmeta->beamsim = "common_convert.C";
105  //leave to user// dkmeta->physics = "bogus";
106 
107  //dkmeta->vintnames.push_back("mytemp_42");
108  //dkmeta->vintnames.push_back("mytemp_ipot");
109  // push entry out to meta-data tree
110  dkmetaTree->Fill();
111 
112  // finish and clean-up
113  treeFile->cd();
114  dk2nuTree->Write();
115  dkmetaTree->Write();
116 
117  treeFile->cd(); // be here so any booked histograms get created inside output file
118  treeFile->mkdir("zzz_diff_hists");
119  treeFile->cd("zzz_diff_hists");
120  histCompare(0,0,"WriteHist");
121 
122  treeFile->Close();
123 
124  delete treeFile; treeFile=0;
125  dk2nuTree=0;
126  dkmetaTree=0;
127 }
128 //____________________________________________________________________________
129 
130 const int kPdgNuE = 12; //
131 const int kPdgAntiNuE = -12; //
132 const int kPdgNuMu = 14; //
133 const int kPdgAntiNuMu = -14; //
134 const int kPdgNuTau = 16; //
135 const int kPdgAntiNuTau = -16; //
136 
137 const int kPdgElectron = 11; //
138 const int kPdgPositron = -11; //
139 const int kPdgMuon = 13; //
140 const int kPdgAntiMuon = -13; //
141 const int kPdgTau = 15; //
142 const int kPdgAntiTau = -15; //
143 
144 const int kPdgUQuark = 2; //
145 const int kPdgAntiUQuark = -2; //
146 const int kPdgDQuark = 1; //
147 const int kPdgAntiDQuark = -1; //
148 const int kPdgSQuark = 3; //
149 const int kPdgAntiSQuark = -3; //
150 const int kPdgCQuark = 4; //
151 const int kPdgAntiCQuark = -4; //
152 const int kPdgBQuark = 5; //
153 const int kPdgAntiBQuark = -5; //
154 const int kPdgTQuark = 6; //
155 const int kPdgAntiTQuark = -6; //
156 
157 const int kPdgDDDiquarkS1 = 1103; // dd, spin = 1
158 const int kPdgUDDiquarkS0 = 2101; // ud, spin = 0
159 const int kPdgUDDiquarkS1 = 2103; // ud, spin = 1
160 const int kPdgUUDiquarkS1 = 2203; // uu, spin = 1
161 const int kPdgSDDiquarkS0 = 3101; // sd, spin = 0
162 const int kPdgSDDiquarkS1 = 3103; // sd, spin = 1
163 const int kPdgSUDiquarkS0 = 3201; // su, spin = 0
164 const int kPdgSUDiquarkS1 = 3203; // su, spin = 1
165 const int kPdgSSDiquarkS1 = 3303; // ss, spin = 1
166 
167 const int kPdgProton = 2212; //
168 const int kPdgAntiProton = -2212; //
169 const int kPdgNeutron = 2112; //
170 const int kPdgAntiNeutron = -2112; //
171 const int kPdgLambda = 3122; // Lambda
172 const int kPdgAntiLambda = -3122; // \bar{Lambda}
173 const int kPdgSigmaP = 3222; // Sigma+
174 const int kPdgSigma0 = 3212; // Sigma0
175 const int kPdgSigmaM = 3112; // Sigma-
176 const int kPdgAntiSigmaP = -3222; // \bar{Sigma+}
177 const int kPdgAntiSigma0 = -3212; // \bar{Sigma0}
178 const int kPdgAntiSigmaM = -3112; // \bar{Sigma-}
179 const int kPdgXi0 = 3322; // Xi0
180 const int kPdgXiM = 3312; // Xi-
181 const int kPdgAntiXi0 = -3322; // \bar{Xi0}
182 const int kPdgAntiXiP = -3312; // \bar{Xi+}
183 const int kPdgOmegaM = 3334; // Omega-
184 const int kPdgAntiOmegaP = -3334; // \bar{Omega+}
185 const int kPdgLambdaPc = 4122; // Lambda+_{c}
186 const int kPdgSigma0c = 4112; // Sigma0_{c}
187 const int kPdgSigmaPc = 4212; // Sigma+_{c}
188 const int kPdgSigmaPPc = 4222; // Sigma++_{c}
189 
190 const int kPdgP33m1232_DeltaM = 1114; // P33(1232) Delta-
191 const int kPdgP33m1232_Delta0 = 2114; // P33(1232) Delta0
192 const int kPdgP33m1232_DeltaP = 2214; // P33(1232) Delta+
193 const int kPdgP33m1232_DeltaPP = 2224; // P33(1232) Delta++
194 const int kPdgS11m1535_N0 = 22112; // S11(1535) N0
195 const int kPdgS11m1535_NP = 22212; // S11(1535) N+
196 const int kPdgD13m1520_N0 = 1214; // D13(1520) N0
197 const int kPdgD13m1520_NP = 2124; // D13(1520) N+
198 const int kPdgS11m1650_N0 = 32112; // S11(1650) N0
199 const int kPdgS11m1650_NP = 32212; // S11(1650) N+
200 const int kPdgD13m1700_N0 = 21214; // D13(1700) N0
201 const int kPdgD13m1700_NP = 22124; // D13(1700) N+
202 const int kPdgD15m1675_N0 = 2116; // D15(1675) N0
203 const int kPdgD15m1675_NP = 2216; // D15(1675) N+
204 const int kPdgS31m1620_DeltaM = 1112; // S31(1620) Delta-
205 const int kPdgS31m1620_Delta0 = 1212; // S31(1620) Delta0
206 const int kPdgS31m1620_DeltaP = 2122; // S31(1620) Delta+
207 const int kPdgS31m1620_DeltaPP = 2222; // S31(1620) Delta++
208 const int kPdgD33m1700_DeltaM = 11114; // D33(1700) Delta-
209 const int kPdgD33m1700_Delta0 = 12114; // D33(1700) Delta0
210 const int kPdgD33m1700_DeltaP = 12214; // D33(1700) Delta+
211 const int kPdgD33m1700_DeltaPP = 12224; // D33(1700) Delta++
212 const int kPdgP11m1440_N0 = 12112; // P11(1440) N0
213 const int kPdgP11m1440_NP = 12212; // P11(1440) N+
214 const int kPdgP13m1720_N0 = 31214; // P13(1720) N0
215 const int kPdgP13m1720_NP = 32124; // P13(1720) N+
216 const int kPdgF15m1680_N0 = 12116; // F15(1680) N0
217 const int kPdgF15m1680_NP = 12216; // F15(1680) N+
218 const int kPdgP31m1910_DeltaM = 21112; // P31(1910) Delta-
219 const int kPdgP31m1910_Delta0 = 21212; // P31(1910) Delta0
220 const int kPdgP31m1910_DeltaP = 22122; // P31(1910) Delta+
221 const int kPdgP31m1910_DeltaPP = 22222; // P31(1910) Delta++
222 const int kPdgP33m1920_DeltaM = 21114; // P33(1920) Delta-
223 const int kPdgP33m1920_Delta0 = 22114; // P33(1920) Delta0
224 const int kPdgP33m1920_DeltaP = 22214; // P33(1920) Delta+
225 const int kPdgP33m1920_DeltaPP = 22224; // P33(1920) Delta++
226 const int kPdgF35m1905_DeltaM = 1116; // F35(1905) Delta-
227 const int kPdgF35m1905_Delta0 = 1216; // F35(1905) Delta0
228 const int kPdgF35m1905_DeltaP = 2126; // F35(1905) Delta+
229 const int kPdgF35m1905_DeltaPP = 2226; // F35(1905) Delta++
230 const int kPdgF37m1950_DeltaM = 1118; // F37(1950) Delta-
231 const int kPdgF37m1950_Delta0 = 2118; // F37(1950) Delta0
232 const int kPdgF37m1950_DeltaP = 2218; // F37(1950) Delta+
233 const int kPdgF37m1950_DeltaPP = 2228; // F37(1950) Delta++
234 const int kPdgP11m1710_N0 = 42112; // P11(1710) N0
235 const int kPdgP11m1710_NP = 42212; // P11(1710) N+
236 
237 const int kPdgPiP = 211; // pi+
238 const int kPdgPiM = -211; // pi-
239 const int kPdgPi0 = 111; // pi0
240 const int kPdgEta = 221; // eta
241 const int kPdgEtaPrm = 331; // eta' (prime)
242 const int kPdgEtac = 441; // eta_{c}
243 const int kPdgEtab = 551; // eta_{b}
244 const int kPdgRhoP = 213; // rho+
245 const int kPdgRhoM = -213; // rho-
246 const int kPdgRho0 = 113; // rho0
247 const int kPdgomega = 223; // omega (the meson, not Omega the baryon)
248 const int kPdgPhi = 333; // phi
249 const int kPdgJpsi = 443; // J/psi
250 const int kPdgY = 553; // Y
251 const int kPdgKP = 321; // K+
252 const int kPdgKM = -321; // K-
253 const int kPdgK0 = 311; // K0
254 const int kPdgAntiK0 = -311; // \bar{K0}
255 const int kPdgK0L = 130; // K0_{long}
256 const int kPdgK0S = 310; // K0_{short}
257 const int kPdgDP = 411; // D+
258 const int kPdgDM = -411; // D-
259 const int kPdgD0 = 421; // D0
260 const int kPdgAntiD0 = -421; // \bar{D0}
261 const int kPdgDPs = 431; // D+_{s}
262 const int kPdgDMs = -431; // D-_{s}
263 
264 const int kPdgGluon = 21; // gluon
265 const int kPdgGamma = 22; // photon
266 const int kPdgZ0 = 23; // Z
267 const int kPdgWP = 24; // W+
268 const int kPdgWM = -24; // W-
269 
270 // Note: PDG codes for nuclear targets can be computed using pdg::IonPdgCode(A,Z)
271 // PDG2006 convention: 10LZZZAAAI
272 // Define names for some commonly used nuclear PDG codes:
273 const int kPdgTgtFreeP = 1000010010;
274 const int kPdgTgtFreeN = 1000000010;
275 const int kPdgTgtDeuterium = 1000010020;
276 const int kPdgTgtC12 = 1000060120;
277 const int kPdgTgtO16 = 1000080160;
278 const int kPdgTgtFe56 = 1000260560;
279 
280 // PDG codes for GENIE special particles
281 const int kPdgHadronicSyst = 2000000001; // dis hadronic system before hadronization
282 const int kPdgHadronicBlob = 2000000002; // unmodelled fraction of the hadronic system
283 const int kPdgBindino = 2000000101; // binding energy subtracted from f/s nucleons
284 const int kPdgCoulobtron = 2000000102; // coulomb energy subtracted from f/s leptons
285 const int kPdgClusterNN = 2000000200; // a nn cluster within a nucleus
286 const int kPdgClusterNP = 2000000201; // a np cluster within a nucleus
287 const int kPdgClusterPP = 2000000202; // a pp cluster within a nucleus
288 
289 // PDG codes for PYTHIA/JETSET special particles
290 const int kPdgCluster = 91;
291 const int kPdgString = 92;
292 const int kPdgIndep = 93;
293 
294 //____________________________________________________________________________
295 int ConvertGeantToPdg(int geant_code, std::string tag="?")
296 {
297  if (geant_code == 3) return kPdgElectron; // 11 / e-
298  if (geant_code == 2) return kPdgPositron; // -11 / e+
299  if (geant_code == 6) return kPdgMuon; // 13 / mu-
300  if (geant_code == 5) return kPdgAntiMuon; // -13 / mu+
301  if (geant_code == 34) return kPdgTau; // 15 / tau-
302  if (geant_code == 33) return kPdgAntiTau; // -15 / tau+
303  if (geant_code == 8) return kPdgPiP; // 211 / pi+
304  if (geant_code == 9) return kPdgPiM; // -211 / pi-
305  if (geant_code == 7) return kPdgPi0; // 111 / pi0
306  if (geant_code == 17) return kPdgEta; // 221 / eta
307  if (geant_code == 11) return kPdgKP; // 321 / K+
308  if (geant_code == 12) return kPdgKM; // -321 / K-
309  if (geant_code == 10) return kPdgK0L; // 130 / K0_{long}
310  if (geant_code == 16) return kPdgK0S; // 310 / K0_{short}
311  if (geant_code == 35) return kPdgDP; // 411 / D+
312  if (geant_code == 36) return kPdgDM; // -411 / D-
313  if (geant_code == 37) return kPdgD0; // 421 / D0
314  if (geant_code == 38) return kPdgAntiD0; // -421 / \bar{D0}
315  if (geant_code == 39) return kPdgDPs; // 431 / D+_{s}
316  if (geant_code == 40) return kPdgDMs; // -431 / D-_{s}
317  if (geant_code == 1) return kPdgGamma; // 22 / photon
318  if (geant_code == 44) return kPdgZ0; // 23 / Z
319  if (geant_code == 42) return kPdgWP; // 24 / W+
320  if (geant_code == 43) return kPdgWM; // -24 / W-
321  if (geant_code == 14) return kPdgProton; // 2212
322  if (geant_code == 15) return kPdgAntiProton; // -2212
323  if (geant_code == 13) return kPdgNeutron; // 2112
324  if (geant_code == 25) return kPdgAntiNeutron; // -2112
325  if (geant_code == 18) return kPdgLambda; // 3122 / Lambda
326  if (geant_code == 26) return kPdgAntiLambda; // -3122 / \bar{Lambda}
327  if (geant_code == 19) return kPdgSigmaP; // 3222 / Sigma+
328  if (geant_code == 20) return kPdgSigma0; // 3212 / Sigma0
329  if (geant_code == 21) return kPdgSigmaM; // 3112 / Sigma-
330  if (geant_code == 29) return kPdgAntiSigmaP; // -3112 / \bar{Sigma+}
331  if (geant_code == 28) return kPdgAntiSigma0; // -3212 / \bar{Sigma0}
332  if (geant_code == 27) return kPdgAntiSigmaM; // -3112 / \bar{Sigma-}
333  if (geant_code == 22) return kPdgXi0; // 3322 / Xi0
334  if (geant_code == 23) return kPdgXiM; // 3312 / Xi-
335  if (geant_code == 30) return kPdgAntiXi0; // -3322 / \bar{Xi0}
336  if (geant_code == 31) return kPdgAntiXiP; // -3312 / \bar{Xi+}
337  if (geant_code == 24) return kPdgOmegaM; // 3334 / Omega-
338  if (geant_code == 32) return kPdgAntiOmegaP; // -3334 / \bar{Omega+}
339 
340  // some rare Geant3 codes that don't really need definitions in PDGCodes.h
341  const int kPdgDeuteron = 1000010020; // pdg::IonPdgCode(2,1);
342  const int kPdgTritium = 1000010030; // pdg::IonPdgCode(3,1);
343  const int kPdgAlpha = 1000020040; // pdg::IonPdgCode(4,2);
344  const int kPdgHe3 = 1000020030; // pdg::IonPdgCode(3,2);
345  if (geant_code == 45) return kPdgDeuteron;
346  if (geant_code == 46) return kPdgTritium;
347  if (geant_code == 47) return kPdgAlpha;
348  if (geant_code == 49) return kPdgHe3;
349 
350  // NuMI ntuples have an odd neutrino "geant" convention
351  if (geant_code == 56) return 14; // nu_mu
352  if (geant_code == 55) return -14; // nu_mu_bar
353  if (geant_code == 53) return 12; // nu_e
354  if (geant_code == 52) return -12; // nu_e_bar
355 
356  if (geant_code != 0 )
357  std::cerr << "## Can not convert geant code: " << geant_code
358  << " to PDG (" << tag << ")" << std::endl;
359  return 0;
360 }
361 //____________________________________________________________________________
362 
363 int Convert5xToPdg(int old_ntype)
364 {
365  // NuMI ntuples have an odd neutrino "geant" convention
366  switch ( old_ntype ) {
367  case 56: return 14; break; // kPdgNuMu; break;
368  case 55: return -14; break; // kPdgAntiNuMu; break;
369  case 53: return 12; break; // kPdgNuE; break;
370  case 52: return -12; break; // kPdgAntiNuE; break;
371  default:
372  std::cerr << "Convert5xToPdg saw ntype " << old_ntype
373  << " -- unknown " << std::endl;
374  assert(0);
375  }
376  return 0;
377 }
378 
379 //____________________________________________________________________________
380 
381 TString ConvertFlukaInteractionCodeToString(int interactionCode)
382 {
383 
384  TString interactionString = "null";
385 
386  if(interactionCode == 100) interactionString = "elastic";
387  if(interactionCode == 101) interactionString = "inelastic";
388  if(interactionCode == 102) interactionString = "decaySecondaries";
389  if(interactionCode == 103) interactionString = "deltaRay";
390  if(interactionCode == 104) interactionString = "pairProduction";
391  if(interactionCode == 105) interactionString = "bremsstrahlung";
392  if(interactionCode == 110) interactionString = "decayProducts";
393  if(interactionCode == 208) interactionString = "bremsstrahlung";
394  if(interactionCode == 210) interactionString = "moller";
395  if(interactionCode == 212) interactionString = "bhabha";
396  if(interactionCode == 214) interactionString = "inFlightAnnihilation";
397  if(interactionCode == 215) interactionString = "annihilationAtRest";
398  if(interactionCode == 217) interactionString = "pairProduction";
399  if(interactionCode == 219) interactionString = "comptonScattering";
400  if(interactionCode == 221) interactionString = "photoelectric";
401  if(interactionCode == 225) interactionString = "rayleighScattering";
402  if(interactionCode == 300) interactionString = "interactionSecndaries";
403  if(interactionCode == 400) interactionString = "deltaRay";
404 
405  return interactionString;
406 }
407 
408 //____________________________________________________________________________
409 
410 TString ConvertVolumeCodeToString(int volumeCode, string volumeFilePath){
411 
412  TString FlukaName;
413  TString GeantName;
414 
415  // Open the Fluka Volumes file
416  ifstream volumeFile;
417  volumeFile.open(volumeFilePath.c_str());
418 
419  if(volumeCode==-2222){
420  FlukaName = "not-defined";
421  return FlukaName;
422  }
423 
424  // Make sure the file opened properly
425  if(!volumeFile.good()) {
426  std::cout << "Error opening Volumes file." << std::endl;
427  FlukaName="ErrorOpeningVolumeFile";
428  return FlukaName;
429  }
430 
431  // Variables for file reading
432  const int maxChar = 1024; // Max number of characters to read from a single line from the input file
433  char volumeLine[maxChar]; // Array to hold a single line from the input file
434  bool volumeFound = false; // Flag whether the requested input volume has been found
435 
436  // Loop over lines from the input file until the volume is found or the end of the file is reached
437  while(!(volumeFound || volumeFile.eof())) {
438  // Get the next line
439  if(volumeFile.getline(volumeLine, maxChar).fail()) {
440  std::cout << "Error accessing the next line. " << std::endl
441  << "It may just be the end of the file, "
442  << "and the EOF bit may not have been set as expected." << std::endl;
443  break;
444  }
445 
446  // Make a std::string from the char array
447  // This must be done BEFORE tokenizing the character array,
448  // as strtok "destroys" the array as it goes
449  std::string stringLine(volumeLine);
450 
451  // Get the first set of contiguous characters that skips leading and contains no spaces
452  // token is the number that labels the volume
453  char* token = strtok(volumeLine, " ");
454  double number = atof(token);
455 
456 
457  if(number==volumeCode) {
458  // If the input volume number is in the current line, set the boolean flag to terminate the while loop
459  volumeFound = true;
460  GeantName = strtok(NULL, " "); // this is Geant volume name
461  FlukaName = strtok(NULL, " "); // this is Fluka volume name
462  }
463  } // end of loop until volume number is found or end of the location file is reached
464 
465  // If the volume was not found, warn the user and exit
466  if(!volumeFound) {
467  std::cout << "Could not find the following volume " <<volumeCode << "." << std::endl;
468  FlukaName="FlukaVolumeNameNotFound";
469  return FlukaName;
470  }
471 
472  //If everything went smoothly, return the FlukaName of the volume
473  return FlukaName;
474 }
475 
476 //____________________________________________________________________________
477 
478 size_t find_loc_index(string match)
479 {
480  // find the index in the array of locations
481  for ( size_t i=0; i < dkmeta->location.size(); ++i ) {
482  //cout << "looking for \"" << match << "\" vs \""
483  // << dkmeta->nameloc[i] << "\"" << endl;
484  if ( dkmeta->location[i].name == match ) return i;
485  }
486  //cout << "failed to find match" << endl;
487  return -1;
488 }
489 //____________________________________________________________________________
490 
491 string construct_outfilename(string infilename)
492 {
493  string ofname = infilename;
494  size_t dot = ofname.find(".root");
495  ofname.insert(dot,".dk2nu");
496  size_t slash = ofname.find_last_of("/");
497  if ( slash != std::string::npos ) {
498  ofname.erase(0,slash+1);
499  }
500  return ofname;
501 }
502 
503 //____________________________________________________________________________
504 
505 double estimate_pots(int highest_potnum)
506 {
507  // looks like low counts are due to "evtno" not including
508  // protons that miss the actual target (hit baffle, etc)
509  // this will vary with beam conditions parameters
510  // so we should round way up, those generating flugg files
511  // aren't going to quantize less than 1000
512  // though 500 would probably be okay, 100 would not.
513  // also can't use _last_ potnum because muons decays don't
514  // have theirs set
515  const Int_t nquant = 1000; // 500; // 100
516  const Double_t rquant = nquant;
517 
518  Int_t estimate = (TMath::FloorNint((highest_potnum-1)/rquant)+1)*nquant;
519  return estimate;
520 }
521 
522 //____________________________________________________________________________
523 bool isCloseEnough(double a, double b, double& difference, double& fdiff,
524  int& nmsg, int mxmsg, string s)
525 {
526  // For float or double, need to determine equality while
527  // accounting for some degree of imprecision
528  // The method applied here is adopted from Theodore C. Belding's
529  // fcmp implementation of Donald Knuth's algorithm.
530  double value[2] = {a,b};
531  //double epsilon = 10.*FLT_EPSILON;
532  double epsilon = 10.*FLT_EPSILON;
533  // if ( isDouble[0] || isDouble[1] ) epsilon = DBL_EPSILON;
534  int exponent;
535  frexp(fabs(value[0]) > fabs(value[1]) ? value[0] : value[1], &exponent);
536  double delta = ldexp(epsilon,exponent);
537  //double
538  difference = value[0] - value[1];
539  if ( difference > delta || difference < -delta ) {
540  double sum = a+b;
541  if ( TMath::Abs(sum) < 1.0e-30 ) sum = 1.0e-30; // protect divide-by-0
542  //double
543  fdiff = TMath::Abs(2.0*(a-b)/(a+b));
544  if ( fdiff > obs_frac_diff_max ) obs_frac_diff_max = fdiff;
545  if ( fdiff < frac_diff_tolerance ) return true; // we'll allow it anyway
546  ++nmsg;
547  if ( nmsg <= mxmsg ) {
548  cout << "## \"" << s << "\" " << a << " " << b
549  //<< " delta " << delta
550  << " diff " << difference
551  << " fdiff " << fdiff << endl;
552  //cout << "## ndecay (dkproc) " << dk2nu->decay.ndecay << endl;
553  if ( nmsg == mxmsg )
554  cout << "## last message for tag \"" << s << "\"" << endl;
555  }
556  return false; // nope
557  }
558  return true; // these are close
559 }
560 
561 bool histCompare(double a, double b, string s)
562 {
563  // booking histograms should end up in the output file
564 
565  static map<std::string,int> nmsg;
566  double diff, frac;
567  const int mxmsg = 10, nbins = 100;
568  static map<std::string,TH1D*> diffmap;
569  static map<std::string,TH1D*> fracmap;
570  /// static bool first = true;
571  //if ( first ) {
572  // first = false;
573  // TH1::SetDefaultSumw2(true);
574  //}
575  if ( s == "WriteHist" ) {
576  map<std::string,TH1D*>::iterator itr;
577  for ( itr=diffmap.begin(); itr != diffmap.end(); ++itr ) (*itr).second->Write();
578  for ( itr=fracmap.begin(); itr != fracmap.end(); ++itr ) (*itr).second->Write();
579  return true;
580  }
581 
582  bool close = isCloseEnough(a,b,diff,frac,nmsg[s],mxmsg,s);
583  TH1D* diffh = diffmap[s];
584  TH1D* frach = fracmap[s];
585  if ( ! diffh ) {
586  string sdh = "difference " + s;
587  diffh = new TH1D(sdh.c_str(),sdh.c_str(),nbins,1,-1); // autoadjust limits by reversed values
588  diffh->Sumw2();
589  diffmap[s] = diffh;
590  }
591  if ( ! frach ) {
592  string sfh = "frac diff " + s;
593  frach = new TH1D(sfh.c_str(),sfh.c_str(),nbins,1,-1); // autoadjust limits by reversed values
594  frach->Sumw2();
595  fracmap[s] = frach;
596  }
597  diffh->Fill(diff);
598  frach->Fill(frac);
599 
600  return close;
601 }
602 
603 //____________________________________________________________________________
604 
const int kPdgF37m1950_DeltaM
const int kPdgP13m1720_N0
const int kPdgGluon
const int kPdgSSDiquarkS1
const int kPdgS31m1620_Delta0
const int kPdgAntiNuTau
const int kPdgF37m1950_DeltaPP
const int kPdgPiP
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
const int kPdgEtaPrm
const int kPdgUQuark
const int kPdgPiM
const int kPdgDPs
const int kPdgRhoM
const int kPdgSUDiquarkS0
const int kPdgPi0
TTree * dkmetaTree
const int kPdgSUDiquarkS1
const int kPdgBindino
const int kPdgS31m1620_DeltaM
const int kPdgAntiProton
const int kPdgAntiCQuark
const int kPdgGamma
const int kPdgP33m1920_Delta0
const int kPdgP31m1910_DeltaP
const int nquant
Definition: getContProf.C:25
double delta
Definition: runWimpSim.h:98
const int kPdgClusterNN
const int kPdgElectron
const int kPdgAntiNuMu
const int kPdgAntiBQuark
const int kPdgomega
const int kPdgClusterNP
const int kPdgCoulobtron
const int kPdgF35m1905_DeltaM
const int kPdgF35m1905_Delta0
const int kPdgAntiMuon
const int kPdgF35m1905_DeltaP
const int kPdgLambda
const int kPdgAntiOmegaP
void readWeightLocations(std::string locfilename, std::vector< bsim::Location > &locations)
OStream cerr
Definition: OStream.cxx:7
const int kPdgPositron
const int kPdgAntiK0
void ConvertInit(std::string locfilename="$(DK2NU)/etc/locations.txt")
bsim::Dk2Nu * dk2nu
const int kPdgAntiDQuark
bsim::DkMeta * dkmeta
const int kPdgD33m1700_DeltaPP
const int kPdgNuE
const int kPdgP31m1910_DeltaPP
const int kPdgAntiNuE
Int_t job
identifying job # (keep files distinct)
Definition: dkmeta.h:86
const int kPdgCQuark
const int kPdgRho0
const int kPdgP33m1232_Delta0
const int kPdgD33m1700_DeltaM
const int kPdgSDDiquarkS0
const int kPdgP33m1920_DeltaM
const int kPdgSigmaM
const int kPdgAntiTQuark
const int kPdgDQuark
double obs_frac_diff_max
const int kPdgAntiSigmaP
const int kPdgK0L
const int kPdgAntiSigma0
const int kPdgY
const int kPdgTQuark
const int kPdgTgtFreeN
const int kPdgXiM
Double_t pots
protons-on-target
Definition: dkmeta.h:87
const int kPdgTgtC12
const int kPdgP11m1440_NP
const int kPdgEtac
const int kPdgBQuark
const int kPdgF37m1950_Delta0
const XML_Char * s
Definition: expat.h:262
const int kPdgClusterPP
const int kPdgSQuark
const int nbins
Definition: cellShifts.C:15
const int kPdgDDDiquarkS1
const int kPdgKM
const int kPdgWM
const int kPdgTgtFe56
const int kPdgTgtFreeP
const int kPdgP13m1720_NP
const int kPdgS11m1650_N0
const int kPdgAntiNeutron
void printWeightLocations(std::vector< bsim::Location > &locations)
print the locations
bool isCloseEnough(double a, double b, double &difference, double &fdiff, int &nmsg, int mxmsg, string s)
const XML_Char int const XML_Char * value
Definition: expat.h:331
const int kPdgS11m1535_NP
const int kPdgTgtO16
const int kPdgS31m1620_DeltaP
const double a
const int kPdgSigmaP
int pots
const int kPdgD15m1675_NP
double estimate_pots(int highest_potnum)
const int kPdgNuTau
const int kPdgS11m1535_N0
const int kPdgOmegaM
const int kPdgKP
const int kPdgS11m1650_NP
const int kPdgK0S
const int kPdgD0
const int kPdgD13m1700_NP
const int kPdgAntiD0
const int kPdgD13m1520_N0
const int kPdgNuMu
const int kPdgXi0
double frac_diff_tolerance
const int kPdgNeutron
const int kPdgP31m1910_Delta0
const int kPdgP33m1920_DeltaPP
double frac(double x)
Fractional part.
const int kPdgP33m1920_DeltaP
const int kPdgUDDiquarkS1
const int kPdgP33m1232_DeltaP
const int kPdgSigmaPc
const int kPdgD13m1700_N0
const int kPdgZ0
int Convert5xToPdg(int old_ntype)
const int kPdgString
const int kPdgD13m1520_NP
TString ConvertFlukaInteractionCodeToString(int interactionCode)
const int kPdgEtab
const int kPdgPhi
TRandom3 * rndm
const int kPdgF35m1905_DeltaPP
int ConvertGeantToPdg(int geant_code, std::string tag="?")
const int kPdgF15m1680_NP
OStream cout
Definition: OStream.cxx:6
const int kPdgWP
const int kPdgSDDiquarkS1
const int kPdgIndep
const int kPdgUDDiquarkS0
TString ConvertVolumeCodeToString(int volumeCode, string volumeFilePath)
const int kPdgD15m1675_N0
const int kPdgD33m1700_Delta0
const int kPdgTau
const int kPdgDM
const int kPdgAntiSQuark
const int kPdgHadronicBlob
const int kPdgProton
double epsilon
int myjob
const int kPdgF15m1680_N0
const int kPdgTgtDeuterium
const int kPdgP11m1710_N0
const int kPdgSigma0c
const hit & b
Definition: hits.cxx:21
const int kPdgJpsi
assert(nhit_max >=nhit_nbins)
const int kPdgDMs
const int kPdgP33m1232_DeltaPP
const int kPdgP33m1232_DeltaM
const int kPdgLambdaPc
const int kPdgF37m1950_DeltaP
TTree * dk2nuTree
const int kPdgSigmaPPc
const int kPdgMuon
const int kPdgD33m1700_DeltaP
size_t find_loc_index(string match)
TFile * treeFile
const int kPdgAntiXiP
const int kPdgS31m1620_DeltaPP
const int kPdgSigma0
bool histCompare(double a, double b, string s)
const int kPdgK0
void ConvertFinish()
const int kPdgAntiTau
const int kPdgCluster
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
::xsd::cxx::tree::token< char, normalized_string > token
Definition: Database.h:156
const int kPdgP31m1910_DeltaM
const int kPdgAntiLambda
const int kPdgP11m1710_NP
const int kPdgRhoP
string construct_outfilename(string infilename)
void ConvertBookNtuple(std::string ofname="test_dk2nu.root")
const int kPdgAntiUQuark
const int kPdgHadronicSyst
const int kPdgAntiSigmaM
const int kPdgEta
procfile close()
const int kPdgUUDiquarkS1
const int kPdgDP
const int kPdgAntiXi0
const int kPdgP11m1440_N0
enum BeamMode string