Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-21 23:14:26

0001 #include "Event.h"
0002 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0003 
0004 //#define DEBUG
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();  // need to reset maps for every event.
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();  // need to reset maps for every event.
0057   }
0058 
0059   void Event::validate() {
0060     // special map needed for sim_val_for_cmssw + set the track scores
0061     if (Config::sim_val_for_cmssw) {
0062       validation_.makeRecoTkToSeedTkMapsDumbCMSSW(*this);
0063       validation_.setTrackScoresDumbCMSSW(*this);
0064     }
0065 
0066     // standard eff/fr/dr validation
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     // special cmssw to mkfit validation
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) {  // fit val for z-phi tuning
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);  // this will be overwritten at the end
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       //sizes are the same as in layerHits_
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     //layerHitMap_ is recreated afterwards
0179 
0180     /*
0181   printf("write %i tracks\n",nt);
0182   for (int it = 0; it<nt; it++) {
0183     printf("track with pT=%5.3f\n",simTracks_[it].pT());
0184     for (int ih=0; ih<simTracks_[it].nTotalHits(); ++ih) {
0185       printf("hit lyr:%2d idx=%i\n", simTracks_[it].getHitLyr(ih), simTracks_[it].getHitIdx(ih));
0186     }
0187   }
0188   printf("write %i layers\n",nl);
0189   for (int il = 0; il<nl; il++) {
0190     printf("write %i hits in layer %i\n",layerHits_[il].size(),il);
0191     for (int ih = 0; ih<layerHits_[il].size(); ih++) {
0192       printf("hit with r=%5.3f x=%5.3f y=%5.3f z=%5.3f\n",layerHits_[il][ih].r(),layerHits_[il][ih].x(),layerHits_[il][ih].y(),layerHits_[il][ih].z());
0193     }
0194   }
0195   */
0196   }
0197 
0198   // #define DUMP_SEEDS
0199   // #define DUMP_SEED_HITS
0200   // #define DUMP_TRACKS
0201   // #define DUMP_TRACK_HITS
0202   // #define DUMP_LAYER_HITS
0203   // #define DUMP_REC_TRACKS
0204   // #define DUMP_REC_TRACK_HITS
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);  //init to 0 by default
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     // HACK TO ONLY SELECT ONE PROBLEMATIC TRACK.
0291     // Note that MC matching gets screwed.
0292     // Works for MC seeding.
0293     //
0294     printf("************** SIM SELECTION HACK IN FORCE ********************\n");
0295     TrackVec x;
0296     x.push_back(simTracks_[3]);
0297     simTracks_.swap(x);
0298     nt = 1;
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     // Returns total number of bytes written.
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     // Returns number of read tracks (negative if actual reading was skipped).
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);  // -2 because data_size counts itself and n_tracks too
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     // Enforce Vxy on all layers, Vz on pixb only.
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     // Sim tracks from cmssw have the following issues:
0499     // - hits are not sorted by layer;
0500     // - there are tracks with too low number of hits, even 0;
0501     // - even with enough hits, there can be too few layers (esp. in endcap);
0502     // - tracks from secondaries can have extremely low pT.
0503     // Possible further checks:
0504     // - make sure enough hits exist in seeding layers.
0505     //
0506     // What is done:
0507     // 1. Hits are sorted by layer;
0508     // 2. Non-findable tracks are marked with Track::Status::not_findable flag.
0509     //
0510     // Returns number of passed simtracks.
0511 
0512     dprintf("Event::clean_cms_simtracks processing %lu simtracks.\n", simTracks_.size());
0513 
0514     int n_acc = 0;
0515     int i = -1;  //wrap in ifdef DEBUG?
0516     for (Track &t : simTracks_) {
0517       i++;
0518 
0519       t.sortHitsByLayer();
0520 
0521       const int lyr_cnt = t.nUniqueLayers();
0522 
0523       //const int lasthit = t.getLastFoundHitPos();
0524       //const float eta = layerHits_[t.getHitLyr(lasthit)][t.getHitIdx(lasthit)].eta();
0525 
0526       if (lyr_cnt < Config::cmsSelMinLayers)  // || Config::TrkInfo.is_transition(eta))
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     //WARNING: Printouts for hits will not make any sense if mkFit is not run with a validation flag such as --quality-val
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;  //FIXME: this speed up prevents transitive masking; check build cost!
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       //#pragma simd /* Vectorization via simd had issues with icc */
0643       for (int tss = ts + 1; tss < ns; tss++) {
0644         const float Pt2 = pt[tss];
0645 
0646         ////// Always require charge consistency. If different charge is assigned, do not remove seed-track
0647         if (charge[tss] != charge[ts])
0648           continue;
0649 
0650         const float thisDPt = std::abs(Pt2 - Pt1);
0651         ////// Require pT consistency between seeds. If dpT is large, do not remove seed-track.
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         ////// Reject tracks within dR-dz elliptical window.
0678         ////// Adaptive thresholds, based on observation that duplicates are more abundant at large pseudo-rapidity and low track pT
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};  //to be stored somewhere common
0723 
0724     //saving seeds by algorithm
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     //saving tracks by algorithm
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   // HitMask handling
0803   //==============================================================================
0804 
0805   void Event::fill_hitmask_bool_vectors(int track_algo, std::vector<std::vector<bool>> &layer_masks) {
0806     // Convert from per-hit uint64_t to per layer bool-vectors for given
0807     // iteration.
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     // Convert from per-hit uint64_t to per layer bool-vectors for a list of
0826     // iterations.
0827     // A hit mask is set if it is set for _all_ listed iterations.
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   // DataFile
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       // File ended, rewind back to beginning
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 }  // end namespace mkfit