Bug Summary

File:libraries/HDDM/DEventSourceHDDM.cc
Location:line 1202, column 13
Description:Value stored to 'mass' is never read

Annotated Source Code

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