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