File indexing completed on 2022-02-21 23:14:26
0001 #include "Event.h"
0002 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0003
0004
0005 #include "RecoTracker/MkFitCore/src/Debug.h"
0006
0007 #ifdef TBB
0008 #include "oneapi/tbb/parallel_for.h"
0009 #endif
0010
0011 #include <memory>
0012
0013 namespace {
0014 std::unique_ptr<mkfit::Validation> dummyValidation(mkfit::Validation::make_validation("dummy", nullptr));
0015 }
0016
0017 namespace mkfit {
0018
0019 std::mutex Event::printmutex;
0020
0021 Event::Event(int evtID, int nLayers) : validation_(*dummyValidation), evtID_(evtID) {
0022 layerHits_.resize(nLayers);
0023 layerHitMasks_.resize(nLayers);
0024 }
0025
0026 Event::Event(Validation &v, int evtID, int nLayers) : validation_(v), evtID_(evtID) {
0027 layerHits_.resize(nLayers);
0028 layerHitMasks_.resize(nLayers);
0029 validation_.resetValidationMaps();
0030 }
0031
0032 void Event::reset(int evtID) {
0033 evtID_ = evtID;
0034
0035 for (auto &&l : layerHits_) {
0036 l.clear();
0037 }
0038 for (auto &&l : layerHitMasks_) {
0039 l.clear();
0040 }
0041
0042 simHitsInfo_.clear();
0043 simTrackStates_.clear();
0044 simTracks_.clear();
0045 simTracksExtra_.clear();
0046 seedTracks_.clear();
0047 seedTracksExtra_.clear();
0048 candidateTracks_.clear();
0049 candidateTracksExtra_.clear();
0050 fitTracks_.clear();
0051 fitTracksExtra_.clear();
0052 cmsswTracks_.clear();
0053 cmsswTracksExtra_.clear();
0054 beamSpot_ = {};
0055
0056 validation_.resetValidationMaps();
0057 }
0058
0059 void Event::validate() {
0060
0061 if (Config::sim_val_for_cmssw) {
0062 validation_.makeRecoTkToSeedTkMapsDumbCMSSW(*this);
0063 validation_.setTrackScoresDumbCMSSW(*this);
0064 }
0065
0066
0067 if (Config::sim_val || Config::sim_val_for_cmssw) {
0068 validation_.setTrackExtras(*this);
0069 validation_.makeSimTkToRecoTksMaps(*this);
0070 validation_.makeSeedTkToRecoTkMaps(*this);
0071 validation_.fillEfficiencyTree(*this);
0072 validation_.fillFakeRateTree(*this);
0073 }
0074
0075
0076 if (Config::cmssw_val) {
0077 validation_.makeCMSSWTkToSeedTkMap(*this);
0078 validation_.makeRecoTkToRecoTkMaps(*this);
0079 validation_.setTrackExtras(*this);
0080 validation_.makeCMSSWTkToRecoTksMaps(*this);
0081 validation_.fillCMSSWEfficiencyTree(*this);
0082 validation_.fillCMSSWFakeRateTree(*this);
0083 }
0084
0085 if (Config::fit_val) {
0086 validation_.fillFitTree(*this);
0087 }
0088 }
0089
0090 void Event::printStats(const TrackVec &trks, TrackExtraVec &trkextras) {
0091 int miss(0), found(0), fp_10(0), fp_20(0), hit8(0), h8_10(0), h8_20(0);
0092
0093 for (auto &&trk : trks) {
0094 auto &&extra = trkextras[trk.label()];
0095 extra.setMCTrackIDInfo(trk, layerHits_, simHitsInfo_, simTracks_, false, true);
0096 if (extra.mcTrackID() < 0) {
0097 ++miss;
0098 } else {
0099 auto &&mctrk = simTracks_[extra.mcTrackID()];
0100 auto pr = trk.pT() / mctrk.pT();
0101 found++;
0102 bool h8 = trk.nFoundHits() >= 8;
0103 bool pt10 = pr > 0.9 && pr < 1.1;
0104 bool pt20 = pr > 0.8 && pr < 1.2;
0105 fp_10 += pt10;
0106 fp_20 += pt20;
0107 hit8 += h8;
0108 h8_10 += h8 && pt10;
0109 h8_20 += h8 && pt20;
0110 }
0111 }
0112 std::cout << "found tracks=" << found << " in pT 10%=" << fp_10 << " in pT 20%=" << fp_20
0113 << " no_mc_assoc=" << miss << std::endl
0114 << " nH >= 8 =" << hit8 << " in pT 10%=" << h8_10 << " in pT 20%=" << h8_20 << std::endl;
0115 }
0116
0117 void Event::write_out(DataFile &data_file) {
0118 FILE *fp = data_file.f_fp;
0119
0120 static std::mutex writemutex;
0121 std::lock_guard<std::mutex> writelock(writemutex);
0122
0123 auto start = ftell(fp);
0124 int evsize = sizeof(int);
0125 fwrite(&evsize, sizeof(int), 1, fp);
0126
0127 evsize += write_tracks(fp, simTracks_);
0128
0129 if (data_file.hasSimTrackStates()) {
0130 int nts = simTrackStates_.size();
0131 fwrite(&nts, sizeof(int), 1, fp);
0132 fwrite(&simTrackStates_[0], sizeof(TrackState), nts, fp);
0133 evsize += sizeof(int) + nts * sizeof(TrackState);
0134 }
0135
0136 int nl = layerHits_.size();
0137 fwrite(&nl, sizeof(int), 1, fp);
0138 evsize += sizeof(int);
0139 for (int il = 0; il < nl; ++il) {
0140 int nh = layerHits_[il].size();
0141 fwrite(&nh, sizeof(int), 1, fp);
0142 fwrite(&layerHits_[il][0], sizeof(Hit), nh, fp);
0143 evsize += sizeof(int) + nh * sizeof(Hit);
0144 }
0145
0146 if (data_file.hasHitIterMasks()) {
0147
0148 for (int il = 0; il < nl; ++il) {
0149 int nh = layerHitMasks_[il].size();
0150 assert(nh == (int)layerHits_[il].size());
0151 fwrite(&layerHitMasks_[il][0], sizeof(uint64_t), nh, fp);
0152 evsize += nh * sizeof(uint64_t);
0153 }
0154 }
0155
0156 int nm = simHitsInfo_.size();
0157 fwrite(&nm, sizeof(int), 1, fp);
0158 fwrite(&simHitsInfo_[0], sizeof(MCHitInfo), nm, fp);
0159 evsize += sizeof(int) + nm * sizeof(MCHitInfo);
0160
0161 if (data_file.hasSeeds()) {
0162 evsize += write_tracks(fp, seedTracks_);
0163 }
0164
0165 if (data_file.hasCmsswTracks()) {
0166 evsize += write_tracks(fp, cmsswTracks_);
0167 }
0168
0169 if (data_file.hasBeamSpot()) {
0170 fwrite(&beamSpot_, sizeof(BeamSpot), 1, fp);
0171 evsize += sizeof(BeamSpot);
0172 }
0173
0174 fseek(fp, start, SEEK_SET);
0175 fwrite(&evsize, sizeof(int), 1, fp);
0176 fseek(fp, 0, SEEK_END);
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196 }
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206 void Event::read_in(DataFile &data_file, FILE *in_fp) {
0207 FILE *fp = in_fp ? in_fp : data_file.f_fp;
0208
0209 data_file.advancePosToNextEvent(fp);
0210
0211 int nt = read_tracks(fp, simTracks_);
0212 Config::nTracks = nt;
0213
0214 if (data_file.hasSimTrackStates()) {
0215 int nts;
0216 fread(&nts, sizeof(int), 1, fp);
0217 simTrackStates_.resize(nts);
0218 fread(&simTrackStates_[0], sizeof(TrackState), nts, fp);
0219 }
0220
0221 int nl;
0222 fread(&nl, sizeof(int), 1, fp);
0223 layerHits_.resize(nl);
0224 layerHitMasks_.resize(nl);
0225 for (int il = 0; il < nl; ++il) {
0226 int nh;
0227 fread(&nh, sizeof(int), 1, fp);
0228 layerHits_[il].resize(nh);
0229 layerHitMasks_[il].resize(nh, 0);
0230 fread(&layerHits_[il][0], sizeof(Hit), nh, fp);
0231 }
0232
0233 if (data_file.hasHitIterMasks()) {
0234 for (int il = 0; il < nl; ++il) {
0235 int nh = layerHits_[il].size();
0236 fread(&layerHitMasks_[il][0], sizeof(uint64_t), nh, fp);
0237 }
0238 }
0239
0240 int nm;
0241 fread(&nm, sizeof(int), 1, fp);
0242 simHitsInfo_.resize(nm);
0243 fread(&simHitsInfo_[0], sizeof(MCHitInfo), nm, fp);
0244
0245 if (data_file.hasSeeds()) {
0246 int ns = read_tracks(fp, seedTracks_, Config::seedInput != cmsswSeeds);
0247 (void)ns;
0248
0249 #ifdef DUMP_SEEDS
0250 printf("Read %i seedtracks (neg value means actual reading was skipped)\n", ns);
0251 for (int it = 0; it < ns; it++) {
0252 const Track &ss = seedTracks_[it];
0253 printf(" %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%i label=%4i algo=%2i\n",
0254 it,
0255 ss.charge(),
0256 ss.pT(),
0257 ss.momEta(),
0258 ss.nFoundHits(),
0259 ss.label(),
0260 (int)ss.algorithm());
0261 #ifdef DUMP_SEED_HITS
0262 for (int ih = 0; ih < seedTracks_[it].nTotalHits(); ++ih) {
0263 int lyr = seedTracks_[it].getHitLyr(ih);
0264 int idx = seedTracks_[it].getHitIdx(ih);
0265 if (idx >= 0) {
0266 const Hit &hit = layerHits_[lyr][idx];
0267 printf(" hit %2d lyr=%3d idx=%4d pos r=%7.3f z=% 8.3f mc_hit=%3d mc_trk=%3d\n",
0268 ih,
0269 lyr,
0270 idx,
0271 layerHits_[lyr][idx].r(),
0272 layerHits_[lyr][idx].z(),
0273 hit.mcHitID(),
0274 hit.mcTrackID(simHitsInfo_));
0275 } else
0276 printf(" hit %2d idx=%i\n", ih, seedTracks_[it].getHitIdx(ih));
0277 }
0278 #endif
0279 }
0280 #endif
0281 }
0282
0283 int nert = -99999;
0284 if (data_file.hasCmsswTracks()) {
0285 nert = read_tracks(fp, cmsswTracks_, !Config::readCmsswTracks);
0286 (void)nert;
0287 }
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301 #ifdef DUMP_TRACKS
0302 printf("Read %i simtracks\n", nt);
0303 for (int it = 0; it < nt; it++) {
0304 const Track &t = simTracks_[it];
0305 printf(" %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%2d label=%4d\n",
0306 it,
0307 t.charge(),
0308 t.pT(),
0309 t.momEta(),
0310 t.nFoundHits(),
0311 t.label());
0312 #ifdef DUMP_TRACK_HITS
0313 for (int ih = 0; ih < t.nTotalHits(); ++ih) {
0314 int lyr = t.getHitLyr(ih);
0315 int idx = t.getHitIdx(ih);
0316 if (idx >= 0) {
0317 const Hit &hit = layerHits_[lyr][idx];
0318 printf(" hit %2d lyr=%2d idx=%3d pos r=%7.3f x=% 8.3f y=% 8.3f z=% 8.3f mc_hit=%3d mc_trk=%3d\n",
0319 ih,
0320 lyr,
0321 idx,
0322 layerHits_[lyr][idx].r(),
0323 layerHits_[lyr][idx].x(),
0324 layerHits_[lyr][idx].y(),
0325 layerHits_[lyr][idx].z(),
0326 hit.mcHitID(),
0327 hit.mcTrackID(simHitsInfo_));
0328 } else
0329 printf(" hit %2d idx=%i\n", ih, t.getHitIdx(ih));
0330 }
0331 #endif
0332 }
0333 #endif
0334 #ifdef DUMP_LAYER_HITS
0335 printf("Read %i layers\n", nl);
0336 int total_hits = 0;
0337 for (int il = 0; il < nl; il++) {
0338 if (layerHits_[il].empty())
0339 continue;
0340
0341 printf("Read %i hits in layer %i\n", (int)layerHits_[il].size(), il);
0342 total_hits += layerHits_[il].size();
0343 for (int ih = 0; ih < (int)layerHits_[il].size(); ih++) {
0344 const Hit &hit = layerHits_[il][ih];
0345 printf(" mcHitID=%5d r=%10g x=%10g y=%10g z=%10g sx=%10.4g sy=%10.4e sz=%10.4e\n",
0346 hit.mcHitID(),
0347 hit.r(),
0348 hit.x(),
0349 hit.y(),
0350 hit.z(),
0351 std::sqrt(hit.exx()),
0352 std::sqrt(hit.eyy()),
0353 std::sqrt(hit.ezz()));
0354 }
0355 }
0356 printf("Total hits in all layers = %d\n", total_hits);
0357 #endif
0358 #ifdef DUMP_REC_TRACKS
0359 printf("Read %i rectracks\n", nert);
0360 for (int it = 0; it < nert; it++) {
0361 const Track &t = cmsswTracks_[it];
0362 printf(" %i with q=%+i pT=%7.3f eta=% 7.3f nHits=%2d label=%4d algo=%2d\n",
0363 it,
0364 t.charge(),
0365 t.pT(),
0366 t.momEta(),
0367 t.nFoundHits(),
0368 t.label(),
0369 (int)t.algorithm());
0370 #ifdef DUMP_REC_TRACK_HITS
0371 for (int ih = 0; ih < t.nTotalHits(); ++ih) {
0372 int lyr = t.getHitLyr(ih);
0373 int idx = t.getHitIdx(ih);
0374 if (idx >= 0) {
0375 const Hit &hit = layerHits_[lyr][idx];
0376 printf(" hit %2d lyr=%2d idx=%3d pos r=%7.3f z=% 8.3f mc_hit=%3d mc_trk=%3d\n",
0377 ih,
0378 lyr,
0379 idx,
0380 hit.r(),
0381 hit.z(),
0382 hit.mcHitID(),
0383 hit.mcTrackID(simHitsInfo_));
0384 } else
0385 printf(" hit %2d idx=%i\n", ih, t.getHitIdx(ih));
0386 }
0387 #endif
0388 }
0389 #endif
0390
0391 if (data_file.hasBeamSpot()) {
0392 fread(&beamSpot_, sizeof(BeamSpot), 1, fp);
0393 }
0394
0395 if (Config::kludgeCmsHitErrors) {
0396 kludge_cms_hit_errors();
0397 }
0398
0399 if (!Config::silent)
0400 printf("Read complete, %d simtracks on file.\n", nt);
0401 }
0402
0403
0404
0405 int Event::write_tracks(FILE *fp, const TrackVec &tracks) {
0406
0407
0408 int n_tracks = tracks.size();
0409 fwrite(&n_tracks, sizeof(int), 1, fp);
0410
0411 auto start = ftell(fp);
0412 int data_size = 2 * sizeof(int) + n_tracks * sizeof(Track);
0413 fwrite(&data_size, sizeof(int), 1, fp);
0414
0415 fwrite(tracks.data(), sizeof(Track), n_tracks, fp);
0416
0417 for (int i = 0; i < n_tracks; ++i) {
0418 fwrite(tracks[i].beginHitsOnTrack(), sizeof(HitOnTrack), tracks[i].nTotalHits(), fp);
0419 data_size += tracks[i].nTotalHits() * sizeof(HitOnTrack);
0420 }
0421
0422 fseek(fp, start, SEEK_SET);
0423 fwrite(&data_size, sizeof(int), 1, fp);
0424 fseek(fp, 0, SEEK_END);
0425
0426 return data_size;
0427 }
0428
0429 int Event::read_tracks(FILE *fp, TrackVec &tracks, bool skip_reading) {
0430
0431
0432 int n_tracks, data_size;
0433 fread(&n_tracks, sizeof(int), 1, fp);
0434 fread(&data_size, sizeof(int), 1, fp);
0435
0436 if (skip_reading) {
0437 fseek(fp, data_size - 2 * sizeof(int), SEEK_CUR);
0438 n_tracks = -n_tracks;
0439 } else {
0440 tracks.resize(n_tracks);
0441
0442 fread(tracks.data(), sizeof(Track), n_tracks, fp);
0443
0444 for (int i = 0; i < n_tracks; ++i) {
0445 tracks[i].resizeHitsForInput();
0446 fread(tracks[i].beginHitsOnTrack_nc(), sizeof(HitOnTrack), tracks[i].nTotalHits(), fp);
0447 }
0448 }
0449
0450 return n_tracks;
0451 }
0452
0453
0454
0455 void Event::setInputFromCMSSW(std::vector<HitVec> hits, TrackVec seeds) {
0456 layerHits_ = std::move(hits);
0457 seedTracks_ = std::move(seeds);
0458 }
0459
0460
0461
0462 void Event::kludge_cms_hit_errors() {
0463
0464
0465 const float Exy = 15 * 1e-4, Vxy = Exy * Exy;
0466 const float Ez = 30 * 1e-4, Vz = Ez * Ez;
0467
0468 int nl = layerHits_.size();
0469
0470 int cnt = 0;
0471
0472 for (int il = 0; il < nl; il++) {
0473 if (layerHits_[il].empty())
0474 continue;
0475
0476 for (Hit &h : layerHits_[il]) {
0477 SVector6 &c = h.error_nc();
0478
0479 float vxy = c[0] + c[2];
0480 if (vxy < Vxy) {
0481 c[0] *= Vxy / vxy;
0482 c[2] *= Vxy / vxy;
0483 ++cnt;
0484 }
0485 if (il < 4 && c[5] < Vz) {
0486 c[5] = Vz;
0487 ++cnt;
0488 }
0489 }
0490 }
0491
0492 printf("Event::kludge_cms_hit_errors processed %d layers, kludged %d entries.\n", nl, cnt);
0493 }
0494
0495
0496
0497 int Event::clean_cms_simtracks() {
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512 dprintf("Event::clean_cms_simtracks processing %lu simtracks.\n", simTracks_.size());
0513
0514 int n_acc = 0;
0515 int i = -1;
0516 for (Track &t : simTracks_) {
0517 i++;
0518
0519 t.sortHitsByLayer();
0520
0521 const int lyr_cnt = t.nUniqueLayers();
0522
0523
0524
0525
0526 if (lyr_cnt < Config::cmsSelMinLayers)
0527 {
0528 dprintf("Rejecting simtrack %d, n_hits=%d, n_layers=%d, pT=%f\n", i, t.nFoundHits(), lyr_cnt, t.pT());
0529 t.setNotFindable();
0530 } else {
0531 dprintf("Accepting simtrack %d, n_hits=%d, n_layers=%d, pT=%f\n", i, t.nFoundHits(), lyr_cnt, t.pT());
0532 ++n_acc;
0533 }
0534 }
0535
0536 return n_acc;
0537 }
0538
0539 void Event::print_tracks(const TrackVec &tracks, bool print_hits) const {
0540 const int nt = tracks.size();
0541
0542
0543 printf("Event::print_tracks printing %d tracks %s hits:\n", nt, (print_hits ? "with" : "without"));
0544 for (int it = 0; it < nt; it++) {
0545 const Track &t = tracks[it];
0546 printf(" %i with q=%+i pT=%7.3f eta=% 7.3f nHits=%2d label=%4d findable=%d score=%7.3f chi2=%7.3f\n",
0547 it,
0548 t.charge(),
0549 t.pT(),
0550 t.momEta(),
0551 t.nFoundHits(),
0552 t.label(),
0553 t.isFindable(),
0554 getScoreCand(t),
0555 t.chi2());
0556
0557 if (print_hits) {
0558 for (int ih = 0; ih < t.nTotalHits(); ++ih) {
0559 int lyr = t.getHitLyr(ih);
0560 int idx = t.getHitIdx(ih);
0561 if (idx >= 0) {
0562 const Hit &hit = layerHits_[lyr][idx];
0563 printf(" hit %2d lyr=%2d idx=%3d pos r=%7.3f z=% 8.3f mc_hit=%3d mc_trk=%3d\n",
0564 ih,
0565 lyr,
0566 idx,
0567 layerHits_[lyr][idx].r(),
0568 layerHits_[lyr][idx].z(),
0569 hit.mcHitID(),
0570 hit.mcTrackID(simHitsInfo_));
0571 } else
0572 printf(" hit %2d lyr=%2d idx=%3d\n", ih, t.getHitLyr(ih), t.getHitIdx(ih));
0573 }
0574 }
0575 }
0576 }
0577
0578 int Event::clean_cms_seedtracks(TrackVec *seed_ptr) {
0579 const float etamax_brl = Config::c_etamax_brl;
0580 const float dpt_common = Config::c_dpt_common;
0581 const float dzmax_brl = Config::c_dzmax_brl;
0582 const float drmax_brl = Config::c_drmax_brl;
0583 const float ptmin_hpt = Config::c_ptmin_hpt;
0584 const float dzmax_hpt = Config::c_dzmax_hpt;
0585 const float drmax_hpt = Config::c_drmax_hpt;
0586 const float dzmax_els = Config::c_dzmax_els;
0587 const float drmax_els = Config::c_drmax_els;
0588
0589 const float dzmax2_inv_brl = 1.f / (dzmax_brl * dzmax_brl);
0590 const float drmax2_inv_brl = 1.f / (drmax_brl * drmax_brl);
0591 const float dzmax2_inv_hpt = 1.f / (dzmax_hpt * dzmax_hpt);
0592 const float drmax2_inv_hpt = 1.f / (drmax_hpt * drmax_hpt);
0593 const float dzmax2_inv_els = 1.f / (dzmax_els * dzmax_els);
0594 const float drmax2_inv_els = 1.f / (drmax_els * drmax_els);
0595
0596 TrackVec &seeds = (seed_ptr != nullptr) ? *seed_ptr : seedTracks_;
0597 const int ns = seeds.size();
0598
0599 TrackVec cleanSeedTracks;
0600 cleanSeedTracks.reserve(ns);
0601 std::vector<bool> writetrack(ns, true);
0602
0603 const float invR1GeV = 1.f / Config::track1GeVradius;
0604
0605 std::vector<int> nHits(ns);
0606 std::vector<int> charge(ns);
0607 std::vector<float> oldPhi(ns);
0608 std::vector<float> pos2(ns);
0609 std::vector<float> eta(ns);
0610 std::vector<float> ctheta(ns);
0611 std::vector<float> invptq(ns);
0612 std::vector<float> pt(ns);
0613 std::vector<float> x(ns);
0614 std::vector<float> y(ns);
0615 std::vector<float> z(ns);
0616
0617 for (int ts = 0; ts < ns; ts++) {
0618 const Track &tk = seeds[ts];
0619 nHits[ts] = tk.nFoundHits();
0620 charge[ts] = tk.charge();
0621 oldPhi[ts] = tk.momPhi();
0622 pos2[ts] = std::pow(tk.x(), 2) + std::pow(tk.y(), 2);
0623 eta[ts] = tk.momEta();
0624 ctheta[ts] = 1.f / std::tan(tk.theta());
0625 invptq[ts] = tk.charge() * tk.invpT();
0626 pt[ts] = tk.pT();
0627 x[ts] = tk.x();
0628 y[ts] = tk.y();
0629 z[ts] = tk.z();
0630 }
0631
0632 for (int ts = 0; ts < ns; ts++) {
0633 if (not writetrack[ts])
0634 continue;
0635
0636 const float oldPhi1 = oldPhi[ts];
0637 const float pos2_first = pos2[ts];
0638 const float Eta1 = eta[ts];
0639 const float Pt1 = pt[ts];
0640 const float invptq_first = invptq[ts];
0641
0642
0643 for (int tss = ts + 1; tss < ns; tss++) {
0644 const float Pt2 = pt[tss];
0645
0646
0647 if (charge[tss] != charge[ts])
0648 continue;
0649
0650 const float thisDPt = std::abs(Pt2 - Pt1);
0651
0652 if (thisDPt > dpt_common * (Pt1))
0653 continue;
0654
0655 const float Eta2 = eta[tss];
0656 const float deta2 = std::pow(Eta1 - Eta2, 2);
0657
0658 const float oldPhi2 = oldPhi[tss];
0659
0660 const float pos2_second = pos2[tss];
0661 const float thisDXYSign05 = pos2_second > pos2_first ? -0.5f : 0.5f;
0662
0663 const float thisDXY = thisDXYSign05 * sqrt(std::pow(x[ts] - x[tss], 2) + std::pow(y[ts] - y[tss], 2));
0664
0665 const float invptq_second = invptq[tss];
0666
0667 const float newPhi1 = oldPhi1 - thisDXY * invR1GeV * invptq_first;
0668 const float newPhi2 = oldPhi2 + thisDXY * invR1GeV * invptq_second;
0669
0670 const float dphi = cdist(std::abs(newPhi1 - newPhi2));
0671
0672 const float dr2 = deta2 + dphi * dphi;
0673
0674 const float thisDZ = z[ts] - z[tss] - thisDXY * (ctheta[ts] + ctheta[tss]);
0675 const float dz2 = thisDZ * thisDZ;
0676
0677
0678
0679 if (std::abs(Eta1) < etamax_brl) {
0680 if (dz2 * dzmax2_inv_brl + dr2 * drmax2_inv_brl < 1.0f)
0681 writetrack[tss] = false;
0682 } else if (Pt1 > ptmin_hpt) {
0683 if (dz2 * dzmax2_inv_hpt + dr2 * drmax2_inv_hpt < 1.0f)
0684 writetrack[tss] = false;
0685 } else {
0686 if (dz2 * dzmax2_inv_els + dr2 * drmax2_inv_els < 1.0f)
0687 writetrack[tss] = false;
0688 }
0689 }
0690
0691 if (writetrack[ts])
0692 cleanSeedTracks.emplace_back(seeds[ts]);
0693 }
0694
0695 seeds.swap(cleanSeedTracks);
0696
0697 #ifdef DEBUG
0698 {
0699 const int ns2 = seeds.size();
0700 printf("Number of CMS seeds before %d --> after %d cleaning\n", ns, ns2);
0701
0702 for (int it = 0; it < ns2; it++) {
0703 const Track &ss = seeds[it];
0704 printf(" %3i q=%+i pT=%7.3f eta=% 7.3f nHits=%i label=% i\n",
0705 it,
0706 ss.charge(),
0707 ss.pT(),
0708 ss.momEta(),
0709 ss.nFoundHits(),
0710 ss.label());
0711 }
0712 }
0713 #endif
0714
0715 return seeds.size();
0716 }
0717
0718 int Event::select_tracks_iter(unsigned int n) {
0719 if (n == 0)
0720 return 1;
0721
0722 unsigned int algorithms[] = {4, 22, 23, 5, 24, 7, 8, 9, 10, 6};
0723
0724
0725 const int ns = seedTracks_.size();
0726
0727 TrackVec cleanSeedTracks;
0728 cleanSeedTracks.reserve(ns);
0729
0730 for (int ts = 0; ts < ns; ts++) {
0731 const Track &tk = seedTracks_[ts];
0732 unsigned int algo = (unsigned int)tk.algorithm();
0733 if (std::find(algorithms, algorithms + n, algo) != algorithms + n)
0734 cleanSeedTracks.emplace_back(seedTracks_[ts]);
0735 }
0736 seedTracks_.swap(cleanSeedTracks);
0737
0738
0739 const int nt = cmsswTracks_.size();
0740
0741 TrackVec cleanTracks;
0742 cleanTracks.reserve(nt);
0743
0744 for (int ts = 0; ts < nt; ts++) {
0745 const Track &tk = cmsswTracks_[ts];
0746 unsigned int algo = (unsigned int)tk.algorithm();
0747 if (std::find(algorithms, algorithms + n, algo) != algorithms + n)
0748 cleanTracks.emplace_back(cmsswTracks_[ts]);
0749 }
0750 cmsswTracks_.swap(cleanTracks);
0751 return cmsswTracks_.size() + seedTracks_.size();
0752 }
0753
0754 int Event::clean_cms_seedtracks_badlabel() {
0755 printf("***\n*** REMOVING SEEDS WITH BAD LABEL. This is a development hack. ***\n***\n");
0756 TrackVec buf;
0757 seedTracks_.swap(buf);
0758 std::copy_if(
0759 buf.begin(), buf.end(), std::back_inserter(seedTracks_), [](const Track &t) { return t.label() >= 0; });
0760 return seedTracks_.size();
0761 }
0762
0763 int Event::use_seeds_from_cmsswtracks() {
0764 int ns = seedTracks_.size();
0765
0766 TrackVec cleanSeedTracks;
0767 cleanSeedTracks.reserve(ns);
0768
0769 for (auto &&cmsswtrack : cmsswTracks_) {
0770 cleanSeedTracks.emplace_back(seedTracks_[cmsswtrack.label()]);
0771 }
0772
0773 seedTracks_.swap(cleanSeedTracks);
0774
0775 return seedTracks_.size();
0776 }
0777
0778 void Event::relabel_bad_seedtracks() {
0779 int newlabel = 0;
0780 for (auto &&track : seedTracks_) {
0781 if (track.label() < 0)
0782 track.setLabel(--newlabel);
0783 }
0784 }
0785
0786 void Event::relabel_cmsswtracks_from_seeds() {
0787 std::map<int, int> cmsswLabelMap;
0788 for (size_t iseed = 0; iseed < seedTracks_.size(); iseed++) {
0789 for (size_t icmssw = 0; icmssw < cmsswTracks_.size(); icmssw++) {
0790 if (cmsswTracks_[icmssw].label() == static_cast<int>(iseed)) {
0791 cmsswLabelMap[icmssw] = seedTracks_[iseed].label();
0792 break;
0793 }
0794 }
0795 }
0796 for (size_t icmssw = 0; icmssw < cmsswTracks_.size(); icmssw++) {
0797 cmsswTracks_[icmssw].setLabel(cmsswLabelMap[icmssw]);
0798 }
0799 }
0800
0801
0802
0803
0804
0805 void Event::fill_hitmask_bool_vectors(int track_algo, std::vector<std::vector<bool>> &layer_masks) {
0806
0807
0808
0809 uint64_t iter_mask = 1 << track_algo;
0810
0811 const int n_lay = (int)layerHits_.size();
0812 layer_masks.resize(n_lay);
0813
0814 for (int l = 0; l < n_lay; ++l) {
0815 const int n_hit = (int)layerHits_[l].size();
0816 layer_masks[l].resize(n_hit);
0817
0818 for (int i = 0; i < n_hit; ++i) {
0819 layer_masks[l][i] = layerHitMasks_[l][i] & iter_mask;
0820 }
0821 }
0822 }
0823
0824 void Event::fill_hitmask_bool_vectors(std::vector<int> &track_algo_vec, std::vector<std::vector<bool>> &layer_masks) {
0825
0826
0827
0828
0829 uint64_t iter_mask = 0;
0830 for (auto ta : track_algo_vec)
0831 iter_mask |= 1 << ta;
0832
0833 const int n_lay = (int)layerHits_.size();
0834 layer_masks.resize(n_lay);
0835
0836 for (int l = 0; l < n_lay; ++l) {
0837 const int n_hit = (int)layerHits_[l].size();
0838 layer_masks[l].resize(n_hit);
0839
0840 for (int i = 0; i < n_hit; ++i) {
0841 uint64_t hitmasks = layerHitMasks_[l][i];
0842 layer_masks[l][i] = ((iter_mask ^ hitmasks) & iter_mask) == 0;
0843 }
0844 }
0845 }
0846
0847
0848
0849
0850
0851 int DataFile::openRead(const std::string &fname, int expected_n_layers) {
0852 constexpr int min_ver = 4;
0853 constexpr int max_ver = 6;
0854
0855 f_fp = fopen(fname.c_str(), "r");
0856 assert(f_fp != 0 && "Opening of input file failed.");
0857
0858 fread(&f_header, sizeof(DataFileHeader), 1, f_fp);
0859
0860 if (f_header.f_magic != 0xBEEF) {
0861 fprintf(stderr, "Incompatible input file (wrong magick).\n");
0862 exit(1);
0863 }
0864 if (f_header.f_format_version < min_ver || f_header.f_format_version > max_ver) {
0865 fprintf(stderr,
0866 "Unsupported file version %d. Supported versions are from %d to %d.\n",
0867 f_header.f_format_version,
0868 min_ver,
0869 max_ver);
0870 exit(1);
0871 }
0872 if (f_header.f_sizeof_track != sizeof(Track)) {
0873 fprintf(stderr,
0874 "sizeof(Track) on file (%d) different from current value (%d).\n",
0875 f_header.f_sizeof_track,
0876 (int)sizeof(Track));
0877 exit(1);
0878 }
0879 if (f_header.f_sizeof_hit != sizeof(Hit)) {
0880 fprintf(stderr,
0881 "sizeof(Hit) on file (%d) different from current value (%d).\n",
0882 f_header.f_sizeof_hit,
0883 (int)sizeof(Hit));
0884 exit(1);
0885 }
0886 if (f_header.f_sizeof_hot != sizeof(HitOnTrack)) {
0887 fprintf(stderr,
0888 "sizeof(HitOnTrack) on file (%d) different from current value (%d).\n",
0889 f_header.f_sizeof_hot,
0890 (int)sizeof(HitOnTrack));
0891 exit(1);
0892 }
0893 if (f_header.f_n_layers != expected_n_layers) {
0894 fprintf(stderr,
0895 "Number of layers on file (%d) is different from current TrackerInfo (%d).\n",
0896 f_header.f_n_layers,
0897 expected_n_layers);
0898 exit(1);
0899 }
0900
0901 printf("Opened file '%s', format version %d, n_layers %d, n_events %d\n",
0902 fname.c_str(),
0903 f_header.f_format_version,
0904 f_header.f_n_layers,
0905 f_header.f_n_events);
0906 if (f_header.f_extra_sections) {
0907 printf(" Extra sections:");
0908 if (f_header.f_extra_sections & ES_SimTrackStates)
0909 printf(" SimTrackStates");
0910 if (f_header.f_extra_sections & ES_Seeds)
0911 printf(" Seeds");
0912 if (f_header.f_extra_sections & ES_CmsswTracks)
0913 printf(" CmsswTracks");
0914 printf("\n");
0915 }
0916
0917 if (Config::seedInput == cmsswSeeds && !hasSeeds()) {
0918 fprintf(stderr, "Reading of CmsswSeeds requested but data not available on file.\n");
0919 exit(1);
0920 }
0921
0922 if (Config::readCmsswTracks && !hasCmsswTracks()) {
0923 fprintf(stderr, "Reading of CmsswTracks requested but data not available on file.\n");
0924 exit(1);
0925 }
0926
0927 return f_header.f_n_events;
0928 }
0929
0930 void DataFile::openWrite(const std::string &fname, int n_layers, int n_ev, int extra_sections) {
0931 f_fp = fopen(fname.c_str(), "w");
0932 f_header.f_n_layers = n_layers;
0933 f_header.f_n_events = n_ev;
0934 f_header.f_extra_sections = extra_sections;
0935
0936 fwrite(&f_header, sizeof(DataFileHeader), 1, f_fp);
0937 }
0938
0939 int DataFile::advancePosToNextEvent(FILE *fp) {
0940 int evsize;
0941
0942 std::lock_guard<std::mutex> readlock(f_next_ev_mutex);
0943
0944 fseek(fp, f_pos, SEEK_SET);
0945 fread(&evsize, sizeof(int), 1, fp);
0946 if (Config::loopOverFile) {
0947
0948 if (feof(fp) != 0) {
0949 f_pos = sizeof(DataFileHeader);
0950 fseek(fp, f_pos, SEEK_SET);
0951 fread(&evsize, sizeof(int), 1, fp);
0952 }
0953 }
0954
0955 f_pos += evsize;
0956
0957 return evsize;
0958 }
0959
0960 void DataFile::skipNEvents(int n_to_skip) {
0961 int evsize;
0962
0963 std::lock_guard<std::mutex> readlock(f_next_ev_mutex);
0964
0965 while (n_to_skip-- > 0) {
0966 fseek(f_fp, f_pos, SEEK_SET);
0967 fread(&evsize, sizeof(int), 1, f_fp);
0968 f_pos += evsize;
0969 }
0970 }
0971
0972 void DataFile::close() {
0973 fclose(f_fp);
0974 f_fp = 0;
0975 f_header = DataFileHeader();
0976 }
0977
0978 void DataFile::CloseWrite(int n_written) {
0979 if (f_header.f_n_events != n_written) {
0980 fseek(f_fp, 0, SEEK_SET);
0981 f_header.f_n_events = n_written;
0982 fwrite(&f_header, sizeof(DataFileHeader), 1, f_fp);
0983 }
0984 close();
0985 }
0986
0987 }
