MkFitCore/standalone/Event.cc

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