EventSummary.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file EventSummary.cxx
3 /// \brief Simple representation of event for LEM use
4 /// \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
9 #include "TCanvas.h"
10 #include "TH2.h"
11 #include "TPad.h"
12 #include "TTree.h"
13 
14 #include <algorithm>
15 #include <cassert>
16 #include <map>
17 
18 #include "pthread.h"
19 
20 #include "Utilities/func/MathUtil.h"
21 
22 namespace lem
23 {
24  double EventSummary::fgDecayRate = 0;
25  double EventSummary::fgChargePower = .5;
26 
27 
28  //......................................................................
29  EventSummary::EventSummary(const std::vector<LiteHit>& hs, int i, std::vector<int> pdgs)
30  : id(i), photL0(0), photL1(0), photE0(0), photE1(0), enrich(false)
31  {
32  nhits = hs.size();
33  hits = new LiteHit[nhits];
34  for(unsigned int i = 0; i < nhits; ++i) hits[i] = hs[i];
35 
36  nParts = pdgs.size();
37  partPdgs = new int[nParts];
38  for(unsigned int i = 0; i < nParts; ++i) partPdgs[i] = pdgs[i];
39 
40  // TODO: the ordering of this with charge decay below makes a difference
41  if(fgChargePower != 1){
42  for(unsigned int i = 0; i < nhits; ++i){
43  hits[i].pecorr = pow(hits[i].pecorr, fgChargePower);
44  }
45  }
46 
47  if(fgDecayRate != 0){
48  for(unsigned int i = 0; i < nhits; ++i){
49  LiteHit& h = hits[i];
50 
51  const int planeDelta = h.Plane()-kVertexPlane;
52  const int cellDelta = h.Cell()-kVertexCell;
53 
54  // TODO: probably this can be optimized...
55  const double dist = util::pythag(planeDelta, cellDelta);
56 
57  const double decay = exp(-fgDecayRate*dist);
58 
59  h.pecorr *= decay;
60  } // end for i
61  } // end if nonzero
62 
63  Finalize();
64  }
65 
66  //......................................................................
68  {
69  *this = evt;
70  }
71 
72  //......................................................................
74  {
75  delete[] hits;
76  delete[] partPdgs;
77  }
78 
79  //......................................................................
81  {
82  // TODO: it would be much nicer for most of this to come from a default
83  // implementation.
84 
85  nhits = evt.nhits;
86  hits = new LiteHit[nhits];
87  for(int i = 0; i < nhits; ++i) hits[i] = evt.hits[i];
88 
89  totalPE = evt.totalPE;
90  trueEVis = evt.trueEVis;
91  y = evt.y;
92  id = evt.id;
93 
94  run = evt.run;
95  subrun = evt.subrun;
96  event = evt.event;
97 
98  ccnc = evt.ccnc;
99  pdg = evt.pdg;
100  origPdg = evt.origPdg;
101  mode = evt.mode;
102 
104  trueVtxCell = evt.trueVtxCell;
106 
107  nParts = evt.nParts;
108  partPdgs = new int[nParts];
109  for(int i = 0; i < nParts; ++i) partPdgs[i] = evt.partPdgs[i];
110 
111  photL0 = evt.photL0;
112  photL1 = evt.photL1;
113  photE0 = evt.photE0;
114  photE1 = evt.photE1;
115 
116  enrich = evt.enrich;
117 
118  return *this;
119  }
120 
121  //......................................................................
122  void EventSummary::FillHists(TH2* h1, TH2* h2, bool flipEven, bool flipOdd, bool idxs) const
123  {
124  TH2* h[2] = {h1, h2};
125 
126  for(unsigned int i = 0; i < nhits; ++i){
127  const int plane = hits[i].Plane();
128  int cell = hits[i].Cell();
129  const int v = plane%2;
130  if(v == 0 && flipEven) cell = 255-cell;
131  if(v == 1 && flipOdd) cell = 255-cell;
132  if(idxs){
133  h[plane%2]->Fill(plane, cell, hits[i].partIdx);
134  }
135  else{
136  h[plane%2]->Fill(plane, cell, hits[i].pecorr);
137  }
138  }
139  h[0]->SetMaximum(std::max(h[0]->GetMaximum(), h[1]->GetMaximum()));
140  h[1]->SetMaximum(std::max(h[0]->GetMaximum(), h[1]->GetMaximum()));
141  }
142 
143  //......................................................................
145  {
146  std::map<int, const char*> partNames;
147  partNames[12] = "#nu_{e}";
148  partNames[14] = "#nu_{#mu}";
149  partNames[16] = "#nu_{#tau}";
150  partNames[-12] = "#bar{#nu}_{e}";
151  partNames[-14] = "#bar{#nu}_{#mu}";
152  partNames[-16] = "#bar{#nu}_{#tau}";
153 
154  return TString::Format("%s %s#rightarrow%s, y=%.2lf, %dPE",
155  ccnc ? "NC" : "CC",
156  partNames[origPdg],
157  partNames[pdg],
158  y,
159  int(totalPE)).Data();
160  }
161 
162  //......................................................................
163  void EventSummary::SaveImage(TString fname, TString title) const
164  {
165  static pthread_mutex_t lock;
166  pthread_mutex_lock(&lock);
167 
168  TCanvas c("a", "b", 1200, 600);
169  c.Divide(2, 1);
170 
171  TH2F h[2] = {TH2F("", ";plane;cell", 128, 32, 160, 128, 64, 192),
172  TH2F("", ";plane;cell", 128, 32, 160, 128, 64, 192)};
173 
174  h[0].SetTitle(Description().c_str());
175 
176  h[1].SetTitle(title);
177 
178  FillHists(&h[0], &h[1]);
179 
180  for(int v = 0; v < 2; ++v){
181  c.cd(v+1);
182  h[v].Draw("colz");
183  }
184 
185  c.cd(0);
186  gPad->Print(fname);
187 
188  pthread_mutex_unlock(&lock);
189  }
190 
191  // Anonymous namespace, no-one else needs to see these
192  namespace{
193  // globals for tree
194  int gRun, gSubrun, gEvent;
195  int gPdg, gOrigPdg, gCCNC, gMode;
196  double gTrueEVis, gY;
197  int gNHits;
198  unsigned short gCellIdx[256*256];
199  float gPECorr[256*256];
200  unsigned char gPartIdx[256*256];
201  int gNParts;
202  int gPartPdgs[256];
203  double gPhotL0, gPhotL1;
204  double gPhotE0, gPhotE1;
205 
206  int gTrueVtxPlane, gTrueVtxCell, gTrueVtxCellOther;
207  }
208 
209  //......................................................................
211  {
212  // Bookkeeping
213  tr->Branch("run", &gRun);
214  tr->Branch("subrun", &gSubrun);
215  tr->Branch("event", &gEvent);
216 
217  // Truth information
218  tr->Branch("pdg", &gPdg);
219  tr->Branch("origpdg", &gOrigPdg);
220  tr->Branch("ccnc", &gCCNC);
221  tr->Branch("mode", &gMode);
222  tr->Branch("trueEVis", &gTrueEVis);
223  tr->Branch("y", &gY);
224 
225  // Length of the arrays
226  tr->Branch("nhits", &gNHits);
227 
228  // Hit positions and charges
229  tr->Branch("cellIdx", &gCellIdx, "cellIdx[nhits]/s");
230  tr->Branch("pecorr", &gPECorr, "pecorr[nhits]/F");
231  tr->Branch("partIdx", &gPartIdx, "partIdx[nhits]/b");
232 
233  // True Vertex
234  tr->Branch("trueVtxPlane", &gTrueVtxPlane);
235  tr->Branch("trueVtxCell", &gTrueVtxCell);
236  tr->Branch("trueVtxCellOther", &gTrueVtxCellOther);
237 
238  // Particles
239  tr->Branch("nParts", &gNParts);
240  tr->Branch("partPdgs", &gPartPdgs, "partPdgs[nParts]/I");
241 
242  // Photons
243  tr->Branch("photL0", &gPhotL0);
244  tr->Branch("photL1", &gPhotL1);
245  tr->Branch("photE0", &gPhotE0);
246  tr->Branch("photE1", &gPhotE1);
247  }
248 
249  //......................................................................
250  void EventSummary::ToTree(TTree* tr) const
251  {
252  gRun = run;
253  gSubrun = subrun;
254  gEvent = event;
255 
256  gPdg = pdg;
257  gOrigPdg = origPdg;
258  gCCNC = ccnc;
259  gMode = mode;
260 
261  gTrueEVis = trueEVis;
262  gY = y;
263 
264  //comparison of constant 65536 with expression of type 'const unsigned short' is always true
265  // so comment out the assert
266  //assert(nhits < 256*256);
267  gNHits = nhits;
268 
269  for(unsigned int n = 0; n < nhits; ++n){
270  gCellIdx[n] = hits[n].cellIdx;
271  gPECorr[n] = hits[n].pecorr;
272  gPartIdx[n] = hits[n].partIdx;
273  }
274 
275  gNParts = nParts;
276  assert(gNParts <= 255);
277  for(int n = 0; n < gNParts; ++n) gPartPdgs[n] = partPdgs[n];
278 
279  gTrueVtxPlane = trueVtxPlane;
280  gTrueVtxCell = trueVtxCell;
281  gTrueVtxCellOther = trueVtxCellOther;
282 
283  gPhotL0 = photL0;
284  gPhotL1 = photL1;
285  gPhotE0 = photE0;
286  gPhotE1 = photE1;
287 
288  tr->Fill();
289  }
290 
291  //......................................................................
293  {
294  // Bookkeeping
295  tr->SetBranchAddress("run", &gRun);
296  tr->SetBranchAddress("subrun", &gSubrun);
297  tr->SetBranchAddress("event", &gEvent);
298 
299  // Truth
300  tr->SetBranchAddress("pdg", &gPdg);
301  tr->SetBranchAddress("origpdg", &gOrigPdg);
302  tr->SetBranchAddress("ccnc", &gCCNC);
303  tr->SetBranchAddress("mode", &gMode);
304  tr->SetBranchAddress("trueEVis", &gTrueEVis);
305  tr->SetBranchAddress("y", &gY);
306 
307  // Hit positions and charges
308  tr->SetBranchAddress("nhits", &gNHits);
309  tr->SetBranchAddress("cellIdx", &gCellIdx);
310  tr->SetBranchAddress("pecorr", &gPECorr);
311  tr->SetBranchAddress("partIdx", &gPartIdx);
312 
313  // Particles
314  tr->SetBranchAddress("nParts", &gNParts);
315  tr->SetBranchAddress("partPdgs", &gPartPdgs);
316 
317  // True vertex
318  tr->SetBranchAddress("trueVtxPlane", &gTrueVtxPlane);
319  tr->SetBranchAddress("trueVtxCell", &gTrueVtxCell);
320  tr->SetBranchAddress("trueVtxCellOther", &gTrueVtxCellOther);
321 
322  // Photons
323  tr->SetBranchAddress("photL0", &gPhotL0);
324  tr->SetBranchAddress("photL1", &gPhotL1);
325  tr->SetBranchAddress("photE0", &gPhotE0);
326  tr->SetBranchAddress("photE1", &gPhotE1);
327  }
328 
329  //......................................................................
330  EventSummary EventSummary::FromTree(TTree* tr, int evtIdx, bool enrich)
331  {
332  tr->GetEntry(evtIdx);
333 
334  std::vector<LiteHit> hs;
335  hs.reserve(gNHits);
336  for(int n = 0; n < gNHits; ++n){
337  lem::LiteHit h(0, 0, gPECorr[n], gPartIdx[n]);
338  h.cellIdx = gCellIdx[n];
339 
340  hs.push_back(h);
341  } // end for n
342 
343  std::vector<int> pdgs(gNParts);
344  for(int n = 0; n < gNParts; ++n) pdgs[n] = gPartPdgs[n];
345 
346  EventSummary ret(hs, -1, pdgs);
347 
348  ret.ccnc = gCCNC;
349  ret.pdg = gPdg;
350  ret.origPdg = gOrigPdg;
351  ret.mode = gMode;
352  ret.trueEVis = gTrueEVis;
353  ret.y = gY;
354 
355  ret.trueVtxPlane = gTrueVtxPlane;
356  ret.trueVtxCell = gTrueVtxCell;
357  ret.trueVtxCellOther = gTrueVtxCellOther;
358 
359  ret.photL0 = gPhotL0;
360  ret.photL1 = gPhotL1;
361  ret.photE0 = gPhotE0;
362  ret.photE1 = gPhotE1;
363 
364  ret.enrich = enrich;
365 
366  return ret;
367  }
368 
369  //......................................................................
371  {
372  totalPE = 0;
373 
374  for(unsigned int i = 0; i < nhits; ++i){
375  LiteHit& h = hits[i];
376  totalPE += h.pecorr;
377  } // end for i
378 
379  // Put in order of cellIdx, might improve locality
380  std::sort(hits, hits+nhits);
381  }
382 
383  //......................................................................
385  {
386  if(nhits < 2) return 0;
387  // Hits are sorted
388  return hits[nhits-1].Plane() - hits[0].Plane();
389  }
390 }
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
unsigned short nhits
Definition: EventSummary.h:53
T max(const caf::Proxy< T > &a, T b)
std::string Description() const
bool enrich
Was this event generated in a filtered-for-pi0 job?
Definition: EventSummary.h:73
Simple representation of event for LEM use.
Definition: EventSummary.h:26
static double fgChargePower
Definition: EventSummary.h:76
constexpr T pow(T x)
Definition: pow.h:75
unsigned char partIdx
Definition: LiteHit.h:26
const int kVertexPlane
Definition: EventSummary.h:22
float pecorr
Definition: LiteHit.h:24
int Cell() const
Definition: LiteHit.h:39
EventSummary & operator=(const EventSummary &evt)
static void InitFromTree(TTree *tr)
double dist
Definition: runWimpSim.h:113
PID
Definition: FillPIDs.h:13
static void InitToTree(TTree *tr)
int Plane() const
Definition: LiteHit.h:38
unsigned short cellIdx
Definition: LiteHit.h:25
void SaveImage(TString fname, TString title="") const
const int kVertexCell
Definition: EventSummary.h:23
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
static EventSummary FromTree(TTree *tr, int evtIdx, bool enrich)
void ToTree(TTree *tr) const
Definition: run.py:1
TH1F * h2
Definition: plot.C:45
TH1F * h1
int Length() const
In planes.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Definition: decay.py:1
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
unsigned char nParts
Definition: EventSummary.h:67
static double fgDecayRate
Definition: EventSummary.h:75
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Simple representation of event for LEM use.
Compressed hit info, basic component of LEM events.
Definition: LiteHit.h:18
void FillHists(TH2 *h1, TH2 *h2, bool flipEven=false, bool flipOdd=false, bool idxs=false) const
int gMode
Definition: gtestDISSF.cxx:53