Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-09 02:02:06 

0001 #include "L1Trigger/L1TTrackMatch/interface/DisplacedVertexProducer.h"
0002 #include "L1TrackUnpacker.h"
0003 
0004 using namespace l1trackunpacker;
0005 
0006 double DisplacedVertexProducer::FloatPtFromBits(const L1TTTrackType& track) const {
0007   ap_uint<14> ptEmulationBits = track.getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1,
0008                                                      TTTrack_TrackWord::TrackBitLocations::kRinvLSB);
0009   ap_ufixed<14, 9> ptEmulation;
0010   ptEmulation.V = (ptEmulationBits.range());
0011   return ptEmulation.to_double();
0012 }
0013 
0014 double DisplacedVertexProducer::FloatEtaFromBits(const L1TTTrackType& track) const {
0015   TTTrack_TrackWord::tanl_t etaBits = track.getTanlWord();
0016   glbeta_intern digieta;
0017   digieta.V = etaBits.range();
0018   return (double)digieta;
0019 }
0020 
0021 double DisplacedVertexProducer::FloatPhiFromBits(const L1TTTrackType& track) const {
0022   int Sector = track.phiSector();
0023   double sector_phi_value = 0;
0024   if (Sector < 5) {
0025     sector_phi_value = 2.0 * M_PI * Sector / 9.0;
0026   } else {
0027     sector_phi_value = (-1.0 * M_PI + M_PI / 9.0 + (Sector - 5) * 2.0 * M_PI / 9.0);
0028   }
0029   glbphi_intern trkphiSector = DoubleToBit(
0030       sector_phi_value, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0031   glbphi_intern local_phiBits = 0;
0032   local_phiBits.V = track.getPhiWord();
0033 
0034   glbphi_intern local_phi =
0035       DoubleToBit(BitToDouble(local_phiBits, TTTrack_TrackWord::TrackBitWidths::kPhiSize, TTTrack_TrackWord::stepPhi0),
0036                   TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit,
0037                   TTTrack_TrackWord::stepPhi0);
0038   glbphi_intern digiphi = local_phi + trkphiSector;
0039   return BitToDouble(
0040       digiphi, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0041 }
0042 
0043 int DisplacedVertexProducer::ChargeFromBits(const L1TTTrackType& track) const {
0044   ap_uint<1> chargeBit = track.getTrackWord()[TTTrack_TrackWord::TrackBitLocations::kRinvMSB];
0045   return 1 - (2 * chargeBit.to_uint());
0046 }
0047 
0048 double convertPtToR(double pt) {
0049   return 100.0 * (1.0 / (0.3 * 3.8)) * pt;  //returns R in cm
0050 }
0051 
0052 bool ComparePtTrack(std::pair<Track_Parameters, edm::Ptr<TrackingParticle>> a,
0053                     std::pair<Track_Parameters, edm::Ptr<TrackingParticle>> b) {
0054   return a.first.pt > b.first.pt;
0055 }
0056 
0057 Double_t dist(Double_t x1, Double_t y1, Double_t x2 = 0, Double_t y2 = 0) {  // Distance between 2 points
0058   return (std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)));
0059 }
0060 
0061 Double_t dist_TPs(Track_Parameters a, Track_Parameters b) {
0062   float x1 = a.x0;   //   Centers of the circles
0063   float y1 = a.y0;   //
0064   float x2 = b.x0;   //
0065   float y2 = b.y0;   //
0066   float R1 = a.rho;  // Radii of the circles
0067   float R2 = b.rho;
0068   float R = dist(x1, y1, x2, y2);  // Distance between centers
0069   if ((R >= fabs(R1 - R2)) && (R <= (R1 + R2))) {
0070     return (0);
0071   } else if (R == 0) {
0072     return (-99999.0);
0073   } else {
0074     return (R - R1 - R2);
0075   }
0076 }
0077 
0078 Int_t calcVertex(Track_Parameters a, Track_Parameters b, Double_t& x_vtx, Double_t& y_vtx, Double_t& z_vtx) {
0079   float x1 = a.x0;   //   Centers of the circles
0080   float y1 = a.y0;   //
0081   float x2 = b.x0;   //
0082   float y2 = b.y0;   //
0083   float R1 = a.rho;  // Radii of the circles
0084   float R2 = b.rho;
0085   float R = dist(x1, y1, x2, y2);  // Distance between centers
0086   if (R == 0)
0087     return -1;
0088   float co1 = (pow(R1, 2) - pow(R2, 2)) / (2 * pow(R, 2));
0089   float radicand = (2 / pow(R, 2)) * (pow(R1, 2) + pow(R2, 2)) - (pow(pow(R1, 2) - pow(R2, 2), 2) / pow(R, 4)) - 1;
0090   float co2 = 0;
0091   if (radicand > 0)
0092     co2 = 0.5 * std::sqrt(radicand);
0093   float ix1_x = 0.5 * (x1 + x2) + co1 * (x2 - x1) + co2 * (y2 - y1);
0094   float ix2_x = 0.5 * (x1 + x2) + co1 * (x2 - x1) - co2 * (y2 - y1);
0095   float ix1_y = 0.5 * (y1 + y2) + co1 * (y2 - y1) + co2 * (x1 - x2);
0096   float ix2_y = 0.5 * (y1 + y2) + co1 * (y2 - y1) - co2 * (x1 - x2);
0097   float ix1_z1 = a.trackZAtVertex(ix1_x, ix1_y);
0098   float ix1_z2 = b.trackZAtVertex(ix1_x, ix1_y);
0099   float ix1_delz = fabs(ix1_z1 - ix1_z2);
0100   float ix2_z1 = a.trackZAtVertex(ix2_x, ix2_y);
0101   float ix2_z2 = b.trackZAtVertex(ix2_x, ix2_y);
0102   float ix2_delz = fabs(ix2_z1 - ix2_z2);
0103   float trk1_POCA[2] = {a.d0 * std::sin(a.phi), -1 * a.d0 * std::cos(a.phi)};
0104   float trk2_POCA[2] = {b.d0 * std::sin(b.phi), -1 * b.d0 * std::cos(b.phi)};
0105   float trk1_ix1_delxy[2] = {ix1_x - trk1_POCA[0], ix1_y - trk1_POCA[1]};
0106   float trk1_ix2_delxy[2] = {ix2_x - trk1_POCA[0], ix2_y - trk1_POCA[1]};
0107   float trk2_ix1_delxy[2] = {ix1_x - trk2_POCA[0], ix1_y - trk2_POCA[1]};
0108   float trk2_ix2_delxy[2] = {ix2_x - trk2_POCA[0], ix2_y - trk2_POCA[1]};
0109   float trk1_traj[2] = {std::cos(a.phi), std::sin(a.phi)};
0110   float trk2_traj[2] = {std::cos(b.phi), std::sin(b.phi)};
0111   bool trk1_ix1_inTraj = ((trk1_ix1_delxy[0] * trk1_traj[0] + trk1_ix1_delxy[1] * trk1_traj[1]) > 0) ? true : false;
0112   bool trk1_ix2_inTraj = ((trk1_ix2_delxy[0] * trk1_traj[0] + trk1_ix2_delxy[1] * trk1_traj[1]) > 0) ? true : false;
0113   bool trk2_ix1_inTraj = ((trk2_ix1_delxy[0] * trk2_traj[0] + trk2_ix1_delxy[1] * trk2_traj[1]) > 0) ? true : false;
0114   bool trk2_ix2_inTraj = ((trk2_ix2_delxy[0] * trk2_traj[0] + trk2_ix2_delxy[1] * trk2_traj[1]) > 0) ? true : false;
0115   if (trk1_ix1_inTraj && trk2_ix1_inTraj && trk1_ix2_inTraj && trk2_ix2_inTraj) {
0116     if (ix1_delz < ix2_delz) {
0117       x_vtx = ix1_x;
0118       y_vtx = ix1_y;
0119       z_vtx = (ix1_z1 + ix1_z2) / 2;
0120       return 0;
0121     } else {
0122       x_vtx = ix2_x;
0123       y_vtx = ix2_y;
0124       z_vtx = (ix2_z1 + ix2_z2) / 2;
0125       return 0;
0126     }
0127   }
0128   if (trk1_ix1_inTraj && trk2_ix1_inTraj) {
0129     x_vtx = ix1_x;
0130     y_vtx = ix1_y;
0131     z_vtx = (ix1_z1 + ix1_z2) / 2;
0132     return 1;
0133   }
0134   if (trk1_ix2_inTraj && trk2_ix2_inTraj) {
0135     x_vtx = ix2_x;
0136     y_vtx = ix2_y;
0137     z_vtx = (ix2_z1 + ix2_z2) / 2;
0138     return 2;
0139   } else {
0140     if (ix1_delz < ix2_delz) {
0141       x_vtx = ix1_x;
0142       y_vtx = ix1_y;
0143       z_vtx = (ix1_z1 + ix1_z2) / 2;
0144       return 3;
0145     } else {
0146       x_vtx = ix2_x;
0147       y_vtx = ix2_y;
0148       z_vtx = (ix2_z1 + ix2_z2) / 2;
0149       return 3;
0150     }
0151   }
0152   return 4;
0153 }
0154 
0155 DisplacedVertexProducer::DisplacedVertexProducer(const edm::ParameterSet& iConfig)
0156     : ttTrackMCTruthToken_(consumes<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>>(
0157           iConfig.getParameter<edm::InputTag>("mcTruthTrackInputTag"))),
0158       trackToken_(consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>(
0159           iConfig.getParameter<edm::InputTag>("l1TracksInputTag"))),
0160       trackGTTToken_(consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>(
0161           iConfig.getParameter<edm::InputTag>("l1TracksGTTInputTag"))),
0162       outputVertexCollectionName_(iConfig.getParameter<std::string>("l1TrackVertexCollectionName")),
0163       model_(iConfig.getParameter<edm::FileInPath>("model")),
0164       runEmulation_(iConfig.getParameter<bool>("runEmulation")),
0165       cutSet_(iConfig.getParameter<edm::ParameterSet>("cutSet")),
0166       chi2rzMax_(cutSet_.getParameter<double>("chi2rzMax")),
0167       promptMVAMin_(cutSet_.getParameter<double>("promptMVAMin")),
0168       ptMin_(cutSet_.getParameter<double>("ptMin")),
0169       etaMax_(cutSet_.getParameter<double>("etaMax")),
0170       dispD0Min_(cutSet_.getParameter<double>("dispD0Min")),
0171       promptMVADispTrackMin_(cutSet_.getParameter<double>("promptMVADispTrackMin")),
0172       overlapEtaMin_(cutSet_.getParameter<double>("overlapEtaMin")),
0173       overlapEtaMax_(cutSet_.getParameter<double>("overlapEtaMax")),
0174       overlapNStubsMin_(cutSet_.getParameter<int>("overlapNStubsMin")),
0175       diskEtaMin_(cutSet_.getParameter<double>("diskEtaMin")),
0176       diskD0Min_(cutSet_.getParameter<double>("diskD0Min")),
0177       barrelD0Min_(cutSet_.getParameter<double>("barrelD0Min")),
0178       RTMin_(cutSet_.getParameter<double>("RTMin")),
0179       RTMax_(cutSet_.getParameter<double>("RTMax")) {
0180   //--- Define EDM output to be written to file (if required)
0181   produces<l1t::DisplacedTrackVertexCollection>(outputVertexCollectionName_);
0182 }
0183 
0184 void DisplacedVertexProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0185   edm::Handle<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTTrackHandle;
0186   iEvent.getByToken(ttTrackMCTruthToken_, MCTruthTTTrackHandle);
0187   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackHandle;
0188   iEvent.getByToken(trackToken_, TTTrackHandle);
0189   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackGTTHandle;
0190   iEvent.getByToken(trackGTTToken_, TTTrackGTTHandle);
0191   std::vector<TTTrack<Ref_Phase2TrackerDigi_>>::const_iterator iterL1Track;
0192   int this_l1track = 0;
0193   std::vector<std::pair<Track_Parameters, edm::Ptr<TrackingParticle>>> selectedTracksWithTruth;
0194 
0195   //track selection loop
0196   for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
0197     edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>> l1track_ptr(TTTrackHandle, this_l1track);
0198     edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>> gtttrack_ptr(TTTrackGTTHandle, this_l1track);
0199     this_l1track++;
0200 
0201     float pt, eta, phi, z0, d0, rho, chi2rphi, chi2rz, bendchi2, MVA1;
0202     int nstub;
0203 
0204     if (runEmulation_) {
0205       pt = FloatPtFromBits(*gtttrack_ptr);
0206       eta = FloatEtaFromBits(*gtttrack_ptr);
0207       phi = FloatPhiFromBits(*gtttrack_ptr);
0208       z0 = gtttrack_ptr->getZ0();  //cm
0209       d0 = gtttrack_ptr->getD0();
0210       int charge = ChargeFromBits(*gtttrack_ptr);
0211       rho = charge * convertPtToR(pt);
0212       chi2rphi = gtttrack_ptr->getChi2RPhi();
0213       chi2rz = gtttrack_ptr->getChi2RZ();
0214       bendchi2 = gtttrack_ptr->getBendChi2();
0215       MVA1 = gtttrack_ptr->getMVAQuality();
0216       nstub = gtttrack_ptr->getNStubs();
0217     } else {
0218       pt = l1track_ptr->momentum().perp();
0219       eta = l1track_ptr->momentum().eta();
0220       phi = l1track_ptr->momentum().phi();
0221       z0 = l1track_ptr->z0();  //cm
0222       float x0 = l1track_ptr->POCA().x();
0223       float y0 = l1track_ptr->POCA().y();
0224       d0 = x0 * std::sin(phi) - y0 * std::cos(phi);
0225       rho = 1 / l1track_ptr->rInv();
0226       chi2rphi = l1track_ptr->chi2XYRed();
0227       chi2rz = l1track_ptr->chi2ZRed();
0228       bendchi2 = l1track_ptr->stubPtConsistency();
0229       MVA1 = l1track_ptr->trkMVA1();
0230       std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
0231           stubRefs = l1track_ptr->getStubRefs();
0232       nstub = (int)stubRefs.size();
0233     }
0234 
0235     if (chi2rz < chi2rzMax_ && MVA1 > promptMVAMin_ && pt > ptMin_ && fabs(eta) < etaMax_) {
0236       if (fabs(d0) > dispD0Min_) {
0237         if (MVA1 <= promptMVADispTrackMin_)
0238           continue;
0239       }
0240       if (fabs(eta) > overlapEtaMin_ && fabs(eta) < overlapEtaMax_) {
0241         if (nstub <= overlapNStubsMin_)
0242           continue;
0243       }
0244       if (fabs(eta) > diskEtaMin_) {
0245         if (fabs(d0) <= diskD0Min_)
0246           continue;
0247       }
0248       if (fabs(eta) <= diskEtaMin_) {
0249         if (fabs(d0) <= barrelD0Min_)
0250           continue;
0251       }
0252 
0253       Track_Parameters track =
0254           Track_Parameters(pt, d0, z0, eta, phi, rho, (this_l1track - 1), nstub, chi2rphi, chi2rz, bendchi2, MVA1);
0255 
0256       edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(l1track_ptr);
0257       selectedTracksWithTruth.push_back(std::make_pair(track, my_tp));
0258     }
0259   }
0260 
0261   //vertex loop
0262   std::unique_ptr<l1t::DisplacedTrackVertexCollection> product(new std::vector<l1t::DisplacedTrackVertex>());
0263   for (int i = 0; i < int(selectedTracksWithTruth.size() - 1); i++) {
0264     for (int j = i + 1; j < int(selectedTracksWithTruth.size()); j++) {
0265       if (dist_TPs(selectedTracksWithTruth[i].first, selectedTracksWithTruth[j].first) != 0)
0266         continue;
0267       Double_t x_dv_trk = -9999.0;
0268       Double_t y_dv_trk = -9999.0;
0269       Double_t z_dv_trk = -9999.0;
0270       edm::Ptr<TrackingParticle> tp_i = selectedTracksWithTruth[i].second;
0271       edm::Ptr<TrackingParticle> tp_j = selectedTracksWithTruth[j].second;
0272       bool isReal = false;
0273       if (!tp_i.isNull() && !tp_j.isNull()) {
0274         bool isHard_i = false;
0275         bool isHard_j = false;
0276         if (!tp_i->genParticles().empty() && !tp_j->genParticles().empty()) {
0277           isHard_i = tp_i->genParticles()[0]->isHardProcess() || tp_i->genParticles()[0]->fromHardProcessFinalState();
0278           isHard_j = tp_j->genParticles()[0]->isHardProcess() || tp_j->genParticles()[0]->fromHardProcessFinalState();
0279         }
0280 
0281         if (tp_i->eventId().event() == 0 && tp_j->eventId().event() == 0 && fabs(tp_i->vx() - tp_j->vx()) < 0.0001 &&
0282             fabs(tp_i->vy() - tp_j->vy()) < 0.0001 && fabs(tp_i->vz() - tp_j->vz()) < 0.0001 && isHard_i && isHard_j &&
0283             ((tp_i->charge() + tp_j->charge()) == 0)) {
0284           isReal = true;
0285         }
0286       }
0287 
0288       int inTraj =
0289           calcVertex(selectedTracksWithTruth[i].first, selectedTracksWithTruth[j].first, x_dv_trk, y_dv_trk, z_dv_trk);
0290       Vertex_Parameters vertex = Vertex_Parameters(
0291           x_dv_trk, y_dv_trk, z_dv_trk, selectedTracksWithTruth[i].first, selectedTracksWithTruth[j].first);
0292 
0293       if (vertex.R_T > RTMax_)
0294         continue;
0295       if (vertex.R_T < RTMin_)
0296         continue;
0297 
0298       l1t::DisplacedTrackVertex outputVertex = l1t::DisplacedTrackVertex(selectedTracksWithTruth[i].first.index,
0299                                                                          selectedTracksWithTruth[j].first.index,
0300                                                                          inTraj,
0301                                                                          vertex.d_T,
0302                                                                          vertex.R_T,
0303                                                                          vertex.cos_T,
0304                                                                          vertex.delta_z,
0305                                                                          vertex.x_dv,
0306                                                                          vertex.y_dv,
0307                                                                          vertex.z_dv,
0308                                                                          vertex.openingAngle,
0309                                                                          vertex.p_mag,
0310                                                                          isReal);
0311 
0312       //Rescaling input features so they all fall within [-20,20]. This reduces bits needed in emulation by 2. See this presentation for more information: https://indico.cern.ch/event/1476881/contributions/6219913/attachments/2964052/5214060/GTT%20Displaced%20Vertexing%20November%208%202024.pdf
0313       float ptRescaling = 0.25;
0314       float deltaZRescaling = 0.125;
0315       if (runEmulation_) {
0316         std::vector<ap_fixed<13, 8, AP_RND_CONV, AP_SAT>> Transformed_features = {
0317             selectedTracksWithTruth[i].first.pt * ptRescaling,
0318             selectedTracksWithTruth[j].first.pt * ptRescaling,
0319             selectedTracksWithTruth[i].first.eta,
0320             selectedTracksWithTruth[j].first.eta,
0321             selectedTracksWithTruth[i].first.phi,
0322             selectedTracksWithTruth[j].first.phi,
0323             selectedTracksWithTruth[i].first.d0,
0324             selectedTracksWithTruth[j].first.d0,
0325             selectedTracksWithTruth[i].first.z0,
0326             selectedTracksWithTruth[j].first.z0,
0327             selectedTracksWithTruth[i].first.chi2rz,
0328             selectedTracksWithTruth[j].first.chi2rz,
0329             selectedTracksWithTruth[i].first.bendchi2,
0330             selectedTracksWithTruth[j].first.bendchi2,
0331             selectedTracksWithTruth[i].first.MVA1,
0332             selectedTracksWithTruth[j].first.MVA1,
0333             vertex.d_T,
0334             vertex.R_T,
0335             vertex.cos_T,
0336             vertex.delta_z * deltaZRescaling};
0337         conifer::BDT<ap_fixed<13, 8, AP_RND_CONV, AP_SAT>, ap_fixed<13, 8, AP_RND_CONV, AP_SAT>> bdt(
0338             this->model_.fullPath());
0339         std::vector<ap_fixed<13, 8, AP_RND_CONV, AP_SAT>> output = bdt.decision_function(Transformed_features);
0340         outputVertex.setScore(output.at(0).to_float());
0341       } else {
0342         std::vector<float> Transformed_features = {float(selectedTracksWithTruth[i].first.pt * ptRescaling),
0343                                                    float(selectedTracksWithTruth[j].first.pt * ptRescaling),
0344                                                    selectedTracksWithTruth[i].first.eta,
0345                                                    selectedTracksWithTruth[j].first.eta,
0346                                                    selectedTracksWithTruth[i].first.phi,
0347                                                    selectedTracksWithTruth[j].first.phi,
0348                                                    selectedTracksWithTruth[i].first.d0,
0349                                                    selectedTracksWithTruth[j].first.d0,
0350                                                    selectedTracksWithTruth[i].first.z0,
0351                                                    selectedTracksWithTruth[j].first.z0,
0352                                                    float(selectedTracksWithTruth[i].first.chi2rz),
0353                                                    float(selectedTracksWithTruth[j].first.chi2rz),
0354                                                    float(selectedTracksWithTruth[i].first.bendchi2),
0355                                                    float(selectedTracksWithTruth[j].first.bendchi2),
0356                                                    float(selectedTracksWithTruth[i].first.MVA1),
0357                                                    float(selectedTracksWithTruth[j].first.MVA1),
0358                                                    vertex.d_T,
0359                                                    vertex.R_T,
0360                                                    vertex.cos_T,
0361                                                    float(vertex.delta_z * deltaZRescaling)};
0362         conifer::BDT<float, float> bdt(this->model_.fullPath());
0363         std::vector<float> output = bdt.decision_function(Transformed_features);
0364         outputVertex.setScore(output.at(0));
0365       }
0366 
0367       product->emplace_back(outputVertex);
0368     }
0369   }
0370 
0371   // //=== Store output
0372   iEvent.put(std::move(product), outputVertexCollectionName_);
0373 }
0374 
0375 DEFINE_FWK_MODULE(DisplacedVertexProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!