Back to home page

Project CMSSW displayed by LXR

 
 

    


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