Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventSourceHDDM.cc
Go to the documentation of this file.
1 // $Id: DEventSourceHDDM.cc 19023 2015-07-14 20:23:27Z beattite $
2 //
3 // Author: David Lawrence June 24, 2004
4 //
5 // changes: Wed Jun 20 17:08:13 EDT 2007 B. Zihlmann
6 // modify TOF section to add several new variables incuding the
7 // GEANT particle type to the Truth hits and the hit and track-hit
8 // list.
9 //
10 // Oct 3, 2012 Yi Qiang: add functions for Cherenkov RICH detector
11 // Oct 11, 2012 Yi Qiang: complete functions for Cherenkov detector
12 // Oct 8, 2013 Yi Qiang: added dedicated object for RICH Truth Hit
13 // July 5, 2014 R.T.Jones: changed over from c to c++ API for hddm
14 // June 22, 2015 J. Stevens: changed RICH -> DIRC and remove CERE
15 // May 7, 2017 R. Dzhygadlo: added DDIRCTruthPmtHit DDIRCTruthBarHit
16 // Oct 20, 2017 A. Somov: Added fields for the DPSHit/DPSCHit
17 //
18 // DEventSourceHDDM methods
19 //
20 
21 #include <iostream>
22 #include <iomanip>
23 #include <cmath>
24 using namespace std;
25 
26 #include <JANA/JFactory_base.h>
27 #include <JANA/JEventLoop.h>
28 #include <JANA/JEvent.h>
29 #include <DANA/DStatusBits.h>
30 
31 #include <JANA/JGeometryXML.h>
32 #include "BCAL/DBCALGeometry.h"
34 
35 #include <DVector2.h>
36 #include <DEventSourceHDDM.h>
37 #include <FDC/DFDCGeometry.h>
38 #include <FCAL/DFCALGeometry.h>
39 #include <FCAL/DFCALHit.h>
40 #include <CCAL/DCCALGeometry.h>
41 #include <CCAL/DCCALHit.h>
42 
43 
44 //------------------------------------------------------------------
45 // Binary predicate used to sort hits
46 //------------------------------------------------------------------
48  public:
49  bool operator()(DMCTrackHit* const &thit1,
50  DMCTrackHit* const &thit2) const {
51  return thit1->z < thit2->z;
52  }
53 };
54 
55 bool MCTrackHitSort_C(DMCTrackHit* const &thit1,
56  DMCTrackHit* const &thit2) {
57  return thit1->z < thit2->z;
58 }
59 
60 
61 //----------------
62 // Constructor
63 //----------------
64 DEventSourceHDDM::DEventSourceHDDM(const char* source_name)
65 : JEventSource(source_name)
66 {
67  /// Constructor for DEventSourceHDDM object
68  ifs = new ifstream(source_name);
69  if (ifs->is_open())
70  fin = new hddm_s::istream(*ifs);
71  else
72  fin = NULL;
73  initialized = false;
74  dapp = NULL;
75  bfield = NULL;
76  geom = NULL;
77 
78  dRunNumber = -1;
79 
80  if( (!gPARMS->Exists("JANA_CALIB_CONTEXT")) && (getenv("JANA_CALIB_CONTEXT")==NULL) ){
81  cout << "============================================================" << endl;
82  cout << " WARNING: JANA_CALIB_CONTEXT not set. " << endl;
83  cout << "You are reading from an HDDM file which is most likely" << endl;
84  cout << "MC data. In most cases, you will want to set this parameter" << endl;
85  cout << "to get proper reconstruction." << endl;
86  cout << "(usually something like \"variation=mc\")" << endl;
87  cout << "============================================================" << endl;
88  }
89 }
90 
91 //----------------
92 // Destructor
93 //----------------
95 {
96  if (fin)
97  delete fin;
98  if (ifs)
99  delete ifs;
100 }
101 
102 //----------------
103 // GetEvent
104 //----------------
105 jerror_t DEventSourceHDDM::GetEvent(JEvent &event)
106 {
107  /// Implementation of JEventSource virtual function
108 
109  if (!fin)
110  return EVENT_SOURCE_NOT_OPEN;
111 
112  // Each open HDDM file takes up about 1M of memory so it's
113  // worthwhile to close it as soon as we can.
114  else if (!ifs->good()) {
115  delete fin;
116  fin = NULL;
117  delete ifs;
118  ifs = NULL;
119  return NO_MORE_EVENTS_IN_SOURCE;
120  }
121 
122  hddm_s::HDDM *record = new hddm_s::HDDM();
123  if (! (*fin >> *record)) {
124  delete fin;
125  fin = NULL;
126  delete ifs;
127  ifs = NULL;
128  return NO_MORE_EVENTS_IN_SOURCE;
129  }
130 
131  ++Nevents_read;
132 
133  int event_number = -1;
134  int run_number = -1;
135 
136  // Get event/run numbers from HDDM
137  hddm_s::PhysicsEvent &pe = record->getPhysicsEvent(0);
138  event_number = pe.getEventNo();
139  run_number = pe.getRunNo();
140 
141  // Copy the reference info into the JEvent object
142  event.SetJEventSource(this);
143  event.SetEventNumber(event_number);
144  event.SetRunNumber(run_number);
145  event.SetRef(record);
146  event.SetStatusBit(kSTATUS_HDDM);
147  event.SetStatusBit(kSTATUS_FROM_FILE);
148  event.SetStatusBit(kSTATUS_PHYSICS_EVENT);
149 
150  return NOERROR;
151 }
152 
153 //----------------
154 // FreeEvent
155 //----------------
156 void DEventSourceHDDM::FreeEvent(JEvent &event)
157 {
158  hddm_s::HDDM *record = (hddm_s::HDDM*)event.GetRef();
159  delete record;
160 }
161 
162 //----------------
163 // GetObjects
164 //----------------
165 jerror_t DEventSourceHDDM::GetObjects(JEvent &event, JFactory_base *factory)
166 {
167  /// This gets called through the virtual method of the
168  /// JEventSource base class. It creates the objects of the type
169  /// on which factory is based. It uses the hddm_s::HDDM* object
170  /// kept in the ref field of the JEvent object passed.
171 
172  // We must have a factory to hold the data
173  if (!factory)
174  throw RESOURCE_UNAVAILABLE;
175 
176  // HDDM doesn't exactly support tagged factories, but the tag
177  // can be used to direct filling of the correct factory.
178  string tag = (factory->Tag()==NULL)? "" : factory->Tag();
179 
180  // The ref field of the JEvent is just the HDDM object pointer.
181  hddm_s::HDDM *record = (hddm_s::HDDM*)event.GetRef();
182  if (!record)
183  throw RESOURCE_UNAVAILABLE;
184 
185  // Get pointer to the B-field object and Geometry object
186  JEventLoop *loop = event.GetJEventLoop();
187  if (initialized == false && loop) {
188  initialized = true;
189  dRunNumber = event.GetRunNumber();
190  dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
191  if (dapp) {
192  jcalib = dapp->GetJCalibration(event.GetRunNumber());
193  // Make sure jcalib is set
194  if (!jcalib) {
195  _DBG_ << "ERROR - no jcalib set!" <<endl;
196  return RESOURCE_UNAVAILABLE;
197  }
198  // Get constants and do basic check on number of elements
199  vector< map<string, float> > tvals;
200  if(jcalib->Get("FDC/strip_calib", tvals))
201  throw JException("Could not load CCDB table: FDC/strip_calib");
202 
203  if (tvals.size() != 192) {
204  _DBG_ << "ERROR - strip calibration vectors are not the right size!"
205  << endl;
206  return VALUE_OUT_OF_RANGE;
207  }
208  map<string,float>::iterator iter;
209  for (iter=tvals[0].begin(); iter!=tvals[0].end(); iter++) {
210  // Copy values into tables. We preserve the order since
211  // that is how it was originally done in hitFDC.c
212  for (unsigned int i=0; i<tvals.size(); i++) {
213  map<string, float> &row = tvals[i];
214  uscale[i]=row["qru"];
215  vscale[i]=row["qrv"];
216  }
217  }
218  }
219  // load BCAL geometry
220  vector<const DBCALGeometry *> BCALGeomVec;
221  loop->Get(BCALGeomVec);
222  if(BCALGeomVec.size() == 0)
223  throw JException("Could not load DBCALGeometry object!");
224  dBCALGeom = BCALGeomVec[0];
225 
226  // load PS geometry
227  vector<const DPSGeometry*> psGeomVect;
228  loop->Get(psGeomVect);
229  if (psGeomVect.size() < 1)
230  return OBJECT_NOT_AVAILABLE;
231  psGeom = psGeomVect[0];
232 
233 
234  }
235 
236  // Warning: This class is not completely thread-safe and can fail if running
237  // running in multithreaded mode over files with events from multiple runs
238  // It is expected that simulated data will rarely contain events from multiple
239  // runs, as this is an intermediate format in the simulation chain, so for
240  // now we just insert a sanity check, and push the problem to the future
241  if(dRunNumber != event.GetRunNumber()) {
242  jerr << endl
243  << "WARNING: DEventSourceHDDM cannot currently handle HDDM files containing" << endl
244  << "events with multiple runs! If you encounter this error message," << endl
245  << "please contact the GlueX Offline Software Group: halld-offline@jlab.org" << endl
246  << endl;
247  exit(-1);
248  }
249 
250  //Get target center
251  //multiple reader threads can access this object: need lock
252  bool locNewRunNumber = false;
253  unsigned int locRunNumber = event.GetRunNumber();
254  LockRead();
255  {
256  locNewRunNumber = (dTargetCenterZMap.find(locRunNumber) == dTargetCenterZMap.end());
257  }
258  UnlockRead();
259  if(locNewRunNumber)
260  {
261  DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
262  DGeometry* locGeometry = dapp->GetDGeometry(loop->GetJEvent().GetRunNumber());
263  double locTargetCenterZ = 0.0;
264  locGeometry->GetTargetZ(locTargetCenterZ);
265 
266  JGeometryXML *jgeom = dynamic_cast<JGeometryXML *>(locGeometry);
267  hddm_s::GeometryList geolist = record->getGeometrys();
268  if (jgeom != 0 && geolist.size() > 0) {
269  std::string md5sim = geolist(0).getMd5simulation();
270  std::string md5smear = geolist(0).getMd5smear();
271  std::string md5recon = jgeom->GetChecksum();
272  geolist(0).setMd5reconstruction(md5recon);
273  if (md5sim != md5smear) {
274  jerr << std::endl
275  << "WARNING: simulation geometry checksum does not match"
276  << " that shown for the mcsmear step."
277  << std::endl;
278  }
279  else if (md5sim != md5recon) {
280  jerr << endl
281  << "WARNING: simulation geometry checksum does not match"
282  << " the geometry being used for reconstruction."
283  << std::endl;
284  }
285  }
286 
287  vector<double> locBeamPeriodVector;
288  if(loop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector))
289  throw runtime_error("Could not load CCDB table: PHOTON_BEAM/RF/beam_period");
290  double locBeamBunchPeriod = locBeamPeriodVector[0];
291 
292  LockRead();
293  {
294  dTargetCenterZMap[locRunNumber] = locTargetCenterZ;
295  dBeamBunchPeriodMap[locRunNumber] = locBeamBunchPeriod;
296  }
297  UnlockRead();
298  }
299 
300  // Get name of data class we're trying to extract
301  string dataClassName = factory->GetDataClassName();
302 
303  if (dataClassName == "DPSHit")
304  return Extract_DPSHit(record,
305  dynamic_cast<JFactory<DPSHit>*>(factory), tag);
306 
307  if (dataClassName == "DPSTruthHit")
308  return Extract_DPSTruthHit(record,
309  dynamic_cast<JFactory<DPSTruthHit>*>(factory), tag);
310 
311  if (dataClassName == "DPSCHit")
312  return Extract_DPSCHit(record,
313  dynamic_cast<JFactory<DPSCHit>*>(factory), tag);
314 
315  if (dataClassName == "DPSCTruthHit")
316  return Extract_DPSCTruthHit(record,
317  dynamic_cast<JFactory<DPSCTruthHit>*>(factory), tag);
318 
319  if (dataClassName == "DRFTime")
320  return Extract_DRFTime(record,
321  dynamic_cast<JFactory<DRFTime>*>(factory), loop);
322 
323  if (dataClassName == "DTAGMHit")
324  return Extract_DTAGMHit(record,
325  dynamic_cast<JFactory<DTAGMHit>*>(factory), tag);
326 
327  if (dataClassName == "DTAGHHit")
328  return Extract_DTAGHHit(record,
329  dynamic_cast<JFactory<DTAGHHit>*>(factory), tag);
330 
331  if (dataClassName == "DMCTrackHit")
332  return Extract_DMCTrackHit(record,
333  dynamic_cast<JFactory<DMCTrackHit>*>(factory), tag);
334 
335  if (dataClassName == "DMCReaction")
336  return Extract_DMCReaction(record,
337  dynamic_cast<JFactory<DMCReaction>*>(factory), tag, loop);
338 
339  if (dataClassName == "DMCThrown")
340  return Extract_DMCThrown(record,
341  dynamic_cast<JFactory<DMCThrown>*>(factory), tag);
342 
343  if (dataClassName == "DBCALTruthShower")
344  return Extract_DBCALTruthShower(record,
345  dynamic_cast<JFactory<DBCALTruthShower>*>(factory), tag);
346 
347  if (dataClassName == "DBCALSiPMSpectrum")
348  return Extract_DBCALSiPMSpectrum(record,
349  dynamic_cast<JFactory<DBCALSiPMSpectrum>*>(factory), tag);
350 
351  if (dataClassName == "DBCALTruthCell")
352  return Extract_DBCALTruthCell(record,
353  dynamic_cast<JFactory<DBCALTruthCell>*>(factory), tag);
354 
355  if (dataClassName == "DBCALSiPMHit")
356  return Extract_DBCALSiPMHit(record,
357  dynamic_cast<JFactory<DBCALSiPMHit>*>(factory), tag);
358 
359  if (dataClassName == "DBCALDigiHit")
360  return Extract_DBCALDigiHit(record,
361  dynamic_cast<JFactory<DBCALDigiHit>*>(factory), tag);
362 
363  if (dataClassName == "DBCALIncidentParticle")
364  return Extract_DBCALIncidentParticle(record,
365  dynamic_cast<JFactory<DBCALIncidentParticle>*>(factory), tag);
366 
367  if (dataClassName == "DBCALTDCDigiHit")
368  return Extract_DBCALTDCDigiHit(record,
369  dynamic_cast<JFactory<DBCALTDCDigiHit>*>(factory), tag);
370 
371  if (dataClassName == "DCDCHit")
372  return Extract_DCDCHit(loop, record,
373  dynamic_cast<JFactory<DCDCHit>*>(factory) , tag);
374 
375  if (dataClassName == "DFDCHit")
376  return Extract_DFDCHit(record,
377  dynamic_cast<JFactory<DFDCHit>*>(factory), tag);
378 
379  if (dataClassName == "DFCALTruthShower")
380  return Extract_DFCALTruthShower(record,
381  dynamic_cast<JFactory<DFCALTruthShower>*>(factory), tag);
382 
383  if (dataClassName == "DFCALHit")
384  return Extract_DFCALHit(record,
385  dynamic_cast<JFactory<DFCALHit>*>(factory), tag,
386  event.GetJEventLoop());
387 
388  if (dataClassName == "DCCALTruthShower")
389  return Extract_DCCALTruthShower(record,
390  dynamic_cast<JFactory<DCCALTruthShower>*>(factory), tag);
391 
392  if (dataClassName == "DCCALHit")
393  return Extract_DCCALHit(record,
394  dynamic_cast<JFactory<DCCALHit>*>(factory), tag,
395  event.GetJEventLoop());
396 
397  if (dataClassName == "DMCTrajectoryPoint" && tag == "")
398  return Extract_DMCTrajectoryPoint(record,
399  dynamic_cast<JFactory<DMCTrajectoryPoint>*>(factory), tag);
400 
401  if (dataClassName == "DTOFTruth")
402  return Extract_DTOFTruth(record,
403  dynamic_cast<JFactory<DTOFTruth>*>(factory), tag);
404 
405  // TOF is a special case: TWO factories are needed at the same time
406  // DTOFHit and DTOFHitMC
407  if (dataClassName == "DTOFHit") {
408  JFactory_base* factory2 = loop->GetFactory("DTOFHitMC", tag.c_str());
409  return Extract_DTOFHit(record,
410  dynamic_cast<JFactory<DTOFHit>*>(factory),
411  dynamic_cast<JFactory<DTOFHitMC>*>(factory2), tag);
412  }
413  if (dataClassName == "DTOFHitMC") {
414  JFactory_base* factory2 = loop->GetFactory("DTOFHit", tag.c_str());
415  return Extract_DTOFHit(record,
416  dynamic_cast<JFactory<DTOFHit>*>(factory2),
417  dynamic_cast<JFactory<DTOFHitMC>*>(factory), tag);
418  }
419 
420  if (dataClassName == "DSCHit")
421  return Extract_DSCHit(record,
422  dynamic_cast<JFactory<DSCHit>*>(factory), tag);
423 
424  if (dataClassName == "DSCTruthHit")
425  return Extract_DSCTruthHit(record,
426  dynamic_cast<JFactory<DSCTruthHit>*>(factory), tag);
427 
428  if (dataClassName == "DFMWPCTruthHit")
429  return Extract_DFMWPCTruthHit(record,
430  dynamic_cast<JFactory<DFMWPCTruthHit>*>(factory), tag);
431 
432  if (dataClassName == "DFMWPCHit")
433  return Extract_DFMWPCHit(record,
434  dynamic_cast<JFactory<DFMWPCHit>*>(factory), tag);
435 
436  if (dataClassName == "DDIRCTruthBarHit")
437  return Extract_DDIRCTruthBarHit(record,
438  dynamic_cast<JFactory<DDIRCTruthBarHit>*>(factory), tag);
439 
440  if (dataClassName == "DDIRCTruthPmtHit")
441  return Extract_DDIRCTruthPmtHit(record,
442  dynamic_cast<JFactory<DDIRCTruthPmtHit>*>(factory), tag);
443 
444  if (dataClassName == "DDIRCPmtHit")
445  return Extract_DDIRCPmtHit(record,
446  dynamic_cast<JFactory<DDIRCPmtHit>*>(factory), tag, event.GetJEventLoop());
447 
448  // extract CereTruth and CereRichHit hits, yqiang Oct 3, 2012
449  // removed CereTruth (merged into MCThrown), added CereHit, yqiang Oct 10 2012
450  if (dataClassName == "DCereHit")
451  return Extract_DCereHit(record,
452  dynamic_cast<JFactory<DCereHit>*>(factory), tag);
453 
454  if (dataClassName == "DTPOLHit")
455  return Extract_DTPOLHit(record,
456  dynamic_cast<JFactory<DTPOLHit>*>(factory), tag);
457 
458  if (dataClassName == "DTPOLTruthHit")
459  return Extract_DTPOLTruthHit(record,
460  dynamic_cast<JFactory<DTPOLTruthHit>*>(factory), tag);
461 
462  return OBJECT_NOT_AVAILABLE;
463 }
464 
465 //------------------
466 // Extract_DRFTime
467 //------------------
468 jerror_t DEventSourceHDDM::Extract_DRFTime(hddm_s::HDDM *record,
469  JFactory<DRFTime> *factory, JEventLoop* locEventLoop)
470 {
471  if (factory==NULL)
472  return OBJECT_NOT_AVAILABLE;
473  string tag = (factory->Tag())? factory->Tag() : "";
474 
475  vector<DRFTime*> locRFTimes;
476 
477  // loop over RF-time records
478  const hddm_s::RFtimeList &rftimes = record->getRFtimes();
479  hddm_s::RFtimeList::iterator iter;
480  for (iter = rftimes.begin(); iter != rftimes.end(); ++iter)
481  {
482  if (iter->getJtag() != tag)
483  continue;
484  DRFTime *locRFTime = new DRFTime;
485  locRFTime->dTime = iter->getTsync();
486  locRFTime->dTimeVariance = 0.0; //SET ME!!
487  locRFTimes.push_back(locRFTime);
488  }
489 
490  if(!locRFTimes.empty())
491  {
492  //found in the file, copy into factory and return
493  factory->CopyTo(locRFTimes);
494  return NOERROR;
495  }
496 
497  //Not found in the file, so either:
498  //Experimental data & it's missing: bail
499  //MC data: generate it
500 
501  vector<const DBeamPhoton*> locMCGENPhotons;
502  locEventLoop->Get(locMCGENPhotons, "MCGEN");
503  if(locMCGENPhotons.empty())
504  return OBJECT_NOT_AVAILABLE; //Experimental data & it's missing: bail
505 
506  //Is MC data. Either:
507  //No tag: return t = 0.0, but true t is 0.0 +/- n*locBeamBunchPeriod: must select the correct beam bunch
508  //TRUTH tag: get exact t from DBeamPhoton tag MCGEN
509 
510  if(tag == "TRUTH")
511  {
512  DRFTime *locRFTime = new DRFTime;
513  locRFTime->dTime = locMCGENPhotons[0]->time();
514  locRFTime->dTimeVariance = 0.0;
515  locRFTimes.push_back(locRFTime);
516  }
517  else
518  {
519  double locBeamBunchPeriod = 0.0;
520  int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
521  LockRead();
522  {
523  locBeamBunchPeriod = dBeamBunchPeriodMap[locRunNumber];
524  }
525  UnlockRead();
526 
527  //start with true RF time, increment/decrement by multiples of locBeamBunchPeriod ns until closest to 0
528  double locTime = locMCGENPhotons[0]->time();
529  int locNumRFBuckets = int(locTime/locBeamBunchPeriod);
530  locTime -= double(locNumRFBuckets)*locBeamBunchPeriod;
531  while(locTime > 0.5*locBeamBunchPeriod)
532  locTime -= locBeamBunchPeriod;
533  while(locTime < -0.5*locBeamBunchPeriod)
534  locTime += locBeamBunchPeriod;
535 
536  DRFTime *locRFTime = new DRFTime;
537  locRFTime->dTime = locTime;
538  locRFTime->dTimeVariance = 0.0;
539  locRFTimes.push_back(locRFTime);
540  }
541 
542  // Copy into factories
543  factory->CopyTo(locRFTimes);
544 
545  return NOERROR;
546 }
547 
548 //------------------
549 // Extract_DMCTrackHit
550 //------------------
551 jerror_t DEventSourceHDDM::Extract_DMCTrackHit(hddm_s::HDDM *record,
552  JFactory<DMCTrackHit> *factory, string tag)
553 {
554  /// Copies the data from the given hddm_s structure. This is called
555  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
556  /// returns OBJECT_NOT_AVAILABLE immediately.
557 
558  if (factory == NULL)
559  return OBJECT_NOT_AVAILABLE;
560  if (tag != "")
561  return OBJECT_NOT_AVAILABLE;
562 
563  // The following routines will create DMCTrackHit objects and add them
564  // to data.
565  vector<DMCTrackHit*> data;
566  GetCDCTruthHits(record, data);
567  GetFDCTruthHits(record, data);
568  GetBCALTruthHits(record, data);
569  GetTOFTruthHits(record, data);
570  GetCherenkovTruthHits(record, data);
571  GetFCALTruthHits(record, data);
572  GetSCTruthHits(record, data);
573 
574  // It has happened that some CDC hits have "nan" for the drift time
575  // in a peculiar event Alex Somov came across. This ultimately caused
576  // a seg. fault in MCTrackHitSort_C. I hate doing this since it
577  // is treating the symptom rather than the cause, but nonetheless,
578  // it patches up the problem for now until there is time to revisit
579  // it later.
580  for (unsigned int i=0; i < data.size(); i++)
581  if (!isfinite(data[i]->z))
582  data[i]->z = -1000.0;
583 
584  // sort hits by z
585  sort(data.begin(), data.end(), MCTrackHitSort_C);
586 
587  // Some systems will use negative phis. Force them all to
588  // be in the 0 to 2pi range
589  for (unsigned int i=0; i < data.size(); i++) {
590  DMCTrackHit *mctrackhit = data[i];
591  if (mctrackhit->phi < 0.0)
592  mctrackhit->phi += 2.0*M_PI;
593  }
594 
595  // Copy into factory
596  factory->CopyTo(data);
597 
598  return NOERROR;
599 }
600 
601 //-------------------
602 // GetCDCTruthHits
603 //-------------------
604 jerror_t DEventSourceHDDM::GetCDCTruthHits(hddm_s::HDDM *record,
605  vector<DMCTrackHit*>& data)
606 {
607  const hddm_s::CdcTruthPointList &points = record->getCdcTruthPoints();
608  hddm_s::CdcTruthPointList::iterator iter;
609  for (iter = points.begin(); iter != points.end(); ++iter) {
610  DMCTrackHit *mctrackhit = new DMCTrackHit;
611  mctrackhit->r = iter->getR();
612  mctrackhit->phi = iter->getPhi();
613  mctrackhit->z = iter->getZ();
614  mctrackhit->track = iter->getTrack();
615  mctrackhit->primary = iter->getPrimary();
616  mctrackhit->ptype = iter->getPtype();
617  mctrackhit->system = SYS_CDC;
618  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
619  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
620  data.push_back(mctrackhit);
621  }
622 
623  return NOERROR;
624 }
625 
626 //-------------------
627 // GetFDCTruthHits
628 //-------------------
629 jerror_t DEventSourceHDDM::GetFDCTruthHits(hddm_s::HDDM *record,
630  vector<DMCTrackHit*>& data)
631 {
632  const hddm_s::FdcTruthPointList &points = record->getFdcTruthPoints();
633  hddm_s::FdcTruthPointList::iterator iter;
634  for (iter = points.begin(); iter != points.end(); ++iter) {
635  float x = iter->getX();
636  float y = iter->getY();
637  DMCTrackHit *mctrackhit = new DMCTrackHit;
638  mctrackhit->r = sqrt(x*x + y*y);
639  mctrackhit->phi = atan2(y,x);
640  mctrackhit->z = iter->getZ();
641  mctrackhit->track = iter->getTrack();
642  mctrackhit->primary = iter->getPrimary();
643  mctrackhit->ptype = iter->getPtype();
644  mctrackhit->system = SYS_FDC;
645  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
646  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
647  data.push_back(mctrackhit);
648  }
649 
650  return NOERROR;
651 }
652 
653 //-------------------
654 // GetBCALTruthHits
655 //-------------------
656 jerror_t DEventSourceHDDM::GetBCALTruthHits(hddm_s::HDDM *record,
657  vector<DMCTrackHit*>& data)
658 {
659  const hddm_s::BcalTruthShowerList &showers = record->getBcalTruthShowers();
660  hddm_s::BcalTruthShowerList::iterator iter;
661  for (iter = showers.begin(); iter != showers.end(); ++iter) {
662  DMCTrackHit *mctrackhit = new DMCTrackHit;
663  mctrackhit->r = iter->getR();
664  mctrackhit->phi = iter->getPhi();
665  mctrackhit->z = iter->getZ();
666  mctrackhit->track = iter->getTrack();
667  mctrackhit->primary = iter->getPrimary();
668  mctrackhit->ptype = iter->getPtype();
669  mctrackhit->system = SYS_BCAL;
670  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
671  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
672  data.push_back(mctrackhit);
673  }
674 
675  return NOERROR;
676 }
677 
678 //-------------------
679 // GetTOFTruthHits
680 //-------------------
681 jerror_t DEventSourceHDDM::GetTOFTruthHits(hddm_s::HDDM *record,
682  vector<DMCTrackHit*>& data)
683 {
684  const hddm_s::FtofTruthPointList &points = record->getFtofTruthPoints();
685  hddm_s::FtofTruthPointList::iterator iter;
686  for (iter = points.begin(); iter != points.end(); ++iter) {
687  float x = iter->getX();
688  float y = iter->getY();
689  DMCTrackHit *mctrackhit = new DMCTrackHit;
690  mctrackhit->r = sqrt(x*x + y*y);
691  mctrackhit->phi = atan2(y,x);
692  mctrackhit->z = iter->getZ();
693  mctrackhit->track = iter->getTrack();
694  mctrackhit->primary = iter->getPrimary();
695  mctrackhit->ptype = iter->getPtype(); // save GEANT particle type
696  mctrackhit->system = SYS_TOF;
697  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
698  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
699  data.push_back(mctrackhit);
700  }
701 
702  return NOERROR;
703 }
704 
705 //-------------------
706 // GetCherenkovTruthHits
707 // modified by yqiang, Oct 10 2012
708 //-------------------
709 jerror_t DEventSourceHDDM::GetCherenkovTruthHits(hddm_s::HDDM *record,
710  vector<DMCTrackHit*>& data)
711 {
712  const hddm_s::CereTruthPointList &points = record->getCereTruthPoints();
713  hddm_s::CereTruthPointList::iterator iter;
714  for (iter = points.begin(); iter != points.end(); ++iter) {
715  float x = iter->getX();
716  float y = iter->getY();
717  DMCTrackHit *mctrackhit = new DMCTrackHit;
718  mctrackhit->r = sqrt(x*x + y*y);
719  mctrackhit->phi = atan2(y,x);
720  mctrackhit->z = iter->getZ();
721  mctrackhit->track = iter->getTrack();
722  mctrackhit->primary = iter->getPrimary();
723  mctrackhit->ptype = iter->getPtype(); // save GEANT particle typ()e
724  mctrackhit->system = SYS_CHERENKOV;
725  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
726  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
727  data.push_back(mctrackhit);
728  }
729 
730  return NOERROR;
731 }
732 
733 //-------------------
734 // GetFCALTruthHits
735 //-------------------
736 jerror_t DEventSourceHDDM::GetFCALTruthHits(hddm_s::HDDM *record,
737  vector<DMCTrackHit*>& data)
738 {
739  const hddm_s::FcalTruthShowerList &points = record->getFcalTruthShowers();
740  hddm_s::FcalTruthShowerList::iterator iter;
741  for (iter = points.begin(); iter != points.end(); ++iter) {
742  float x = iter->getX();
743  float y = iter->getY();
744  DMCTrackHit *mctrackhit = new DMCTrackHit;
745  mctrackhit->r = sqrt(x*x + y*y);
746  mctrackhit->phi = atan2(y,x);
747  mctrackhit->z = iter->getZ();
748  mctrackhit->track = iter->getTrack();
749  mctrackhit->primary = iter->getPrimary();
750  mctrackhit->ptype = iter->getPtype();
751  mctrackhit->system = SYS_FCAL;
752  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
753  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
754  data.push_back(mctrackhit);
755  }
756 
757  return NOERROR;
758 }
759 
760 //-------------------
761 // GetCCALTruthHits
762 //-------------------
763 jerror_t DEventSourceHDDM::GetCCALTruthHits(hddm_s::HDDM *record,
764  vector<DMCTrackHit*>& data)
765 {
766  const hddm_s::CcalTruthShowerList &points = record->getCcalTruthShowers();
767  hddm_s::CcalTruthShowerList::iterator iter;
768  for (iter = points.begin(); iter != points.end(); ++iter) {
769  float x = iter->getX();
770  float y = iter->getY();
771  DMCTrackHit *mctrackhit = new DMCTrackHit;
772  mctrackhit->r = sqrt(x*x + y*y);
773  mctrackhit->phi = atan2(y,x);
774  mctrackhit->z = iter->getZ();
775  mctrackhit->track = iter->getTrack();
776  mctrackhit->primary = iter->getPrimary();
777  mctrackhit->ptype = iter->getPtype();
778  mctrackhit->system = SYS_CCAL;
779  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
780  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
781  data.push_back(mctrackhit);
782  }
783 
784  return NOERROR;
785 }
786 
787 
788 //-------------------
789 // GetSCTruthHits
790 //-------------------
791 jerror_t DEventSourceHDDM::GetSCTruthHits(hddm_s::HDDM *record,
792  vector<DMCTrackHit*>& data)
793 {
794  const hddm_s::StcTruthPointList &points = record->getStcTruthPoints();
795  hddm_s::StcTruthPointList::iterator iter;
796  for (iter = points.begin(); iter != points.end(); ++iter) {
797  DMCTrackHit *mctrackhit = new DMCTrackHit;
798  mctrackhit->r = iter->getR();
799  mctrackhit->phi = iter->getPhi();
800  mctrackhit->z = iter->getZ();
801  mctrackhit->track = iter->getTrack();
802  mctrackhit->primary = iter->getPrimary();
803  mctrackhit->ptype = iter->getPtype(); // save GEANT particle type
804  mctrackhit->system = SYS_START;
805  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
806  mctrackhit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
807  data.push_back(mctrackhit);
808  }
809 
810  return NOERROR;
811 }
812 
813 //------------------
814 // Extract_DBCALSiPMHit
815 //------------------
816 jerror_t DEventSourceHDDM::Extract_DBCALSiPMHit(hddm_s::HDDM *record,
817  JFactory<DBCALSiPMHit> *factory, string tag)
818 {
819  /// Copies the data from the given hddm_s structure. This is called
820  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
821  /// returns OBJECT_NOT_AVAILABLE immediately.
822 
823  if (factory == NULL)
824  return OBJECT_NOT_AVAILABLE;
825  if (tag != "")
826  return OBJECT_NOT_AVAILABLE;
827 
828  vector<DBCALSiPMHit*> data;
829 
830  const hddm_s::BcalSiPMUpHitList &uphits = record->getBcalSiPMUpHits();
831  hddm_s::BcalSiPMUpHitList::iterator uiter;
832  for (uiter = uphits.begin(); uiter != uphits.end(); ++uiter) {
833  DBCALSiPMHit *response = new DBCALSiPMHit;
834  response->module = uiter->getModule();
835  response->layer = uiter->getLayer();
836  response->sector = uiter->getSector();
837  response->E = uiter->getE();
838  response->t = uiter->getT();
839  response->end = DBCALGeometry::kUpstream;
840  response->cellId = dBCALGeom->cellId(uiter->getModule(),
841  uiter->getLayer(),
842  uiter->getSector());
843  data.push_back(response);
844  }
845 
846  const hddm_s::BcalSiPMDownHitList &downhits = record->getBcalSiPMDownHits();
847  hddm_s::BcalSiPMDownHitList::iterator diter;
848  for (diter = downhits.begin(); diter != downhits.end(); ++diter) {
849  DBCALSiPMHit *response = new DBCALSiPMHit;
850  response->module = diter->getModule();
851  response->layer = diter->getLayer();
852  response->sector = diter->getSector();
853  response->E = diter->getE();
854  response->t = diter->getT();
855  response->end = DBCALGeometry::kDownstream;
856  response->cellId = dBCALGeom->cellId(diter->getModule(),
857  diter->getLayer(),
858  diter->getSector());
859  data.push_back(response);
860  }
861 
862  // Copy into factory
863  factory->CopyTo(data);
864 
865  return NOERROR;
866 }
867 
868 //------------------
869 // Extract_DBCALDigiHit
870 //------------------
871 jerror_t DEventSourceHDDM::Extract_DBCALDigiHit(hddm_s::HDDM *record,
872  JFactory<DBCALDigiHit> *factory, string tag)
873 {
874  /// Copies the data from the given hddm_s structure. This is called
875  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
876  /// returns OBJECT_NOT_AVAILABLE immediately.
877 
878  if (factory == NULL)
879  return OBJECT_NOT_AVAILABLE;
880  if (tag != "")
881  return OBJECT_NOT_AVAILABLE;
882 
883  vector<DBCALDigiHit*> data;
884 
885  const hddm_s::BcalfADCDigiHitList &digihits = record->getBcalfADCDigiHits();
886  hddm_s::BcalfADCDigiHitList::iterator iter;
887  for (iter = digihits.begin(); iter != digihits.end(); ++iter) {
888  DBCALDigiHit *response = new DBCALDigiHit;
889  response->module = iter->getModule();
890  response->layer = iter->getLayer();
891  response->sector = iter->getSector();
892  response->pulse_integral = (uint32_t)iter->getPulse_integral();
893  response->pulse_peak = 0;
894  if(iter->getBcalfADCPeaks().size() > 0) {
895  response->pulse_peak = iter->getBcalfADCPeak().getPeakAmp();
896  }
897  response->pulse_time = (uint32_t)iter->getPulse_time();
898  response->pedestal = 1;
899  response->QF = 1;
900  response->nsamples_integral = 1;
901  response->nsamples_pedestal = 1;
902  response->datasource = 3;
903  response->end = (iter->getEnd() == 0)? DBCALGeometry::kUpstream :
905  data.push_back(response);
906  }
907 
908  // Copy into factory
909  factory->CopyTo(data);
910 
911  return NOERROR;
912 }
913 
914 //------------------
915 // Extract_DBCALIncidentParticle
916 //------------------
918  JFactory<DBCALIncidentParticle> *factory,
919  string tag)
920 {
921  /// Copies the data from the given hddm_s structure. This is called
922  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
923  /// returns OBJECT_NOT_AVAILABLE immediately.
924 
925  if (factory == NULL)
926  return OBJECT_NOT_AVAILABLE;
927  if (tag != "")
928  return OBJECT_NOT_AVAILABLE;
929 
930  vector<DBCALIncidentParticle*> data;
931 
932  const hddm_s::BcalTruthIncidentParticleList &plist =
933  record->getBcalTruthIncidentParticles();
934  hddm_s::BcalTruthIncidentParticleList::iterator iter;
935  for (iter = plist.begin(); iter != plist.end(); ++iter) {
937  part->ptype = iter->getPtype();
938  part->px = iter->getPx();
939  part->py = iter->getPy();
940  part->pz = iter->getPz();
941  part->x = iter->getX();
942  part->y = iter->getY();
943  part->z = iter->getZ();
944  data.push_back(part);
945  }
946 
947  factory->CopyTo(data);
948 
949  return NOERROR;
950 }
951 
952 //------------------
953 // Extract_DBCALSiPMSpectrum
954 //------------------
955 jerror_t DEventSourceHDDM::Extract_DBCALSiPMSpectrum(hddm_s::HDDM *record,
956  JFactory<DBCALSiPMSpectrum> *factory,
957  string tag)
958 {
959  /// Copies the data from the given hddm_s structure. This is called
960  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
961  /// returns OBJECT_NOT_AVAILABLE immediately.
962 
963  if (factory == NULL)
964  return OBJECT_NOT_AVAILABLE;
965  if (tag != "" && tag != "TRUTH")
966  return OBJECT_NOT_AVAILABLE;
967 
968  vector<DBCALSiPMSpectrum*> data;
969 
970  const hddm_s::BcalSiPMSpectrumList &specs = record->getBcalSiPMSpectrums();
971  hddm_s::BcalSiPMSpectrumList::iterator iter;
972  for (iter = specs.begin(); iter != specs.end(); ++iter) {
973  DBCALSiPMSpectrum *dana_spectrum = new DBCALSiPMSpectrum;
974  dana_spectrum->module = iter->getModule();
975  dana_spectrum->layer = iter->getLayer();
976  dana_spectrum->sector = iter->getSector();
977  dana_spectrum->end = (iter->getEnd()==0)?
980  if (tag == "")
981  dana_spectrum->incident_id = 0;
982  else if (tag == "TRUTH") {
983  const hddm_s::BcalSiPMTruthList &truths = iter->getBcalSiPMTruths();
984  if (truths.size() > 0)
985  dana_spectrum->incident_id = truths.begin()->getIncident_id();
986  else
987  dana_spectrum->incident_id = 0;
988  }
989 
990  double t = iter->getTstart();
991  double bin_width = iter->getBin_width();
992  stringstream ss(iter->getVals());
993 
994  // Extract values and use them to fill histo
995  string entry;
996  while (ss >> entry) {
997  if (entry[0] == 'X') {
998  // get rid of the X, the rest of the entry is the number of zeroes to add
999  stringstream sss(entry.substr(1));
1000  int num_zeros;
1001  sss >> num_zeros;
1002 
1003  for(int i=0; i<num_zeros; i++) {
1004  dana_spectrum->spectrum.Fill(t, 0.0);
1005  t += bin_width;
1006  }
1007  } else {
1008  stringstream sss(entry);
1009  double dE;
1010  sss >> dE;
1011  dana_spectrum->spectrum.Fill(t, dE);
1012  t += bin_width;
1013  }
1014  }
1015 
1016  data.push_back(dana_spectrum);
1017  }
1018 
1019  // Copy into factory
1020  factory->CopyTo(data);
1021 
1022  return NOERROR;
1023 }
1024 
1025 //------------------
1026 // Extract_DBCALTDCDigiHit
1027 //------------------
1028 jerror_t DEventSourceHDDM::Extract_DBCALTDCDigiHit(hddm_s::HDDM *record,
1029  JFactory<DBCALTDCDigiHit> *factory, string tag)
1030 {
1031  /// Copies the data from the given hddm_s structure. This is called
1032  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1033  /// returns OBJECT_NOT_AVAILABLE immediately.
1034 
1035  if (factory == NULL)
1036  return OBJECT_NOT_AVAILABLE;
1037  if (tag != "")
1038  return OBJECT_NOT_AVAILABLE;
1039 
1040  vector<DBCALTDCDigiHit*> data;
1041 
1042  const hddm_s::BcalTDCDigiHitList &hits = record->getBcalTDCDigiHits();
1043  hddm_s::BcalTDCDigiHitList::iterator iter;
1044  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1045  DBCALTDCDigiHit *bcaltdchit = new DBCALTDCDigiHit;
1046  bcaltdchit->module = iter->getModule();
1047  bcaltdchit->layer = iter->getLayer();
1048  bcaltdchit->sector = iter->getSector();
1049  bcaltdchit->end = (iter->getEnd() == 0)? DBCALGeometry::kUpstream :
1051  bcaltdchit->time = (uint32_t)iter->getTime();
1052  data.push_back(bcaltdchit);
1053  }
1054 
1055  // Copy into factory
1056  factory->CopyTo(data);
1057 
1058  return NOERROR;
1059 }
1060 
1061 //------------------
1062 // Extract_DMCReaction
1063 //------------------
1064 jerror_t DEventSourceHDDM::Extract_DMCReaction(hddm_s::HDDM *record,
1065  JFactory<DMCReaction> *factory, string tag,
1066  JEventLoop *loop)
1067 {
1068  /// Copies the data from the given hddm_s structure. This is called
1069  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1070  /// returns OBJECT_NOT_AVAILABLE immediately.
1071 
1072  if (tag != "")
1073  return OBJECT_NOT_AVAILABLE;
1074 
1075  if (factory == NULL)
1076  return OBJECT_NOT_AVAILABLE;
1077 
1078  double locTargetCenterZ = 0.0;
1079  int locRunNumber = loop->GetJEvent().GetRunNumber();
1080  LockRead();
1081  {
1082  locTargetCenterZ = dTargetCenterZMap[locRunNumber];
1083  }
1084  UnlockRead();
1085  DVector3 locPosition(0.0, 0.0, locTargetCenterZ);
1086 
1087  vector<DMCReaction*> dmcreactions;
1088 
1089  const hddm_s::ReactionList &reacts = record->getReactions();
1090  hddm_s::ReactionList::iterator iter;
1091  for (iter = reacts.begin(); iter != reacts.end(); ++iter) {
1092  DMCReaction *mcreaction = new DMCReaction;
1093  dmcreactions.push_back(mcreaction);
1094  mcreaction->type = iter->getType();
1095  mcreaction->weight = iter->getWeight();
1096  hddm_s::Origin &origin = iter->getVertex().getOrigin();
1097  double torig = origin.getT();
1098  double zorig = origin.getVz();
1099 
1100  const hddm_s::BeamList &beams = record->getBeams();
1101  if (beams.size() > 0) {
1102  hddm_s::Beam &beam = iter->getBeam();
1103  DVector3 mom(beam.getMomentum().getPx(),
1104  beam.getMomentum().getPy(),
1105  beam.getMomentum().getPz());
1106  mcreaction->beam.setPosition(locPosition);
1107  mcreaction->beam.setMomentum(mom);
1108  mcreaction->beam.setPID(Gamma);
1109  mcreaction->target.setPID(IDTrack(mcreaction->beam.charge(),
1110  mcreaction->beam.mass()));
1111  mcreaction->beam.setTime(torig - (zorig - locTargetCenterZ)/29.9792458);
1112  }
1113  else {
1114  // fake values for DMCReaction
1115  mcreaction->beam.setPosition(locPosition);
1116  mcreaction->beam.setPID(Gamma);
1117  }
1118 
1119  const hddm_s::TargetList &targets = record->getTargets();
1120  if (targets.size() > 0) {
1121  hddm_s::Target &target = iter->getTarget();
1122  DKinematicData target_kd;
1123  DVector3 mom(target.getMomentum().getPx(),
1124  target.getMomentum().getPy(),
1125  target.getMomentum().getPz());
1126  mcreaction->target.setPosition(locPosition);
1127  mcreaction->target.setMomentum(mom);
1128  mcreaction->target.setPID(IDTrack(target.getProperties().getCharge(),
1129  target.getProperties().getMass()));
1130  mcreaction->target.setTime(torig - (zorig - locTargetCenterZ)/29.9792458);
1131  }
1132  else {
1133  // fake values for DMCReaction
1134  mcreaction->target.setPosition(locPosition);
1135  }
1136  }
1137 
1138  // Copy into factories
1139  //_DBG_<<"Creating "<<dmcreactions.size()<<" DMCReaction objects"<<endl;
1140 
1141  factory->CopyTo(dmcreactions);
1142 
1143  return NOERROR;
1144 }
1145 
1146 //------------------
1147 // Extract_DMCThrown
1148 //------------------
1149 jerror_t DEventSourceHDDM::Extract_DMCThrown(hddm_s::HDDM *record,
1150  JFactory<DMCThrown> *factory, string tag)
1151 {
1152  /// Copies the data from the given hddm_s structure. This is called
1153  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1154  /// returns OBJECT_NOT_AVAILABLE immediately.
1155 
1156  if (factory == NULL)
1157  return OBJECT_NOT_AVAILABLE;
1158  if (tag != "")
1159  return OBJECT_NOT_AVAILABLE;
1160 
1161  vector<DMCThrown*> data;
1162 
1163  const hddm_s::VertexList &verts = record->getVertices();
1164  hddm_s::VertexList::iterator iter;
1165  for (iter = verts.begin(); iter != verts.end(); ++iter) {
1166  const hddm_s::OriginList &origs = iter->getOrigins();
1167  const hddm_s::ProductList &prods = iter->getProducts();
1168  double vertex[4] = {0., 0., 0., 0.};
1169  if (origs.size() > 0) {
1170  vertex[0] = iter->getOrigin().getT();
1171  vertex[1] = iter->getOrigin().getVx();
1172  vertex[2] = iter->getOrigin().getVy();
1173  vertex[3] = iter->getOrigin().getVz();
1174  }
1175  hddm_s::ProductList::iterator piter;
1176  for (piter = prods.begin(); piter != prods.end(); ++piter) {
1177  double E = piter->getMomentum().getE();
1178  double px = piter->getMomentum().getPx();
1179  double py = piter->getMomentum().getPy();
1180  double pz = piter->getMomentum().getPz();
1181  double mass = sqrt(E*E - (px*px + py*py + pz*pz));
1182  if (!isfinite(mass))
1183  mass = 0.0;
1184  DMCThrown *mcthrown = new DMCThrown;
1185  mcthrown->type = piter->getType();
1186  mcthrown->myid = piter->getId();
1187  mcthrown->parentid = piter->getParentid();
1188  mcthrown->mech = piter->getMech();
1189  mcthrown->pdgtype = piter->getPdgtype();
1190  mcthrown->setPID((Particle_t)mcthrown->type);
1191  mcthrown->setMomentum(DVector3(px, py, pz));
1192  mcthrown->setPosition(DVector3(vertex[1], vertex[2], vertex[3]));
1193  mcthrown->setTime(vertex[0]);
1194  data.push_back(mcthrown);
1195  }
1196  }
1197 
1198  // Copy into factory
1199  factory->CopyTo(data);
1200 
1201  return NOERROR;
1202 }
1203 
1204 //------------------
1205 // Extract_DCDCHit
1206 //------------------
1207 jerror_t DEventSourceHDDM::Extract_DCDCHit(JEventLoop* locEventLoop, hddm_s::HDDM *record,
1208  JFactory<DCDCHit> *factory, string tag)
1209 {
1210  /// Copies the data from the given hddm_s structure. This is called
1211  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1212  /// returns OBJECT_NOT_AVAILABLE immediately.
1213 
1214  if (factory == NULL)
1215  return OBJECT_NOT_AVAILABLE;
1216 
1217  // Since we are writing out CDC hits with the new "Calib" tag by default
1218  // assume that is what we are reading in, so that we don't support the
1219  // default tag anymore
1220  // sdobbs -- 3/13/2018
1221  //if (tag != "" && tag != "TRUTH" && tag != "Calib")
1222  if (tag != "TRUTH" && tag != "Calib")
1223  return OBJECT_NOT_AVAILABLE;
1224 
1225  vector<DCDCHit*> data;
1226 
1227  if ( tag == "" || tag == "Calib" ) {
1228  vector<const DCDCHit*> locTruthHits;
1229  locEventLoop->Get(locTruthHits, "TRUTH");
1230 
1231  //pre-sort truth hits
1232  map<pair<int, int>, vector<const DCDCHit*>> locTruthHitMap; //key pair: ring, straw
1233  for(auto& locTruthHit : locTruthHits)
1234  locTruthHitMap[std::make_pair(locTruthHit->ring, locTruthHit->straw)].push_back(locTruthHit);
1235 
1236  const hddm_s::CdcStrawHitList &hits = record->getCdcStrawHits();
1237  hddm_s::CdcStrawHitList::iterator iter;
1238  int locIndex = 0;
1239  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1240  DCDCHit *hit = new DCDCHit;
1241  hit->ring = iter->getRing();
1242  hit->straw = iter->getStraw();
1243  hit->q = iter->getQ();
1244  hit->t = iter->getT();
1245  if(iter->getCdcDigihits().size() > 0) {
1246  hit->amp = iter->getCdcDigihit().getPeakAmp();
1247  }
1248  else{
1249  // for generated events (not folded-in background events) for which we
1250  // have no digi hits return q
1251  hit->amp=hit->q;
1252  }
1253  hit->QF = 0;
1254  if(iter->getCdcHitQFs().size() > 0) {
1255  hit->QF = iter->getCdcHitQF().getQF();
1256  }
1257  hit->d = 0.; // initialize to zero to avoid any NaN
1258  hit->itrack = 0; // track information is in TRUTH tag
1259  hit->ptype = 0; // ditto
1260 
1261  //match hit between truth & recon
1262  auto& locPotentialTruthHits = locTruthHitMap[std::make_pair(hit->ring, hit->straw)];
1263  double locBestDeltaT = 9.9E99;
1264  const DCDCHit* locBestTruthHit = nullptr;
1265  for(auto& locTruthHit : locPotentialTruthHits)
1266  {
1267  auto locDeltaT = fabs(hit->t - locTruthHit->t);
1268  if(locDeltaT >= locBestDeltaT)
1269  continue;
1270  locBestDeltaT = locDeltaT;
1271  locBestTruthHit = locTruthHit;
1272  }
1273  if(locBestTruthHit != nullptr)
1274  hit->AddAssociatedObject(locBestTruthHit);
1275 
1276  data.push_back(hit);
1277  ++locIndex;
1278  }
1279  }
1280  else if (tag == "TRUTH") {
1281  const hddm_s::CdcStrawTruthHitList &thits = record->getCdcStrawTruthHits();
1282  hddm_s::CdcStrawTruthHitList::iterator iter;
1283  for (iter = thits.begin(); iter != thits.end(); ++iter) {
1284  DCDCHit *hit = new DCDCHit;
1285  hit->ring = iter->getRing();
1286  hit->straw = iter->getStraw();
1287  hit->q = iter->getQ();
1288  hit->t = iter->getT();
1289  hit->d = iter->getD();
1290  hit->itrack = iter->getItrack();
1291  hit->ptype = iter->getPtype();
1292  data.push_back(hit);
1293  }
1294  }
1295 
1296  // Copy into factory
1297  factory->CopyTo(data);
1298 
1299  return NOERROR;
1300 }
1301 
1302 
1303 //------------------
1304 // Extract_DFDCHit
1305 //------------------
1306 jerror_t DEventSourceHDDM::Extract_DFDCHit(hddm_s::HDDM *record,
1307  JFactory<DFDCHit> *factory, string tag)
1308 {
1309  /// Copies the data from the given hddm_s structure. This is called
1310  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1311  /// returns OBJECT_NOT_AVAILABLE immediately.
1312 
1313  if (factory == NULL)
1314  return OBJECT_NOT_AVAILABLE;
1315  if (tag != "" && tag != "TRUTH" && tag != "CALIB")
1316  return OBJECT_NOT_AVAILABLE;
1317 
1318  vector<DFDCHit*> data;
1319 
1320  if (tag == "") {
1321  const hddm_s::FdcAnodeHitList &ahits = record->getFdcAnodeHits();
1322  hddm_s::FdcAnodeHitList::iterator ahiter;
1323  for (ahiter = ahits.begin(); ahiter != ahits.end(); ++ahiter) {
1324  DFDCHit* newHit = new DFDCHit();
1325  newHit->layer = ahiter->getLayer();
1326  newHit->module = ahiter->getModule();
1327  newHit->element = ahiter->getWire();
1328  newHit->q = ahiter->getDE();
1329  newHit->pulse_height = 0.; // not measured
1330  newHit->t = ahiter->getT();
1331  newHit->d = 0.; // initialize to zero to avoid any NaN
1332  newHit->itrack = 0; // track information is in TRUTH tag
1333  newHit->ptype = 0; // ditto
1334  newHit->plane = 2;
1335  newHit->type = 0;
1336  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1337  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1338  newHit->r = DFDCGeometry::getWireR(newHit);
1339  data.push_back(newHit);
1340  }
1341 
1342  // Ditto for the cathodes.
1343  const hddm_s::FdcCathodeHitList &chits = record->getFdcCathodeHits();
1344  hddm_s::FdcCathodeHitList::iterator chiter;
1345  for (chiter = chits.begin(); chiter != chits.end(); ++chiter) {
1346  DFDCHit* newHit = new DFDCHit();
1347  newHit->layer = chiter->getLayer();
1348  newHit->module = chiter->getModule();
1349  newHit->element = chiter->getStrip();
1350  if (newHit->element > 1000)
1351  newHit->element -= 1000;
1352  newHit->plane = chiter->getPlane();
1353  newHit->q = chiter->getQ();
1354  newHit->pulse_height = newHit->q;
1355  if(chiter->getFdcDigihits().size() > 0) {
1356  newHit->pulse_height = chiter->getFdcDigihit().getPeakAmp();
1357  }
1358  newHit->t = chiter->getT();
1359  newHit->d = 0.; // initialize to zero to avoid any NaN
1360  newHit->itrack = 0; // track information is in TRUTH tag
1361  newHit->ptype = 0; // ditto
1362  newHit->type = 1;
1363  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1364  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1365  newHit->r = DFDCGeometry::getStripR(newHit);
1366  data.push_back(newHit);
1367  }
1368  }
1369 
1370  else if (tag == "TRUTH"){
1371  const hddm_s::FdcAnodeTruthHitList &aths = record->getFdcAnodeTruthHits();
1372  hddm_s::FdcAnodeTruthHitList::iterator atiter;
1373  for (atiter = aths.begin(); atiter != aths.end(); ++atiter) {
1374  DFDCHit* newHit = new DFDCHit();
1375  newHit->layer = atiter->getLayer();
1376  newHit->module = atiter->getModule();
1377  newHit->element = atiter->getWire();
1378  newHit->q = atiter->getDE();
1379  newHit->pulse_height=0.; // not measured
1380  newHit->t = atiter->getT();
1381  newHit->d = atiter->getD();
1382  newHit->itrack = atiter->getItrack();
1383  newHit->ptype = atiter->getPtype();
1384  newHit->plane = 2;
1385  newHit->type = 0;
1386  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1387  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1388  newHit->r = DFDCGeometry::getWireR(newHit);
1389  data.push_back(newHit);
1390  }
1391 
1392  // Ditto for the cathodes.
1393  const hddm_s::FdcCathodeTruthHitList &cths =
1394  record->getFdcCathodeTruthHits();
1395  hddm_s::FdcCathodeTruthHitList::iterator ctiter;
1396  for (ctiter = cths.begin(); ctiter != cths.end(); ++ctiter) {
1397  DFDCHit* newHit = new DFDCHit();
1398  newHit->layer = ctiter->getLayer();
1399  newHit->module = ctiter->getModule();
1400  newHit->element = ctiter->getStrip();
1401  if (newHit->element > 1000)
1402  newHit->element -= 1000;
1403  newHit->plane = ctiter->getPlane();
1404  newHit->q = ctiter->getQ();
1405  newHit->pulse_height = newHit->q;
1406  newHit->t = ctiter->getT();
1407  newHit->d = 0.; // initialize to zero to avoid any NaN
1408  newHit->itrack = ctiter->getItrack();
1409  newHit->ptype = ctiter->getPtype();
1410  newHit->type = 1;
1411  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1412  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1413  newHit->r = DFDCGeometry::getStripR(newHit);
1414  data.push_back(newHit);
1415  }
1416  }
1417 
1418  else if (tag == "CALIB") {
1419  // Deal with the wires
1420  const hddm_s::FdcAnodeHitList &ahits = record->getFdcAnodeHits();
1421  hddm_s::FdcAnodeHitList::iterator ahiter;
1422  for (ahiter = ahits.begin(); ahiter != ahits.end(); ++ahiter) {
1423  DFDCHit* newHit = new DFDCHit();
1424  newHit->layer = ahiter->getLayer();
1425  newHit->module = ahiter->getModule();
1426  newHit->element = ahiter->getWire();
1427  newHit->q = ahiter->getDE();
1428  newHit->t = ahiter->getT();
1429  newHit->d = 0.; // initialize to zero to avoid any NaN
1430  newHit->itrack = 0; // track information is in TRUTH tag
1431  newHit->ptype = 0; // ditto
1432  newHit->plane = 2;
1433  newHit->type = 0;
1434  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1435  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1436  newHit->r = DFDCGeometry::getWireR(newHit);
1437  data.push_back(newHit);
1438  }
1439 
1440  // Ditto for the cathodes.
1441  const hddm_s::FdcCathodeHitList &chits = record->getFdcCathodeHits();
1442  hddm_s::FdcCathodeHitList::iterator chiter;
1443  for (chiter = chits.begin(); chiter != chits.end(); ++chiter) {
1444  DFDCHit* newHit = new DFDCHit();
1445  newHit->layer = chiter->getLayer();
1446  newHit->module = chiter->getModule();
1447  newHit->element = chiter->getStrip();
1448  if (newHit->element > 1000)
1449  newHit->element-=1000;
1450  newHit->plane = chiter->getPlane();
1451  if (newHit->plane == 1) // v
1452  newHit->q = chiter->getQ()*vscale[newHit->element-1];
1453  else // u
1454  newHit->q = chiter->getQ()*uscale[newHit->element-1];
1455  newHit->t = chiter->getT();
1456  newHit->d = 0.; // initialize to zero to avoid any NaN
1457  newHit->itrack = 0; // track information is in TRUTH tag
1458  newHit->ptype = 0; // ditto
1459  newHit->type = 1;
1460  newHit->gPlane = DFDCGeometry::gPlane(newHit);
1461  newHit->gLayer = DFDCGeometry::gLayer(newHit);
1462  newHit->r = DFDCGeometry::getStripR(newHit);
1463  data.push_back(newHit);
1464  }
1465  }
1466 
1467  // Copy into factory
1468  factory->CopyTo(data);
1469 
1470  return NOERROR;
1471 }
1472 
1473 //------------------
1474 // Extract_DBCALTruthShower
1475 //------------------
1476 jerror_t DEventSourceHDDM::Extract_DBCALTruthShower(hddm_s::HDDM *record,
1477  JFactory<DBCALTruthShower> *factory,
1478  string tag)
1479 {
1480  /// Copies the data from the given hddm_s structure. This is called
1481  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1482  /// returns OBJECT_NOT_AVAILABLE immediately.
1483 
1484  if (factory == NULL)
1485  return OBJECT_NOT_AVAILABLE;
1486  if (tag != "")
1487  return OBJECT_NOT_AVAILABLE;
1488 
1489  vector<DBCALTruthShower*> data;
1490 
1491  const hddm_s::BcalTruthShowerList &shows = record->getBcalTruthShowers();
1492  hddm_s::BcalTruthShowerList::iterator iter;
1493  for (iter = shows.begin(); iter != shows.end(); ++iter) {
1494  DBCALTruthShower *bcaltruth = new DBCALTruthShower;
1495  bcaltruth->track = iter->getTrack();
1496  bcaltruth->ptype = iter->getPtype();
1497  bcaltruth->primary = (iter->getPrimary())? 1 : 0;
1498  bcaltruth->phi = iter->getPhi();
1499  bcaltruth->r = iter->getR();
1500  bcaltruth->z = iter->getZ();
1501  bcaltruth->t = iter->getT();
1502  bcaltruth->E = iter->getE();
1503  bcaltruth->px = iter->getPx();
1504  bcaltruth->py = iter->getPy();
1505  bcaltruth->pz = iter->getPz();
1506  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
1507  bcaltruth->itrack = (ids.size())? ids.begin()->getItrack() : 0;
1508  data.push_back(bcaltruth);
1509  }
1510 
1511  // Copy into factory
1512  factory->CopyTo(data);
1513 
1514  return NOERROR;
1515 }
1516 
1517 //------------------
1518 // Extract_DBCALTruthCell
1519 //------------------
1520 jerror_t DEventSourceHDDM::Extract_DBCALTruthCell(hddm_s::HDDM *record,
1521  JFactory<DBCALTruthCell> *factory,
1522  string tag)
1523 {
1524  /// Copies the data from the given hddm_s structure. This is called
1525  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1526  /// returns OBJECT_NOT_AVAILABLE immediately.
1527 
1528  if (factory == NULL)
1529  return OBJECT_NOT_AVAILABLE;
1530  if (tag != "")
1531  return OBJECT_NOT_AVAILABLE;
1532 
1533  vector<DBCALTruthCell*> data;
1534 
1535  const hddm_s::BcalTruthHitList &hits = record->getBcalTruthHits();
1536  hddm_s::BcalTruthHitList::iterator hiter;
1537  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1538  DBCALTruthCell *truthcell = new DBCALTruthCell();
1539  truthcell->module = hiter->getModule();
1540  truthcell->layer = hiter->getLayer();
1541  truthcell->sector = hiter->getSector();
1542  truthcell->E = hiter->getE();
1543  truthcell->t = hiter->getT();
1544  truthcell->zLocal = hiter->getZLocal();
1545  data.push_back(truthcell);
1546  }
1547 
1548  // Copy into factory
1549  factory->CopyTo(data);
1550 
1551  return NOERROR;
1552 }
1553 
1554 //------------------
1555 // Extract_DFCALTruthShower
1556 //------------------
1557 jerror_t DEventSourceHDDM::Extract_DFCALTruthShower(hddm_s::HDDM *record,
1558  JFactory<DFCALTruthShower> *factory,
1559  string tag)
1560 {
1561  /// Copies the data from the given hddm_s structure. This is called
1562  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1563  /// returns OBJECT_NOT_AVAILABLE immediately.
1564 
1565  if (factory == NULL)
1566  return OBJECT_NOT_AVAILABLE;
1567  if (tag != "")
1568  return OBJECT_NOT_AVAILABLE;
1569 
1570  vector<DFCALTruthShower*> data;
1571  JObject::oid_t id=1;
1572 
1573  const hddm_s::FcalTruthShowerList &shows = record->getFcalTruthShowers();
1574  hddm_s::FcalTruthShowerList::iterator iter;
1575  for (iter = shows.begin(); iter != shows.end(); ++iter) {
1576  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
1577  int itrack = (ids.size())? ids.begin()->getItrack() : 0;
1578  DFCALTruthShower *dfcaltruthshower = new DFCALTruthShower(
1579  id++,
1580  iter->getX(),
1581  iter->getY(),
1582  iter->getZ(),
1583  iter->getPx(),
1584  iter->getPy(),
1585  iter->getPz(),
1586  iter->getE(),
1587  iter->getT(),
1588  iter->getPrimary(),
1589  iter->getTrack(),
1590  iter->getPtype(),
1591  itrack
1592  );
1593  data.push_back(dfcaltruthshower);
1594  }
1595 
1596  // Copy into factory
1597  factory->CopyTo(data);
1598 
1599  return NOERROR;
1600 }
1601 
1602 //------------------
1603 // Extract_DFCALHit
1604 //------------------
1605 jerror_t DEventSourceHDDM::Extract_DFCALHit(hddm_s::HDDM *record,
1606  JFactory<DFCALHit> *factory, string tag,
1607  JEventLoop* eventLoop)
1608 {
1609  /// Copies the data from the given hddm_s structure. This is called
1610  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1611  /// returs OBJECT_NOT_AVAILABLE immediately.
1612 
1613  if (factory == NULL)
1614  return OBJECT_NOT_AVAILABLE;
1615  if (tag != "" && tag != "TRUTH")
1616  return OBJECT_NOT_AVAILABLE;
1617 
1618  // extract the FCAL Geometry (for isBlockActive() and positionOnFace())
1619  vector<const DFCALGeometry*> fcalGeomVect;
1620  eventLoop->Get( fcalGeomVect );
1621  if (fcalGeomVect.size() < 1)
1622  return OBJECT_NOT_AVAILABLE;
1623  const DFCALGeometry& fcalGeom = *(fcalGeomVect[0]);
1624 
1625  vector<DFCALHit*> data;
1626 
1627  if (tag == "") {
1628  const hddm_s::FcalHitList &hits = record->getFcalHits();
1629  hddm_s::FcalHitList::iterator iter;
1630  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1631  int row = iter->getRow();
1632  int column = iter->getColumn();
1633 
1634  // Filter out non-physical blocks here
1635  if (!fcalGeom.isBlockActive(row, column))
1636  continue;
1637 
1638  // Get position of blocks on front face. (This should really come from
1639  // hdgeant directly so the positions can be shifted in mcsmear.)
1640  DVector2 pos = fcalGeom.positionOnFace(row, column);
1641 
1642  DFCALHit *mchit = new DFCALHit();
1643  mchit->row = row;
1644  mchit->column = column;
1645  mchit->x = pos.X();
1646  mchit->y = pos.Y();
1647  mchit->E = iter->getE();
1648  mchit->t = iter->getT();
1649  mchit->intOverPeak = 6.;
1650  if(iter->getFcalDigihits().size() > 0) {
1651  mchit->intOverPeak = iter->getFcalDigihit().getIntegralOverPeak();
1652  }
1653  data.push_back(mchit);
1654  }
1655  }
1656  else if (tag == "TRUTH") {
1657  const hddm_s::FcalTruthHitList &hits = record->getFcalTruthHits();
1658  hddm_s::FcalTruthHitList::iterator iter;
1659  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1660  int row = iter->getRow();
1661  int column = iter->getColumn();
1662 
1663  // Filter out non-physical blocks here
1664  if (!fcalGeom.isBlockActive(row, column))
1665  continue;
1666 
1667  // Get position of blocks on front face. (This should really come from
1668  // hdgeant directly so the positions can be shifted in mcsmear.)
1669  DVector2 pos = fcalGeom.positionOnFace(row, column);
1670 
1671  DFCALHit *mchit = new DFCALHit();
1672  mchit->row = row;
1673  mchit->column = column;
1674  mchit->x = pos.X();
1675  mchit->y = pos.Y();
1676  mchit->E = iter->getE();
1677  mchit->t = iter->getT();
1678  mchit->intOverPeak = 6.;
1679  data.push_back(mchit);
1680  }
1681  }
1682 
1683  // Copy into factory
1684  factory->CopyTo(data);
1685 
1686  return NOERROR;
1687 }
1688 
1689 //------------------
1690 // Extract_DCCALTruthShower
1691 //------------------
1692 jerror_t DEventSourceHDDM::Extract_DCCALTruthShower(hddm_s::HDDM *record,
1693  JFactory<DCCALTruthShower> *factory,
1694  string tag)
1695 {
1696  /// Copies the data from the given hddm_s structure. This is called
1697  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1698  /// returns OBJECT_NOT_AVAILABLE immediately.
1699 
1700  if (factory == NULL)
1701  return OBJECT_NOT_AVAILABLE;
1702  if (tag != "")
1703  return OBJECT_NOT_AVAILABLE;
1704 
1705  vector<DCCALTruthShower*> data;
1706  JObject::oid_t id=1;
1707 
1708  const hddm_s::CcalTruthShowerList &shows = record->getCcalTruthShowers();
1709  hddm_s::CcalTruthShowerList::iterator iter;
1710  for (iter = shows.begin(); iter != shows.end(); ++iter) {
1711  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
1712  int itrack = (ids.size())? ids.begin()->getItrack() : 0;
1713  DCCALTruthShower *dccaltruthshower = new DCCALTruthShower(
1714  id++,
1715  iter->getX(),
1716  iter->getY(),
1717  iter->getZ(),
1718  iter->getPx(),
1719  iter->getPy(),
1720  iter->getPz(),
1721  iter->getE(),
1722  iter->getT(),
1723  iter->getPrimary(),
1724  iter->getTrack(),
1725  iter->getPtype(),
1726  itrack
1727  );
1728  data.push_back(dccaltruthshower);
1729  }
1730 
1731  // Copy into factory
1732  factory->CopyTo(data);
1733 
1734  return NOERROR;
1735 }
1736 
1737 //------------------
1738 // Extract_DCCALHit
1739 //------------------
1740 jerror_t DEventSourceHDDM::Extract_DCCALHit(hddm_s::HDDM *record,
1741  JFactory<DCCALHit> *factory, string tag,
1742  JEventLoop* eventLoop)
1743 {
1744  /// Copies the data from the given hddm_s structure. This is called
1745  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1746  /// returs OBJECT_NOT_AVAILABLE immediately.
1747 
1748  if (factory == NULL)
1749  return OBJECT_NOT_AVAILABLE;
1750  if (tag != "" && tag != "TRUTH")
1751  return OBJECT_NOT_AVAILABLE;
1752 
1753  // extract the CCAL Geometry (for isBlockActive() and positionOnFace())
1754  vector<const DCCALGeometry*> ccalGeomVect;
1755  eventLoop->Get( ccalGeomVect );
1756  if (ccalGeomVect.size() < 1)
1757  return OBJECT_NOT_AVAILABLE;
1758  const DCCALGeometry& ccalGeom = *(ccalGeomVect[0]);
1759 
1760  vector<DCCALHit*> data;
1761  int hitId = 0;
1762 
1763  if (tag == "") {
1764  const hddm_s::CcalHitList &hits = record->getCcalHits();
1765  hddm_s::CcalHitList::iterator iter;
1766  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1767  int row = iter->getRow();
1768  int column = iter->getColumn();
1769 
1770  // Filter out non-physical blocks here
1771  if (!ccalGeom.isBlockActive(row, column))
1772  continue;
1773 
1774  // Get position of blocks on front face. (This should really come from
1775  // hdgeant directly so the poisitions can be shifted in mcsmear.)
1776  DVector2 pos = ccalGeom.positionOnFace(row, column);
1777 
1778  DCCALHit *mchit = new DCCALHit();
1779  mchit->row = row;
1780  mchit->column = column;
1781  mchit->x = pos.X();
1782  mchit->y = pos.Y();
1783  mchit->E = iter->getE();
1784  mchit->t = iter->getT();
1785  mchit->id = hitId++;
1786  data.push_back(mchit);
1787  }
1788  }
1789 
1790  else if (tag == "TRUTH") {
1791  const hddm_s::CcalTruthHitList &hits = record->getCcalTruthHits();
1792  hddm_s::CcalTruthHitList::iterator iter;
1793  for (iter = hits.begin(); iter != hits.end(); ++iter) {
1794  int row = iter->getRow();
1795  int column = iter->getColumn();
1796 
1797  // Filter out non-physical blocks here
1798  if (!ccalGeom.isBlockActive(row, column))
1799  continue;
1800 
1801  // Get position of blocks on front face. (This should really come from
1802  // hdgeant directly so the poisitions can be shifted in mcsmear.)
1803  DVector2 pos = ccalGeom.positionOnFace(row, column);
1804 
1805  DCCALHit *mchit = new DCCALHit();
1806  mchit->row = row;
1807  mchit->column = column;
1808  mchit->x = pos.X();
1809  mchit->y = pos.Y();
1810  mchit->E = iter->getE();
1811  mchit->t = iter->getT();
1812  mchit->id = hitId++;
1813  data.push_back(mchit);
1814  }
1815  }
1816 
1817  // Copy into factory
1818  factory->CopyTo(data);
1819 
1820  return NOERROR;
1821 }
1822 
1823 //------------------
1824 // Extract_DMCTrajectoryPoint
1825 //------------------
1826 jerror_t DEventSourceHDDM::Extract_DMCTrajectoryPoint(hddm_s::HDDM *record,
1827  JFactory<DMCTrajectoryPoint> *factory,
1828  string tag)
1829 {
1830  /// Copies the data from the given hddm_s structure. This is called
1831  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1832  /// returns OBJECT_NOT_AVAILABLE immediately.
1833 
1834  if (factory == NULL)
1835  return OBJECT_NOT_AVAILABLE;
1836  if (tag != "")
1837  return OBJECT_NOT_AVAILABLE;
1838 
1839  vector<DMCTrajectoryPoint*> data;
1840 
1841  const hddm_s::McTrajectoryPointList &pts = record->getMcTrajectoryPoints();
1842  hddm_s::McTrajectoryPointList::iterator iter;
1843  for (iter = pts.begin(); iter != pts.end(); ++iter) {
1845  p->x = iter->getX();
1846  p->y = iter->getY();
1847  p->z = iter->getZ();
1848  p->t = iter->getT();
1849  p->px = iter->getPx();
1850  p->py = iter->getPy();
1851  p->pz = iter->getPz();
1852  p->E = iter->getE();
1853  p->dE = iter->getDE();
1854  p->primary_track = iter->getPrimary_track();
1855  p->track = iter->getTrack();
1856  p->part = iter->getPart();
1857  p->radlen = iter->getRadlen();
1858  p->step = iter->getStep();
1859  p->mech = iter->getMech();
1860  data.push_back(p);
1861  }
1862 
1863  // Copy into factory
1864  factory->CopyTo(data);
1865 
1866  return NOERROR;
1867 }
1868 
1869 //------------------
1870 // Extract_DTOFTruth
1871 //------------------
1872 jerror_t DEventSourceHDDM::Extract_DTOFTruth(hddm_s::HDDM *record,
1873  JFactory<DTOFTruth>* factory, string tag)
1874 {
1875  /// Copies the data from the given hddm_s structure. This is called
1876  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1877  /// returns OBJECT_NOT_AVAILABLE immediately.
1878 
1879  if (factory == NULL)
1880  return OBJECT_NOT_AVAILABLE;
1881  if (tag != "")
1882  return OBJECT_NOT_AVAILABLE;
1883 
1884  vector<DTOFTruth*> data;
1885 
1886  const hddm_s::FtofTruthPointList &points = record->getFtofTruthPoints();
1887  hddm_s::FtofTruthPointList::iterator iter;
1888  for (iter = points.begin(); iter != points.end(); ++iter) {
1889  DTOFTruth *toftruth = new DTOFTruth;
1890  toftruth->primary = iter->getPrimary();
1891  toftruth->track = iter->getTrack();
1892  toftruth->x = iter->getX();
1893  toftruth->y = iter->getY();
1894  toftruth->z = iter->getZ();
1895  toftruth->t = iter->getT();
1896  toftruth->px = iter->getPx();
1897  toftruth->py = iter->getPy();
1898  toftruth->pz = iter->getPz();
1899  toftruth->E = iter->getE();
1900  toftruth->ptype = iter->getPtype();
1901  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
1902  toftruth->itrack = (ids.size())? ids.begin()->getItrack() : 0;
1903  data.push_back(toftruth);
1904  }
1905 
1906  // Copy into factory
1907  factory->CopyTo(data);
1908 
1909  return NOERROR;
1910 }
1911 
1912 //------------------
1913 // Extract_DTOFHit
1914 //------------------
1915 jerror_t DEventSourceHDDM::Extract_DTOFHit( hddm_s::HDDM *record,
1916  JFactory<DTOFHit>* factory,
1917  JFactory<DTOFHitMC> *factoryMC,
1918  string tag)
1919 {
1920  /// Copies the data from the given hddm_s structure. This is called
1921  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
1922  /// returns OBJECT_NOT_AVAILABLE immediately.
1923 
1924  if (factory == NULL)
1925  return OBJECT_NOT_AVAILABLE;
1926  if (tag != "" && tag != "TRUTH")
1927  return OBJECT_NOT_AVAILABLE;
1928 
1929  vector<DTOFHit*> data;
1930  vector<DTOFHitMC*> dataMC;
1931 
1932  const hddm_s::FtofCounterList &ctrs = record->getFtofCounters();
1933  hddm_s::FtofCounterList::iterator iter;
1934  for (iter = ctrs.begin(); iter != ctrs.end(); ++iter) {
1935  if (tag == "") {
1936  vector<DTOFHit*> north_hits;
1937  vector<DTOFHit*> south_hits;
1938 
1939  // Loop over north AND south hits
1940  const hddm_s::FtofHitList &hits = iter->getFtofHits();
1941  hddm_s::FtofHitList::iterator hiter;
1942  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1943  DTOFHit *tofhit = new DTOFHit;
1944  tofhit->bar = hiter->getBar();
1945  tofhit->plane = hiter->getPlane();
1946  tofhit->end = hiter->getEnd();
1947  tofhit->dE = hiter->getDE();
1948  tofhit->Amp = 0.;
1949  if(hiter->getFtofDigihits().size() > 0) {
1950  tofhit->Amp = hiter->getFtofDigihit().getPeakAmp();
1951  }
1952  tofhit->t = hiter->getT();
1953  tofhit->t_TDC = tofhit->t;
1954  tofhit->t_fADC= tofhit->t;
1955  tofhit->has_TDC=true;
1956  tofhit->has_fADC=true;
1957  data.push_back(tofhit);
1958  if (tofhit->end == 0)
1959  north_hits.push_back(tofhit);
1960  else
1961  south_hits.push_back(tofhit);
1962  }
1963 
1964  // return truth hits in a different factory
1965  const hddm_s::FtofTruthHitList &truths = iter->getFtofTruthHits();
1966  hddm_s::FtofTruthHitList::iterator titer;
1967  unsigned int north_mchits = 0;
1968  unsigned int south_mchits = 0;
1969  for (titer = truths.begin(); titer != truths.end(); ++titer) {
1970  DTOFHitMC *tofmchit = new DTOFHitMC;
1971  tofmchit->bar = titer->getBar();
1972  tofmchit->plane = titer->getPlane();
1973  tofmchit->end = titer->getEnd();
1974  tofmchit->itrack = titer->getFtofTruthExtra(0).getItrack();
1975  tofmchit->ptype = titer->getFtofTruthExtra(0).getPtype();
1976  tofmchit->dist = titer->getFtofTruthExtra(0).getDist();
1977  tofmchit->x = titer->getFtofTruthExtra(0).getX();
1978  tofmchit->y = titer->getFtofTruthExtra(0).getY();
1979  tofmchit->z = titer->getFtofTruthExtra(0).getZ();
1980  tofmchit->px = titer->getFtofTruthExtra(0).getPx();
1981  tofmchit->py = titer->getFtofTruthExtra(0).getPy();
1982  tofmchit->pz = titer->getFtofTruthExtra(0).getPz();
1983  tofmchit->E = titer->getFtofTruthExtra(0).getE();
1984  dataMC.push_back(tofmchit);
1985 
1986  // best-guess at tofhit-tofMChit association, not exact
1987  if (tofmchit->end == 0) {
1988  if (north_mchits < north_hits.size()) {
1989  north_hits[north_mchits]->AddAssociatedObject(tofmchit);
1990  }
1991  north_mchits++;
1992  }
1993  else {
1994  if (south_mchits < south_hits.size()) {
1995  south_hits[south_mchits]->AddAssociatedObject(tofmchit);
1996  }
1997  south_mchits++;
1998  }
1999  }
2000  }
2001 
2002  else if (tag == "TRUTH") {
2003  const hddm_s::FtofTruthHitList &truths = iter->getFtofTruthHits();
2004  hddm_s::FtofTruthHitList::iterator titer;
2005  for (titer = truths.begin(); titer != truths.end(); ++titer) {
2006  DTOFHit *tofhit = new DTOFHit;
2007  tofhit->bar = titer->getBar();
2008  tofhit->plane = titer->getPlane();
2009  tofhit->end = titer->getEnd();
2010  tofhit->dE = titer->getDE();
2011  tofhit->t = titer->getT();
2012  tofhit->t_fADC= tofhit->t;
2013  tofhit->t_TDC = tofhit->t;
2014  tofhit->has_TDC=true;
2015  tofhit->has_fADC=true;
2016  data.push_back(tofhit);
2017 
2018  DTOFHitMC *tofmchit = new DTOFHitMC;
2019  tofmchit->bar = tofhit->bar;
2020  tofmchit->plane = tofhit->plane;
2021  tofmchit->end = tofhit->end;
2022  tofmchit->itrack = titer->getFtofTruthExtra().getItrack();
2023  tofmchit->ptype = titer->getFtofTruthExtra().getPtype();
2024  tofmchit->dist = titer->getFtofTruthExtra().getDist();
2025  tofmchit->x = titer->getFtofTruthExtra().getX();
2026  tofmchit->y = titer->getFtofTruthExtra().getY();
2027  tofmchit->z = titer->getFtofTruthExtra().getZ();
2028  tofmchit->px = titer->getFtofTruthExtra().getPx();
2029  tofmchit->py = titer->getFtofTruthExtra().getPy();
2030  tofmchit->pz = titer->getFtofTruthExtra().getPz();
2031  tofmchit->E = titer->getFtofTruthExtra().getE();
2032  dataMC.push_back(tofmchit);
2033  tofhit->AddAssociatedObject(tofmchit);
2034  }
2035  }
2036  }
2037 
2038  // Copy into factory
2039  factory->CopyTo(data);
2040  factoryMC->CopyTo(dataMC);
2041 
2042  return NOERROR;
2043 }
2044 
2045 //------------------
2046 // Extract_DSCHit
2047 //------------------
2048 jerror_t DEventSourceHDDM::Extract_DSCHit(hddm_s::HDDM *record,
2049  JFactory<DSCHit>* factory, string tag)
2050 {
2051  /// Copies the data from the given hddm_s structure. This is called
2052  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2053  /// returns OBJECT_NOT_AVAILABLE immediately.
2054 
2055  if (factory == NULL)
2056  return OBJECT_NOT_AVAILABLE;
2057  if (tag != "" && tag != "TRUTH")
2058  return OBJECT_NOT_AVAILABLE;
2059 
2060  vector<DSCHit*> data;
2061 
2062  if (tag == "") {
2063  const hddm_s::StcHitList &hits = record->getStcHits();
2064  hddm_s::StcHitList::iterator iter;
2065  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2066  DSCHit *hit = new DSCHit;
2067  hit->sector = iter->getSector();
2068  hit->dE = iter->getDE();
2069  hit->t = iter->getT();
2070  hit->t_TDC=hit->t;
2071  hit->t_fADC=hit->t;
2072  hit->pulse_height = 0.;
2073  if(iter->getStcDigihits().size() > 0) {
2074  hit->pulse_height = iter->getStcDigihit().getPeakAmp();
2075  }
2076  hit->has_TDC=true;
2077  hit->has_fADC=true;
2078  data.push_back(hit);
2079  }
2080  }
2081  else if (tag == "TRUTH") {
2082  const hddm_s::StcTruthHitList &hits = record->getStcTruthHits();
2083  hddm_s::StcTruthHitList::iterator iter;
2084  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2085  DSCHit *hit = new DSCHit;
2086  hit->sector = iter->getSector();
2087  hit->dE = iter->getDE();
2088  hit->t = iter->getT();
2089  hit->t_TDC=hit->t;
2090  hit->t_fADC=hit->t;
2091  hit->has_TDC=true;
2092  hit->has_fADC=true;
2093  data.push_back(hit);
2094  }
2095  }
2096 
2097  // Copy into factory
2098  factory->CopyTo(data);
2099 
2100  return NOERROR;
2101 }
2102 
2103 //------------------
2104 // Extract_DSCTruthHit
2105 //------------------
2106 jerror_t DEventSourceHDDM::Extract_DSCTruthHit(hddm_s::HDDM *record,
2107  JFactory<DSCTruthHit>* factory, string tag)
2108 {
2109  /// Copies the data from the given hddm_s structure. This is called
2110  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2111  /// returns OBJECT_NOT_AVAILABLE immediately.
2112 
2113  if (factory == NULL)
2114  return OBJECT_NOT_AVAILABLE;
2115  if (tag != "")
2116  return OBJECT_NOT_AVAILABLE;
2117 
2118  vector<DSCTruthHit*> data;
2119 
2120  const hddm_s::StcTruthPointList &points = record->getStcTruthPoints();
2121  hddm_s::StcTruthPointList::iterator iter;
2122  for (iter = points.begin(); iter != points.end(); ++iter) {
2123  DSCTruthHit *hit = new DSCTruthHit;
2124  hit->dEdx = iter->getDEdx();
2125  hit->phi = iter->getPhi();
2126  hit->primary = iter->getPrimary();
2127  hit->ptype = iter->getPtype();
2128  hit->r = iter->getR();
2129  hit->t = iter->getT();
2130  hit->z = iter->getZ();
2131  hit->track = iter->getTrack();
2132  hit->sector = iter->getSector();
2133  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
2134  hit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
2135  data.push_back(hit);
2136  }
2137 
2138  // Copy into factory
2139  factory->CopyTo(data);
2140 
2141  return NOERROR;
2142 }
2143 
2144 //------------------
2145 // Extract_DTrackTimeBased
2146 //------------------
2147 jerror_t DEventSourceHDDM::Extract_DTrackTimeBased(hddm_s::HDDM *record,
2148  JFactory<DTrackTimeBased> *factory,
2149  string tag, int32_t runnumber, JEventLoop* locEventLoop)
2150 {
2151  // Note: Since this is a reconstructed factory, we want to generally return OBJECT_NOT_AVAILABLE
2152  // rather than NOERROR. The reason being that the caller interprets "NOERROR" to mean "yes I
2153  // usually can provide objects of that type, but this event has none." This will cause it to
2154  // skip any attempt at reconstruction. On the other hand, a value of "OBJECT_NOT_AVAILABLE" tells
2155  // it "I cannot provide those type of objects for this event."
2156 
2157  if (factory == NULL)
2158  return OBJECT_NOT_AVAILABLE;
2159  if (tag != "")
2160  return OBJECT_NOT_AVAILABLE;
2161 
2162  vector<DTrackTimeBased*> data;
2163  vector<DReferenceTrajectory*> rts;
2164 
2165  const hddm_s::TracktimebasedList &ttbs = record->getTracktimebaseds();
2166  hddm_s::TracktimebasedList::iterator iter;
2167  for (iter = ttbs.begin(); iter != ttbs.end(); ++iter) {
2168  DVector3 mom(iter->getMomentum().getPx(),
2169  iter->getMomentum().getPy(),
2170  iter->getMomentum().getPz());
2171  DVector3 pos(iter->getOrigin().getVx(),
2172  iter->getOrigin().getVy(),
2173  iter->getOrigin().getVz());
2175  track->setMomentum(mom);
2176  track->setPosition(pos);
2177  track->setPID(IDTrack(iter->getProperties().getCharge(),
2178  iter->getProperties().getMass()));
2179  track->chisq = iter->getChisq();
2180  track->Ndof = iter->getNdof();
2181  track->FOM = iter->getFOM();
2182  track->candidateid = iter->getCandidateid();
2183  track->id = iter->getId();
2184 
2185  // Reconstitute errorMatrix
2186  auto locCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
2187  locCovarianceMatrix->ResizeTo(7, 7);
2188  string str_vals = iter->getErrorMatrix().getVals();
2189  StringToTMatrixFSym(str_vals, locCovarianceMatrix.get(),
2190  iter->getErrorMatrix().getNrows(),
2191  iter->getErrorMatrix().getNcols());
2192  track->setErrorMatrix(locCovarianceMatrix);
2193 
2194  // Reconstitute TrackingErrorMatrix
2195  str_vals = iter->getTrackingErrorMatrix().getVals();
2196  auto TrackingErrorMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
2197  TrackingErrorMatrix->ResizeTo(5, 5);
2198  StringToTMatrixFSym(str_vals, TrackingErrorMatrix.get(),
2199  iter->getTrackingErrorMatrix().getNrows(),
2200  iter->getTrackingErrorMatrix().getNcols());
2201  track->setTrackingErrorMatrix(TrackingErrorMatrix);
2202 
2203  data.push_back(track);
2204  }
2205 
2206  // Copy into factory
2207  if (ttbs.size() > 0){
2208  factory->CopyTo(data);
2209 
2210  // If the event had a s_Tracktimebased_t pointer, then report
2211  // back that we read them in from the file. Otherwise, report
2212  // OBJECT_NOT_AVAILABLE
2213  return NOERROR;
2214  }
2215 
2216  // If we get to here then there was not even a placeholder in the HDDM file.
2217  // Return OBJECT_NOT_AVAILABLE to indicate reconstruction should be tried.
2218  return OBJECT_NOT_AVAILABLE;
2219 }
2220 
2221 
2222 //-------------------------------
2223 // StringToTMatrixFSym
2224 //-------------------------------
2225 string DEventSourceHDDM::StringToTMatrixFSym(string &str_vals, TMatrixFSym* mat,
2226  int Nrows, int Ncols)
2227 {
2228  /// This is the inverse of the DMatrixDSymToString method in the
2229  /// danahddm plugin.
2230 
2231  // Convert the given string into a symmetric matrix
2232  mat->ResizeTo(Nrows, Ncols);
2233  stringstream ss(str_vals);
2234  for (int irow=0; irow<mat->GetNrows(); irow++) {
2235  for (int icol=irow; icol<mat->GetNcols(); icol++) {
2236  ss >> (*mat)[irow][icol];
2237  (*mat)[icol][irow] = (*mat)[irow][icol];
2238  }
2239  }
2240 
2241  return ss.str();
2242 }
2243 
2244 //------------------
2245 // Extract_DTAGMHit
2246 //------------------
2247 jerror_t DEventSourceHDDM::Extract_DTAGMHit(hddm_s::HDDM *record,
2248  JFactory<DTAGMHit>* factory,
2249  string tag)
2250 {
2251  /// Copies the data from the given hddm_s structure. This is called
2252  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2253  /// returns OBJECT_NOT_AVAILABLE immediately.
2254 
2255  if (factory == NULL)
2256  return OBJECT_NOT_AVAILABLE;
2257  if (tag != "" && tag != "TRUTH")
2258  return OBJECT_NOT_AVAILABLE;
2259 
2260  vector<DTAGMHit*> data;
2261 
2262  // loop over microChannel/taggerHit records
2263  const hddm_s::MicroChannelList &tags = record->getMicroChannels();
2264  hddm_s::MicroChannelList::iterator iter;
2265  for (iter = tags.begin(); iter != tags.end(); ++iter) {
2266  if (tag == "") {
2267  const hddm_s::TaggerHitList &hits = iter->getTaggerHits();
2268  hddm_s::TaggerHitList::iterator hiter;
2269  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
2270  DTAGMHit *taghit = new DTAGMHit();
2271  taghit->E = hiter->getE();
2272  taghit->t = hiter->getT();
2273  taghit->npix_fadc = hiter->getNpe();
2274  taghit->time_fadc = hiter->getTADC();
2275  taghit->column = hiter->getColumn();
2276  taghit->row = hiter->getRow();
2277  taghit->has_fADC = true;
2278  taghit->has_TDC = true;
2279  data.push_back(taghit);
2280  }
2281  }
2282  else if (tag == "TRUTH") {
2283  const hddm_s::TaggerTruthHitList &hits = iter->getTaggerTruthHits();
2284  hddm_s::TaggerTruthHitList::iterator hiter;
2285  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
2286  DTAGMHit *taghit = new DTAGMHit();
2287  taghit->E = hiter->getE();
2288  taghit->t = hiter->getT();
2289  taghit->npix_fadc = hiter->getDE() * 1e5; // ~1e5 pixels/GeV
2290  taghit->time_fadc = hiter->getT();
2291  taghit->column = hiter->getColumn();
2292  taghit->row = hiter->getRow();
2293  taghit->has_fADC = true;
2294  taghit->has_TDC = true;
2295  taghit->bg = hiter->getBg();
2296  data.push_back(taghit);
2297  }
2298  }
2299  }
2300 
2301  // Copy into factory
2302  factory->CopyTo(data);
2303 
2304  return NOERROR;
2305 }
2306 
2307 //------------------
2308 // Extract_DTAGHHit
2309 //------------------
2310 jerror_t DEventSourceHDDM::Extract_DTAGHHit( hddm_s::HDDM *record,
2311  JFactory<DTAGHHit>* factory,
2312  string tag)
2313 {
2314  /// Copies the data from the given hddm_s structure. This is called
2315  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2316  /// returns OBJECT_NOT_AVAILABLE immediately.
2317 
2318  if (factory == NULL)
2319  return OBJECT_NOT_AVAILABLE;
2320  if (tag != "" && tag != "TRUTH")
2321  return OBJECT_NOT_AVAILABLE;
2322 
2323  vector<DTAGHHit*> data;
2324 
2325  // loop over hodoChannel/taggerHit records
2326  const hddm_s::HodoChannelList &tags = record->getHodoChannels();
2327  hddm_s::HodoChannelList::iterator iter;
2328  for (iter = tags.begin(); iter != tags.end(); ++iter) {
2329  if (tag == "") {
2330  const hddm_s::TaggerHitList &hits = iter->getTaggerHits();
2331  hddm_s::TaggerHitList::iterator hiter;
2332  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
2333  DTAGHHit *taghit = new DTAGHHit();
2334  taghit->E = hiter->getE();
2335  taghit->t = hiter->getT();
2336  taghit->npe_fadc = hiter->getNpe();
2337  taghit->time_fadc = hiter->getTADC();
2338  taghit->counter_id = hiter->getCounterId();
2339  taghit->has_fADC = true;
2340  taghit->has_TDC = true;
2341  data.push_back(taghit);
2342  }
2343  }
2344  else if (tag == "TRUTH") {
2345  const hddm_s::TaggerTruthHitList &hits = iter->getTaggerTruthHits();
2346  hddm_s::TaggerTruthHitList::iterator hiter;
2347  for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
2348  DTAGHHit *taghit = new DTAGHHit();
2349  taghit->E = hiter->getE();
2350  taghit->t = hiter->getT();
2351  taghit->npe_fadc = hiter->getDE() * 5e5; // ~5e5 pe/GeV
2352  taghit->time_fadc = hiter->getT();
2353  taghit->counter_id = hiter->getCounterId();
2354  taghit->has_fADC = true;
2355  taghit->has_TDC = true;
2356  taghit->bg = hiter->getBg();
2357 
2358  data.push_back(taghit);
2359  }
2360  }
2361  }
2362 
2363  // Copy into factory
2364  factory->CopyTo(data);
2365 
2366  return NOERROR;
2367 }
2368 
2369 //------------------
2370 // Extract_DPSHit
2371 //------------------
2372 jerror_t DEventSourceHDDM::Extract_DPSHit(hddm_s::HDDM *record,
2373  JFactory<DPSHit>* factory, string tag)
2374 {
2375  /// Copies the data from the given hddm_s structure. This is called
2376  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2377  /// returns OBJECT_NOT_AVAILABLE immediately.
2378 
2379  if (factory == NULL)
2380  return OBJECT_NOT_AVAILABLE;
2381  if (tag != "" && tag != "TRUTH")
2382  return OBJECT_NOT_AVAILABLE;
2383 
2384  vector<DPSHit*> data;
2385 
2386  if (tag == "") {
2387  const hddm_s::PsHitList &hits = record->getPsHits();
2388  hddm_s::PsHitList::iterator iter;
2389  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2390  DPSHit *hit = new DPSHit;
2391  if(iter->getArm() == 0)
2392  hit->arm = DPSGeometry::Arm::kNorth;
2393  else
2394  hit->arm = DPSGeometry::Arm::kSouth;
2395  hit->column = iter->getColumn();
2396  double npix_fadc = iter->getDE()*0.5e5; // 100 pixels in 2 MeV
2397  hit->npix_fadc = npix_fadc;
2398  hit->t = iter->getT();
2399 
2400  hit->E = 0.5*(psGeom->getElow(hit->arm,hit->column) + psGeom->getEhigh(hit->arm,hit->column));
2401  hit->pulse_peak = npix_fadc*21; // 1 pixel 21 fadc counts
2402  hit->integral = npix_fadc*21*5.1; // integral/peak = 5.1
2403  data.push_back(hit);
2404  }
2405  }
2406  else if (tag == "TRUTH") {
2407  const hddm_s::PsTruthHitList &hits = record->getPsTruthHits();
2408  hddm_s::PsTruthHitList::iterator iter;
2409  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2410  DPSHit *hit = new DPSHit;
2411  if(iter->getArm() == 0)
2412  hit->arm = DPSGeometry::Arm::kNorth;
2413  else
2414  hit->arm = DPSGeometry::Arm::kSouth;
2415  hit->column = iter->getColumn();
2416  hit->npix_fadc = iter->getDE() * 1e5; // ~1e5 pixels/GeV
2417  hit->t = iter->getT();
2418  data.push_back(hit);
2419  }
2420  }
2421 
2422  // Copy into factory
2423  factory->CopyTo(data);
2424 
2425  return NOERROR;
2426 }
2427 
2428 //------------------
2429 // Extract_DPSTruthHit
2430 //------------------
2431 jerror_t DEventSourceHDDM::Extract_DPSTruthHit(hddm_s::HDDM *record,
2432  JFactory<DPSTruthHit>* factory, string tag)
2433 {
2434  /// Copies the data from the given hddm_s structure. This is called
2435  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2436  /// returns OBJECT_NOT_AVAILABLE immediately.
2437 
2438  if (factory == NULL)
2439  return OBJECT_NOT_AVAILABLE;
2440  if (tag != "")
2441  return OBJECT_NOT_AVAILABLE;
2442 
2443  vector<DPSTruthHit*> data;
2444 
2445  const hddm_s::PsTruthPointList &points = record->getPsTruthPoints();
2446  hddm_s::PsTruthPointList::iterator iter;
2447  for (iter = points.begin(); iter != points.end(); ++iter) {
2448  DPSTruthHit *hit = new DPSTruthHit;
2449  hit->dEdx = iter->getDEdx();
2450  hit->primary = iter->getPrimary();
2451  hit->ptype = iter->getPtype();
2452  hit->t = iter->getT();
2453  hit->x = iter->getX();
2454  hit->y = iter->getY();
2455  hit->z = iter->getZ();
2456  hit->track = iter->getTrack();
2457  hit->column = iter->getColumn();
2458  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
2459  hit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
2460  data.push_back(hit);
2461  }
2462 
2463  // Copy into factory
2464  factory->CopyTo(data);
2465 
2466  return NOERROR;
2467 }
2468 
2469 //------------------
2470 // Extract_DPSCHit
2471 //------------------
2472 jerror_t DEventSourceHDDM::Extract_DPSCHit(hddm_s::HDDM *record,
2473  JFactory<DPSCHit>* factory, string tag)
2474 {
2475  /// Copies the data from the given hddm_s structure. This is called
2476  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2477  /// returns OBJECT_NOT_AVAILABLE immediately.
2478 
2479  if (factory == NULL)
2480  return OBJECT_NOT_AVAILABLE;
2481  if (tag != "" && tag != "TRUTH")
2482  return OBJECT_NOT_AVAILABLE;
2483 
2484  vector<DPSCHit*> data;
2485 
2486  if (tag == "") {
2487  const hddm_s::PscHitList &hits = record->getPscHits();
2488  hddm_s::PscHitList::iterator iter;
2489  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2490  DPSCHit *hit = new DPSCHit;
2491  if(iter->getArm() == 0)
2492  hit->arm = DPSGeometry::Arm::kNorth;
2493  else
2494  hit->arm = DPSGeometry::Arm::kSouth;
2495  hit->module = iter->getModule();
2496 
2497  double npe_fadc = iter->getDE()*2.5e5;
2498  hit->npe_fadc = npe_fadc;
2499  hit->pulse_peak = npe_fadc*0.4; // 1000 pe - 400 fadc count
2500  hit->integral = npe_fadc*0.4*3; // integral/peak = 3.
2501 
2502  hit->t = iter->getT();
2503  hit->time_tdc = iter->getT();
2504  hit->time_fadc = iter->getT();
2505 
2506  hit->has_fADC = true;
2507  hit->has_TDC = true;
2508 
2509  data.push_back(hit);
2510  }
2511  }
2512  else if (tag == "TRUTH") {
2513  const hddm_s::PscTruthHitList &hits = record->getPscTruthHits();
2514  hddm_s::PscTruthHitList::iterator iter;
2515  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2516  DPSCHit *hit = new DPSCHit;
2517  if(iter->getArm() == 0)
2518  hit->arm = DPSGeometry::Arm::kNorth;
2519  else
2520  hit->arm = DPSGeometry::Arm::kSouth;
2521  hit->module = iter->getModule();
2522  hit->npe_fadc = iter->getDE() * 5e5; // ~5e5 pe/GeV
2523  hit->t = iter->getT();
2524  data.push_back(hit);
2525  }
2526  }
2527 
2528  // Copy into factory
2529  factory->CopyTo(data);
2530 
2531  return NOERROR;
2532 }
2533 
2534 //------------------
2535 // Extract_DPSCTruthHit
2536 //------------------
2537 jerror_t DEventSourceHDDM::Extract_DPSCTruthHit(hddm_s::HDDM *record,
2538  JFactory<DPSCTruthHit>* factory, string tag)
2539 {
2540  /// Copies the data from the given hddm_s structure. This is called
2541  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2542  /// returns OBJECT_NOT_AVAILABLE immediately.
2543 
2544  if (factory == NULL)
2545  return OBJECT_NOT_AVAILABLE;
2546  if (tag != "")
2547  return OBJECT_NOT_AVAILABLE;
2548 
2549  vector<DPSCTruthHit*> data;
2550 
2551  const hddm_s::PscTruthPointList &points = record->getPscTruthPoints();
2552  hddm_s::PscTruthPointList::iterator iter;
2553  for (iter = points.begin(); iter != points.end(); ++iter) {
2554  DPSCTruthHit *hit = new DPSCTruthHit;
2555  hit->dEdx = iter->getDEdx();
2556  hit->primary = iter->getPrimary();
2557  hit->ptype = iter->getPtype();
2558  hit->t = iter->getT();
2559  hit->x = iter->getX();
2560  hit->y = iter->getY();
2561  hit->z = iter->getZ();
2562  hit->track = iter->getTrack();
2563  hit->column = iter->getModule();
2564  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
2565  hit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
2566  data.push_back(hit);
2567  }
2568 
2569  // Copy into factory
2570  factory->CopyTo(data);
2571 
2572  return NOERROR;
2573 }
2574 
2575 //------------------
2576 // Etract_DTPOLHit
2577 //------------------
2578 jerror_t DEventSourceHDDM::Extract_DTPOLHit(hddm_s::HDDM *record,
2579  JFactory<DTPOLHit>* factory, string tag)
2580 {
2581  if (factory == NULL)
2582  return OBJECT_NOT_AVAILABLE;
2583  if (tag != "")
2584  return OBJECT_NOT_AVAILABLE;
2585 
2586  vector<DTPOLHit*> data;
2587 
2588  if (tag == "")
2589  {
2590  const hddm_s::TpolHitList &hits = record->getTpolHits();
2591  hddm_s::TpolHitList::iterator iter;
2592  for (iter = hits.begin(); iter != hits.end(); ++iter)
2593  {
2594  DTPOLHit *hit = new DTPOLHit;
2595  hit->sector = iter->getSector();
2596  hit->ring = iter->getRing();
2597  hit->dE = iter->getDE();
2598  hit->t = iter->getT();
2599 
2600  data.push_back(hit);
2601  }
2602  }
2603  else if (tag == "Truth")
2604  {
2605  const hddm_s::TpolTruthHitList &truthHits = record->getTpolTruthHits();
2606  hddm_s::TpolTruthHitList::iterator iter;
2607  for (iter = truthHits.begin(); iter != truthHits.end(); ++iter)
2608  {
2609  DTPOLHit *hit = new DTPOLHit;
2610  hit->sector = iter->getSector();
2611  hit->t = iter->getT();
2612 
2613  data.push_back(hit);
2614  }
2615  }
2616 
2617  factory->CopyTo(data);
2618 
2619  return NOERROR;
2620 }
2621 
2622 //------------------------
2623 // Extract_DTPOLTruthHit
2624 //------------------------
2625 jerror_t DEventSourceHDDM::Extract_DTPOLTruthHit(hddm_s::HDDM *record, JFactory<DTPOLTruthHit>* factory, string tag)
2626 {
2627  if (factory == NULL)
2628  return OBJECT_NOT_AVAILABLE;
2629  if (tag != "")
2630  return OBJECT_NOT_AVAILABLE;
2631 
2632  vector<DTPOLTruthHit*> data;
2633 
2634  const hddm_s::TpolTruthPointList &points = record->getTpolTruthPoints();
2635  hddm_s::TpolTruthPointList::iterator iter;
2636  for (iter = points.begin(); iter != points.end(); ++iter)
2637  {
2638  DTPOLTruthHit *hit = new DTPOLTruthHit;
2639  hit->dEdx = iter->getDEdx();
2640  hit->primary = iter->getPrimary();
2641  hit->track = iter->getTrack();
2642  hit->ptype = iter->getPtype();
2643  hit->r = iter->getR();
2644  hit->phi = iter->getPhi();
2645  hit->t = iter->getT();
2646  const hddm_s::TrackIDList &ids = iter->getTrackIDs();
2647  hit->itrack = (ids.size())? ids.begin()->getItrack() : 0;
2648 
2649  data.push_back(hit);
2650  }
2651 
2652  factory->CopyTo(data);
2653 
2654  return NOERROR;
2655 }
2656 
2657 Particle_t DEventSourceHDDM::IDTrack(float locCharge, float locMass) const
2658 {
2659  float locMassTolerance = 0.010;
2660  if (locCharge > 0.1) // Positive particles
2661  {
2662  if (fabs(locMass - ParticleMass(Proton)) < locMassTolerance) return Proton;
2663  if (fabs(locMass - ParticleMass(PiPlus)) < locMassTolerance) return PiPlus;
2664  if (fabs(locMass - ParticleMass(KPlus)) < locMassTolerance) return KPlus;
2665  if (fabs(locMass - ParticleMass(Positron)) < locMassTolerance) return Positron;
2666  if (fabs(locMass - ParticleMass(MuonPlus)) < locMassTolerance) return MuonPlus;
2667  }
2668  else if(locCharge < -0.1) // Negative particles
2669  {
2670  if (fabs(locMass - ParticleMass(PiMinus)) < locMassTolerance) return PiMinus;
2671  if (fabs(locMass - ParticleMass(KMinus)) < locMassTolerance) return KMinus;
2672  if (fabs(locMass - ParticleMass(MuonMinus)) < locMassTolerance) return MuonMinus;
2673  if (fabs(locMass - ParticleMass(Electron)) < locMassTolerance) return Electron;
2674  if (fabs(locMass - ParticleMass(AntiProton)) < locMassTolerance) return AntiProton;
2675  }
2676  else //Neutral Track
2677  {
2678  if (fabs(locMass - ParticleMass(Gamma)) < locMassTolerance) return Gamma;
2679  if (fabs(locMass - ParticleMass(Neutron)) < locMassTolerance) return Neutron;
2680  }
2681  return Unknown;
2682 }
2683 
2684 //------------------
2685 // Extract_DFMWPCTruthHit
2686 //------------------
2687 jerror_t DEventSourceHDDM::Extract_DFMWPCTruthHit(hddm_s::HDDM *record, JFactory<DFMWPCTruthHit> *factory, string tag)
2688 {
2689  /// Copies the data from the given hddm_s record. This is called
2690  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2691  /// returns OBJECT_NOT_AVAILABLE immediately.
2692 
2693  if (factory == NULL) return OBJECT_NOT_AVAILABLE;
2694  if (tag != "") return OBJECT_NOT_AVAILABLE;
2695 
2696  vector<DFMWPCTruthHit*> data;
2697 
2698  const hddm_s::FmwpcTruthHitList &points = record->getFmwpcTruthHits();
2699  hddm_s::FmwpcTruthHitList::iterator iter;
2700  for (iter = points.begin(); iter != points.end(); ++iter) {
2701  DFMWPCTruthHit *hit = new DFMWPCTruthHit;
2702  hit->layer = iter->getLayer();
2703  hit->wire = iter->getWire();
2704  hit->dE = iter->getDE();
2705  hit->dx = iter->getDx();
2706  hit->t = iter->getT();
2707  data.push_back(hit);
2708  }
2709 
2710  // Copy into factory
2711  factory->CopyTo(data);
2712 
2713  return NOERROR;
2714 }
2715 
2716 //------------------
2717 // Extract_DFMWPCHit
2718 //------------------
2719 jerror_t DEventSourceHDDM::Extract_DFMWPCHit(hddm_s::HDDM *record, JFactory<DFMWPCHit> *factory, string tag)
2720 {
2721  /// Copies the data from the given hddm_s record. This is called
2722  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2723  /// returns OBJECT_NOT_AVAILABLE immediately.
2724 
2725  if (factory == NULL) return OBJECT_NOT_AVAILABLE;
2726  if (tag != "") return OBJECT_NOT_AVAILABLE;
2727 
2728  vector<DFMWPCHit*> data;
2729 
2730  const hddm_s::FmwpcHitList &points = record->getFmwpcHits();
2731  hddm_s::FmwpcHitList::iterator iter;
2732  for (iter = points.begin(); iter != points.end(); ++iter) {
2733  DFMWPCHit *hit = new DFMWPCHit;
2734  hit->layer = iter->getLayer();
2735  hit->wire = iter->getWire();
2736  hit->dE = iter->getDE();
2737  hit->t = iter->getT();
2738  data.push_back(hit);
2739  }
2740 
2741  // Copy into factory
2742  factory->CopyTo(data);
2743 
2744  return NOERROR;
2745 }
2746 
2747 //------------------
2748 // Extract_DDIRCPmtHit
2749 //------------------
2750 jerror_t DEventSourceHDDM::Extract_DDIRCPmtHit(hddm_s::HDDM *record,
2751  JFactory<DDIRCPmtHit> *factory, string tag,
2752  JEventLoop* eventLoop)
2753 {
2754  /// Copies the data from the given hddm_s structure. This is called
2755  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2756  /// returs OBJECT_NOT_AVAILABLE immediately.
2757 
2758  if (factory == NULL)
2759  return OBJECT_NOT_AVAILABLE;
2760  if (tag != "")
2761  return OBJECT_NOT_AVAILABLE;
2762 
2763  vector<DDIRCPmtHit*> data;
2764 
2765  if (tag == "") {
2766  vector<const DDIRCTruthPmtHit*> locDIRCTruthPmtHit;
2767  eventLoop->Get(locDIRCTruthPmtHit);
2768 
2769  const hddm_s::DircPmtHitList &hits = record->getDircPmtHits();
2770  hddm_s::DircPmtHitList::iterator iter;
2771  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2772  double time = iter->getT();
2773  int channel = iter->getCh();
2774 
2775  DDIRCPmtHit *hit = new DDIRCPmtHit();
2776  hit->t = time;
2777  hit->ch = channel;
2778 
2779  for (auto& iterTruth : locDIRCTruthPmtHit) { //.begin(); iterTruth != locDIRCTruthPmtHit.end(); ++iterTruth) {
2780 
2781  // must match channel and time
2782  if(channel == iterTruth->ch && fabs(time-iterTruth->t) < 5.0) {
2783 
2784  hit->AddAssociatedObject(iterTruth);
2785 
2786  break;
2787  }
2788  }
2789 
2790  data.push_back(hit);
2791  }
2792  }
2793 
2794  // Copy into factory
2795  factory->CopyTo(data);
2796 
2797  return NOERROR;
2798 }
2799 
2800 //------------------
2801 // Extract_DCereHit
2802 // added by yqiang Oct 11, 2012
2803 //------------------
2804 jerror_t DEventSourceHDDM::Extract_DCereHit(hddm_s::HDDM *record,
2805  JFactory<DCereHit>* factory, string tag)
2806 {
2807  if (factory == NULL)
2808  return OBJECT_NOT_AVAILABLE;
2809  if (tag != "" && tag != "TRUTH")
2810  return OBJECT_NOT_AVAILABLE;
2811 
2812  vector<DCereHit*> data;
2813 
2814  if (tag == "") {
2815  const hddm_s::CereHitList &hits = record->getCereHits();
2816  hddm_s::CereHitList::iterator iter;
2817  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2818  DCereHit *hit = new DCereHit;
2819  hit->sector = iter->getSector();
2820  hit->pe = iter->getPe();
2821  hit->t = iter->getT();
2822  data.push_back(hit);
2823  }
2824  }
2825  else if (tag == "TRUTH") {
2826  const hddm_s::CereTruthHitList &hits = record->getCereTruthHits();
2827  hddm_s::CereTruthHitList::iterator iter;
2828  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2829  DCereHit *hit = new DCereHit;
2830  hit->sector = iter->getSector();
2831  hit->pe = iter->getPe();
2832  hit->t = iter->getT();
2833  data.push_back(hit);
2834  }
2835  }
2836 
2837  // copy into factory
2838  factory->CopyTo(data);
2839 
2840  return NOERROR;
2841 }
2842 
2843 //------------------
2844 // Extract_DDIRCTruthBarHit
2845 //------------------
2846 jerror_t DEventSourceHDDM::Extract_DDIRCTruthBarHit(hddm_s::HDDM *record,
2847  JFactory<DDIRCTruthBarHit>* factory, string tag)
2848 {
2849  /// Copies the data from the given hddm_s structure. This is called
2850  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2851  /// returns OBJECT_NOT_AVAILABLE immediately.
2852 
2853  if (factory == NULL)
2854  return OBJECT_NOT_AVAILABLE;
2855  if (tag != "")
2856  return OBJECT_NOT_AVAILABLE;
2857 
2858  vector<DDIRCTruthBarHit*> data;
2859 
2860  const hddm_s::DircTruthBarHitList &hits = record->getDircTruthBarHits();
2861  hddm_s::DircTruthBarHitList::iterator iter;
2862  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2864  hit->x = iter->getX();
2865  hit->y = iter->getY();
2866  hit->z = iter->getZ();
2867  hit->px = iter->getPx();
2868  hit->py = iter->getPy();
2869  hit->pz = iter->getPz();
2870  hit->t = iter->getT();
2871  hit->E = iter->getE();
2872  hit->pdg = iter->getPdg();
2873  hit->bar = iter->getBar();
2874  hit->track = iter->getTrack();
2875  data.push_back(hit);
2876  }
2877 
2878  // Copy into factory
2879  factory->CopyTo(data);
2880 
2881  return NOERROR;
2882 }
2883 
2884 //------------------
2885 // Extract_DDIRCTruthPmtHit
2886 //------------------
2887 jerror_t DEventSourceHDDM::Extract_DDIRCTruthPmtHit(hddm_s::HDDM *record,
2888  JFactory<DDIRCTruthPmtHit>* factory, string tag)
2889 {
2890  /// Copies the data from the given hddm_s structure. This is called
2891  /// from JEventSourceHDDM::GetObjects. If factory is NULL, this
2892  /// returns OBJECT_NOT_AVAILABLE immediately.
2893 
2894  if (factory == NULL)
2895  return OBJECT_NOT_AVAILABLE;
2896  if (tag != "")
2897  return OBJECT_NOT_AVAILABLE;
2898 
2899  vector<DDIRCTruthPmtHit*> data;
2900 
2901 
2902  const hddm_s::DircTruthPmtHitList &hits = record->getDircTruthPmtHits();
2903  hddm_s::DircTruthPmtHitList::iterator iter;
2904  for (iter = hits.begin(); iter != hits.end(); ++iter) {
2906  hit->x = iter->getX();
2907  hit->y = iter->getY();
2908  hit->z = iter->getZ();
2909  hit->t = iter->getT();
2910  hit->E = iter->getE();
2911  hit->ch = iter->getCh();
2912  hit->key_bar = iter->getKey_bar();
2913  hddm_s::DircTruthPmtHitExtraList &hitextras = iter->getDircTruthPmtHitExtras();
2914  if(hitextras.size() > 0) {
2915  hit->t_fixed = hitextras(0).getT_fixed();
2916  hit->path = hitextras(0).getPath();
2917  hit->refl = hitextras(0).getRefl();
2918  hit->bbrefl = hitextras(0).getBbrefl();
2919  }
2920  data.push_back(hit);
2921  }
2922 
2923  // Copy into factory
2924  factory->CopyTo(data);
2925 
2926  return NOERROR;
2927 }
float q
Definition: DFDCHit.h:29
float t
Definition: DCCALHit.h:27
int sector
Definition: DTPOLHit.h:11
jerror_t Extract_DMCThrown(hddm_s::HDDM *record, JFactory< DMCThrown > *factory, string tag)
DBCALGeometry::End end
Definition: DBCALSiPMHit.h:35
Definition: track.h:16
double t
Definition: DTAGHHit.h:19
float x
Definition: DCCALHit.h:24
void setMomentum(const DVector3 &aMomentum)
DVector2 positionOnFace(int row, int column) const
double E
Definition: DTAGHHit.h:18
void setTime(double locTime)
int itrack
Definition: DCDCHit.h:25
DKinematicData target
Definition: DMCReaction.h:24
jerror_t GetTOFTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
hddm_s::istream * fin
double E
Definition: DTAGMHit.h:18
double npe_fadc
Definition: DTAGHHit.h:25
jerror_t Extract_DTAGHHit(hddm_s::HDDM *record, JFactory< DTAGHHit > *factory, string tag)
double integral
Definition: DPSHit.h:23
int module
Definition: DPSCHit.h:20
int row
Definition: DCCALHit.h:22
jerror_t Extract_DBCALTruthCell(hddm_s::HDDM *record, JFactory< DBCALTruthCell > *factory, string tag)
jerror_t GetObjects(JEvent &event, JFactory_base *factory)
Definition: DPSHit.h:15
int bg
Definition: DTAGHHit.h:27
int end
Definition: DTOFHit.h:23
const int Ncols
float pulse_height
Definition: DSCHit.h:23
const DMagneticFieldMap * bfield
float z
Definition: DTOFTruth.h:23
double pulse_peak
Definition: DPSCHit.h:23
jerror_t Extract_DTOFTruth(hddm_s::HDDM *record, JFactory< DTOFTruth > *factory, string tag)
JCalibration * jcalib
double npix_fadc
Definition: DPSHit.h:25
jerror_t Extract_DDIRCTruthBarHit(hddm_s::HDDM *record, JFactory< DDIRCTruthBarHit > *factory, string tag)
int type
Definition: DFDCHit.h:44
TVector2 DVector2
Definition: DVector2.h:9
int column
Definition: DCCALHit.h:23
float x
Definition: DTOFTruth.h:23
double integral
Definition: DPSCHit.h:22
bool has_TDC
Definition: DSCHit.h:25
TVector3 DVector3
Definition: DVector3.h:14
int plane
Definition: DTOFHit.h:21
bool has_fADC
Definition: DTOFHit.h:29
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
DBCALGeometry::End end
static int gPlane(const DFDCHit *h)
DFDCGeometry::gPlane(): Compute the global plane (1-74) number based on module, layer, and plane.
Definition: DFDCGeometry.h:67
double dTimeVariance
Definition: DRFTime.h:25
char string[256]
int row
Definition: DTAGMHit.h:20
int track
This is the unique number that GEANT has assigned the particle.
bool has_TDC
Definition: DTAGHHit.h:26
float E
Definition: DTOFTruth.h:26
#define y
float amp
Definition: DCDCHit.h:21
int ptype
Definition: DTOFHitMC.h:24
double dE
Definition: DTPOLHit.h:17
jerror_t Extract_DPSCTruthHit(hddm_s::HDDM *record, JFactory< DPSCTruthHit > *factory, string tag)
int primary
Definition: DTOFTruth.h:22
double dTime
Definition: DRFTime.h:24
DPSGeometry::Arm arm
Definition: DPSCHit.h:19
int gLayer
Definition: DFDCHit.h:28
float pulse_height
Definition: DFDCHit.h:30
int itrack
This is the index within the MCThrown structure of this track.
int bg
Definition: DTAGMHit.h:28
const DGeometry * geom
double E
Definition: DPSHit.h:21
jerror_t Extract_DBCALTDCDigiHit(hddm_s::HDDM *record, JFactory< DBCALTDCDigiHit > *factory, string tag)
jerror_t Extract_DTPOLHit(hddm_s::HDDM *record, JFactory< DTPOLHit > *factory, string tag)
double time_fadc
Definition: DTAGMHit.h:25
int ring
Definition: DTPOLHit.h:13
int sector
Definition: DCereHit.h:20
uint32_t nsamples_integral
number of samples used in integral
Definition: DBCALDigiHit.h:34
float y
Definition: DTOFHitMC.h:28
jerror_t Extract_DFDCHit(hddm_s::HDDM *record, JFactory< DFDCHit > *factory, string tag)
double getEhigh(int arm, int column) const
Definition: DPSGeometry.cc:48
float t
Definition: DTOFHit.h:28
int ptype
Definition: DFDCHit.h:48
double t
Definition: DTAGMHit.h:19
int layer
Definition: DFDCHit.h:23
jerror_t Extract_DMCTrackHit(hddm_s::HDDM *record, JFactory< DMCTrackHit > *factory, string tag)
jerror_t Extract_DDIRCTruthPmtHit(hddm_s::HDDM *record, JFactory< DDIRCTruthPmtHit > *factory, string tag)
int pdgtype
PDG particle type (not used by GEANT)
Definition: DMCThrown.h:21
double t
Definition: DPSCHit.h:21
Definition: GlueX.h:27
DVector2 positionOnFace(int row, int column) const
int ptype
This is the particle ID number.
DApplication * dapp
bool has_TDC
Definition: DPSCHit.h:27
Definition: GlueX.h:17
Definition: GlueX.h:19
float y
Definition: DTOFTruth.h:23
int sector
Definition: DSCHit.h:18
DEventSourceHDDM(const char *source_name)
jerror_t GetFCALTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
DetectorSystem_t system
particle type
Definition: DMCTrackHit.h:25
float pz
Definition: DTOFHitMC.h:32
DKinematicData beam
Definition: DMCReaction.h:25
DBCALGeometry::End end
Definition: DBCALDigiHit.h:28
bool has_fADC
Definition: DPSCHit.h:27
int track
Definition: DTOFTruth.h:20
float t_TDC
Definition: DTOFHit.h:27
uint32_t datasource
0=window raw data, 1=old(pre-Fall16) firmware, 2=Df250PulseData, 3=MC
Definition: DBCALDigiHit.h:37
uint32_t pulse_peak
identified pulse height as returned by FPGA algorithm
Definition: DBCALDigiHit.h:30
double weight
Definition: DMCReaction.h:23
jerror_t Extract_DMCTrajectoryPoint(hddm_s::HDDM *record, JFactory< DMCTrajectoryPoint > *factory, string tag)
jerror_t Extract_DFMWPCHit(hddm_s::HDDM *record, JFactory< DFMWPCHit > *factory, string tag)
jerror_t GetCCALTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
float Amp
Definition: DTOFHit.h:25
jerror_t Extract_DTOFHit(hddm_s::HDDM *record, JFactory< DTOFHit > *factory, JFactory< DTOFHitMC > *factoryMC, string tag)
int counter_id
Definition: DTAGHHit.h:20
bool has_fADC
Definition: DTAGMHit.h:27
jerror_t Extract_DSCTruthHit(hddm_s::HDDM *record, JFactory< DSCTruthHit > *factory, string tag)
Definition: DSCHit.h:14
jerror_t GetBCALTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
float t
Definition: DCDCHit.h:22
uint32_t pulse_integral
identified pulse integral as returned by FPGA algorithm
Definition: DBCALDigiHit.h:29
float dE
Definition: DFMWPCHit.h:19
float t
Definition: DFMWPCHit.h:20
static int gLayer(const DFDCHit *h)
DFDCGeometry::gLayer(): Compute the global layer (detection layer 1-24) number based on module and la...
Definition: DFDCGeometry.h:59
map< unsigned int, double > dBeamBunchPeriodMap
int mech
production mechanism of this partcle (generator specific)
Definition: DMCThrown.h:24
float z
coordinates of hit in cm and rad
Definition: DMCTrackHit.h:20
jerror_t Extract_DBCALTruthShower(hddm_s::HDDM *record, JFactory< DBCALTruthShower > *factory, string tag)
jerror_t Extract_DCCALHit(hddm_s::HDDM *record, JFactory< DCCALHit > *factory, string tag, JEventLoop *eventLoop)
Definition: GlueX.h:20
static float getStripR(const DFDCHit *h)
DFDCGeometry::getStripR(): Return coordinate in U or V space of a strip.
Definition: DFDCGeometry.h:84
Definition: GlueX.h:18
double npe_fadc
Definition: DPSCHit.h:26
float dist
Definition: DTOFHitMC.h:26
float x
Definition: DTOFHitMC.h:27
int ptype
Definition: DCDCHit.h:26
int row
Definition: DFCALHit.h:23
uint32_t nsamples_pedestal
number of samples used in pedestal
Definition: DBCALDigiHit.h:35
int parentid
id of parent of this particle from original generator
Definition: DMCThrown.h:23
jerror_t Extract_DCCALTruthShower(hddm_s::HDDM *record, JFactory< DCCALTruthShower > *factory, string tag)
DGeometry * GetDGeometry(unsigned int run_number)
Definition: GlueX.h:22
float dEdx
Definition: DPSTruthHit.h:13
bool primary
Definition: DSCTruthHit.h:19
int end
Definition: DTOFHitMC.h:23
const DPSGeometry * psGeom
float py
Definition: DTOFHitMC.h:31
jerror_t GetCherenkovTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
float px
Definition: DTOFTruth.h:24
int plane
Definition: DFDCHit.h:26
File: DTOFHitMC.h Created: Wed Jan 19 14:22:41 EST 2011 Creator: B. Zihlmann Purpose: Container class...
Definition: DTOFHitMC.h:16
#define _DBG_
Definition: HDEVIO.h:12
int ring
Definition: DCDCHit.h:18
int layer
Definition: DFMWPCHit.h:17
bool has_TDC
Definition: DTAGMHit.h:27
float d
Definition: DFDCHit.h:35
const DBCALGeometry * dBCALGeom
map< unsigned int, double > dTargetCenterZMap
TFile * Target
Definition: root_merge.cc:40
DBCALGeometry::End end
jerror_t Extract_DMCReaction(hddm_s::HDDM *record, JFactory< DMCReaction > *factory, string tag, JEventLoop *loop)
int wire
Definition: DFMWPCHit.h:18
uint32_t QF
Quality Factor from FPGA algorithms.
Definition: DBCALDigiHit.h:33
float px
Definition: DTOFHitMC.h:30
float py
Definition: DTOFTruth.h:24
float E
Definition: DTOFHitMC.h:33
jerror_t Extract_DTrackTimeBased(hddm_s::HDDM *record, JFactory< DTrackTimeBased > *factory, string tag, int32_t runnumber, JEventLoop *locEventLoop)
void FreeEvent(JEvent &event)
double charge(void) const
bool has_TDC
Definition: DTOFHit.h:30
int element
Definition: DFDCHit.h:25
float Fill(float x, float weight=1.0)
Definition: DHistogram.h:173
double npix_fadc
Definition: DTAGMHit.h:26
jerror_t Extract_DPSTruthHit(hddm_s::HDDM *record, JFactory< DPSTruthHit > *factory, string tag)
bool primary
Definition: DPSTruthHit.h:14
float t_fADC
Definition: DSCHit.h:22
float t_fADC
Definition: DTOFHit.h:26
float t_TDC
Definition: DSCHit.h:21
TH1D * py[NCHANNELS]
string StringToTMatrixFSym(string &str_vals, TMatrixFSym *mat, int Nrows, int Ncols)
jerror_t Extract_DBCALSiPMSpectrum(hddm_s::HDDM *record, JFactory< DBCALSiPMSpectrum > *factory, string tag)
float y
Definition: DFCALHit.h:26
void setPID(Particle_t locPID)
jerror_t Extract_DTAGMHit(hddm_s::HDDM *record, JFactory< DTAGMHit > *factory, string tag)
float q
Definition: DCDCHit.h:20
static double ParticleMass(Particle_t p)
jerror_t Extract_DPSCHit(hddm_s::HDDM *record, JFactory< DPSCHit > *factory, string tag)
double sqrt(double)
float E
Definition: DCCALHit.h:26
virtual ~DEventSourceHDDM()
static thread_local shared_ptr< DResourcePool< TMatrixFSym > > dResourcePool_TMatrixFSym
jerror_t Extract_DFCALTruthShower(hddm_s::HDDM *record, JFactory< DFCALTruthShower > *factory, string tag)
double time_tdc
Definition: DPSCHit.h:24
float d
Definition: DCDCHit.h:23
int itrack
Definition: DFDCHit.h:47
jerror_t Extract_DBCALSiPMHit(hddm_s::HDDM *record, JFactory< DBCALSiPMHit > *factory, string tag)
uint32_t pedestal
pedestal info used by FPGA (if any)
Definition: DBCALDigiHit.h:32
int primary
primary track=1 not primary track=0
Definition: DMCTrackHit.h:23
float dEdx
Definition: DSCTruthHit.h:18
jerror_t GetFDCTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
jerror_t Extract_DTPOLTruthHit(hddm_s::HDDM *record, JFactory< DTPOLTruthHit > *factory, string tag)
float pz
Definition: DTOFTruth.h:24
float t
Definition: DTOFTruth.h:25
float t
Definition: DFDCHit.h:32
jerror_t Extract_DRFTime(hddm_s::HDDM *record, JFactory< DRFTime > *factory, JEventLoop *locEventLoop)
int sss[NCHANNELS][NOXbins][NOYbins]
Definition: SC_PTC.C:75
bool has_fADC
Definition: DTAGHHit.h:26
float t
Definition: DCereHit.h:22
int bar
Definition: DTOFHit.h:22
int itrack
Definition: DTOFTruth.h:21
bool MCTrackHitSort_C(DMCTrackHit *const &thit1, DMCTrackHit *const &thit2)
double pulse_peak
Definition: DPSHit.h:24
jerror_t Extract_DFMWPCTruthHit(hddm_s::HDDM *record, JFactory< DFMWPCTruthHit > *factory, string tag)
float dE
Definition: DSCHit.h:19
int itrack
MC track index.
Definition: DMCTrackHit.h:22
double time_fadc
Definition: DTAGHHit.h:24
jerror_t GetCDCTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
std::ifstream * ifs
jerror_t GetEvent(JEvent &event)
bool has_fADC
Definition: DSCHit.h:24
int bar
Definition: DTOFHitMC.h:22
uint32_t pulse_time
identified pulse time as returned by FPGA algorithm
Definition: DBCALDigiHit.h:31
int cellId(int module, int layer, int sector) const
these functions are about encoding/decoding module/layer/sector info in a cellId
float t
Definition: DSCHit.h:20
double t
Definition: DPSHit.h:22
File: DTOFHit.h Created: Tue Jan 18 16:15:26 EST 2011 Creator: B. Zihlmann Purpose: Container class t...
Definition: DTOFHit.h:16
int module
Definition: DFDCHit.h:24
int column
Definition: DTAGMHit.h:21
float x
Definition: DFCALHit.h:25
float E
Definition: DFCALHit.h:27
jerror_t Extract_DCDCHit(JEventLoop *locEventLoop, hddm_s::HDDM *record, JFactory< DCDCHit > *factory, string tag)
DPSGeometry::Arm arm
Definition: DPSHit.h:19
double t
Definition: DTPOLHit.h:18
jerror_t Extract_DPSHit(hddm_s::HDDM *record, JFactory< DPSHit > *factory, string tag)
static float getWireR(const DFDCHit *h)
DFDCGeometry::getWireR(): Return X coordinate of a wire.
Definition: DFDCGeometry.h:75
void setPosition(const DVector3 &aPosition)
int column
Definition: DPSHit.h:20
bool isBlockActive(int row, int column) const
int type
GEANT particle ID.
Definition: DMCThrown.h:20
float t
Definition: DFCALHit.h:28
int column
Definition: DFCALHit.h:24
int track
Track number.
Definition: DMCTrackHit.h:21
bool isBlockActive(int row, int column) const
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
jerror_t Extract_DDIRCPmtHit(hddm_s::HDDM *record, JFactory< DDIRCPmtHit > *factory, string tag, JEventLoop *eventLoop)
Particle_t IDTrack(float locCharge, float locMass) const
jerror_t Extract_DSCHit(hddm_s::HDDM *record, JFactory< DSCHit > *factory, string tag)
int straw
Definition: DCDCHit.h:19
int plane
Definition: DTOFHitMC.h:21
float intOverPeak
Definition: DFCALHit.h:29
int itrack
Definition: DTOFHitMC.h:25
int QF
Definition: DCDCHit.h:24
jerror_t Extract_DCereHit(hddm_s::HDDM *record, JFactory< DCereHit > *factory, string tag)
jerror_t Extract_DBCALIncidentParticle(hddm_s::HDDM *record, JFactory< DBCALIncidentParticle > *factory, string tag)
int myid
id of this particle from original generator
Definition: DMCThrown.h:22
float dE
Definition: DTOFHit.h:24
int ptype
Definition: DTOFTruth.h:27
jerror_t Extract_DBCALDigiHit(hddm_s::HDDM *record, JFactory< DBCALDigiHit > *factory, string tag)
int gPlane
Definition: DFDCHit.h:27
float y
Definition: DCCALHit.h:25
float pe
Definition: DCereHit.h:21
jerror_t Extract_DFCALHit(hddm_s::HDDM *record, JFactory< DFCALHit > *factory, string tag, JEventLoop *eventLoop)
double time_fadc
Definition: DPSCHit.h:25
bool operator()(DMCTrackHit *const &thit1, DMCTrackHit *const &thit2) const
float r
Definition: DFDCHit.h:33
double mass(void) const
class DFDCHit: definition for a basic FDC hit data type.
Definition: DFDCHit.h:20
float z
Definition: DTOFHitMC.h:29
double getElow(int arm, int column) const
Definition: DPSGeometry.cc:40
jerror_t GetSCTruthHits(hddm_s::HDDM *record, vector< DMCTrackHit * > &data)
Particle_t
Definition: particleType.h:12