Bug Summary

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

Annotated Source Code

1//
2// Author: Richard Jones June 29, 2012
3//
4//
5// DEventSourceREST methods
6//
7
8#include <iostream>
9#include <iomanip>
10#include <fstream>
11#include <climits>
12
13#include <JANA/JFactory_base.h>
14#include <JANA/JEventLoop.h>
15#include <JANA/JEvent.h>
16#include <DANA/DStatusBits.h>
17
18#include <DVector2.h>
19#include <DEventSourceREST.h>
20
21//----------------
22// Constructor
23//----------------
24DEventSourceREST::DEventSourceREST(const char* source_name)
25 : JEventSource(source_name)
26{
27 /// Constructor for DEventSourceREST object
28 ifs = new ifstream(source_name);
29 ifs->get();
30 ifs->unget();
31 if (ifs->rdbuf()->in_avail() > 30) {
32 class nonstd_streambuf: public std::streambuf {
33 public: char *pub_gptr() {return gptr();}
34 };
35 void *buf = (void*)ifs->rdbuf();
36 std::stringstream sbuf(((nonstd_streambuf*)buf)->pub_gptr());
37 std::string head;
38 std::getline(sbuf, head);
39 std::string expected = " class=\"r\" ";
40 if (head.find(expected) == head.npos) {
41 std::string msg("Unexpected header found in input REST stream: ");
42 throw std::runtime_error(msg + head + source_name);
43 }
44 }
45
46 fin = new hddm_r::istream(*ifs);
47
48 PRUNE_DUPLICATE_TRACKS = true;
49 gPARMS->SetDefaultParameter("REST:PRUNE_DUPLICATE_TRACKS", PRUNE_DUPLICATE_TRACKS,
50 "Turn on/off cleaning up multiple tracks with the same hypothesis from the same candidate. Set to \"0\" to turn off (it's on by default)");
51
52 RECO_DIRC_CALC_LUT = true;
53 gPARMS->SetDefaultParameter("REST:DIRC_CALC_LUT", RECO_DIRC_CALC_LUT, "Turn on/off DIRC LUT reconstruction");
54
55 dDIRCMaxChannels = 108*64;
56
57 // any other initialization which needs to happen
58 dBCALShowerFactory = nullptr;
59 dFCALShowerFactory = nullptr;
60
61 USE_CCDB_BCAL_COVARIANCE = false;
62 gPARMS->SetDefaultParameter("REST:USE_CCDB_BCAL_COVARIANCE", USE_CCDB_BCAL_COVARIANCE,
63 "Load REST BCAL Shower covariance matrices from CCDB instead of the file.");
64 USE_CCDB_FCAL_COVARIANCE = false;
65 gPARMS->SetDefaultParameter("REST:USE_CCDB_FCAL_COVARIANCE", USE_CCDB_FCAL_COVARIANCE,
66 "Load REST BFAL Shower covariance matrices from CCDB instead of the file.");
67
68 gPARMS->SetDefaultParameter("REST:JANACALIBCONTEXT", REST_JANA_CALIB_CONTEXT);
69 calib_generator = new JCalibrationGeneratorCCDB; // keep this around in case we need to use it
70}
71
72//----------------
73// Destructor
74//----------------
75DEventSourceREST::~DEventSourceREST()
76{
77 if (fin) {
78 delete fin;
79 }
80 if (ifs) {
81 delete ifs;
82 }
83
84 for(auto &entry : dJCalib_olds)
85 delete entry.second;
86 for(auto &entry : dTAGHGeoms)
87 delete entry.second;
88 for(auto &entry : dTAGMGeoms)
89 delete entry.second;
90}
91
92//----------------
93// GetEvent
94//----------------
95jerror_t DEventSourceREST::GetEvent(JEvent &event)
96{
97 /// Implementation of JEventSource virtual function
98
99 if (!fin) {
100 return EVENT_SOURCE_NOT_OPEN;
101 }
102
103 // Each open hddm file takes up about 1M of memory so it's
104 // worthwhile to close it as soon as we can.
105 if (ifs->eof()) {
106 delete fin;
107 fin = NULL__null;
108 delete ifs;
109 ifs = NULL__null;
110
111 return NO_MORE_EVENTS_IN_SOURCE;
112 }
113
114#if HDDM_SETPOSITION_EXAMPLE
115 static std::ifstream fevlist("events.list");
116 static int events_to_go = 0;
117 if (events_to_go-- == 0 && fevlist.good()) {
118 uint64_t start;
119 uint32_t offset, status;
120 fevlist >> start >> offset >> status >> events_to_go;
121 if (fevlist.good())
122 fin->setPosition(hddm_r::streamposition(start, offset, status));
123 }
124#endif
125
126#if HDDM_GETPOSITION_EXAMPLE
127 hddm_r::streamposition pos(fin->getPosition());
128 // Later on below, if this event passes all of your selection cuts
129 // then you might want to write the event position to output, as in
130 // std::cout << "interesting event found at "
131 // << pos.start << "," << pos.offset << "," << pos.status
132 // << std::endl;
133#endif
134
135 hddm_r::HDDM *record = new hddm_r::HDDM();
136 try{
137 while (record->getReconstructedPhysicsEvents().size() == 0) {
138 if (! (*fin >> *record)) {
139 delete fin;
140 fin = NULL__null;
141 delete ifs;
142 ifs = NULL__null;
143 return NO_MORE_EVENTS_IN_SOURCE;
144 }
145 }
146 }catch(std::runtime_error &e){
147 cerr << "Exception caught while trying to read REST file!" << endl;
148 cerr << e.what() << endl;
149 _DBG__std::cerr<<"libraries/HDDM/DEventSourceREST.cc"<<
":"<<149<<std::endl;
150 return NO_MORE_EVENTS_IN_SOURCE;
151 }
152
153 // Copy the reference info into the JEvent object
154 while (true) {
155 hddm_r::ReconstructedPhysicsEvent &re
156 = record->getReconstructedPhysicsEvent();
157 int runno = re.getRunNo();
158 int eventno = re.getEventNo();
159 if (runno == 0 && eventno == 0) {
160 // found a comment record, print comment strings and continue
161 const hddm_r::CommentList &comments = re.getComments();
162 hddm_r::CommentList::iterator iter;
163 for (iter = comments.begin(); iter != comments.end(); ++iter) {
164 std::cout << " | " << iter->getText() << std::endl;
165 }
166
167 //set version string
168 const hddm_r::DataVersionStringList& locVersionStrings = re.getDataVersionStrings();
169 hddm_r::DataVersionStringList::iterator Versioniter;
170 for (Versioniter = locVersionStrings.begin(); Versioniter != locVersionStrings.end(); ++Versioniter) {
171 string HDDM_DATA_VERSION_STRING = Versioniter->getText();
172 if(gPARMS->Exists("REST:DATAVERSIONSTRING"))
173 gPARMS->SetParameter("REST:DATAVERSIONSTRING", HDDM_DATA_VERSION_STRING);
174 else
175 gPARMS->SetDefaultParameter("REST:DATAVERSIONSTRING", HDDM_DATA_VERSION_STRING);
176 break;
177 }
178
179 // set REST calib context - use this to load calibration constants that were used
180 // to create the REST files, if needed, but let this be overridden by command-line options
181 if( REST_JANA_CALIB_CONTEXT == "" ) {
182 const hddm_r::CcdbContextList& locContextStrings = re.getCcdbContexts();
183 hddm_r::CcdbContextList::iterator Contextiter;
184 for (Contextiter = locContextStrings.begin(); Contextiter != locContextStrings.end(); ++Contextiter) {
185 REST_JANA_CALIB_CONTEXT = Contextiter->getText();
186 jout << " REST file next CCDB context = " << REST_JANA_CALIB_CONTEXT << endl; // DEBUG?
187 }
188 }
189
190 record->clear();
191 while (record->getReconstructedPhysicsEvents().size() == 0) {
192 if (! (*fin >> *record)) {
193 delete fin;
194 fin = NULL__null;
195 delete ifs;
196 ifs = NULL__null;
197 return NO_MORE_EVENTS_IN_SOURCE;
198 }
199 }
200
201 continue;
202 }
203 event.SetEventNumber(re.getEventNo());
204 event.SetRunNumber(re.getRunNo());
205 event.SetJEventSource(this);
206 event.SetRef(record);
207 event.SetStatusBit(kSTATUS_REST);
208 event.SetStatusBit(kSTATUS_FROM_FILE);
209 event.SetStatusBit(kSTATUS_PHYSICS_EVENT);
210
211 ++Nevents_read;
212 break;
213 }
214
215 return NOERROR;
216}
217
218//----------------
219// FreeEvent
220//----------------
221void DEventSourceREST::FreeEvent(JEvent &event)
222{
223 hddm_r::HDDM *record = (hddm_r::HDDM*)event.GetRef();
224 delete record;
225}
226
227//----------------
228// GetObjects
229//----------------
230jerror_t DEventSourceREST::GetObjects(JEvent &event, JFactory_base *factory)
231{
232 /// This gets called through the virtual method of the
233 /// JEventSource base class. It creates the objects of the type
234 /// on which factory is based. It uses the HDDM_Element object
235 /// kept in the ref field of the JEvent object passed.
236
237 // We must have a factory to hold the data
238 if (!factory) {
239 throw RESOURCE_UNAVAILABLE;
240 }
241
242 // The ref field of the JEvent is just the HDDM record pointer.
243 hddm_r::HDDM *record = (hddm_r::HDDM*)event.GetRef();
244 if (!record) {
245 throw RESOURCE_UNAVAILABLE;
246 }
247
248 JEventLoop* locEventLoop = event.GetJEventLoop();
249 string dataClassName = factory->GetDataClassName();
250
251 DApplication* dapp = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
252 JCalibration *jcalib = dapp->GetJCalibration(event.GetRunNumber());
253
254
255 //Get target center
256 //multiple reader threads can access this object: need lock
257 bool locNewRunNumber = false;
258 unsigned int locRunNumber = event.GetRunNumber();
259 LockRead();
260 {
261 locNewRunNumber = (dTargetCenterZMap.find(locRunNumber) == dTargetCenterZMap.end());
262 }
263 UnlockRead();
264 if(locNewRunNumber)
265 {
266 DGeometry* locGeometry = dapp->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
267 double locTargetCenterZ = 0.0;
268 locGeometry->GetTargetZ(locTargetCenterZ);
269
270 map<string, double> beam_vals;
271 if (locEventLoop->GetCalib("PHOTON_BEAM/beam_spot",beam_vals))
272 throw JException("Could not load CCDB table: PHOTON_BEAM/beam_spot");
273
274 vector<double> locBeamPeriodVector;
275 if(locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector))
276 throw JException("Could not load CCDB table: PHOTON_BEAM/RF/beam_period");
277 double locBeamBunchPeriod = locBeamPeriodVector[0];
278
279 vector< vector <int> > locDIRCChannelStatus;
280 vector<int> new_dirc_status(dDIRCMaxChannels);
281 locDIRCChannelStatus.push_back(new_dirc_status);
282 locDIRCChannelStatus.push_back(new_dirc_status);
283 if(RECO_DIRC_CALC_LUT) { // get DIRC channel status from DB
284 if (locEventLoop->GetCalib("/DIRC/North/channel_status", locDIRCChannelStatus[0]))
285 jout << "Error loading /DIRC/North/channel_status !" << endl;
286 if (locEventLoop->GetCalib("/DIRC/South/channel_status", locDIRCChannelStatus[1]))
287 jout << "Error loading /DIRC/South/channel_status !" << endl;
288 }
289
290
291 LockRead();
292 {
293 dTargetCenterZMap[locRunNumber] = locTargetCenterZ;
294 dBeamBunchPeriodMap[locRunNumber] = locBeamBunchPeriod;
295 dDIRCChannelStatusMap[locRunNumber] = locDIRCChannelStatus;
296
297 dBeamCenterMap[locRunNumber].Set(beam_vals["x"],
298 beam_vals["y"]);
299 dBeamDirMap[locRunNumber].Set(beam_vals["dxdz"],
300 beam_vals["dydz"]);
301 dBeamZ0Map[locRunNumber]=beam_vals["z"];
302
303 jout << "Run " << locRunNumber << " beam spot:"
304 << " x=" << dBeamCenterMap[locRunNumber].X()
305 << " y=" << dBeamCenterMap[locRunNumber].Y()
306 << " z=" << dBeamZ0Map[locRunNumber]
307 << " dx/dz=" << dBeamDirMap[locRunNumber].X()
308 << " dy/dz=" << dBeamDirMap[locRunNumber].Y()
309 << endl;
310
311 // tagger related configs for reverse mapping tagger energy to counter number
312 if( REST_JANA_CALIB_CONTEXT != "" ) {
313 JCalibration *jcalib_old = calib_generator->MakeJCalibration(jcalib->GetURL(), locRunNumber, REST_JANA_CALIB_CONTEXT );
314 dTAGHGeoms[locRunNumber] = new DTAGHGeometry(jcalib_old, locRunNumber);
315 dTAGMGeoms[locRunNumber] = new DTAGMGeometry(jcalib_old, locRunNumber);
316 dJCalib_olds[locRunNumber] = jcalib_old;
317 }
318 }
319 UnlockRead();
320
321 // do multiple things to limit the number of locks
322 // make sure that we have a handle to the FCAL shower factory
323 if(USE_CCDB_FCAL_COVARIANCE) {
324 if(dFCALShowerFactory==nullptr) {
325 dFCALShowerFactory = static_cast<DFCALShower_factory*>(locEventLoop->GetFactory("DFCALShower"));
326 if(dFCALShowerFactory==nullptr)
327 throw JException("Couldn't find DFCALShower_factory???");
328 }
329 dFCALShowerFactory->LoadCovarianceLookupTables(locEventLoop);
330 }
331
332 // same with BCAL
333 if(USE_CCDB_BCAL_COVARIANCE) {
334 if(dBCALShowerFactory==nullptr) {
335 dBCALShowerFactory = static_cast<DBCALShower_factory_IU*>(locEventLoop->GetFactory("DBCALShower", "IU"));
336 if(dBCALShowerFactory==nullptr)
337 throw JException("Couldn't find DBCALShower_factory???");
338 }
339 dBCALShowerFactory->LoadCovarianceLookupTables(locEventLoop);
340 }
341
342 }
343
344 if (dataClassName =="DMCReaction") {
345 return Extract_DMCReaction(record,
346 dynamic_cast<JFactory<DMCReaction>*>(factory), locEventLoop);
347 }
348 if (dataClassName =="DRFTime") {
349 return Extract_DRFTime(record,
350 dynamic_cast<JFactory<DRFTime>*>(factory), locEventLoop);
351 }
352 if (dataClassName =="DBeamPhoton") {
353 return Extract_DBeamPhoton(record,
354 dynamic_cast<JFactory<DBeamPhoton>*>(factory),
355 locEventLoop);
356 }
357 if (dataClassName =="DMCThrown") {
358 return Extract_DMCThrown(record,
359 dynamic_cast<JFactory<DMCThrown>*>(factory));
360 }
361 if (dataClassName =="DTOFPoint") {
362 return Extract_DTOFPoint(record,
363 dynamic_cast<JFactory<DTOFPoint>*>(factory));
364 }
365 if (dataClassName =="DCTOFPoint") {
366 return Extract_DCTOFPoint(record,
367 dynamic_cast<JFactory<DCTOFPoint>*>(factory));
368 }
369 if (dataClassName =="DSCHit") {
370 return Extract_DSCHit(record,
371 dynamic_cast<JFactory<DSCHit>*>(factory));
372 }
373 if (dataClassName =="DFCALShower") {
374 return Extract_DFCALShower(record,
375 dynamic_cast<JFactory<DFCALShower>*>(factory));
376 }
377 if (dataClassName =="DBCALShower") {
378 return Extract_DBCALShower(record,
379 dynamic_cast<JFactory<DBCALShower>*>(factory));
380 }
381 if (dataClassName =="DCCALShower") {
382 return Extract_DCCALShower(record,
383 dynamic_cast<JFactory<DCCALShower>*>(factory));
384 }
385 if (dataClassName =="DTrackTimeBased") {
386 return Extract_DTrackTimeBased(record,
387 dynamic_cast<JFactory<DTrackTimeBased>*>(factory), locEventLoop);
388 }
389 if (dataClassName =="DTrigger") {
390 return Extract_DTrigger(record,
391 dynamic_cast<JFactory<DTrigger>*>(factory));
392 }
393 if (dataClassName =="DDIRCPmtHit") {
394 return Extract_DDIRCPmtHit(record,
395 dynamic_cast<JFactory<DDIRCPmtHit>*>(factory), locEventLoop);
396 }
397 if (dataClassName =="DFMWPCHit") {
398 return Extract_DFMWPCHit(record,
399 dynamic_cast<JFactory<DFMWPCHit>*>(factory), locEventLoop);
400 }
401 if (dataClassName =="DDetectorMatches") {
402 return Extract_DDetectorMatches(locEventLoop, record,
403 dynamic_cast<JFactory<DDetectorMatches>*>(factory));
404 }
405 if (dataClassName =="DEventHitStatistics") {
406 return Extract_DEventHitStatistics(record,
407 dynamic_cast<JFactory<DEventHitStatistics>*>(factory));
408 }
409
410 return OBJECT_NOT_AVAILABLE;
411}
412
413//------------------
414// Extract_DMCReaction
415//------------------
416jerror_t DEventSourceREST::Extract_DMCReaction(hddm_r::HDDM *record,
417 JFactory<DMCReaction> *factory, JEventLoop* locEventLoop)
418{
419 /// Copies the data from the Reaction hddm class. This is called
420 /// from JEventSourceREST::GetObjects. If factory is NULL, this
421 /// returns OBJECT_NOT_AVAILABLE immediately.
422
423 if (factory==NULL__null) {
424 return OBJECT_NOT_AVAILABLE;
425 }
426 std::string tag = (factory->Tag())? factory->Tag() : "";
427
428 double locTargetCenterZ = 0.0;
429 int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
430 LockRead();
431 {
432 locTargetCenterZ = dTargetCenterZMap[locRunNumber];
433 }
434 UnlockRead();
435 DVector3 locPosition(0.0, 0.0, locTargetCenterZ);
436
437 vector<DMCReaction*> dmcreactions;
438
439 // loop over reaction records
440 const hddm_r::ReactionList &reactions = record->getReactions();
441 hddm_r::ReactionList::iterator iter;
442 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
443 if (iter->getJtag() != tag) {
444 continue;
445 }
446 DMCReaction *mcreaction = new DMCReaction;
447 dmcreactions.push_back(mcreaction);
448 mcreaction->type = iter->getType();
449 mcreaction->weight = iter->getWeight();
450 double Ebeam = iter->getEbeam();
451
452 hddm_r::Origin &origin = iter->getVertex().getOrigin();
453 double torig = origin.getT();
454 double zorig = origin.getVz();
455
456 mcreaction->beam.setPosition(locPosition);
457 mcreaction->beam.setMomentum(DVector3(0.0, 0.0, Ebeam));
458 mcreaction->beam.setTime(torig - (zorig - locTargetCenterZ)/29.9792458);
459 mcreaction->beam.setPID(Gamma);
460
461 mcreaction->target.setPosition(locPosition);
462 Particle_t ttype = iter->getTargetType();
463 mcreaction->target.setPID((Particle_t)ttype);
464 mcreaction->target.setTime(torig - (zorig - locTargetCenterZ)/29.9792458);
465 }
466
467 // Copy into factories
468 factory->CopyTo(dmcreactions);
469
470 return NOERROR;
471}
472
473//------------------
474// Extract_DRFTime
475//------------------
476jerror_t DEventSourceREST::Extract_DRFTime(hddm_r::HDDM *record,
477 JFactory<DRFTime> *factory, JEventLoop* locEventLoop)
478{
479 if (factory==NULL__null)
480 return OBJECT_NOT_AVAILABLE;
481 string tag = (factory->Tag())? factory->Tag() : "";
482
483 vector<DRFTime*> locRFTimes;
484
485 // loop over RF-time records
486 const hddm_r::RFtimeList &rftimes = record->getRFtimes();
487 hddm_r::RFtimeList::iterator iter;
488 for (iter = rftimes.begin(); iter != rftimes.end(); ++iter)
489 {
490 if (iter->getJtag() != tag)
491 continue;
492 DRFTime *locRFTime = new DRFTime;
493 locRFTime->dTime = iter->getTsync();
494 locRFTime->dTimeVariance = 0.0; //SET ME!!
495 locRFTimes.push_back(locRFTime);
496 }
497
498 if(!locRFTimes.empty())
499 {
500 //found in the file, copy into factory and return
501 factory->CopyTo(locRFTimes);
502 return NOERROR;
503 }
504
505 //Not found in the file, so either:
506 //Experimental data & it's missing: bail
507 //MC data: generate it
508
509 vector<const DBeamPhoton*> locMCGENPhotons;
510 locEventLoop->Get(locMCGENPhotons, "MCGEN");
511 if(locMCGENPhotons.empty())
512 return OBJECT_NOT_AVAILABLE; //Experimental data & it's missing: bail
513
514 //Is MC data. Either:
515 //No tag: return photon_time propagated by +/- n*locBeamBunchPeriod to get close to 0.0
516 //TRUTH tag: get exact t from DBeamPhoton tag MCGEN
517
518 if(tag == "TRUTH")
519 {
520 DRFTime *locRFTime = new DRFTime;
521 locRFTime->dTime = locMCGENPhotons[0]->time();
522 locRFTime->dTimeVariance = 0.0;
523 locRFTimes.push_back(locRFTime);
524 }
525 else
526 {
527 double locBeamBunchPeriod = 0.0;
528 int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
529 LockRead();
530 {
531 locBeamBunchPeriod = dBeamBunchPeriodMap[locRunNumber];
532 }
533 UnlockRead();
534
535 //start with true RF time, increment/decrement by multiples of locBeamBunchPeriod ns until closest to 0
536 double locTime = locMCGENPhotons[0]->time();
537 int locNumRFBuckets = int(locTime/locBeamBunchPeriod);
538 locTime -= double(locNumRFBuckets)*locBeamBunchPeriod;
539 while(locTime > 0.5*locBeamBunchPeriod)
540 locTime -= locBeamBunchPeriod;
541 while(locTime < -0.5*locBeamBunchPeriod)
542 locTime += locBeamBunchPeriod;
543 DRFTime *locRFTime = new DRFTime;
544 locRFTime->dTime = locTime;
545 locRFTime->dTimeVariance = 0.0;
546 locRFTimes.push_back(locRFTime);
547 }
548
549 // Copy into factories
550 factory->CopyTo(locRFTimes);
551
552 return NOERROR;
553}
554
555//------------------
556// Extract_DBeamPhoton
557//------------------
558jerror_t DEventSourceREST::Extract_DBeamPhoton(hddm_r::HDDM *record,
559 JFactory<DBeamPhoton> *factory,
560 JEventLoop *eventLoop)
561{
562 /// This is called from JEventSourceREST::GetObjects. If factory is NULL,
563 /// return OBJECT_NOT_AVAILABLE immediately. If factory tag="MCGEN" then
564 /// copy the beam photon data from the Reaction hddm class.
565
566 if (factory==NULL__null)
567 return OBJECT_NOT_AVAILABLE;
568 string tag = (factory->Tag())? factory->Tag() : "";
569
570 vector<DBeamPhoton*> dbeam_photons;
571
572 // extract the TAGH geometry
573 vector<const DTAGHGeometry*> taghGeomVect;
574 eventLoop->Get(taghGeomVect);
575 if (taghGeomVect.empty())
576 return OBJECT_NOT_AVAILABLE;
577 const DTAGHGeometry* taghGeom = taghGeomVect[0];
578
579 // extract the TAGM geometry
580 vector<const DTAGMGeometry*> tagmGeomVect;
581 eventLoop->Get(tagmGeomVect);
582 if (tagmGeomVect.empty())
583 return OBJECT_NOT_AVAILABLE;
584 const DTAGMGeometry* tagmGeom = tagmGeomVect[0];
585
586 if(tag == "MCGEN")
587 {
588 vector<const DMCReaction*> dmcreactions;
589 eventLoop->Get(dmcreactions);
590
591 for(size_t loc_i = 0; loc_i < dmcreactions.size(); ++loc_i)
592 {
593 DBeamPhoton *beamphoton = new DBeamPhoton;
594 *(DKinematicData*)beamphoton = dmcreactions[loc_i]->beam;
595 if(tagmGeom->E_to_column(beamphoton->energy(), beamphoton->dCounter))
596 beamphoton->dSystem = SYS_TAGM;
597 else if(taghGeom->E_to_counter(beamphoton->energy(), beamphoton->dCounter))
598 beamphoton->dSystem = SYS_TAGH;
599 else
600 beamphoton->dSystem = SYS_NULL;
601 dbeam_photons.push_back(beamphoton);
602 }
603
604 // Copy into factories
605 factory->CopyTo(dbeam_photons);
606
607 return NOERROR;
608 }
609
610 double locTargetCenterZ = 0.0;
611 int locRunNumber = eventLoop->GetJEvent().GetRunNumber();
612 LockRead();
613 {
614 locTargetCenterZ = dTargetCenterZMap[locRunNumber];
615 }
616 UnlockRead();
617
618 DVector3 pos(0.0, 0.0, locTargetCenterZ);
619
620 //now get the objects
621 const hddm_r::TagmBeamPhotonList &locTagmBeamPhotonList = record->getTagmBeamPhotons();
622 hddm_r::TagmBeamPhotonList::iterator locTAGMiter;
623 for(locTAGMiter = locTagmBeamPhotonList.begin(); locTAGMiter != locTagmBeamPhotonList.end(); ++locTAGMiter)
624 {
625 if (locTAGMiter->getJtag() != tag)
626 continue;
627
628 // load the counter number (if it exists) and set the energy based on the counter
629 unsigned int column = 0;
630 hddm_r::TagmChannelList &locTagmChannelList = locTAGMiter->getTagmChannels();
631 if (locTagmChannelList.size() > 0) {
632 // it's easy if the column is already set
633 column = locTagmChannelList().getColumn();
634 } else {
635 // if the TAGM column isn't saved in the REST file, then we do one of two things
636 // 1) if there's no special CCDB context associated with the file, we can just
637 // reverse engineer the counter, assuming the latest CCDB
638 // 2) If there is a special CCDB context specified, then use that instead
639 if (dJCalib_olds[locRunNumber] == nullptr) {
640 if (!tagmGeom->E_to_column(locTAGMiter->getE(), column)) {
641 column = 0;
642 }
643 } else {
644 if (!dTAGMGeoms[locRunNumber]->E_to_column(locTAGMiter->getE(), column)) {
645 column = 0;
646 }
647 }
648
649 if(column == 0)
650 std::cerr << "Error in DEventSourceREST - tagger microscope could not look up column for energy "
651 << locTAGMiter->getE() << std::endl;
652 }
653
654 // sometimes the simulation will set photons that miss the tagger counters to have a column of zero - skip these
655 if(column == 0) {
656 continue;
657 }
658
659 DBeamPhoton* gamma = new DBeamPhoton();
660
661 double Elo_tagm = tagmGeom->getElow(column);
662 double Ehi_tagm = tagmGeom->getEhigh(column);
663 double Ebeam_tagm = (Elo_tagm + Ehi_tagm)/2.;
664
665 // read the rest of the data from the REST file
666 DVector3 mom(0.0, 0.0, Ebeam_tagm);
667 gamma->setPID(Gamma);
668 gamma->setMomentum(mom);
669 gamma->setPosition(pos);
670 gamma->setTime(locTAGMiter->getT());
671 gamma->dSystem = SYS_TAGM;
672 gamma->dCounter = column;
673
674 auto locCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
675 locCovarianceMatrix->ResizeTo(7, 7);
676 locCovarianceMatrix->Zero();
677 gamma->setErrorMatrix(locCovarianceMatrix);
678
679 dbeam_photons.push_back(gamma);
680 }
681
682 const hddm_r::TaghBeamPhotonList &locTaghBeamPhotonList = record->getTaghBeamPhotons();
683 hddm_r::TaghBeamPhotonList::iterator locTAGHiter;
684 for(locTAGHiter = locTaghBeamPhotonList.begin(); locTAGHiter != locTaghBeamPhotonList.end(); ++locTAGHiter)
685 {
686 if (locTAGHiter->getJtag() != tag)
687 continue;
688
689 // load the counter number (if it exists) and set the energy based on the counter
690 unsigned int counter = 0;
691 hddm_r::TaghChannelList &locTaghChannelList = locTAGHiter->getTaghChannels();
692 if (locTaghChannelList.size() > 0) {
693 // it's easy if the column is already set
694 counter = locTaghChannelList().getCounter();
695 } else {
696 // if the TAGH column isn't saved in the REST file, then we do one of two things
697 // 1) if there's no special CCDB context associated with the file, we can just
698 // reverse engineer the counter, assuming the latest CCDB
699 // 2) If there is a special CCDB context specified, then use that instead
700 if (dJCalib_olds[locRunNumber] == nullptr) {
701 if (!taghGeom->E_to_counter(locTAGHiter->getE(), counter)) {
702 counter = 0;
703 }
704 } else {
705 if (!dTAGHGeoms[locRunNumber]->E_to_counter(locTAGHiter->getE(), counter)) {
706 counter = 0;
707 }
708 }
709
710 if(counter == 0)
711 std::cerr << "Error in DEventSourceREST - tagger hodoscope could not look up counter for energy "
712 << locTAGHiter->getE() << std::endl;
713 }
714
715 // sometimes the simulation will set photons that miss the tagger counters to have a column of zero - skip these
716 if(counter == 0) {
717 continue;
718 }
719
720 DBeamPhoton* gamma = new DBeamPhoton();
721
722 double Elo_tagh = taghGeom->getElow(counter);
723 double Ehi_tagh = taghGeom->getEhigh(counter);
724 double Ebeam_tagh = (Elo_tagh + Ehi_tagh)/2.;
725
726 DVector3 mom(0.0, 0.0, Ebeam_tagh);
727 gamma->setPID(Gamma);
728 gamma->setMomentum(mom);
729 gamma->setPosition(pos);
730 gamma->setTime(locTAGHiter->getT());
731 gamma->dSystem = SYS_TAGH;
732 gamma->dCounter = counter;
733
734 auto locCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
735 locCovarianceMatrix->ResizeTo(7, 7);
736 locCovarianceMatrix->Zero();
737 gamma->setErrorMatrix(locCovarianceMatrix);
738
739 dbeam_photons.push_back(gamma);
740 }
741
742 if((tag == "TAGGEDMCGEN") && dbeam_photons.empty())
743 return OBJECT_NOT_AVAILABLE; //EITHER: didn't hit a tagger counter //OR: old MC data (pre-saving TAGGEDMCGEN): try using TAGGEDMCGEN factory
744
745 // Copy into factories
746 factory->CopyTo(dbeam_photons);
747
748 return NOERROR;
749}
750
751//------------------
752// Extract_DMCThrown
753//------------------
754jerror_t DEventSourceREST::Extract_DMCThrown(hddm_r::HDDM *record,
755 JFactory<DMCThrown> *factory)
756{
757 /// Copies the data from the hddm vertex records. This is called
758 /// from JEventSourceREST::GetObjects. If factory is NULL, this
759 /// returns OBJECT_NOT_AVAILABLE immediately.
760
761 if (factory==NULL__null) {
762 return OBJECT_NOT_AVAILABLE;
763 }
764 string tag = (factory->Tag())? factory->Tag() : "";
765
766 vector<DMCThrown*> data;
767
768 // loop over vertex records
769 hddm_r::VertexList vertices = record->getVertices();
770 hddm_r::VertexList::iterator iter;
771 for (iter = vertices.begin(); iter != vertices.end(); ++iter) {
772 if (iter->getJtag() != tag) {
773 continue;
774 }
775 const hddm_r::Origin &orig = iter->getOrigin();
776 double vx = orig.getVx();
777 double vy = orig.getVy();
778 double vz = orig.getVz();
779 double vt = orig.getT();
780 const hddm_r::ProductList &products = iter->getProducts();
781 hddm_r::ProductList::iterator piter;
782 for (piter = products.begin(); piter != products.end(); ++piter) {
783 double E = piter->getMomentum().getE();
784 double px = piter->getMomentum().getPx();
785 double py = piter->getMomentum().getPy();
786 double pz = piter->getMomentum().getPz();
787 double mass = sqrt(E*E - (px*px + py*py + pz*pz));
788 if (!isfinite(mass)) {
789 mass = 0.0;
Value stored to 'mass' is never read
790 }
791 DMCThrown *mcthrown = new DMCThrown;
792 int pdgtype = piter->getPdgtype();
793 Particle_t ptype = PDGtoPType(pdgtype);
794 mcthrown->type = ptype;
795 mcthrown->pdgtype = pdgtype;
796 mcthrown->myid = piter->getId();
797 mcthrown->parentid = piter->getParentId();
798 mcthrown->mech = 0;
799 mcthrown->setPID(ptype);
800 mcthrown->setMomentum(DVector3(px, py, pz));
801 mcthrown->setPosition(DVector3(vx, vy, vz));
802 mcthrown->setTime(vt);
803 data.push_back(mcthrown);
804 }
805 }
806
807 // Copy into factory
808 factory->CopyTo(data);
809
810 return NOERROR;
811}
812
813//------------------
814// Extract_DTOFPoint
815//------------------
816jerror_t DEventSourceREST::Extract_DTOFPoint(hddm_r::HDDM *record,
817 JFactory<DTOFPoint>* factory)
818{
819 /// Copies the data from the tofPoint hddm record. This is called
820 /// from JEventSourceREST::GetObjects. If factory is NULL, this
821 /// returns OBJECT_NOT_AVAILABLE immediately.
822
823 if (factory==NULL__null) {
824 return OBJECT_NOT_AVAILABLE;
825 }
826 string tag = (factory->Tag())? factory->Tag() : "";
827
828 vector<DTOFPoint*> data;
829
830 // loop over tofPoint records
831 const hddm_r::TofPointList &tofs = record->getTofPoints();
832 hddm_r::TofPointList::iterator iter;
833 for (iter = tofs.begin(); iter != tofs.end(); ++iter) {
834 if (iter->getJtag() != tag) {
835 continue;
836 }
837 DTOFPoint *tofpoint = new DTOFPoint();
838 tofpoint->pos = DVector3(iter->getX(),iter->getY(),iter->getZ());
839 tofpoint->t = iter->getT();
840 tofpoint->dE = iter->getDE();
841 tofpoint->tErr = iter->getTerr();
842
843 //Status
844 const hddm_r::TofStatusList& locTofStatusList = iter->getTofStatuses();
845 hddm_r::TofStatusList::iterator locStatusIterator = locTofStatusList.begin();
846 if(locStatusIterator == locTofStatusList.end())
847 {
848 tofpoint->dHorizontalBar = 0;
849 tofpoint->dVerticalBar = 0;
850 tofpoint->dHorizontalBarStatus = 3;
851 tofpoint->dVerticalBarStatus = 3;
852 }
853 else //should only be 1
854 {
855 for(; locStatusIterator != locTofStatusList.end(); ++locStatusIterator)
856 {
857 int locStatus = locStatusIterator->getStatus(); //horizontal_bar + 45*vertical_bar + 45*45*horizontal_status + 45*45*4*vertical_status
858 tofpoint->dVerticalBarStatus = locStatus/(45*45*4);
859 locStatus %= 45*45*4; //Assume compiler optimizes multiplication
860 tofpoint->dHorizontalBarStatus = locStatus/(45*45);
861 locStatus %= 45*45;
862 tofpoint->dVerticalBar = locStatus/45;
863 tofpoint->dHorizontalBar = locStatus % 45;
864 }
865 }
866 // Energy deposition
867 const hddm_r::TofEnergyDepositionList& locTofEnergyDepositionList = iter->getTofEnergyDepositions();
868 hddm_r::TofEnergyDepositionList::iterator locEnergyDepositionIterator = locTofEnergyDepositionList.begin();
869 if(locEnergyDepositionIterator == locTofEnergyDepositionList.end())
870 {
871 tofpoint->dE1 = 0.;
872 tofpoint->dE2 = 0.;
873 }
874 else //should only be 1
875 {
876 for(; locEnergyDepositionIterator != locTofEnergyDepositionList.end(); ++locEnergyDepositionIterator)
877 {
878 tofpoint->dE1 = locEnergyDepositionIterator->getDE1();
879 tofpoint->dE2 = locEnergyDepositionIterator->getDE2();
880 }
881 }
882
883 data.push_back(tofpoint);
884 }
885
886 // Copy into factory
887 factory->CopyTo(data);
888
889 return NOERROR;
890}
891
892//------------------
893// Extract_DCTOFPoint
894//------------------
895jerror_t DEventSourceREST::Extract_DCTOFPoint(hddm_r::HDDM *record,
896 JFactory<DCTOFPoint>* factory)
897{
898 /// Copies the data from the ctofPoint hddm record. This is called
899 /// from JEventSourceREST::GetObjects. If factory is NULL, this
900 /// returns OBJECT_NOT_AVAILABLE immediately.
901
902 if (factory==NULL__null) {
903 return OBJECT_NOT_AVAILABLE;
904 }
905 string tag = (factory->Tag())? factory->Tag() : "";
906
907 vector<DCTOFPoint*> data;
908
909 // loop over ctofPoint records
910 const hddm_r::CtofPointList &ctofs = record->getCtofPoints();
911 hddm_r::CtofPointList::iterator iter;
912 for (iter = ctofs.begin(); iter != ctofs.end(); ++iter) {
913 if (iter->getJtag() != tag) {
914 continue;
915 }
916 DCTOFPoint *ctofpoint = new DCTOFPoint();
917 ctofpoint->bar = iter->getBar();
918 ctofpoint->pos = DVector3(iter->getX(),iter->getY(),iter->getZ());
919 ctofpoint->t = iter->getT();
920 ctofpoint->dE = iter->getDE();
921
922 data.push_back(ctofpoint);
923 }
924
925 // Copy into factory
926 factory->CopyTo(data);
927
928 return NOERROR;
929}
930
931
932//------------------
933// Extract_DSCHit
934//------------------
935jerror_t DEventSourceREST::Extract_DSCHit(hddm_r::HDDM *record,
936 JFactory<DSCHit>* factory)
937{
938 /// Copies the data from the startHit hddm record. This is called
939 /// from JEventSourceREST::GetObjects. If factory is NULL, this
940 /// returns OBJECT_NOT_AVAILABLE immediately.
941
942 if (factory==NULL__null) {
943 return OBJECT_NOT_AVAILABLE;
944 }
945 string tag = (factory->Tag())? factory->Tag() : "";
946
947 vector<DSCHit*> data;
948
949 // loop over startHit records
950 const hddm_r::StartHitList &starts = record->getStartHits();
951 hddm_r::StartHitList::iterator iter;
952 for (iter = starts.begin(); iter != starts.end(); ++iter) {
953 if (iter->getJtag() != tag) {
954 continue;
955 }
956 DSCHit *start = new DSCHit();
957 start->sector = iter->getSector();
958 start->dE = iter->getDE();
959 start->t = iter->getT();
960 data.push_back(start);
961 }
962
963 // Copy into factory
964 factory->CopyTo(data);
965
966 return NOERROR;
967}
968
969//-----------------------
970// Extract_DTrigger
971//-----------------------
972jerror_t DEventSourceREST::Extract_DTrigger(hddm_r::HDDM *record, JFactory<DTrigger>* factory)
973{
974 /// Copies the data from the trigger hddm record. This is
975 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
976 /// returns OBJECT_NOT_AVAILABLE immediately.
977
978 if (factory==NULL__null)
979 return OBJECT_NOT_AVAILABLE;
980 string tag = (factory->Tag())? factory->Tag() : "";
981
982 vector<DTrigger*> data;
983
984 // loop over trigger records
985 const hddm_r::TriggerList &triggers = record->getTriggers();
986 hddm_r::TriggerList::iterator iter;
987 for (iter = triggers.begin(); iter != triggers.end(); ++iter)
988 {
989 if (iter->getJtag() != tag)
990 continue;
991
992 DTrigger *locTrigger = new DTrigger();
993 locTrigger->Set_L1TriggerBits(Convert_SignedIntToUnsigned(iter->getL1_trig_bits()));
994 locTrigger->Set_L1FrontPanelTriggerBits(Convert_SignedIntToUnsigned(iter->getL1_fp_trig_bits()));
995
996 const hddm_r::TriggerEnergySumsList& locTriggerEnergySumsList = iter->getTriggerEnergySumses();
997 hddm_r::TriggerEnergySumsList::iterator locTriggerEnergySumsIterator = locTriggerEnergySumsList.begin();
998 if(locTriggerEnergySumsIterator == locTriggerEnergySumsList.end()) {
999 locTrigger->Set_GTP_BCALEnergy(0);
1000 locTrigger->Set_GTP_FCALEnergy(0);
1001 } else { //should only be 1
1002 for(; locTriggerEnergySumsIterator != locTriggerEnergySumsList.end(); ++locTriggerEnergySumsIterator) {
1003 locTrigger->Set_GTP_BCALEnergy(locTriggerEnergySumsIterator->getBCALEnergySum());
1004 locTrigger->Set_GTP_FCALEnergy(locTriggerEnergySumsIterator->getFCALEnergySum());
1005 }
1006 }
1007
1008 data.push_back(locTrigger);
1009 }
1010
1011 // Copy into factory
1012 factory->CopyTo(data);
1013
1014 return NOERROR;
1015}
1016
1017//-----------------------
1018// Extract_DFCALShower
1019//-----------------------
1020jerror_t DEventSourceREST::Extract_DFCALShower(hddm_r::HDDM *record,
1021 JFactory<DFCALShower>* factory)
1022{
1023 /// Copies the data from the fcalShower hddm record. This is
1024 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1025 /// returns OBJECT_NOT_AVAILABLE immediately.
1026
1027 if (factory==NULL__null) {
1028 return OBJECT_NOT_AVAILABLE;
1029 }
1030 string tag = (factory->Tag())? factory->Tag() : "";
1031
1032 vector<DFCALShower*> data;
1033
1034 // loop over fcal shower records
1035 const hddm_r::FcalShowerList &showers =
1036 record->getFcalShowers();
1037 hddm_r::FcalShowerList::iterator iter;
1038 for (iter = showers.begin(); iter != showers.end(); ++iter) {
1039 if (iter->getJtag() != tag)
1040 continue;
1041
1042 DFCALShower *shower = new DFCALShower();
1043 shower->setPosition(DVector3(iter->getX(),iter->getY(),iter->getZ()));
1044 shower->setEnergy(iter->getE());
1045 shower->setTime(iter->getT());
1046
1047 if(USE_CCDB_FCAL_COVARIANCE) {
1048 dFCALShowerFactory->FillCovarianceMatrix(shower);
1049 } else {
1050 TMatrixFSym covariance(5);
1051 covariance(0,0) = iter->getEerr()*iter->getEerr();
1052 covariance(1,1) = iter->getXerr()*iter->getXerr();
1053 covariance(2,2) = iter->getYerr()*iter->getYerr();
1054 covariance(3,3) = iter->getZerr()*iter->getZerr();
1055 covariance(4,4) = iter->getTerr()*iter->getTerr();
1056 covariance(1,2) = covariance(2,1) = iter->getXycorr()*iter->getXerr()*iter->getYerr();
1057 covariance(1,3) = covariance(3,1) = iter->getXzcorr()*iter->getXerr()*iter->getZerr();
1058 covariance(2,3) = covariance(3,2) = iter->getYzcorr()*iter->getYerr()*iter->getZerr();
1059 covariance(0,3) = covariance(3,0) = iter->getEzcorr()*iter->getEerr()*iter->getZerr();
1060 covariance(3,4) = covariance(4,3) = iter->getTzcorr()*iter->getTerr()*iter->getZerr();
1061
1062 // further correlations (an extension of REST format, so code is different.)
1063 const hddm_r::FcalCorrelationsList& locFcalCorrelationsList = iter->getFcalCorrelationses();
1064 hddm_r::FcalCorrelationsList::iterator locFcalCorrelationsIterator = locFcalCorrelationsList.begin();
1065 if(locFcalCorrelationsIterator != locFcalCorrelationsList.end()) {
1066 covariance(0,4) = covariance(4,0) = locFcalCorrelationsIterator->getEtcorr()*iter->getEerr()*iter->getTerr();
1067 covariance(0,1) = covariance(1,0) = locFcalCorrelationsIterator->getExcorr()*iter->getEerr()*iter->getXerr();
1068 covariance(0,2) = covariance(2,0) = locFcalCorrelationsIterator->getEycorr()*iter->getEerr()*iter->getYerr();
1069 covariance(1,4) = covariance(4,1) = locFcalCorrelationsIterator->getTxcorr()*iter->getTerr()*iter->getXerr();
1070 covariance(2,4) = covariance(4,2) = locFcalCorrelationsIterator->getTycorr()*iter->getTerr()*iter->getYerr();
1071 }
1072 shower->ExyztCovariance = covariance;
1073 }
1074
1075 // MVA classifier output - this information is being calculated in DNeutralShower now!
1076 //const hddm_r::FcalShowerClassificationList& locFcalShowerClassificationList = iter->getFcalShowerClassifications();
1077 //hddm_r::FcalShowerClassificationList::iterator locFcalShowerClassificationIterator = locFcalShowerClassificationList.begin();
1078 //if(locFcalShowerClassificationIterator != locFcalShowerClassificationList.end()) {
1079 // shower->setClassifierOutput(locFcalShowerClassificationIterator->getClassifierOuput());
1080 //}
1081
1082 // shower shape and other parameters. used e.g. as input to MVA classifier
1083 const hddm_r::FcalShowerPropertiesList& locFcalShowerPropertiesList = iter->getFcalShowerPropertiesList();
1084 hddm_r::FcalShowerPropertiesList::iterator locFcalShowerPropertiesIterator = locFcalShowerPropertiesList.begin();
1085 if(locFcalShowerPropertiesIterator != locFcalShowerPropertiesList.end()) {
1086 shower->setDocaTrack(locFcalShowerPropertiesIterator->getDocaTrack());
1087 shower->setTimeTrack(locFcalShowerPropertiesIterator->getTimeTrack());
1088 shower->setSumU(locFcalShowerPropertiesIterator->getSumU());
1089 shower->setSumV(locFcalShowerPropertiesIterator->getSumV());
1090 shower->setE1E9(locFcalShowerPropertiesIterator->getE1E9());
1091 shower->setE9E25(locFcalShowerPropertiesIterator->getE9E25());
1092 }
1093
1094 const hddm_r::FcalShowerNBlocksList& locFcalShowerNBlocksList = iter->getFcalShowerNBlockses();
1095 hddm_r::FcalShowerNBlocksList::iterator locFcalShowerNBlocksIterator = locFcalShowerNBlocksList.begin();
1096 if(locFcalShowerNBlocksIterator != locFcalShowerNBlocksList.end()) {
1097 shower->setNumBlocks(locFcalShowerNBlocksIterator->getNumBlocks());
1098 }
1099 data.push_back(shower);
1100 }
1101
1102 // Copy into factory
1103 factory->CopyTo(data);
1104
1105 return NOERROR;
1106}
1107
1108//-----------------------
1109// Extract_DBCALShower
1110//-----------------------
1111jerror_t DEventSourceREST::Extract_DBCALShower(hddm_r::HDDM *record,
1112 JFactory<DBCALShower>* factory)
1113{
1114 /// Copies the data from the bcalShower hddm record. This is
1115 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1116 /// returns OBJECT_NOT_AVAILABLE immediately.
1117
1118 if (factory==NULL__null) {
1119 return OBJECT_NOT_AVAILABLE;
1120 }
1121 string tag = (factory->Tag())? factory->Tag() : "";
1122
1123 vector<DBCALShower*> data;
1124
1125 // loop over bcal shower records
1126 const hddm_r::BcalShowerList &showers =
1127 record->getBcalShowers();
1128 hddm_r::BcalShowerList::iterator iter;
1129 for (iter = showers.begin(); iter != showers.end(); ++iter) {
1130 if (iter->getJtag() != tag)
1131 continue;
1132
1133 DBCALShower *shower = new DBCALShower();
1134 shower->E = iter->getE();
1135 shower->E_raw = -1;
1136 shower->x = iter->getX();
1137 shower->y = iter->getY();
1138 shower->z = iter->getZ();
1139 shower->t = iter->getT();
1140 shower->Q = 0; // Fix this to zero for now, can add to REST if it's ever used in higher-level analyses
1141
1142 if(USE_CCDB_BCAL_COVARIANCE) {
1143 dBCALShowerFactory->FillCovarianceMatrix(shower);
1144 } else {
1145 TMatrixFSym covariance(5);
1146 covariance(0,0) = iter->getEerr()*iter->getEerr();
1147 covariance(1,1) = iter->getXerr()*iter->getXerr();
1148 covariance(2,2) = iter->getYerr()*iter->getYerr();
1149 covariance(3,3) = iter->getZerr()*iter->getZerr();
1150 covariance(4,4) = iter->getTerr()*iter->getTerr();
1151 covariance(1,2) = covariance(2,1) = iter->getXycorr()*iter->getXerr()*iter->getYerr();
1152 covariance(1,3) = covariance(3,1) = iter->getXzcorr()*iter->getXerr()*iter->getZerr();
1153 covariance(2,3) = covariance(3,2) = iter->getYzcorr()*iter->getYerr()*iter->getZerr();
1154 covariance(0,3) = covariance(3,0) = iter->getEzcorr()*iter->getEerr()*iter->getZerr();
1155 covariance(3,4) = covariance(4,3) = iter->getTzcorr()*iter->getTerr()*iter->getZerr();
1156
1157 // further correlations (an extension of REST format, so code is different.)
1158 const hddm_r::BcalCorrelationsList& locBcalCorrelationsList = iter->getBcalCorrelationses();
1159 hddm_r::BcalCorrelationsList::iterator locBcalCorrelationsIterator = locBcalCorrelationsList.begin();
1160 if(locBcalCorrelationsIterator != locBcalCorrelationsList.end()) {
1161 covariance(0,4) = covariance(4,0) = locBcalCorrelationsIterator->getEtcorr()*iter->getEerr()*iter->getTerr();
1162 covariance(0,1) = covariance(1,0) = locBcalCorrelationsIterator->getExcorr()*iter->getEerr()*iter->getXerr();
1163 covariance(0,2) = covariance(2,0) = locBcalCorrelationsIterator->getEycorr()*iter->getEerr()*iter->getYerr();
1164 covariance(1,4) = covariance(4,1) = locBcalCorrelationsIterator->getTxcorr()*iter->getTerr()*iter->getXerr();
1165 covariance(2,4) = covariance(4,2) = locBcalCorrelationsIterator->getTycorr()*iter->getTerr()*iter->getYerr();
1166 }
1167 shower->ExyztCovariance = covariance;
1168 }
1169
1170 // preshower
1171 const hddm_r::PreshowerList& locPreShowerList = iter->getPreshowers();
1172 hddm_r::PreshowerList::iterator locPreShowerIterator = locPreShowerList.begin();
1173 if(locPreShowerIterator == locPreShowerList.end())
1174 shower->E_preshower = 0.0;
1175 else //should only be 1
1176 {
1177 for(; locPreShowerIterator != locPreShowerList.end(); ++locPreShowerIterator)
1178 shower->E_preshower = locPreShowerIterator->getPreshowerE();
1179 }
1180
1181 // width
1182 const hddm_r::WidthList& locWidthList = iter->getWidths();
1183 hddm_r::WidthList::iterator locWidthIterator = locWidthList.begin();
1184 if(locWidthIterator == locWidthList.end()) {
1185 shower->sigLong = -1.;
1186 shower->sigTrans = -1.;
1187 shower->sigTheta = -1.;
1188 }
1189 else //should only be 1
1190 {
1191 for(; locWidthIterator != locWidthList.end(); ++locWidthIterator) {
1192 shower->sigLong = locWidthIterator->getSigLong();
1193 shower->sigTrans = locWidthIterator->getSigTrans();
1194 shower->sigTheta = locWidthIterator->getSigTheta();
1195 }
1196 }
1197
1198 const hddm_r::BcalClusterList& locBcalClusterList = iter->getBcalClusters();
1199 hddm_r::BcalClusterList::iterator locBcalClusterIterator = locBcalClusterList.begin();
1200 if(locBcalClusterIterator == locBcalClusterList.end())
1201 shower->N_cell = -1;
1202 else //should only be 1
1203 {
1204 for(; locBcalClusterIterator != locBcalClusterList.end(); ++locBcalClusterIterator)
1205 shower->N_cell = locBcalClusterIterator->getNcell();
1206 }
1207
1208 const hddm_r::BcalLayersList& locBcalLayersList = iter->getBcalLayerses();
1209 hddm_r::BcalLayersList::iterator locBcalLayersIterator = locBcalLayersList.begin();
1210 if(locBcalLayersIterator == locBcalLayersList.end()) {
1211 shower->E_L2 = 0.;
1212 shower->E_L3 = 0.;
1213 shower->E_L4 = 0.;
1214 shower->rmsTime = -1;
1215 }
1216 else //should only be 1
1217 {
1218 for(; locBcalLayersIterator != locBcalLayersList.end(); ++locBcalLayersIterator) {
1219 shower->rmsTime = locBcalLayersIterator->getRmsTime();
1220 shower->E_L2 = locBcalLayersIterator->getE_L2();
1221 shower->E_L3 = locBcalLayersIterator->getE_L3();
1222 shower->E_L4 = locBcalLayersIterator->getE_L4();
1223 }
1224 }
1225
1226 data.push_back(shower);
1227 }
1228
1229 // Copy into factory
1230 factory->CopyTo(data);
1231
1232 return NOERROR;
1233}
1234
1235//-----------------------
1236// Extract_DCCALShower
1237//-----------------------
1238jerror_t DEventSourceREST::Extract_DCCALShower(hddm_r::HDDM *record,
1239 JFactory<DCCALShower>* factory)
1240{
1241 /// Copies the data from the ccalShower hddm record. This is
1242 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1243 /// returns OBJECT_NOT_AVAILABLE immediately.
1244
1245 if (factory==NULL__null) {
1246 return OBJECT_NOT_AVAILABLE;
1247 }
1248 string tag = (factory->Tag())? factory->Tag() : "";
1249
1250 vector<DCCALShower*> data;
1251
1252 // loop over ccal shower records
1253 const hddm_r::CcalShowerList &showers =
1254 record->getCcalShowers();
1255 hddm_r::CcalShowerList::iterator iter;
1256 for (iter = showers.begin(); iter != showers.end(); ++iter) {
1257 if (iter->getJtag() != tag)
1258 continue;
1259
1260 DCCALShower *shower = new DCCALShower();
1261 shower->E = iter->getE();
1262 shower->x = iter->getX();
1263 shower->y = iter->getY();
1264 shower->z = iter->getZ();
1265 shower->time = iter->getT();
1266 shower->sigma_t = iter->getTerr();
1267 shower->sigma_E = iter->getEerr();
1268 shower->Emax = iter->getEmax();
1269 shower->x1 = iter->getX1();
1270 shower->y1 = iter->getY1();
1271 shower->chi2 = iter->getChi2();
1272
1273 shower->type = iter->getType();
1274 shower->dime = iter->getDime();
1275 shower->id = iter->getId();
1276 shower->idmax = iter->getIdmax();
1277
1278 data.push_back(shower);
1279 }
1280
1281 // Copy into factory
1282 factory->CopyTo(data);
1283
1284 return NOERROR;
1285}
1286
1287//--------------------------------
1288// Extract_DTrackTimeBased
1289//--------------------------------
1290jerror_t DEventSourceREST::Extract_DTrackTimeBased(hddm_r::HDDM *record,
1291 JFactory<DTrackTimeBased>* factory, JEventLoop* locEventLoop)
1292{
1293 /// Copies the data from the chargedTrack hddm record. This is
1294 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1295 /// returns OBJECT_NOT_AVAILABLE immediately.
1296
1297 if (factory==NULL__null) {
1298 return OBJECT_NOT_AVAILABLE;
1299 }
1300
1301 int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
1302 DVector2 locBeamCenter,locBeamDir;
1303 double locBeamZ0=0.;
1304 LockRead();
1305 {
1306 locBeamCenter = dBeamCenterMap[locRunNumber];
1307 locBeamDir = dBeamDirMap[locRunNumber];
1308 locBeamZ0 = dBeamZ0Map[locRunNumber];
1309 }
1310 UnlockRead();
1311
1312 string tag = (factory->Tag())? factory->Tag() : "";
1313
1314 vector<DTrackTimeBased*> data;
1315
1316 // loop over chargedTrack records
1317 const hddm_r::ChargedTrackList &tracks = record->getChargedTracks();
1318 hddm_r::ChargedTrackList::iterator iter;
1319 for (iter = tracks.begin(); iter != tracks.end(); ++iter) {
1320 if (iter->getJtag() != tag) {
1321 continue;
1322 }
1323 DTrackTimeBased *tra = new DTrackTimeBased();
1324 tra->trackid = 0;
1325 tra->candidateid = iter->getCandidateId();
1326 Particle_t ptype = iter->getPtype();
1327 tra->setPID(ptype);
1328
1329 const hddm_r::TrackFit &fit = iter->getTrackFit();
1330 tra->Ndof = fit.getNdof();
1331 tra->chisq = fit.getChisq();
1332 tra->FOM = TMath::Prob(tra->chisq, tra->Ndof);
1333 tra->setTime(fit.getT0());
1334 DVector3 track_pos(fit.getX0(),fit.getY0(),fit.getZ0());
1335 DVector3 track_mom(fit.getPx(),fit.getPy(),fit.getPz());
1336 tra->setPosition(track_pos);
1337 tra->setMomentum(track_mom);
1338
1339 auto loc5x5ErrorMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
1340 loc5x5ErrorMatrix->ResizeTo(5, 5);
1341 (*loc5x5ErrorMatrix)(0,0) = fit.getE11();
1342 (*loc5x5ErrorMatrix)(0,1) = (*loc5x5ErrorMatrix)(1,0) = fit.getE12();
1343 (*loc5x5ErrorMatrix)(0,2) = (*loc5x5ErrorMatrix)(2,0) = fit.getE13();
1344 (*loc5x5ErrorMatrix)(0,3) = (*loc5x5ErrorMatrix)(3,0) = fit.getE14();
1345 (*loc5x5ErrorMatrix)(0,4) = (*loc5x5ErrorMatrix)(4,0) = fit.getE15();
1346 (*loc5x5ErrorMatrix)(1,1) = fit.getE22();
1347 (*loc5x5ErrorMatrix)(1,2) = (*loc5x5ErrorMatrix)(2,1) = fit.getE23();
1348 (*loc5x5ErrorMatrix)(1,3) = (*loc5x5ErrorMatrix)(3,1) = fit.getE24();
1349 (*loc5x5ErrorMatrix)(1,4) = (*loc5x5ErrorMatrix)(4,1) = fit.getE25();
1350 (*loc5x5ErrorMatrix)(2,2) = fit.getE33();
1351 (*loc5x5ErrorMatrix)(2,3) = (*loc5x5ErrorMatrix)(3,2) = fit.getE34();
1352 (*loc5x5ErrorMatrix)(2,4) = (*loc5x5ErrorMatrix)(4,2) = fit.getE35();
1353 (*loc5x5ErrorMatrix)(3,3) = fit.getE44();
1354 (*loc5x5ErrorMatrix)(3,4) = (*loc5x5ErrorMatrix)(4,3) = fit.getE45();
1355 (*loc5x5ErrorMatrix)(4,4) = fit.getE55();
1356 tra->setTrackingErrorMatrix(loc5x5ErrorMatrix);
1357
1358 // Convert from cartesian coordinates to the 5x1 state vector corresponding to the tracking error matrix.
1359 double vect[5];
1360 DVector2 beam_pos=locBeamCenter+(track_pos.Z()-locBeamZ0)*locBeamDir;
1361 DVector2 diff(track_pos.X()-beam_pos.X(),track_pos.Y()-beam_pos.Y());
1362 vect[2]=tan(M_PI_21.57079632679489661923 - track_mom.Theta());
1363 vect[1]=track_mom.Phi();
1364 double sinphi=sin(vect[1]);
1365 double cosphi=cos(vect[1]);
1366 vect[0]=tra->charge()/track_mom.Perp();
1367 vect[4]=track_pos.Z();
1368 vect[3]=diff.Mod();
1369
1370 if ((diff.X() > 0 && sinphi>0) || (diff.Y() <0 && cosphi>0) || (diff.Y() >0 && cosphi<0) || (diff.X() <0 && sinphi<0))
1371 vect[3] *= -1.;
1372 tra->setTrackingStateVector(vect[0], vect[1], vect[2], vect[3], vect[4]);
1373
1374 // Set the 7x7 covariance matrix.
1375 auto loc7x7ErrorMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
1376 loc7x7ErrorMatrix->ResizeTo(7, 7);
1377 Get7x7ErrorMatrix(tra->mass(), vect, loc5x5ErrorMatrix.get(), loc7x7ErrorMatrix.get());
1378 tra->setErrorMatrix(loc7x7ErrorMatrix);
1379 (*loc7x7ErrorMatrix)(6, 6) = fit.getT0err()*fit.getT0err();
1380
1381 // Track parameters at exit of tracking volume
1382 const hddm_r::ExitParamsList& locExitParamsList = iter->getExitParamses();
1383 hddm_r::ExitParamsList::iterator locExitParamsIterator = locExitParamsList.begin();
1384 if (locExitParamsIterator!=locExitParamsList.end()){
1385 // Create the extrapolation vector
1386 vector<DTrackFitter::Extrapolation_t>myvector;
1387 tra->extrapolations.emplace(SYS_NULL,myvector);
1388
1389 for(; locExitParamsIterator != locExitParamsList.end(); ++locExitParamsIterator){
1390 DVector3 pos(locExitParamsIterator->getX1(),
1391 locExitParamsIterator->getY1(),
1392 locExitParamsIterator->getZ1());
1393 DVector3 mom(locExitParamsIterator->getPx1(),
1394 locExitParamsIterator->getPy1(),
1395 locExitParamsIterator->getPz1());
1396 tra->extrapolations[SYS_NULL].push_back(DTrackFitter::Extrapolation_t(pos,mom,locExitParamsIterator->getT1(),0.));
1397 }
1398 }
1399 // Hit layers
1400 const hddm_r::ExpectedhitsList& locExpectedhitsList = iter->getExpectedhitses();
1401 hddm_r::ExpectedhitsList::iterator locExpectedhitsIterator = locExpectedhitsList.begin();
1402 if(locExpectedhitsIterator == locExpectedhitsList.end())
1403 {
1404 tra->potential_cdc_hits_on_track = 0;
1405 tra->potential_fdc_hits_on_track = 0;
1406 tra->measured_cdc_hits_on_track = 0;
1407 tra->measured_fdc_hits_on_track = 0;
1408 //tra->cdc_hit_usage.total_hits = 0;
1409 //tra->fdc_hit_usage.total_hits = 0;
1410 }
1411 else //should only be 1
1412 {
1413 for(; locExpectedhitsIterator != locExpectedhitsList.end(); ++locExpectedhitsIterator)
1414 {
1415 tra->potential_cdc_hits_on_track = locExpectedhitsIterator->getExpectedCDChits();
1416 tra->potential_fdc_hits_on_track = locExpectedhitsIterator->getExpectedFDChits();
1417 tra->measured_cdc_hits_on_track = locExpectedhitsIterator->getMeasuredCDChits();
1418 tra->measured_fdc_hits_on_track = locExpectedhitsIterator->getMeasuredFDChits();
1419 //tra->cdc_hit_usage.total_hits = locExpectedhitsIterator->getMeasuredCDChits();
1420 //tra->fdc_hit_usage.total_hits = locExpectedhitsIterator->getMeasuredFDChits();
1421 }
1422 }
1423
1424 // Expected number of hits
1425 const hddm_r::HitlayersList& locHitlayersList = iter->getHitlayerses();
1426 hddm_r::HitlayersList::iterator locHitlayersIterator = locHitlayersList.begin();
1427 if(locHitlayersIterator == locHitlayersList.end())
1428 {
1429 tra->dCDCRings = 0;
1430 tra->dFDCPlanes = 0;
1431 }
1432 else //should only be 1
1433 {
1434 for(; locHitlayersIterator != locHitlayersList.end(); ++locHitlayersIterator)
1435 {
1436 tra->dCDCRings = locHitlayersIterator->getCDCrings();
1437 tra->dFDCPlanes = locHitlayersIterator->getFDCplanes();
1438 }
1439 }
1440
1441 // MC match hit info
1442 const hddm_r::McmatchList& locMCMatchesList = iter->getMcmatchs();
1443 hddm_r::McmatchList::iterator locMcmatchIterator = locMCMatchesList.begin();
1444 if(locMcmatchIterator == locMCMatchesList.end())
1445 {
1446 tra->dCDCRings = 0;
1447 tra->dFDCPlanes = 0;
1448 }
1449 else //should only be 1
1450 {
1451 for(; locMcmatchIterator != locMCMatchesList.end(); ++locMcmatchIterator)
1452 {
1453 tra->dMCThrownMatchMyID = locMcmatchIterator->getIthrown();
1454 tra->dNumHitsMatchedToThrown = locMcmatchIterator->getNumhitsmatch();
1455 }
1456 }
1457
1458 // add the drift chamber dE/dx information
1459 const hddm_r::DEdxDCList &el = iter->getDEdxDCs();
1460 hddm_r::DEdxDCList::iterator diter = el.begin();
1461 tra->ddx_CDC_trunc.clear();
1462 tra->ddx_CDC_amp_trunc.clear();
1463 tra->ddEdx_CDC_trunc.clear();
1464 tra->ddEdx_CDC_amp_trunc.clear();
1465 tra->ddx_FDC_trunc.clear();
1466 tra->ddx_FDC_amp_trunc.clear();
1467 tra->ddEdx_FDC_trunc.clear();
1468 tra->ddEdx_FDC_amp_trunc.clear();
1469 if (diter != el.end()) {
1470 tra->dNumHitsUsedFordEdx_FDC = diter->getNsampleFDC();
1471 tra->dNumHitsUsedFordEdx_CDC = diter->getNsampleCDC();
1472 tra->ddEdx_FDC = diter->getDEdxFDC();
1473 tra->ddEdx_CDC = diter->getDEdxCDC();
1474 tra->ddx_FDC = diter->getDxFDC();
1475 tra->ddx_CDC = diter->getDxCDC();
1476 const hddm_r::CDCAmpdEdxList &el2 = diter->getCDCAmpdEdxs();
1477 hddm_r::CDCAmpdEdxList::iterator diter2 = el2.begin();
1478 if (diter2 != el2.end()){
1479 tra->ddx_CDC_amp= diter2->getDxCDCAmp();
1480 tra->ddEdx_CDC_amp = diter2->getDEdxCDCAmp();
1481 }
1482 else{
1483 tra->ddx_CDC_amp=tra->ddx_CDC;
1484 tra->ddEdx_CDC_amp=tra->ddEdx_CDC;
1485 }
1486 const hddm_r::CDCdEdxTruncList &cdctruncs = diter->getCDCdEdxTruncs();
1487 hddm_r::CDCdEdxTruncList::iterator itcdc;
1488 for (itcdc = cdctruncs.begin(); itcdc != cdctruncs.end(); ++itcdc) {
1489 int ntrunc = itcdc->getNtrunc();
1490 for (int s=(int)tra->ddx_CDC_trunc.size(); s <= ntrunc; ++s) {
1491 tra->ddx_CDC_trunc.push_back(0);
1492 tra->ddx_CDC_amp_trunc.push_back(0);
1493 tra->ddEdx_CDC_trunc.push_back(0);
1494 tra->ddEdx_CDC_amp_trunc.push_back(0);
1495 }
1496 tra->ddx_CDC_trunc[ntrunc] = itcdc->getDx();
1497 tra->ddx_CDC_amp_trunc[ntrunc] = itcdc->getDxAmp();
1498 tra->ddEdx_CDC_trunc[ntrunc] = itcdc->getDEdx();
1499 tra->ddEdx_CDC_amp_trunc[ntrunc] = itcdc->getDEdxAmp();
1500 }
1501 const hddm_r::FDCdEdxTruncList &fdctruncs = diter->getFDCdEdxTruncs();
1502 hddm_r::FDCdEdxTruncList::iterator itfdc;
1503 for (itfdc = fdctruncs.begin(); itfdc != fdctruncs.end(); ++itfdc) {
1504 int ntrunc = itfdc->getNtrunc();
1505 for (int s=(int)tra->ddx_FDC_trunc.size(); s <= ntrunc; ++s) {
1506 tra->ddx_FDC_trunc.push_back(0);
1507 tra->ddx_FDC_amp_trunc.push_back(0);
1508 tra->ddEdx_FDC_trunc.push_back(0);
1509 tra->ddEdx_FDC_amp_trunc.push_back(0);
1510 }
1511 tra->ddx_FDC_trunc[ntrunc] = itfdc->getDx();
1512 tra->ddx_FDC_amp_trunc[ntrunc] = itfdc->getDxAmp();
1513 tra->ddEdx_FDC_trunc[ntrunc] = itfdc->getDEdx();
1514 tra->ddEdx_FDC_amp_trunc[ntrunc] = itfdc->getDEdxAmp();
1515 }
1516 }
1517 else {
1518 tra->dNumHitsUsedFordEdx_FDC = 0;
1519 tra->dNumHitsUsedFordEdx_CDC = 0;
1520 tra->ddEdx_FDC = 0.0;
1521 tra->ddEdx_CDC = 0.0;
1522 tra->ddx_FDC = 0.0;
1523 tra->ddx_CDC = 0.0;
1524 tra->ddEdx_CDC_amp = 0.0;
1525 tra->ddx_CDC_amp = 0.0;
1526 }
1527
1528 data.push_back(tra);
1529 }
1530
1531 if( PRUNE_DUPLICATE_TRACKS && (data.size() > 1) ) {
1532 vector< int > indices_to_erase;
1533
1534 for( unsigned int i=0; i<data.size()-1; i++ ) {
1535 for( unsigned int j=i+1; j<data.size(); j++ ) {
1536 if(find(indices_to_erase.begin(), indices_to_erase.end(), j) != indices_to_erase.end())
1537 continue;
1538
1539 // look through the remaining tracks for duplicates
1540 // (1) if there is a track with the same candidate/PID and worse chi^2, reject that track
1541 // (2) if there is a track with the same candidate/PID and better chi^2, reject this track
1542 if( (data[i]->candidateid == data[j]->candidateid)
1543 && (data[i]->PID() == data[j]->PID()) ) { // is a duplicate track
1544 if(data[i]->chisq < data[j]->chisq) {
1545 indices_to_erase.push_back(j);
1546 } else {
1547 indices_to_erase.push_back(i);
1548 }
1549 }
1550 }
1551 }
1552
1553 // create the new set of tracks
1554 vector<DTrackTimeBased*> new_data;
1555 for( unsigned int i=0; i<data.size(); i++ ) {
1556 if(find(indices_to_erase.begin(), indices_to_erase.end(), i) != indices_to_erase.end())
1557 continue;
1558
1559 new_data.push_back(data[i]);
1560 }
1561 data = new_data; // replace the set of tracks with the pruned one
1562 }
1563
1564 // Copy into factory
1565 factory->CopyTo(data);
1566
1567 return NOERROR;
1568}
1569
1570//--------------------------------
1571// Extract_DDetectorMatches
1572//--------------------------------
1573jerror_t DEventSourceREST::Extract_DDetectorMatches(JEventLoop* locEventLoop, hddm_r::HDDM *record, JFactory<DDetectorMatches>* factory)
1574{
1575 /// Copies the data from the detectorMatches hddm record. This is
1576 /// called from JEventSourceREST::GetObjects. If factory is NULL, this
1577 /// returns OBJECT_NOT_AVAILABLE immediately.
1578
1579 if(factory==NULL__null)
1580 return OBJECT_NOT_AVAILABLE;
1581
1582 string tag = (factory->Tag())? factory->Tag() : "";
1583 vector<DDetectorMatches*> data;
1584
1585 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
1586 locEventLoop->Get(locTrackTimeBasedVector);
1587
1588 vector<const DSCHit*> locSCHits;
1589 locEventLoop->Get(locSCHits);
1590
1591 vector<const DTOFPoint*> locTOFPoints;
1592 locEventLoop->Get(locTOFPoints);
1593
1594 vector<const DCTOFPoint*> locCTOFPoints;
1595 locEventLoop->Get(locCTOFPoints);
1596
1597 vector<const DBCALShower*> locBCALShowers;
1598 locEventLoop->Get(locBCALShowers);
1599
1600 vector<const DFCALShower*> locFCALShowers;
1601 locEventLoop->Get(locFCALShowers);
1602
1603 const DParticleID* locParticleID = NULL__null;
1604 vector<const DDIRCPmtHit*> locDIRCHits;
1605 vector<const DDIRCTruthBarHit*> locDIRCBarHits;
1606 if(RECO_DIRC_CALC_LUT) {
1607 locEventLoop->GetSingle(locParticleID);
1608 locEventLoop->Get(locDIRCHits);
1609 locEventLoop->Get(locDIRCBarHits);
1610 }
1611
1612 const hddm_r::DetectorMatchesList &detectormatches = record->getDetectorMatcheses();
1613
1614 // loop over chargedTrack records
1615 hddm_r::DetectorMatchesList::iterator iter;
1616 for(iter = detectormatches.begin(); iter != detectormatches.end(); ++iter)
1617 {
1618 if(iter->getJtag() != tag)
1619 continue;
1620
1621 DDetectorMatches *locDetectorMatches = new DDetectorMatches();
1622
1623 const hddm_r::DircMatchParamsList &dircList = iter->getDircMatchParamses();
1624 hddm_r::DircMatchParamsList::iterator dircIter = dircList.begin();
1625 const hddm_r::DircMatchHitList &dircMatchHitList = iter->getDircMatchHits();
1626
1627 for(; dircIter != dircList.end(); ++dircIter)
1628 {
1629 size_t locTrackIndex = dircIter->getTrack();
1630 if(locTrackIndex >= locTrackTimeBasedVector.size()) continue;
1631
1632 auto locTrackTimeBased = locTrackTimeBasedVector[locTrackIndex];
1633 if( !locTrackTimeBased ) continue;
1634
1635 auto locDIRCMatchParams = std::make_shared<DDIRCMatchParams>();
1636 map<shared_ptr<const DDIRCMatchParams> ,vector<const DDIRCPmtHit*> > locDIRCTrackMatchParams;
1637 locDetectorMatches->Get_DIRCTrackMatchParamsMap(locDIRCTrackMatchParams);
1638
1639 if(RECO_DIRC_CALC_LUT) {
1640 TVector3 locProjPos(dircIter->getX(),dircIter->getY(),dircIter->getZ());
1641 TVector3 locProjMom(dircIter->getPx(),dircIter->getPy(),dircIter->getPz());
1642 double locFlightTime = dircIter->getT();
1643
1644 if( locParticleID->Get_DIRCLut()->CalcLUT(locProjPos, locProjMom, locDIRCHits, locFlightTime, locTrackTimeBased->mass(), locDIRCMatchParams, locDIRCBarHits, locDIRCTrackMatchParams) )
1645 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], std::const_pointer_cast<const DDIRCMatchParams>(locDIRCMatchParams));
1646 }
1647 else {
1648 // add hits to match list
1649 hddm_r::DircMatchHitList::iterator dircMatchHitIter = dircMatchHitList.begin();
1650 for(; dircMatchHitIter != dircMatchHitList.end(); ++dircMatchHitIter) {
1651 size_t locMatchHitTrackIndex = dircMatchHitIter->getTrack();
1652 if(locMatchHitTrackIndex == locTrackIndex) {
1653 size_t locMatchHitIndex = dircMatchHitIter->getHit();
1654 locDIRCTrackMatchParams[locDIRCMatchParams].push_back(locDIRCHits[locMatchHitIndex]);
1655 }
1656 }
1657
1658 locDIRCMatchParams->dExtrapolatedPos = DVector3(dircIter->getX(),dircIter->getY(),dircIter->getZ());
1659 locDIRCMatchParams->dExtrapolatedMom = DVector3(dircIter->getPx(),dircIter->getPy(),dircIter->getPz());
1660 locDIRCMatchParams->dExtrapolatedTime = dircIter->getT();
1661 locDIRCMatchParams->dExpectedThetaC = dircIter->getExpectthetac();
1662 locDIRCMatchParams->dThetaC = dircIter->getThetac();
1663 locDIRCMatchParams->dDeltaT = dircIter->getDeltat();
1664 locDIRCMatchParams->dLikelihoodElectron = dircIter->getLele();
1665 locDIRCMatchParams->dLikelihoodPion = dircIter->getLpi();
1666 locDIRCMatchParams->dLikelihoodKaon = dircIter->getLk();
1667 locDIRCMatchParams->dLikelihoodProton = dircIter->getLp();
1668 locDIRCMatchParams->dNPhotons = dircIter->getNphotons();
1669 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], std::const_pointer_cast<const DDIRCMatchParams>(locDIRCMatchParams));
1670 }
1671 }
1672
1673 const hddm_r::BcalMatchParamsList &bcalList = iter->getBcalMatchParamses();
1674 hddm_r::BcalMatchParamsList::iterator bcalIter = bcalList.begin();
1675 for(; bcalIter != bcalList.end(); ++bcalIter)
1676 {
1677 size_t locShowerIndex = bcalIter->getShower();
1678 size_t locTrackIndex = bcalIter->getTrack();
1679
1680 auto locShowerMatchParams = std::make_shared<DBCALShowerMatchParams>();
1681 locShowerMatchParams->dBCALShower = locBCALShowers[locShowerIndex];
1682 locShowerMatchParams->dx = bcalIter->getDx();
1683 locShowerMatchParams->dFlightTime = bcalIter->getTflight();
1684 locShowerMatchParams->dFlightTimeVariance = bcalIter->getTflightvar();
1685 locShowerMatchParams->dPathLength = bcalIter->getPathlength();
1686 locShowerMatchParams->dDeltaPhiToShower = bcalIter->getDeltaphi();
1687 locShowerMatchParams->dDeltaZToShower = bcalIter->getDeltaz();
1688
1689 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locBCALShowers[locShowerIndex], std::const_pointer_cast<const DBCALShowerMatchParams>(locShowerMatchParams));
1690 }
1691
1692 const hddm_r::FcalMatchParamsList &fcalList = iter->getFcalMatchParamses();
1693 hddm_r::FcalMatchParamsList::iterator fcalIter = fcalList.begin();
1694 for(; fcalIter != fcalList.end(); ++fcalIter)
1695 {
1696 size_t locShowerIndex = fcalIter->getShower();
1697 size_t locTrackIndex = fcalIter->getTrack();
1698
1699 auto locShowerMatchParams = std::make_shared<DFCALShowerMatchParams>();
1700 locShowerMatchParams->dFCALShower = locFCALShowers[locShowerIndex];
1701 locShowerMatchParams->dx = fcalIter->getDx();
1702 locShowerMatchParams->dFlightTime = fcalIter->getTflight();
1703 locShowerMatchParams->dFlightTimeVariance = fcalIter->getTflightvar();
1704 locShowerMatchParams->dPathLength = fcalIter->getPathlength();
1705 locShowerMatchParams->dDOCAToShower = fcalIter->getDoca();
1706 locShowerMatchParams->dEcenter=0.;
1707 locShowerMatchParams->dE3x3=0.;
1708 locShowerMatchParams->dE5x5=0.;
1709 const hddm_r::FcalEnergyParamsList &fcalEnergyList = fcalIter->getFcalEnergyParamses();
1710 hddm_r::FcalEnergyParamsList::iterator fcalEnergyIter = fcalEnergyList.begin();
1711 for(; fcalEnergyIter != fcalEnergyList.end(); ++fcalEnergyIter){
1712 locShowerMatchParams->dEcenter=fcalEnergyIter->getEcenter();
1713 locShowerMatchParams->dE3x3=fcalEnergyIter->getE3x3();
1714 locShowerMatchParams->dE5x5=fcalEnergyIter->getE5x5();
1715 }
1716 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locFCALShowers[locShowerIndex], std::const_pointer_cast<const DFCALShowerMatchParams>(locShowerMatchParams));
1717 }
1718
1719 const hddm_r::FcalSingleHitMatchParamsList &fcalSingleHitList = iter->getFcalSingleHitMatchParamses();
1720 hddm_r::FcalSingleHitMatchParamsList::iterator fcalSingleHitIter = fcalSingleHitList.begin();
1721 for(; fcalSingleHitIter != fcalSingleHitList.end(); ++fcalSingleHitIter)
1722 {
1723 size_t locTrackIndex = fcalSingleHitIter->getTrack();
1724
1725 auto locSingleHitMatchParams = std::make_shared<DFCALSingleHitMatchParams>();
1726 locSingleHitMatchParams->dEHit = fcalSingleHitIter->getEhit();
1727 locSingleHitMatchParams->dTHit = fcalSingleHitIter->getThit();
1728 locSingleHitMatchParams->dx = fcalSingleHitIter->getDx();
1729 locSingleHitMatchParams->dFlightTime = fcalSingleHitIter->getTflight();
1730 locSingleHitMatchParams->dFlightTimeVariance = fcalSingleHitIter->getTflightvar();
1731 locSingleHitMatchParams->dPathLength = fcalSingleHitIter->getPathlength();
1732 locSingleHitMatchParams->dDOCAToHit = fcalSingleHitIter->getDoca();
1733
1734 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex],std::const_pointer_cast<const DFCALSingleHitMatchParams>(locSingleHitMatchParams));
1735 }
1736
1737
1738 const hddm_r::ScMatchParamsList &scList = iter->getScMatchParamses();
1739 hddm_r::ScMatchParamsList::iterator scIter = scList.begin();
1740 for(; scIter != scList.end(); ++scIter)
1741 {
1742 size_t locHitIndex = scIter->getHit();
1743 size_t locTrackIndex = scIter->getTrack();
1744
1745 auto locSCHitMatchParams = std::make_shared<DSCHitMatchParams>();
1746 locSCHitMatchParams->dSCHit = locSCHits[locHitIndex];
1747 locSCHitMatchParams->dEdx = scIter->getDEdx();
1748 locSCHitMatchParams->dHitTime = scIter->getThit();
1749 locSCHitMatchParams->dHitTimeVariance = scIter->getThitvar();
1750 locSCHitMatchParams->dHitEnergy = scIter->getEhit();
1751 locSCHitMatchParams->dFlightTime = scIter->getTflight();
1752 locSCHitMatchParams->dFlightTimeVariance = scIter->getTflightvar();
1753 locSCHitMatchParams->dPathLength = scIter->getPathlength();
1754 locSCHitMatchParams->dDeltaPhiToHit = scIter->getDeltaphi();
1755
1756 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locSCHits[locHitIndex], std::const_pointer_cast<const DSCHitMatchParams>(locSCHitMatchParams));
1757 }
1758
1759 const hddm_r::TofMatchParamsList &tofList = iter->getTofMatchParamses();
1760 hddm_r::TofMatchParamsList::iterator tofIter = tofList.begin();
1761 for(; tofIter != tofList.end(); ++tofIter)
1762 {
1763 size_t locHitIndex = tofIter->getHit();
1764 size_t locTrackIndex = tofIter->getTrack();
1765
1766 auto locTOFHitMatchParams = std::make_shared<DTOFHitMatchParams>();
1767 locTOFHitMatchParams->dTOFPoint = locTOFPoints[locHitIndex];
1768
1769 locTOFHitMatchParams->dHitTime = tofIter->getThit();
1770 locTOFHitMatchParams->dHitTimeVariance = tofIter->getThitvar();
1771 locTOFHitMatchParams->dHitEnergy = tofIter->getEhit();
1772
1773 locTOFHitMatchParams->dEdx = tofIter->getDEdx();
1774 locTOFHitMatchParams->dFlightTime = tofIter->getTflight();
1775 locTOFHitMatchParams->dFlightTimeVariance = tofIter->getTflightvar();
1776 locTOFHitMatchParams->dPathLength = tofIter->getPathlength();
1777 locTOFHitMatchParams->dDeltaXToHit = tofIter->getDeltax();
1778 locTOFHitMatchParams->dDeltaYToHit = tofIter->getDeltay();
1779
1780 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locTOFPoints[locHitIndex], std::const_pointer_cast<const DTOFHitMatchParams>(locTOFHitMatchParams));
1781
1782 // dE/dx per plane
1783 const hddm_r::TofDedxList& locTofDedxList = tofIter->getTofDedxs();
1784 hddm_r::TofDedxList::iterator locTofDedxIterator = locTofDedxList.begin();
1785 if(locTofDedxIterator == locTofDedxList.end())
1786 {
1787 locTOFHitMatchParams->dEdx1 = 0.;
1788 locTOFHitMatchParams->dEdx2 = 0.;
1789 }
1790 else //should only be 1
1791 {
1792 for(; locTofDedxIterator != locTofDedxList.end(); ++locTofDedxIterator)
1793 {
1794 locTOFHitMatchParams->dEdx1 = locTofDedxIterator->getDEdx1();
1795 locTOFHitMatchParams->dEdx2 = locTofDedxIterator->getDEdx2();
1796
1797 // check if already have average dE/dx
1798 if(locTOFHitMatchParams->dEdx > 0) continue;
1799
1800 // average dE/dx is missing: take average if hits in both planes, otherwise use single plane value
1801 if(locTOFHitMatchParams->dEdx1>0 && locTOFHitMatchParams->dEdx2>0)
1802 locTOFHitMatchParams->dEdx = (locTOFHitMatchParams->dEdx1 + locTOFHitMatchParams->dEdx2) / 2.0;
1803 else if(locTOFHitMatchParams->dEdx1>0)
1804 locTOFHitMatchParams->dEdx = locTOFHitMatchParams->dEdx1;
1805 else if(locTOFHitMatchParams->dEdx2>0)
1806 locTOFHitMatchParams->dEdx = locTOFHitMatchParams->dEdx2;
1807 }
1808 }
1809 }
1810
1811 // Extract track matching data for FMPWCs
1812 const hddm_r::FmwpcMatchParamsList &fmwpcList = iter->getFmwpcMatchParamses();
1813 hddm_r::FmwpcMatchParamsList::iterator fmwpcIter = fmwpcList.begin();
1814 for(; fmwpcIter != fmwpcList.end(); ++fmwpcIter)
1815 {
1816 size_t locTrackIndex = fmwpcIter->getTrack();
1817 const hddm_r::FmwpcDataList &fmwpcDataList = fmwpcIter->getFmwpcDatas();
1818 hddm_r::FmwpcDataList::iterator fmwpcDataIter = fmwpcDataList.begin();
1819
1820 auto locFMWPCMatchParams = std::make_shared<DFMWPCMatchParams>();
1821 for(; fmwpcDataIter != fmwpcDataList.end(); ++fmwpcDataIter){
1822 locFMWPCMatchParams->dLayers.push_back(fmwpcDataIter->getLayer());
1823 locFMWPCMatchParams->dNhits.push_back(fmwpcDataIter->getNhits());
1824 locFMWPCMatchParams->dDists.push_back(fmwpcDataIter->getDist());
1825 locFMWPCMatchParams->dClosestWires.push_back(fmwpcDataIter->getClosestwire());
1826 }
1827 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], std::const_pointer_cast<const DFMWPCMatchParams>(locFMWPCMatchParams));
1828 }
1829
1830 // Extract track matching data for CTOF
1831 const hddm_r::CtofMatchParamsList &ctofList = iter->getCtofMatchParamses();
1832 hddm_r::CtofMatchParamsList::iterator ctofIter = ctofList.begin();
1833 for(; ctofIter != ctofList.end(); ++ctofIter)
1834 {
1835 size_t locHitIndex = ctofIter->getHit();
1836 size_t locTrackIndex = ctofIter->getTrack();
1837
1838 auto locCTOFHitMatchParams = std::make_shared<DCTOFHitMatchParams>();
1839 locCTOFHitMatchParams->dCTOFPoint = locCTOFPoints[locHitIndex];
1840 locCTOFHitMatchParams->dEdx = ctofIter->getDEdx();
1841 locCTOFHitMatchParams->dFlightTime = ctofIter->getTflight();
1842 locCTOFHitMatchParams->dDeltaXToHit = ctofIter->getDeltax();
1843 locCTOFHitMatchParams->dDeltaYToHit = ctofIter->getDeltay();
1844
1845 locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locCTOFPoints[locHitIndex], std::const_pointer_cast<const DCTOFHitMatchParams>(locCTOFHitMatchParams));
1846 }
1847
1848 const hddm_r::BcalDOCAtoTrackList &bcaldocaList = iter->getBcalDOCAtoTracks();
1849 hddm_r::BcalDOCAtoTrackList::iterator bcaldocaIter = bcaldocaList.begin();
1850 for(; bcaldocaIter != bcaldocaList.end(); ++bcaldocaIter)
1851 {
1852 size_t locShowerIndex = bcaldocaIter->getShower();
1853 double locDeltaPhi = bcaldocaIter->getDeltaphi();
1854 double locDeltaZ = bcaldocaIter->getDeltaz();
1855 locDetectorMatches->Set_DistanceToNearestTrack(locBCALShowers[locShowerIndex], locDeltaPhi, locDeltaZ);
1856 }
1857
1858 const hddm_r::FcalDOCAtoTrackList &fcaldocaList = iter->getFcalDOCAtoTracks();
1859 hddm_r::FcalDOCAtoTrackList::iterator fcaldocaIter = fcaldocaList.begin();
1860 for(; fcaldocaIter != fcaldocaList.end(); ++fcaldocaIter)
1861 {
1862 size_t locShowerIndex = fcaldocaIter->getShower();
1863 double locDOCA = fcaldocaIter->getDoca();
1864 locDetectorMatches->Set_DistanceToNearestTrack(locFCALShowers[locShowerIndex], locDOCA);
1865 }
1866
1867 const hddm_r::TflightPCorrelationList &correlationList = iter->getTflightPCorrelations();
1868 hddm_r::TflightPCorrelationList::iterator correlationIter = correlationList.begin();
1869 for(; correlationIter != correlationList.end(); ++correlationIter)
1870 {
1871 size_t locTrackIndex = correlationIter->getTrack();
1872 DetectorSystem_t locDetectorSystem = (DetectorSystem_t)correlationIter->getSystem();
1873 double locCorrelation = correlationIter->getCorrelation();
1874 locDetectorMatches->Set_FlightTimePCorrelation(locTrackTimeBasedVector[locTrackIndex], locDetectorSystem, locCorrelation);
1875 }
1876
1877 data.push_back(locDetectorMatches);
1878 }
1879
1880 // Copy data to factory
1881 factory->CopyTo(data);
1882
1883 return NOERROR;
1884}
1885
1886// Transform the 5x5 tracking error matrix into a 7x7 error matrix in cartesian
1887// coordinates.
1888// This was copied and transformed from DKinFit.cc
1889void DEventSourceREST::Get7x7ErrorMatrix(double mass, const double vec[5], const TMatrixFSym* C5x5, TMatrixFSym* loc7x7ErrorMatrix)
1890{
1891 TMatrixF J(7,5);
1892
1893 // State vector
1894 double q_over_pt=vec[0];
1895 double phi=vec[1];
1896 double tanl=vec[2];
1897 double D=vec[3];
1898
1899 double pt=1./fabs(q_over_pt);
1900 double pt_sq=pt*pt;
1901 double cosphi=cos(phi);
1902 double sinphi=sin(phi);
1903 double q=(q_over_pt>0)?1.:-1.;
1904
1905 J(0, 0)=-q*pt_sq*cosphi;
1906 J(0, 1)=-pt*sinphi;
1907
1908 J(1, 0)=-q*pt_sq*sinphi;
1909 J(1, 1)=pt*cosphi;
1910
1911 J(2, 0)=-q*pt_sq*tanl;
1912 J(2, 2)=pt;
1913
1914 J(3, 1)=-D*cosphi;
1915 J(3, 3)=-sinphi;
1916
1917 J(4, 1)=-D*sinphi;
1918 J(4, 3)=cosphi;
1919
1920 J(5, 4)=1.;
1921
1922 // C'= JCJ^T
1923 TMatrixFSym locTempMatrix(*C5x5);
1924 *loc7x7ErrorMatrix=locTempMatrix.Similarity(J);
1925}
1926
1927uint32_t DEventSourceREST::Convert_SignedIntToUnsigned(int32_t locSignedInt) const
1928{
1929 //Convert uint32_t to int32_t
1930 //Reverse scheme (from DEventWriterREST):
1931 //If is >= 0, then the int32_t is the same as the uint32_t (last bit not set)
1932 //If is the minimum int: bit 32 is 1, and all of the other bits are zero
1933 //Else, bit 32 is 1, then the uint32_t is -1 * int32_t, + add the last bit
1934 if(locSignedInt >= 0)
1935 return uint32_t(locSignedInt); //bit 32 is zero
1936 else if(locSignedInt == numeric_limits<int32_t>::min()) // -(2^31)
1937 return uint32_t(0x80000000); //bit 32 is 1, all others are 0
1938 return uint32_t(-1*locSignedInt) + uint32_t(0x80000000); //bit 32 is 1, all others are negative of signed int (which was negative)
1939}
1940
1941//-----------------------
1942// Extract_DDIRCPmtHit
1943//-----------------------
1944jerror_t DEventSourceREST::Extract_DDIRCPmtHit(hddm_r::HDDM *record,
1945 JFactory<DDIRCPmtHit>* factory, JEventLoop* locEventLoop)
1946{
1947 /// Copies the data from the fcalShower hddm record. This is
1948 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1949 /// returns OBJECT_NOT_AVAILABLE immediately.
1950
1951 if (factory==NULL__null) {
1952 return OBJECT_NOT_AVAILABLE;
1953 }
1954 string tag = (factory->Tag())? factory->Tag() : "";
1955
1956 vector<DDIRCPmtHit*> data;
1957
1958 // loop over fcal shower records
1959 const hddm_r::DircHitList &hits =
1960 record->getDircHits();
1961 hddm_r::DircHitList::iterator iter;
1962 for (iter = hits.begin(); iter != hits.end(); ++iter) {
1963 if (iter->getJtag() != tag)
1964 continue;
1965
1966 // throw away hits from bad or noisy channels (after REST reconstruction)
1967 int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
1968 int box = (iter->getCh() < dDIRCMaxChannels) ? 1 : 0;
1969 int channel = iter->getCh() % dDIRCMaxChannels;
1970 dirc_status_state status = static_cast<dirc_status_state>(dDIRCChannelStatusMap[locRunNumber][box][channel]);
1971 if ( (status==BAD) || (status==NOISY) ) {
1972 continue;
1973 }
1974
1975 DDIRCPmtHit *hit = new DDIRCPmtHit();
1976 hit->setChannel(iter->getCh());
1977 hit->setTime(iter->getT());
1978 hit->setTOT(iter->getTot());
1979
1980 data.push_back(hit);
1981 }
1982
1983 // Copy into factory
1984 factory->CopyTo(data);
1985
1986 return NOERROR;
1987}
1988
1989//-----------------------
1990// Extract_DFMWPCHit
1991//-----------------------
1992jerror_t DEventSourceREST::Extract_DFMWPCHit(hddm_r::HDDM *record,
1993 JFactory<DFMWPCHit>* factory, JEventLoop* locEventLoop)
1994{
1995 /// Copies the data from the fmwpc hit hddm record. This is
1996 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
1997 /// returns OBJECT_NOT_AVAILABLE immediately.
1998
1999 if (factory==NULL__null) {
2000 return OBJECT_NOT_AVAILABLE;
2001 }
2002 string tag = (factory->Tag())? factory->Tag() : "";
2003
2004 vector<DFMWPCHit*> data;
2005
2006 // loop over fmwpc hit records
2007 const hddm_r::FmwpcHitList &hits =
2008 record->getFmwpcHits();
2009 hddm_r::FmwpcHitList::iterator iter;
2010 for (iter = hits.begin(); iter != hits.end(); ++iter) {
2011 if (iter->getJtag() != tag)
2012 continue;
2013
2014 DFMWPCHit *hit = new DFMWPCHit();
2015 hit->layer = iter->getLayer();
2016 hit->wire = iter->getWire();
2017 hit->q = iter->getQ();
2018 hit->amp = iter->getAmp();
2019 hit->t = iter->getT();
2020 hit->QF = iter->getQf();
2021 hit->ped = iter->getPed();
2022
2023 data.push_back(hit);
2024 }
2025
2026 // Copy into factory
2027 factory->CopyTo(data);
2028
2029 return NOERROR;
2030}
2031
2032//----------------------------
2033// Extract_DEventHitStatistics
2034//----------------------------
2035jerror_t DEventSourceREST::Extract_DEventHitStatistics(hddm_r::HDDM *record,
2036 JFactory<DEventHitStatistics>* factory)
2037{
2038 /// Copies the data from the hitStatistics hddm record. This is
2039 /// call from JEventSourceREST::GetObjects. If factory is NULL, this
2040 /// returns OBJECT_NOT_AVAILABLE immediately.
2041
2042 if (factory==NULL__null) {
2043 return OBJECT_NOT_AVAILABLE;
2044 }
2045 string tag = (factory->Tag())? factory->Tag() : "";
2046
2047 vector<DEventHitStatistics*> data;
2048
2049 hddm_r::HitStatisticsList slist = (hddm_r::HitStatisticsList)record->getHitStatisticses();
2050 if (slist.size() != 0 && slist().getJtag() == tag) {
2051 DEventHitStatistics *stats = new DEventHitStatistics;
2052 hddm_r::StartCountersList starts = slist().getStartCounterses();
2053 stats->start_counters = (starts.size() > 0)? starts().getCount() : 0;
2054 hddm_r::CdcStrawsList straws = slist().getCdcStrawses();
2055 stats->cdc_straws = (straws.size() > 0)? straws().getCount() : 0;
2056 hddm_r::FdcPseudosList pseudos = slist().getFdcPseudoses();
2057 stats->fdc_pseudos = (pseudos.size() > 0)? pseudos().getCount() : 0;
2058 hddm_r::BcalCellsList cells = slist().getBcalCellses();
2059 stats->bcal_cells = (cells.size() > 0)? cells().getCount() : 0;
2060 hddm_r::FcalBlocksList blocks = slist().getFcalBlockses();
2061 stats->fcal_blocks = (blocks.size() > 0)? blocks().getCount() : 0;
2062 hddm_r::CcalBlocksList bloccs = slist().getCcalBlockses();
2063 stats->ccal_blocks = (bloccs.size() > 0)? bloccs().getCount() : 0;
2064 hddm_r::TofPaddlesList paddles = slist().getTofPaddleses();
2065 stats->tof_paddles = (paddles.size() > 0)? paddles().getCount() : 0;
2066 hddm_r::DircPMTsList pmts = slist().getDircPMTses();
2067 stats->dirc_PMTs = (pmts.size() > 0)? pmts().getCount() : 0;
2068
2069 data.push_back(stats);
2070 }
2071
2072 // Copy into factory
2073 factory->CopyTo(data);
2074
2075 return NOERROR;
2076}