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 kPdgAntiSQuark
const int kPdgAntiUQuark
const int kPdgF37m1950_DeltaP
const int kPdgUDDiquarkS1
const int kPdgAntiD0
const int kPdgSDDiquarkS0
const int kPdgAntiBQuark
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
const int kPdgWP
const int kPdgHadronicSyst
const int kPdgP33m1920_DeltaP
const int kPdgProton
TString ConvertVolumeCodeToString(int volumeCode, string volumeFilePath)
const int kPdgomega
const int kPdgClusterNN
const int kPdgP33m1232_DeltaP
const int kPdgP33m1920_DeltaPP
const int kPdgNuMu
const int kPdgF15m1680_N0
TTree * dk2nuTree
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const int kPdgAntiK0
const int kPdgDM
const int nquant
Definition: getContProf.C:25
const int kPdgMuon
const int kPdgSigmaPc
const int kPdgTgtO16
const int kPdgP11m1710_NP
const int kPdgP11m1710_N0
const int kPdgP33m1232_DeltaPP
const int kPdgSUDiquarkS0
double delta
Definition: runWimpSim.h:98
const int kPdgS11m1535_N0
const int kPdgS31m1620_DeltaPP
const int kPdgBQuark
const int kPdgPi0
const int kPdgAntiLambda
const int kPdgDQuark
int pots
const int kPdgAntiXi0
const int kPdgDMs
const int kPdgAntiNuMu
const int kPdgSUDiquarkS1
const int kPdgBindino
const int kPdgSigmaPPc
const int kPdgDP
TTree * dkmetaTree
T ldexp(const T &a, int b)
Definition: ldexp.hpp:19
const int kPdgF15m1680_NP
const int kPdgP33m1232_Delta0
const int kPdgAntiDQuark
const int kPdgF37m1950_DeltaPP
const int kPdgString
const int kPdgOmegaM
const int kPdgK0L
const int kPdgS11m1650_N0
const int kPdgRhoP
const int kPdgCQuark
const int kPdgGluon
const int kPdgAntiCQuark
const int kPdgTQuark
double estimate_pots(int highest_potnum)
void readWeightLocations(std::string locfilename, std::vector< bsim::Location > &locations)
OStream cerr
Definition: OStream.cxx:7
const int kPdgSigmaP
int myjob
const int kPdgAntiTau
const int kPdgAntiXiP
int ConvertGeantToPdg(int geant_code, std::string tag="?")
const int kPdgAntiSigma0
const int kPdgF35m1905_DeltaM
string construct_outfilename(string infilename)
const int kPdgPhi
const int kPdgPiM
Int_t job
identifying job # (keep files distinct)
Definition: dkmeta.h:86
const int kPdgClusterPP
bsim::Dk2Nu * dk2nu
const int kPdgCoulobtron
const int kPdgSigmaM
const int kPdgK0S
TRandom3 * rndm
const int kPdgTgtDeuterium
const int kPdgZ0
const int kPdgD13m1520_NP
const int kPdgS11m1650_NP
Double_t pots
protons-on-target
Definition: dkmeta.h:87
const int kPdgTgtFreeP
const int kPdgP13m1720_NP
const int kPdgSigma0
const XML_Char * s
Definition: expat.h:262
const int kPdgF37m1950_Delta0
const int kPdgEtab
void ConvertFinish()
const int nbins
Definition: cellShifts.C:15
const int kPdgPositron
double frac_diff_tolerance
const int kPdgP31m1910_DeltaP
const int kPdgD13m1520_N0
const int kPdgF35m1905_Delta0
TFile * treeFile
void printWeightLocations(std::vector< bsim::Location > &locations)
print the locations
const int kPdgD13m1700_NP
const XML_Char int const XML_Char * value
Definition: expat.h:331
const int kPdgP33m1232_DeltaM
const int kPdgAntiNeutron
const int kPdgAntiSigmaM
const int kPdgHadronicBlob
const int kPdgUUDiquarkS1
const double a
const int kPdgAntiNuTau
const int kPdgCluster
const int kPdgTau
const int kPdgSigma0c
const int kPdgTgtFreeN
const int kPdgEta
const int kPdgS31m1620_Delta0
const int kPdgY
void ConvertBookNtuple(std::string ofname="test_dk2nu.root")
double frac(double x)
Fractional part.
double obs_frac_diff_max
const int kPdgRho0
const int kPdgAntiMuon
const int kPdgAntiProton
const int kPdgIndep
const int kPdgP11m1440_NP
const int kPdgUQuark
OStream cout
Definition: OStream.cxx:6
const int kPdgK0
const int kPdgLambdaPc
const int kPdgWM
const int kPdgF35m1905_DeltaPP
bool isCloseEnough(double a, double b, double &difference, double &fdiff, int &nmsg, int mxmsg, string s)
const int kPdgRhoM
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
const int kPdgGamma
const int kPdgTgtC12
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const int kPdgElectron
const int kPdgEtaPrm
bool histCompare(double a, double b, string s)
bsim::DkMeta * dkmeta
const int kPdgD33m1700_Delta0
const int kPdgKM
const int kPdgD33m1700_DeltaPP
const int kPdgSQuark
const int kPdgAntiNuE
const int kPdgSSDiquarkS1
const int kPdgClusterNP
const int kPdgXi0
int Convert5xToPdg(int old_ntype)
double epsilon
size_t find_loc_index(string match)
const int kPdgNeutron
const hit & b
Definition: hits.cxx:21
const int kPdgP33m1920_DeltaM
assert(nhit_max >=nhit_nbins)
const int kPdgS11m1535_NP
const int kPdgF37m1950_DeltaM
const int kPdgTgtFe56
const int kPdgF35m1905_DeltaP
const int kPdgDPs
const int kPdgP31m1910_Delta0
const int kPdgNuTau
const int kPdgP31m1910_DeltaM
const int kPdgS31m1620_DeltaM
const int kPdgP33m1920_Delta0
const int kPdgNuE
const int kPdgAntiTQuark
const int kPdgS31m1620_DeltaP
const int kPdgD15m1675_N0
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
const int kPdgDDDiquarkS1
const int kPdgKP
::xsd::cxx::tree::token< char, normalized_string > token
Definition: Database.h:156
TString ConvertFlukaInteractionCodeToString(int interactionCode)
const int kPdgJpsi
const int kPdgEtac
const int kPdgP13m1720_N0
const int kPdgD15m1675_NP
const int kPdgD13m1700_N0
const int kPdgSDDiquarkS1
const int kPdgXiM
const int kPdgLambda
const int kPdgAntiSigmaP
procfile close()
const int kPdgD33m1700_DeltaP
const int kPdgUDDiquarkS0
const int kPdgP11m1440_N0
const int kPdgP31m1910_DeltaPP
const int kPdgD33m1700_DeltaM
void ConvertInit(std::string locfilename="$(DK2NU)/etc/locations.txt")
const int kPdgPiP
const int kPdgAntiOmegaP
const int kPdgD0