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