Library.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file Library.cxx
3 /// \brief Collection of events for matching
4 /// \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "LEM/func/Library.h"
8 
9 #include "LEM/func/OnDemand.h"
10 #include "LEM/func/Util.h"
11 
12 #include "Utilities/func/MathUtil.h"
13 
14 #include "TChain.h"
15 
16 #include <iostream>
17 
18 #include <sys/mman.h>
19 
20 namespace lem
21 {
22  // Private to this file
23  namespace{
24  char* kLibAddr = (char*)0x00000100CAFE0000;
25  }
26 
27  //......................................................................
29  {
30  if(evt.IsSig()){
31  // Scale up to keep as much signal as possible in sample, weighted
32  // correctly internally. The missing prefactor is (ss2th13*ssth23)
33  return util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
34  }
35  if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 14){
36  // CC mu
37  return 1-util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
38  }
39  if(evt.ccnc == 0 && abs(evt.pdg) == 12 && abs(evt.origPdg) == 12){
40  // Beam nue
41  // return 1-.1*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
42  // Want library matching to ignore nues
43  return 0;
44  }
45  if(evt.ccnc == 1 && evt.pdg == evt.origPdg){
46  // NC
47  return 1;
48  }
49  if(evt.ccnc == 1 && evt.pdg != evt.origPdg){
50  // Swapped NC. We play clever reweighting games so this is OK
51  return 1;
52  }
53  if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 12){
54  // CC mu appearance, take ss2th13=0.1, maximal th23
55  return .1*.5*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
56  }
57  // Don't know what this is
58  return 0;
59  }
60 
61  //......................................................................
63  bool touchAll)
64  {
65  Library* lib = new Library;
66 
67  lib->fMappingSize = MMapFileAtAddress(libPath+"/map.bin", kLibAddr,
68  touchAll, "/tmp/lem_preload_lib.lock");
69 
70  lib->fSummary = *((LibrarySummary*)kLibAddr);
71 
72  lib->fLibSize = lib->fSummary.N;
73 
74  lib->fLib = (MatchableEvent*)((char*)kLibAddr+sizeof(LibrarySummary));
75 
76  lib->fLibDownsample = (MatchableEvent*)((char*)lib->fLib+lib->fLibSize*sizeof(MatchableEvent));
77 
78  return lib;
79  }
80 
81  //......................................................................
83  {
84  CreateOnDemandMapping((libPath+"/map.bin").c_str(), kLibAddr);
85 
86  Library* lib = new Library;
87 
88  lib->fSummary = *((LibrarySummary*)kLibAddr);
89 
90  lib->fLibSize = lib->fSummary.N;
91 
92  lib->fLib = (MatchableEvent*)((char*)kLibAddr+sizeof(LibrarySummary));
93 
94  lib->fLibDownsample = (MatchableEvent*)((char*)lib->fLib+lib->fLibSize*sizeof(MatchableEvent));
95 
96  return lib;
97  }
98 
99  //......................................................................
101  int splitLibIndex,
102  int splitLibFactor,
103  int downsampleFactor)
104  {
105  // TODO: this order won't match what heads is expecting
106 
107  TChain ch("makelib/tree");
108  ch.Add((libPath+"/lem_hist_sig.root").c_str());
109  ch.Add((libPath+"/lem_hist_bkg.root").c_str());
110  const int Nunenrich = ch.GetEntries();
111  ch.Add((libPath+"/lem_hist_enrich.root").c_str());
112 
114 
115  const int startIdx = (splitLibIndex*ch.GetEntries())/splitLibFactor;
116  const int endIdx = ((splitLibIndex+1)*ch.GetEntries())/splitLibFactor;
117 
118  const int N = endIdx-startIdx;
119 
120  assert(N > 0);
121 
122  std::cout << "Loading libraries..." << std::endl;
123  int nSig = 0, nBkg = 0;
124  double totSig = 0, totBkg = 0, totEnrich = 0;
125  int nTrueNC = 0, nSwapNC = 0;
126 
127  Library* lib = new Library;
128  lib->fMappingSize = 0;
129 
130  lib->fLibSize = N;
131  lib->fLib = new MatchableEvent[N];
132  lib->fLibDownsample = new MatchableEvent[N];
133 
134  for(int n = startIdx; n < endIdx; ++n){
135  if((n-startIdx)%10000 == 0)
136  std::cout << (n-startIdx)/10000 << "/" << N/10000 << std::endl;
137 
138  const bool isEnrich = n >= Nunenrich;
139  EventSummary sum = EventSummary::FromTree(&ch, n, isEnrich);
140  sum.id = n;
141 
142  // This has to happen here, and not in MakePID where the other
143  // reweighting happens, because it feeds into the totals for the
144  // different event classes that we write out.
145  const double P = isEnrich ? 1 : SimpleSurvivalProb(sum);
146  // Screws up indexing, don't do it. FindMatches can cope
147  // if(P == 0) continue; // beam nue
148 
150  evt.weight = P;
151 
152  if(!isEnrich && sum.ccnc == 1){
153  if(sum.pdg == sum.origPdg) ++nTrueNC; else ++nSwapNC;
154  }
155 
156  lib->fLib[n-startIdx] = evt;
157 
158  lib->fLibDownsample[n-startIdx] = evt.Downsampled(downsampleFactor);
159 
160  if(isEnrich){
161  ++totEnrich;
162  }
163  else{
164  if(evt.IsSig()){
165  ++nSig;
166  totSig += P;
167  }
168  else{
169  ++nBkg;
170  totBkg += P;
171  }
172  }
173  } // end for n
174 
175  lib->fSummary.nTrueNC = nTrueNC;
176  lib->fSummary.nSwapNC = nSwapNC;
177  lib->fSummary.totBkg = totBkg;
178  lib->fSummary.totSig = totSig;
179  lib->fSummary.totEnrich = totEnrich;
180 
181  std::cout << "Loaded " << lib->fSummary << std::endl;
182 
183  return lib;
184  }
185 
186  //......................................................................
188  {
189  if(fMappingSize != 0){
190  // munmap(kLibAddr, fMappingSize);
191  }
192  else{
193  // Was loaded from trees
194  delete[] fLib;
195  delete[] fLibDownsample;
196  }
197  }
198 } // namespace
long fMappingSize
Definition: Library.h:49
Simple representation of event for LEM use.
Definition: EventSummary.h:26
Attach some information used in matching to an EventSummary.
Collection of events for matching.
Collection of events for matching.
Definition: Library.h:18
void abs(TH1 *hist)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
static void InitFromTree(TTree *tr)
LibrarySummary fSummary
Definition: Library.h:44
Definition: __init__.py:1
MatchableEvent * fLib
Definition: Library.h:46
PID
Definition: FillPIDs.h:14
#define P(a, b, c, d, e, x)
long MMapFileAtAddress(const std::string &fname, void *addr, bool touchAll, const std::string &lockName)
Definition: Util.cxx:18
bool IsSig() const
Definition: EventSummary.h:36
string lib
Definition: cvnie.py:9
int evt
static EventSummary FromTree(TTree *tr, int evtIdx, bool enrich)
MatchableEvent Downsampled(int factor) const
Details of the library LEM matches were made against.
OStream cout
Definition: OStream.cxx:6
static double SimpleSurvivalProb(const EventSummary &evt)
Definition: Library.cxx:28
Utilities for mapping in library files on-demand.
T sin(T number)
Definition: d0nt_math.hpp:132
static Library * FromMMapOnDemand(const std::string &libPath)
Definition: Library.cxx:82
static Library * FromTrees(const std::string &libPath, int splitLibFactor, int splitLibIndex, int downsampleFactor)
Definition: Library.cxx:100
assert(nhit_max >=nhit_nbins)
MatchableEvent * fLibDownsample
Definition: Library.h:47
static Library * FromMMap(const std::string &libPath, bool touchAll)
Definition: Library.cxx:62
unsigned int fLibSize
Definition: Library.h:45
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
void CreateOnDemandMapping(const std::string &fname, const char *base)
Definition: OnDemand.cxx:123
enum BeamMode string