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