MCCheater_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MCCheater.cxx
3 // \brief
4 // \version $Id: MCCheater.cxx,v 1.16 2012-09-20 22:01:30 greenc Exp $
5 // \author $Author: greenc $
6 ////////////////////////////////////////////////////////////////////////
7 // C++ includes
8 #include <functional>
9 #include <cmath>
10 #include <iostream>
11 #include <map>
12 #include <string>
13 #include <vector>
14 
15 // ROOT includes
16 #include "TArray.h"
17 #include "TArrayD.h"
18 #include "TH1.h"
19 #include "TH2.h"
20 #include "TMatrixD.h"
21 #include "TObject.h"
22 
23 // ART includes
27 #include "art_root_io/TFileDirectory.h"
28 #include "art_root_io/TFileService.h"
33 #include "fhiclcpp/ParameterSet.h"
35 
36 // NOvA includes
37 #include "Geometry/Geometry.h"
39 #include "MCCheater/BackTracker.h"
40 #include "MCCheater/SimHit.h"
41 #include "MCCheater/SimTrack.h"
42 #include "Simulation/FLSHit.h"
43 #include "Simulation/FLSHitList.h"
44 #include "Simulation/Particle.h"
46 
47 class TH1F;
48 class TH2F;
49 
50 namespace cheat
51 {
52 
53  class MCCheater : public art::EDProducer {
54 
55  public:
56 
57  explicit MCCheater(fhicl::ParameterSet const &pset); // the explicit tag tells the compiler that GENIEGen is not an art::ParameterSet
58  virtual ~MCCheater();
59 
60  // Optional, read/write access to event
61  void produce(art::Event& evt);
62 
63  // Optional if you want to be able to configure from event display, for example
64 // void reconfigure(std::istream &is,
65 // std::ostream &os,
66 // fhicl::ParameterSet pset);
67 
68  // Optional use if you have histograms, ntuples, etc you want around for every event
69  void beginJob();
70 
71  private:
72 
73  //double prod[6][6];
74  //double tran[7][7];
75  int ntest;
76 
79 // std::string fSimHitListLabel;
80 // std::string fSimTrackListLabel;
81 
82  // histograms
83  TH1F* fMuE;
84  TH1F* fElE;
85  TH1F* felL;
86  TH1F* fmuL;
87  TH1F* ftotalE;
88  TH1F* ftotEXV;
89  TH1F* ftotEYV;
90  TH1F* fdcosX;
91  TH1F* fdcosY;
92  TH1F* fdcosZ;
93  //TH1F* fdelE;
94  TH1F* fsimhitE;
95  //TH1F* fEpar;
96  TH1F* fFLSE;
97  TH1F* fnumsh;
98  TH1F* fNumFLS;
99  TH1F* ftlen;
100  //TH1F* fshlen;
101  TH1F* fthreshE;
102  TH1F* fnplanes;
103  //TH2F* fVertx;
104  //TH2F* fFidVert;
105  //TH1F* fzvtx;
106  TH1F* fparte;
107  TH1F* fnpart;
108  TH1F* fntrajpts;
109  //TH1F* fXRange;
110  //TH1F* fYRange;
111  //TH1F* fZRange;
112  //TH1F* fnshptr;
113  TH1F* fcorx;
114  TH1F* fcory;
115  //TH2F* fstsxy;
116  //TH2F* fstexy;
117  TH1F* fclhit;
118  TH1F* fzstart;
119 
120  };
121 ////////////////////////////////////////////////////////////////////////
123  : EDProducer(pset)
124  , fMCTruthListLabel (pset.get< std::string >("MCTruthListLabel" ))
125  , fFLSHitListLabel (pset.get< std::string >("FLSHitListLabel" ))
126  {
127  // tell module what its making
128  produces< std::vector<cheat::SimHit> >();
129  produces< std::vector<cheat::SimTrack> >();
130  // produces< std::vector<cheat::SimVertex> >();
131  // produces< std::vector<cheat::SimCluster> >();
132 
133  // Be noisy to demonstrate what's happening
134  mf::LogInfo("MCCheater") << " MCCheater::MCCheater()\n";
135 
136  // print this to keep track of version number
137  mf::LogInfo("MCCheater") << " ****************************************\n"
138  << " * MCCheater Version Number 4 *\n"
139  << " ****************************************\n";
140 
141  }
142 
143  //..............................................................
145  {
146  //======================================================================
147  // Clean up any memory allocated by your module
148  //======================================================================
149  }
150 
151  //......................................................................
153  {
154 
155  //
156  // Book histograms, ntuples, initialize counts etc., etc., ...
157  //
159 
160  fMuE = tfs->make<TH1F>("fMuE", ";True Muon Energy (GeV)", 50, 0., 10.);
161  fElE = tfs->make<TH1F>("fElE", ";True Electron Energy (GeV)", 50, 0., 10.);
162  felL = tfs->make<TH1F>("felL", ";length(cm)", 60, 0., 1200.);
163  fmuL = tfs->make<TH1F>("fmuL", ";length(cm)", 60, 0., 1200.);
164  ftotalE = tfs->make<TH1F>("ftotalE", ";E of SimHit (MeV)", 100, 0., 50.);
165  ftotEXV = tfs->make<TH1F>("ftotEXV", ";E of Xview SimHit (MeV)", 100, 0., 50.);
166  ftotEYV = tfs->make<TH1F>("ftotEYV", ";E of Yview SimHit (MeV)", 100, 0., 50.);
167  fdcosX = tfs->make<TH1F>("fdcosX", ";Dir cos X", 100, -1., 1.);
168  fdcosY = tfs->make<TH1F>("fdcosY", ";Dir cos Y", 100, -1., 1.);
169  fdcosZ = tfs->make<TH1F>("fdcosZ", ";Dir cos Z", 100, -1., 1.);
170  fthreshE = tfs->make<TH1F>("fthreshE", ";E of threshold SimHit (MeV)", 100, 0., 50.);
171  fsimhitE = tfs->make<TH1F>("fsimhitE", ";Total SimHit Energy per Event (GeV)", 200, 0., 20.);
172  // fshlen = tfs->make<TH1F>("fshlen", "Shower Length (cm)", 250, 0., 1000.);
173  ftlen = tfs->make<TH1F>("ftlen", ";Track Length (cm)", 250, 0., 500.);
174  // fEpar = tfs->make<TH1F>("fEpar", ";SimHit E per Event Tr 0 (GeV)", 250, 0., 5.);
175  fFLSE = tfs->make<TH1F>("fFLSE", ";Total FLSHit Energy per Event (GeV)", 200, 0., 20.);
176  fnumsh = tfs->make<TH1F>("fnumsh", ";Number of SimHits per Event", 50, -0.5, 999.5);
177  fnplanes = tfs->make<TH1F>("fnplanes", ";Number of Planes in Track", 100, 0., 200.);
178  // fVertx = tfs->make<TH2F>("fVertx", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
179  // fstsxy = tfs->make<TH2F>("fstsxy", ";SimTrack start y vs x", 150, -150., 150., 240, -220., 220.);
180  // fstexy = tfs->make<TH2F>("fstexy", ";SimTrack end y vs x", 150, -150., 150., 240, -220., 220.);
181  // fFidVert = tfs->make<TH2F>("fFidfVert", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
182  // fzvtx = tfs->make<TH1F>("fzvtx", ";Zvtx (cm)", 350, 0., 7000.);
183  fNumFLS = tfs->make<TH1F>("fNumFLS", ";Number of FLSHits in Event", 100, -0.5, 1999.5);
184  fnpart = tfs->make<TH1F>("fnpart", ";Number of Particles in Event", 50, -0.5, 49.5);
185  fparte = tfs->make<TH1F>("fparte", ";Total Energy of Particles in Event (GeV)", 200, 0., 20.);
186  fntrajpts = tfs->make<TH1F>("fntrajpts", ";Number of Trajectory Points per Particle", 50, 0., 50.);
187  // fXRange = tfs->make<TH1F>("fXRange", ";Range of SimHit X Coords", 120, -210., 210.);
188  // fYRange = tfs->make<TH1F>("fYRange", ";Range of SimHit Y Coords", 120, -210., 210.);
189  // fZRange = tfs->make<TH1F>("fZRange", ";Range of SimHit Z Coords", 100, 0., 1000.);
190  // fnshptr = tfs->make<TH1F>("fnshptr", ";Number of SimHits per Track", 50, 0., 100.);
191  fcorx = tfs->make<TH1F>("fcorx", ";X Coord of xview plane", 120, -810., 810.);
192  fcory = tfs->make<TH1F>("fcory", ";Y Coord of yview plane", 120, -810., 810.);
193  fclhit = tfs->make<TH1F>("fclhit", ";Penetration Depth with SimHits", 125, 0., 250.);
194  fzstart = tfs->make<TH1F>("fzstart", ";Plane of particle start", 100, -1.0, 199.0);
195 
196  }
197 
199  {
200  //======================================================================
201  // Called for every event. "Reco" implies that you are adding
202  // information to the event
203  //======================================================================
204 
205  mf::LogInfo("MCCheater") << "MCCheater::produce() " << " New Reco Event\n";
206 
207  // increment counter for M matrix
208  ntest++;
209 
210  /*
211  // get MCTruth info - should loop over this list since there may be more
212  // than 1 interaction, a neutrino event and a cosmic event for example
213  // other than that, don't need this list for the products of the cheater
214  art::Handle< std::vector<simb::MCTruth> > mclist;
215  evt.getByLabel(fMCTruthListLabel, mclist);
216  if (mclist->empty())
217  {
218  mf::LogError("MCCheater") << "No MCTruth List - bummer\n";
219  return;
220  }
221 
222  // for now, only have 1 MCTruth per event, so only need mclist[0]
223  // art::PtrVector<simb::MCTruth> mctruth; don't need this either?
224  // for (unsigned int i=0; i<mclist->size(); ++i)
225  // {
226  art::Ptr<simb::MCTruth> mct(mclist, 0);
227  std::cout << " Number of Particles from MCTruth = " << mct->NParticles() << "\n";
228  // }
229 
230  mf::LogInfo("MCCheater") << "Number of MCTruth lists = "<< mclist->size() << "\n";
231 
232 
233  // fill a histogram with the size of each MCTruth
234  // for now, there is only 1 MCTruth list per MC event, so only need to use i=0
235 
236  for(unsigned int i=0; i < mclist->size(); ++i)
237  {
238  std::cout << " Number of particles from MCTruth = " << mctruth[i]->NParticles() << "\n";
239  }
240  */
241 
242  // PARTICLE LIST NEEDED FOR CHEATER
243  // get Particle info - needed for Sim Objects
245  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
246 
247  // First list needed for Cheater - the particle list
248  // using the particle list, make some plots,
249  // make a map linking trackID with particles
250  // and a map linking trackID with its particleID - this is all that's needed for cheater
251  std::map<int, const sim::Particle*> trpartmap; // dont really need this map
252  std::map<int, int> trpidmap; // TrackID ParticleID map
253  double PartE = 0.;
254  unsigned int NPart = 0;
255  int npar = 0;
256  int pcount = 0;
257  for (sim::ParticleNavigator::const_iterator j = pnav.begin(); j != pnav.end(); ++j){
258 
259  const sim::Particle p = *((*j).second);
260  NPart += 1;
261  npar += 1;
262 
263  int pid = p.PdgCode(); // what be this particle
264  // energy of particle
265  PartE += p.E();
266  unsigned int ntrajpts = p.NumberTrajectoryPoints();
267  fntrajpts->Fill(ntrajpts);
268  // for now, just check length of particle
269  // double plength = sqrt((enx-stx)*(enx-stx)+(eny-sty)*(eny-sty)+(enz-stz)*(enz-stz));
270  if (abs(pid)==13) fMuE->Fill(p.E()); // true muon energy for some reason?
271  if (abs(pid)==11) fElE->Fill(p.E()); // true electron energy
272  int tid = p.TrackId();
273  // reduce number of particles with some obvious? cuts, no neutrinos, pi0s, photons, nuclei
274  if (abs(pid)==12 || abs(pid)==14) continue; // no neutrinos is only cut so far
275  pcount++;
276  trpartmap[tid] = &p;
277  trpidmap[tid]=pid;
278  } // end loop over particles in this list
279  mf::LogInfo("MCCheater") << "Final particle count after cuts in particle list = " << pcount << "\n";
280  fparte->Fill(PartE);
281  fnpart->Fill(NPart);
282 
283  // mf::LogInfo("MCCheater") << "TrackID Particle Map size " << trpartmap.size() << "\n";
284  // mf::LogInfo("MCCheater") << "TrackID ParticleID Map size " << trpidmap.size() << "\n";
285 
286 
287  // FLSHITS ALSO NEEDED FOR CHEATER
288  // need FLSHits also - get FLSHit info
290  evt.getByLabel(fFLSHitListLabel, hitlist);
291  if (hitlist->empty()){
292  mf::LogError("MCCheater") << "No FLSHit Lists - bummer\n";
293  return;
294  }
295  // mf::LogInfo("MCCheater") << "Number of FLSHit hitlists = " << hitlist->size() << "\n";
296 
297  // Make maps if needed - these maps not necessary for the Cheater
298  // std::map<unsigned int, std::pair<int, int> > tidmgmmap; // map linking track ID to its mother and grandmother
299  std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidpcpmap;
300  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcptidmap; // map linking pair plane cell to vector of TrackIDs
301  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcpmotmap; // map linking pair plane cell to vector of TrackID mothers
302  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcppidmap; // map linking pair plane cell to vector of particleIDs
303  // std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcpendmap; // map linking pair plane cell to vector of Edeps
304  // std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcptimmap; // map linking pair plane cell to vector of times
305 
306  // FLSHits are formed for each particle that deposits energy in a cell, so there
307  // can be more than 1 FLSHit in a particular cell corresponding to each particle
308  // that contributed to the energy. And the cell IDs in FLSHits are not unique, when
309  // combined with the PlaneId, then a unique channel ID for the whole detector can be determined.
310  int nFLShits = 0;
311  double flshitE = 0.;
312  // geo::Geometry& geo = util::EDMUtils::GetGeometry(evt);
314  // Make a map linking FLSHits to TrackID
315  std::map<int, std::vector<const sim::FLSHit*> > trhitmap; // map linking TrackID to FLSHits
316  // can define a map with a pair of keys and a value, so can make a map with plane and cell as key and
317  // the value a vector of eg, FLSHits
318  std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> > pcpflsmap; // map linking plane cell pair to FLShits
319 
320  for (unsigned int i=0; i<hitlist->size(); ++i) // loop over FLSHit Lists - 1 per process interaction
321  {
322 
323  art::Ptr<sim::FLSHitList> hlist(hitlist, i);
324  // mf::LogInfo("MCCheater") << "Num of FLSHits in list = " << hlist->fHits.size() << "\n";
325 
326  const std::vector<sim::FLSHit>& hits = hlist->fHits;
327 
328  nFLShits += hits.size();
329  // loop over FLSHits in this list
330  for (unsigned int j=0; j<hits.size(); ++j)
331  {
332  int tid = hits[j].GetTrackID();
333  if (tid<0) continue; // just protection - shouldn't be here this is incoming particle (neutrino)
334  // here, loop over particle list to get mother of this trackID
335  int tmoth = -2; // set mother to -2 for now
336  int tgmoth = -2; // set grandmother to -2 for now
337  int tpid = 0;
338  // need to use only the particle id from particle list for this track - FLSHit pid can be different one
339  for (sim::ParticleNavigator::const_iterator k = pnav.begin(); k != pnav.end(); ++k)
340  {
341  const sim::Particle part = *((*k).second);
342  if (part.TrackId()==tid)
343  {
344  tpid = part.PdgCode();
345  tmoth = part.Mother();
346  }
347  // have to find grandma myself -
348  tgmoth = pnav.EveId(part.TrackId());
349  }
350  int smoth = tmoth;
351  if (smoth > -2) smoth = 1;
352  if (tmoth==-2) mf::LogInfo("MCCheater") << " THIS TrackID NOT IN PARTICLE LIST = " << tid << " " << hits[i].GetPDG() <<"\n";
353  // make a pair out of tmoth and tgmoth
354  std::pair<int, int> mgmpair(tmoth, tgmoth);
355  // if (tpid>6000) continue;
356  unsigned int planeid = hits[j].GetPlaneID();
357  unsigned int cellid = hits[j].GetCellID();
358  // define the key as a plane cell pair here
359  std::pair<unsigned int, unsigned int> pcpair(planeid, cellid);
360  // this is pair of trackID and particleID
361  std::pair<unsigned int, unsigned int> tppair(tid, tpid);
362  double xyzloc[3];
363  double xyzwor[3];
364  xyzloc[0] = 0.5*(hits[j].GetEntryX()+hits[j].GetExitX());
365  xyzloc[1] = 0.5*(hits[j].GetEntryY()+hits[j].GetExitY());
366  xyzloc[2] = 0.5*(hits[j].GetEntryZ()+hits[j].GetExitZ());
367  fgeo->Plane(planeid)->Cell(cellid)->LocalToWorld(xyzloc,xyzwor);
368  // if (h.GetEdep()>0.001) trhitmap[tid].push_back(& h);
369  // no selection - fill maps with everything
370  trhitmap[tid].push_back(& hits[j]);
371  pcpflsmap[pcpair].push_back(& hits[j]);
372  double end = hits[j].GetEdep();
373  flshitE += end;
374  // this is time, but don't use it yet
375  // double tim = hits[j].fT;
376  // geo::CellUniqueId cid = hits[j].fId;
377  // fill maps which will be used to make SimHits, so don't make any cuts yet
378  // if tpid is 22, a photon, don't use
379 
380  if (end>0.000000000001)
381  {
382  tidpcpmap[tid].push_back(pcpair);
383  // tidmgmmap[tid].push_back(mgmpair);
384  // tidmgmmap[tid]=mgmpair;
385  // pcptidmap[pcpair].push_back(tid);
386  // pcpmotmap[pcpair].push_back(tmoth);
387  // pcppidmap[pcpair].push_back(tpid);
388  // pcpendmap[pcpair].push_back(end);
389  // pcptimmap[pcpair].push_back(tim);
390  }
391 
392  } // end of FLSHit hits loop
393  } // end of FLSHit lists loop
394  fFLSE->Fill(flshitE); // histogram of total FLSHit E
395  fNumFLS->Fill(nFLShits); // histogram of number of FLSHits
396  /*
397  std::cout << " TrackID FLSHit Map size " << trhitmap.size() << "\n";
398  std::cout << " TrackID PlaneCell pair Map size " << tidpcpmap.size() << "\n";
399  std::cout << " TrackID MotherGrandmother pair Map size = " << tidmgmmap.size() << "\n";
400  std::cout << " PlaneCell pair FLSHit Map size " << pcpflsmap.size() << "\n";
401  std::cout << " PlaneCell pair TrackID Map size " << pcptidmap.size() << "\n";
402  std::cout << " PlaneCell pair TID Mother Map size " << pcpmotmap.size() << "\n";
403  std::cout << " PlaneCell pair ParticleID Map size " << pcppidmap.size() << "\n";
404  std::cout << " PlaneCell pair Edep Map size " << pcpendmap.size() << "\n";
405  std::cout << " PlaneCell pair Tdep Map size " << pcptimmap.size() << "\n";
406  */
407  // test some maps
408 
409  // mf::LogInfo("MCCheater") << "Testing PlaneCell pair FLSHit map\n";
410  /*
411  for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
412  {
413  std::pair<unsigned int, unsigned int> planecell = itm->first;
414  std::vector<const sim::FLSHit*> harr = itm->second;
415  std::cout <<" Number of FLSHits = " << harr.size() << " For Plane and Cell = " << planecell.first << " " << planecell.second << "\n";
416  }
417  */
418 
419  // NEED THIS MAP FOR SIMTRACKS
420  // find unique planecells for trackIDs, so can tell whether a track is longer than n cells
421  // make a map linking each trackID with its unique planecells and the size of the value
422  // planecell array is the number of unique planecells on this track
423  std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidupcpmap;
424  for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidpcpmap.begin(); itm != tidpcpmap.end(); ++itm)
425  {
426  unsigned int trackid = itm->first;
427  std::vector<std::pair<unsigned int, unsigned int> > pcarray = itm->second;
428 
429  // eliminate duplicates, resize array with erase
430  // stl unique moves duplicates to the end of the array, doesn't eliminate
431  // erase deletes the items after the new location from unique
432  // this is part of that magic of C++ (R Rechenmacher)
433  std::sort(pcarray.begin(), pcarray.end());
434  std::vector<std::pair<unsigned int, unsigned int> >::iterator endloc;
435  endloc = std::unique(pcarray.begin(), pcarray.end());
436  pcarray.erase(endloc, pcarray.end());
437  for (unsigned int i=0; i<pcarray.size(); ++i)
438  {
439  tidupcpmap[trackid].push_back(pcarray[i]);
440  }
441  }
442 
443 
444  // MAKE SIMHITS HERE
445  // make SimHits using the pc pair, FLSHit map
446  double Etotev = 0.;
447  unsigned int NumSH = 0;
448  std::vector<cheat::SimHit> simhits;
449  for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
450  {
451  // cheat::SimHit* hit = new cheat::SimHit();
453  std::pair<unsigned int, unsigned int> pcp = itm->first; // plane-cell pair of this hit
454  std::vector<const sim::FLSHit*> fh = itm->second; // get the FLSHits in this pc pair
455  unsigned int hsize = fh.size(); // number of FLSHits in this cell
456 
457  // can set plane cell and view for this SimHit now
458  hit.SetPlane(pcp.first); // set plane of this hit
459  unsigned int planen = pcp.first; // can use this later to set fiducial cut on fd
460  hit.SetCell(pcp.second); // set cell of this hit
461  unsigned int celln = pcp.second; // can use this later for fiducial cuts on fd
462  // next 2 lines set view for this plane
463  geo::View_t view = fgeo->Plane(pcp.first)->View();
464  hit.SetView(view); // set view of this hit
465  uint32_t ichan = 1; // have to find this somewhere?
466  hit.SetChannel(ichan);
467 
468  // make the arrays of E, T, PID, and TrID
469  std::vector<double> earr;
470  std::vector<double> carr;
471  std::vector<int> tarr;
472  std::vector<int> parr;
473 
474  std::vector<std::pair<int, int> > tparr; // trackID, PID pair
475  std::vector<std::pair<int, double> > tearr;
476  std::vector<std::pair<int, double> > tcarr;
477  double ene = 0.;
478  double coorx = 0.;
479  double coory = 0.;
480  double coorz = 0.;
481 
482  double_t xxyyzz[3];
483  fgeo->Plane(planen)->Cell(celln)->GetCenter(xxyyzz);
484  coorz = xxyyzz[2];
485  if (view==geo::kX)
486  {
487  coorx = xxyyzz[0];
488  fcorx->Fill(coorx);
489  // need to get y from flshit
490  } else if (view==geo::kY)
491  {
492  coory = xxyyzz[1];
493  fcory->Fill(coory);
494  // need to get x from flshit
495  }
496 
497  for (unsigned int i = 0; i<hsize; ++i)
498  {
499  earr.push_back(fh[i]->GetEdep()); //energy of this hit
500  std::pair<int, double> tep(fh[i]->GetTrackID(), fh[i]->GetEdep()); // pair of trackID, hitE
501  tearr.push_back(tep);
502  carr.push_back(fh[i]->GetEntryT()); // time of this hit
503  std::pair<int, double> ttp(fh[i]->GetTrackID(), fh[i]->GetEntryT()); // pair of trackID, time
504  tcarr.push_back(ttp);
505 
506  tarr.push_back(fh[i]->GetTrackID()); // trackID of this hit
507  // get particle ID from TrackID ParticleID map, not from FLSHit since its different!
508  for (std::map<int, int>::iterator itm = trpidmap.begin(); itm != trpidmap.end(); ++itm)
509  {
510  int trid = itm->first;
511  if (trid == fh[i]->GetTrackID())
512  {
513  parr.push_back(itm->second); // particle ID of this trackID
514  std::pair<int, int> tpp(fh[i]->GetTrackID(), itm->second); // pair of trackID, particleID
515  tparr.push_back(tpp);
516  }
517  }
518  ene += fh[i]->GetEdep(); // energy sum of all contributions to this hit
519 
520  // this uses first FLS contribution to get the non-view coordinate for a SimHit
521  // replace this with an array of the actual values for each contribution
522  if (i==0)
523  {
524  double loc0[3];
525  double wor0[3];
526  double zexo = fh[0]->GetExitZ();
527  double zeno = fh[0]->GetEntryZ();
528  loc0[2] = 0.5*(zexo+zeno);
529  double xexo = fh[0]->GetExitX();
530  double xeno = fh[0]->GetEntryX();
531  loc0[0] = 0.5*(xexo+xeno);
532  double yexo = fh[0]->GetExitY();
533  double yeno = fh[0]->GetEntryY();
534  loc0[1] = 0.5*(yexo+yeno);
535  // change from local to global
536  fgeo->Plane(planen)->Cell(celln)->LocalToWorld(loc0,wor0);
537  if (view==geo::kX) coory = wor0[1];
538  if (view==geo::kY) coorx = wor0[0];
539  }
540  }
541  hit.SetTotE(ene); // set E of this hit
542  ftotalE->Fill(ene * 1000.); // fill this hist also in units of MeV
543  hit.SetPE(ene*5000.); // took a guess to set the number of photoelectrons
544  // hit.SetPECorr(ene*5000.); // took another guess to seet the number of corrected photoelectrons
545  hit.SetCoorX(coorx);
546  hit.SetCoorY(coory);
547  hit.SetCoorZ(coorz);
548 
549  // now need to finish definition based on number of FLSHits in this cell
550  if (hsize == 1) // there is only 1 hit in this cell, can make SimHit directly
551  {
552  hit.SetNP(1); // just 1 here
553  hit.SetEnD(earr);
554  hit.SetTim(carr);
555  hit.SetTrID(tarr);
556  hit.SetPID(parr);
557  }
558  else // size greater than 1 combine if duplicates or mother-daughters
559  {
560  // without sorting, can do with remove - then can adjust other arrays as well
561  std::vector<int>::iterator startloc(tarr.begin());
562  std::vector<int>::iterator endloc(tarr.end());
563  int value(*startloc);
564 
565  while (++startloc < endloc)
566  {
567  endloc = std::remove(startloc, endloc, value);
568  value = *startloc;
569  }
570  tarr.erase(endloc,tarr.end());
571  if (tarr.size()==1)
572  {
573  hit.SetNP(1);
574  hit.SetTrID(tarr);
575  std::vector<int> mparr (1,parr[0]);
576  hit.SetPID(mparr);
577  std::vector<double> mcarr (1,carr[0]);
578  hit.SetTim(mcarr);
579  double esum = 0.;
580  for (unsigned int i=0; i<hsize; ++i)
581  {
582  esum += earr[i];
583  }
584  std::vector<double> mearr (1,esum);
585  hit.SetEnD(mearr);
586  }
587  else // tarr size greater than 1
588  {
589  // if the size of the tarr array is the same as the size of tpp, for now, all hits are
590  // valid and unique and the arrays should be used to make a SimHit
591  if (hsize==tarr.size())
592  {
593  hit.SetNP(hsize);
594  hit.SetPID(parr);
595  hit.SetEnD(earr);
596  hit.SetTrID(tarr);
597  hit.SetTim(carr);
598  }
599  else
600  {
601  // tarr size less than hsize so must analyze
602  // can set the size to tarr size
603  hit.SetNP(tarr.size());
604  hit.SetTrID(tarr);
605 
606  // use pairs to re-generate PID, Energy, time arrays
607  std::vector<int> pvec(tarr.size());
608  std::vector<double> evec(tarr.size());
609  std::vector<double> tvec(tarr.size());
610  for (unsigned int i=0; i<tarr.size(); ++i)
611  {
612  for (unsigned int j=0; j<tparr.size(); ++j)
613  {
614  double lowt = 9999999999.;
615  if (tparr[j].first == tarr[i])
616  {
617  pvec[i] = tparr[j].second;
618  evec[i] += tearr[j].second;
619  if (tcarr[j].second > 0. && tcarr[j].second < lowt)
620  {
621  tvec[i] = tcarr[j].second;
622  lowt = tcarr[j].second;
623  }
624  }
625  }
626  } // end of tarr loop
627  hit.SetPID(pvec);
628  hit.SetEnD(evec);
629  hit.SetTim(tvec);
630  }
631  }
632  }
633 
634  // can put threshold and timing cuts here if wanted before making hit list
635  // 0.0005 is half MeV threshold cut
636  if( hit.fTotE > 0.00000005) // for now, just use a threshold cut
637  {
638  // finally, put this hit into list
639  simhits.push_back(hit);
640  Etotev += hit.fTotE;
641  NumSH++;
642  }
643  }
644 
645  fsimhitE->Fill(Etotev); //fill histogram of total SimHit E in event
646  fnumsh->Fill(NumSH); // fill histogram of number of SimHits in Event
647  std::cout << " Total E (GeV) in SimHits = " << Etotev << "\n";
648 
649  // check some SimHit stuff
650  std::cout << "Size of SimHit List = " << simhits.size() << "\n";
651  // loop over all SimHits to get Info
652  for (unsigned int i=0; i<simhits.size(); ++i)
653  {
654  cheat::SimHit h=simhits[i];
655  // std::cout <<" Num of particles for "<< i << " th hit = " << h.fNP << "\n"
656  // <<" Plane and Cell of this hit = " << h.fPlane << " " << h.fCell << "\n"
657  // <<" View of this hit = " << h.fView << "\n"
658  // <<" Total E of this hit = "<< h.fTotE << "\n";
659  // fill some histograms
660  fthreshE->Fill(h.fTotE * 1000.); // this is energy per hit in MeV
661  if (h.View()==geo::kX) ftotEXV->Fill(h.fTotE * 1000.); // this is energy per hit in MeV per X view
662  if (h.View()==geo::kY) ftotEYV->Fill(h.fTotE * 1000.); // this is energy per hti in MeV per Y view
663  // for (unsigned int j=0; j<h.fNP; ++j)
664  // {
665  // std::cout <<" Particle Edep, ID, TrID, Time = "<< h.fEnD[j] << " " << h.fPID[j] << " " << h.fTrID[j] << " " << h.fTim[j] << "\n";
666  // }
667  }
668 
669  // MAKE SIMTRACKS HERE
670  // to make SimTracks, use TrackID unique PlaneCell pair map - if a TrackID has more than
671  // 2 unique planecell pairs, both views represented, its a SimTrack
672  // will need a map linking trackID to SimHits for SimTracks
673  // std::map<unsigned int, art::Ptr<cheat::SimHit> > simtrhitmap;
674  std::map<unsigned int, std::vector<cheat::SimHit*> > simtrhitmap;
675  std::vector<cheat::SimTrack> simtracks; // initialize SimTracks
676  for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidupcpmap.begin(); itm != tidupcpmap.end(); ++itm)
677  {
678  unsigned int tid = itm->first;
679  unsigned int upcsize = itm->second.size();
680  // if only 1 or 2 unique pc pairs, this is not reconstructable as a track
681  // also eventually probably want to add view requirement for both views
682  if (upcsize>2)
683  {
684  // this is a SimTrack
686  // first, get start and end coords from true particle length using trackID particle map
687  int trpid = 0;
688  // make some variables to hold the trajectory information
689  TVector3 scoord; // Track start point
690  bool foundStart = false;
691  unsigned int startpoint = 0;
692  TVector3 ecoord; // Track end point
693  bool foundEnd = false;
694  unsigned int endpoint = 0;
695  std::vector<TLorentzVector> trpoints; // vector holding the particle trajectory points between the start and end pont.
696 
697 
698  std::map<int, const sim::Particle*>::const_iterator it;
699  for (it = trpartmap.begin(); it != trpartmap.end(); ++it)
700  {
701  unsigned int trid = it->first;
702  const sim::Particle* p=it->second;
703  if (trid==tid)
704  {
705 
706  trpid = p->PdgCode();
707 
708  // Get the starting position of the track as the first trajectory point inside the detector volume
709  unsigned int trajPoint = 0;
710  while(!foundStart && trajPoint<p->NumberTrajectoryPoints()){
711  TLorentzVector start = p->Position(trajPoint);
712  if(TMath::Abs(start.X())<fgeo->DetHalfWidth() && TMath::Abs(start.Y())<fgeo->DetHalfHeight() && start.Z()>0.0 && start.Z()<fgeo->DetLength()){
713  scoord.SetXYZ(start.X(),start.Y(),start.Z());
714  foundStart = true;
715  startpoint = trajPoint;
716  // std::cout<<"Found the start point: "<<startpoint<<std::endl;
717  }
718  else{
719  ++trajPoint;
720  }
721  }
722  //initiate endpoint to same as startpoint
723  endpoint = startpoint;
724  trajPoint = p->NumberTrajectoryPoints()-1;
725  // Get the endpoint as the last trajectory point that is inside the detector volume
726  while(!foundEnd && trajPoint>=startpoint){
727  TLorentzVector end = p->Position(trajPoint);
728  if(TMath::Abs(end.X())<fgeo->DetHalfWidth() && TMath::Abs(end.Y())<fgeo->DetHalfHeight() && end.Z()>0.0 && end.Z()<fgeo->DetLength()){
729  ecoord.SetXYZ(end.X(),end.Y(),end.Z());
730  foundEnd = true;
731  endpoint = trajPoint;
732  // std::cout<<"Found the end point: "<<endpoint<<std::endl;
733  }
734  else{
735  --trajPoint;
736  }
737  }
738  // Save all the trajectory points after start point up to the end point.
739  if(foundStart && foundEnd && endpoint>startpoint){
740  for(unsigned int trpoint = startpoint+1; trpoint <= endpoint;++trpoint){
741  // std::cout<<"Start point: "<<startpoint<<" End point: "<<endpoint<<std::endl;
742  // std::cout<<"Putting trajectory point into track list: "<<trpoint<<std::endl;
743  trpoints.push_back(p->Position(trpoint));
744  }
745  }
746 
747  }
748  }
749  // now only make a simtrack object if the track start and stop where found
750  if(foundStart && foundEnd && endpoint>=startpoint){
752 
753  // calculate direction cosines from start and end of this MC track
754  strack.SetPID(trpid);
755  strack.SetTrackID(tid);
756  strack.SetStart(scoord);
757  for(unsigned int trpoint = 0; trpoint < trpoints.size();++trpoint){
758  TVector3 pos(trpoints[trpoint].X(),trpoints[trpoint].Y(),trpoints[trpoint].Z());
759  strack.AppendTrajectoryPoint(pos);
760  // Get initial track direction
761  if(trpoint == 0){
762  // calculate direction cosines from the first and second trajectory points of this MC track
763 
764  double dist = sqrt((ecoord.X()-trpoints[trpoint].X())*(ecoord.X()-trpoints[trpoint].X())+
765  (ecoord.Y()-trpoints[trpoint].Y())*(ecoord.Y()-trpoints[trpoint].Y())+
766  (ecoord.Z()-trpoints[trpoint].Z())*(ecoord.Z()-trpoints[trpoint].Z()));
767  TVector3 dir((ecoord.X()-trpoints[trpoint].X())/dist,(ecoord.Y()-trpoints[trpoint].Y())/dist,(ecoord.Z()-trpoints[trpoint].Z())/dist);
768  strack.SetDir(dir);
769 
770  }
771  }
772  // add to SimTracks
773  simtracks.push_back(strack);
774  }
775 
776  // loop over SimHits, adding to the TrackID, SimHit map
777  // simtrhitmap[tid].resize(0); // Jon suggested this to fix map?
778  for (unsigned int i=0; i<simhits.size(); ++i)
779  {
780  bool ismatch = false;
781  // const cheat::SimHit h=simhits[i];
782  cheat::SimHit h=simhits[i];
783  for (unsigned int j=0; j<h.fNP; ++j)
784  {
785  if (h.fTrID[j]==(int) tid) ismatch = true;
786  }
787  if (ismatch) simtrhitmap[tid].push_back(& h);
788  }
789  } // end of upc greater than 2 loop
790  } // end of loop over trackID planecell pair map
791 
792  // write SimTracks to event
793  std::unique_ptr<std::vector<cheat::SimTrack> > simtrackcol(new std::vector<cheat::SimTrack>(simtracks) );
794  evt.put(std::move(simtrackcol));
795 
796  /*
797  // ********************** Diagnostics for MCCheater ***********************
798  // setup a test for single particles - if trajectory along z axis, test variable is true
799  bool zdir = false;
800  // report some SimTrack stuff
801  std::cout << "Size of SimTrack List = " << simtracks.size() << "\n";
802  double tlensum = 0.;
803  // define min and max coords for NDOS
804  double xmin = -130.1;
805  double xmax = 130.1;
806  double ymin = -195.1;
807  double ymax = 195.1;
808  double zmin = 0.;
809  double zmax = 1200.;
810  for (unsigned int i=0; i<simtracks.size(); ++i)
811  {
812  fzstart->Fill(simtracks[i].StartPos()[2]/6.0);
813  double xlen = simtracks[i].EndPos()[0]-simtracks[i].StartPos()[0];
814  double ylen = simtracks[i].EndPos()[1]-simtracks[i].StartPos()[1];
815  double zlen = simtracks[i].EndPos()[2]-simtracks[i].StartPos()[2];
816  fnplanes->Fill(zlen/6.);
817  double tlength = sqrt(xlen*xlen+ylen*ylen+zlen*zlen);
818  tlensum += tlength;
819  fdcosX->Fill(simtracks[i].Dir()[0]);
820  fdcosY->Fill(simtracks[i].Dir()[1]);
821  fdcosZ->Fill(simtracks[i].Dir()[2]);
822  if (simtracks[i].Dir()[2]>0.5) zdir = true;
823  if (zdir && simtracks[i].TrackID()==0) ftlen->Fill(tlength);
824  if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==13) fmuL->Fill(tlength);
825  if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==11) felL->Fill(tlength);
826  std::cout << " Track Number " << i << "\n";
827  std::cout << " Track ID = " << simtracks[i].TrackID() << " Pdg ID = " << simtracks[i].PID() << "\n";
828  std::cout << " Start Coord = " << simtracks[i].StartPos()[0] << " " << simtracks[i].StartPos()[1] << " " << simtracks[i].StartPos()[2] << "\n";
829  std::cout << " End Coord = " << simtracks[i].EndPos()[0] << " " << simtracks[i].EndPos()[1] << " " << simtracks[i].EndPos()[2] << "\n";
830  std::cout << " Direction = " << simtracks[i].Dir()[0] << " " << simtracks[i].Dir()[1] << " " << simtracks[i].Dir()[2] << "\n";
831  // find point of intersection with plane of NDOS
832  // first, use top plane
833  double ttop = (195.-simtracks[i].StartPos()[1])/ylen;
834  double xint = simtracks[i].StartPos()[0]+ttop*(xlen);
835  double yint = simtracks[i].StartPos()[1]+ttop*(ylen);
836  double zint = simtracks[i].StartPos()[2]+ttop*(zlen);
837  std::cout << " TP ttrack = " << ttop << " Coords = " << xint << " " << yint << " " << zint << "\n";
838  if (xint>xmin && xint<xmax && yint>ymin && yint<ymax && zint>zmin && zint<zmax)
839  {
840  // loop over all simhits to find closest to this point
841  double sdist = 999.;
842  for (unsigned int j=0; j<simhits.size(); ++j)
843  {
844  cheat::SimHit h=simhits[j];
845  // does this hit contain contribution from this track?
846  bool htid = false;
847  for (unsigned int k=0; k<h.fNP; ++k)
848  {
849  if (h.fTrID[k]==simtracks[i].TrackID()) htid = true;
850  }
851  if (htid)
852  {
853  double hdist = sqrt((h.fCoorX-xint)*(h.fCoorX-xint)+(h.fCoorY-yint)*(h.fCoorY-yint)+(h.fCoorZ-zint)*(h.fCoorZ-zint));
854  // std::cout << " Hit distance = " << hdist << "\n";
855  if (hdist<sdist) sdist = hdist;
856  } // satisfies trackID
857  } // end of SimHit loop
858  std::cout << " Closest distance to edge = " << sdist << " for TrackID = " << simtracks[i].TrackID() << "\n";
859  fclhit->Fill(sdist);
860  }
861  }
862  */
863 
864  /*
865  double E0SHTot = 0.;
866  if (zdir)
867  {
868  for (unsigned int i=0; i<simhits.size(); ++i)
869  {
870  // cheat::SimHit hit=simhits[i];
871  E0SHTot += simhits[i].fTotE;
872  }
873  fEpar->Fill(E0SHTot); // histogram of forward-going energy sum
874  }
875  */
876 
877  // check map of SimTrack TrackID SimHit
878  // this map doesn't work?
879  /*
880  for (std::map<unsigned int, std::vector<cheat::SimHit*> >::iterator itm = simtrhitmap.begin(); itm != simtrhitmap.end(); ++itm)
881  {
882  std::vector<cheat::SimHit*> sharr = itm->second;
883  mf::LogInfo("MCCheater") <<" Number of SimHits = " << sharr.size() << " For SimTrackID = " << itm->first << "\n";
884  fnshptr->Fill(sharr.size());
885  double stc[3];
886  double enc[3];
887  stc[1] = -210.;
888  enc[1] = 210.;
889 
890  for (unsigned int i=0; i<sharr.size(); ++i)
891  {
892  cheat::SimHit* h=sharr[i];
893  double xcor = h->fCoorX;
894  double ycor = h->fCoorY;
895  double zcor = h->fCoorZ;
896  // histogram these coordinates
897  fXRange->Fill(xcor);
898  fYRange->Fill(ycor);
899  fZRange->Fill(zcor);
900 
901  if (ycor > stc[1])
902  {
903  stc[1] = ycor;
904  stc[2] = zcor;
905  stc[0] = xcor;
906  continue;
907  } else if (ycor < stc[1])
908  {
909  if (ycor < enc[1])
910  {
911  enc[1] = ycor;
912  enc[2] = zcor;
913  enc[0] = xcor;
914  continue;
915  }
916  }
917  }
918  fstsxy->Fill(stc[0],stc[1]);
919  fstexy->Fill(enc[0],enc[1]);
920  }
921  */
922 
923  // write SimHits to event
924  std::unique_ptr<std::vector<cheat::SimHit> > simhitcol(new std::vector<cheat::SimHit>(simhits) );
925  std::cout << " Wrote SimHits to event of size = " << simhits.size() << "\n";
926  evt.put(std::move(simhitcol));
927 
928  return;
929  }
930 
932 } // end namespace cheat
933 ////////////////////////////////////////////////////////////////////////
double E(const int i=0) const
Definition: MCParticle.h:232
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
virtual void SetStart(TVector3 start)
Definition: Track.cxx:168
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
void produce(art::Event &evt)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
set< int >::iterator it
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:748
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
int Mother() const
Definition: MCParticle.h:212
void SetPID(int pid)
Definition: SimTrack.h:28
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double DetLength() const
void SetChannel(uint32_t iChan)
Definition: RawDigit.h:99
void abs(TH1 *hist)
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetTotE(double tote)
Definition: SimHit.h:50
void SetTrackID(int trackid)
Definition: SimTrack.h:29
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
static const unsigned int npar
int TrackId() const
Definition: MCParticle.h:209
Float_t Y
Definition: plot.C:38
void SetView(geo::View_t view)
Definition: CellHit.h:54
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
double dist
Definition: runWimpSim.h:113
Float_t Z
Definition: plot.C:38
void SetCell(unsigned short cell)
Definition: CellHit.h:52
TString part[npart]
Definition: Style.C:32
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void hits()
Definition: readHits.C:15
void SetPE(float pe)
Definition: CellHit.h:55
MCCheater(fhicl::ParameterSet const &pset)
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetNP(unsigned short np)
Definition: SimHit.h:51
void SetCoorZ(double coorz)
Definition: SimHit.h:49
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
Float_t strack
Definition: plot.C:35
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
void SetPID(std::vector< int > PID)
Definition: SimHit.h:53
void SetTim(std::vector< double > Tim)
Definition: SimHit.h:55
int evt
const double j
Definition: BetheBloch.cxx:29
std::string fMCTruthListLabel
double DetHalfHeight() const
void SetTrID(std::vector< int > TrID)
Definition: SimHit.h:52
OStream cout
Definition: OStream.cxx:6
void SetCoorX(double coorx)
Definition: SimHit.h:47
void SetEnD(std::vector< double > EnD)
Definition: SimHit.h:54
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:17
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
std::string fFLSHitListLabel
Definition: structs.h:12
double fTotE
Total Energy dep from all particles.
Definition: SimHit.h:65
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
Float_t X
Definition: plot.C:38
Definition: fwd.h:29
void SetCoorY(double coory)
Definition: SimHit.h:48
unsigned short fNP
size of vectors
Definition: SimHit.h:66
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< int > fTrID
vector of sim track IDs
Definition: SimHit.h:67
int EveId(const int trackID) const
enum BeamMode string