Bug Summary

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