TrackFindingTracklet/src/HybridFit.cc

0001 #include "L1Trigger/TrackFindingTracklet/interface/HybridFit.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/L1TStub.h"
0005
0006 #ifdef USEHYBRID
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DataFormats/Math/interface/deltaPhi.h"
0009
0010 using namespace std;
0011 using namespace trklet;
0012
0013 HybridFit::HybridFit(unsigned int iSector, Settings const& settings, Globals* globals) : settings_(settings) {
0014   iSector_ = iSector;
0015   globals_ = globals;
0016 }
0017
0018 void HybridFit::Fit(Tracklet* tracklet, std::vector<const Stub*>& trackstublist) {
0019   if (settings_.fakefit()) {
0020     vector<const L1TStub*> l1stubsFromFitTrack;
0021     for (unsigned int k = 0; k < trackstublist.size(); k++) {
0022       const L1TStub* L1stub = trackstublist[k]->l1tstub();
0023       l1stubsFromFitTrack.push_back(L1stub);
0024     }
0025     tracklet->setFitPars(tracklet->rinvapprox(),
0026                          tracklet->phi0approx(),
0027                          tracklet->d0approx(),
0028                          tracklet->tapprox(),
0029                          tracklet->z0approx(),
0030                          0.,
0031                          0.,
0032                          tracklet->rinv(),
0033                          tracklet->phi0(),
0034                          tracklet->d0(),
0035                          tracklet->t(),
0036                          tracklet->z0(),
0037                          0.,
0038                          0.,
0039                          tracklet->fpgarinv().value(),
0040                          tracklet->fpgaphi0().value(),
0041                          tracklet->fpgad0().value(),
0042                          tracklet->fpgat().value(),
0043                          tracklet->fpgaz0().value(),
0044                          0,
0045                          0,
0046                          0,
0047                          l1stubsFromFitTrack);
0048     return;
0049   }
0050
0051   std::vector<tmtt::Stub*> TMTTstubs;
0052   std::map<unsigned int, const L1TStub*> L1StubIndices;
0053   unsigned int L1stubID = 0;
0054
0055   if (globals_->tmttSettings() == nullptr) {
0056     if (settings_.printDebugKF())
0057       edm::LogVerbatim("L1track") << "Creating TMTT::Settings in HybridFit::Fit";
0058     globals_->tmttSettings() = make_unique<tmtt::Settings>();
0059     globals_->tmttSettings()->setMagneticField(settings_.bfield());
0060   }
0061
0062   const tmtt::Settings& TMTTsettings = *globals_->tmttSettings();
0063
0064   int kf_phi_sec = iSector_;
0065
0066   for (unsigned int k = 0; k < trackstublist.size(); k++) {
0067     const L1TStub* L1stubptr = trackstublist[k]->l1tstub();
0068
0069     double kfphi = L1stubptr->phi();
0070     double kfr = L1stubptr->r();
0071     double kfz = L1stubptr->z();
0072     double kfbend = L1stubptr->bend();
0073     bool psmodule = L1stubptr->isPSmodule();
0074     unsigned int iphi = L1stubptr->iphi();
0075     double alpha = L1stubptr->alpha(settings_.stripPitch(psmodule));
0076     bool isTilted = L1stubptr->isTilted();
0077
0078     bool isBarrel = trackstublist[k]->layerdisk() < N_LAYER;
0079     int kflayer;
0080
0081     if (isBarrel) {  // Barrel-specific
0082       kflayer = L1stubptr->layer() + 1;
0083       if (settings_.printDebugKF())
0084         edm::LogVerbatim("L1track") << "Will create layer stub with : ";
0085     } else {  // Disk-specific
0086       kflayer = abs(L1stubptr->disk());
0087       if (kfz > 0) {
0088         kflayer += 10;
0089       } else {
0090         kflayer += 20;
0091       }
0092       if (settings_.printDebugKF())
0093         edm::LogVerbatim("L1track") << "Will create disk stub with : ";
0094     }
0095
0096     float stripPitch = settings_.stripPitch(psmodule);
0097     float stripLength = settings_.stripLength(psmodule);
0098     unsigned int nStrips = settings_.nStrips(psmodule);
0099
0100     if (settings_.printDebugKF()) {
0101       edm::LogVerbatim("L1track") << kfphi << " " << kfr << " " << kfz << " " << kfbend << " " << kflayer << " "
0102                                   << isBarrel << " " << psmodule << " " << isTilted << " \n"
0103                                   << stripPitch << " " << stripLength << " " << nStrips;
0104     }
0105
0106     unsigned int uniqueStubIndex = 1000 * L1stubID + L1stubptr->allStubIndex();
0107     tmtt::Stub* TMTTstubptr = new tmtt::Stub(&TMTTsettings,
0108                                              uniqueStubIndex,
0109                                              kfphi,
0110                                              kfr,
0111                                              kfz,
0112                                              kfbend,
0113                                              iphi,
0114                                              -alpha,
0115                                              kflayer,
0116                                              kf_phi_sec,
0117                                              psmodule,
0118                                              isBarrel,
0119                                              isTilted,
0120                                              stripPitch,
0121                                              stripLength,
0122                                              nStrips);
0123     TMTTstubs.push_back(TMTTstubptr);
0124     L1StubIndices[uniqueStubIndex] = L1stubptr;
0125     L1stubID++;
0126   }
0127
0128   if (settings_.printDebugKF()) {
0129     edm::LogVerbatim("L1track") << "Made TMTTstubs: trackstublist.size() = " << trackstublist.size();
0130   }
0131
0132   double kfrinv = tracklet->rinvapprox();
0133   double kfphi0 = tracklet->phi0approx();
0134   double kfz0 = tracklet->z0approx();
0135   double kft = tracklet->tapprox();
0136   double kfd0 = tracklet->d0approx();
0137
0138   if (settings_.printDebugKF()) {
0139     edm::LogVerbatim("L1track") << "tracklet phi0 = " << kfphi0 << "\n"
0140                                 << "iSector = " << iSector_ << "\n"
0141                                 << "dphisectorHG = " << settings_.dphisectorHG();
0142   }
0143
0144   // KF wants global phi0, not phi0 measured with respect to lower edge of sector (Tracklet convention).
0145   kfphi0 = reco::reduceRange(kfphi0 + iSector_ * settings_.dphisector() - 0.5 * settings_.dphisectorHG());
0146
0147   std::pair<float, float> helixrphi(kfrinv / (0.01 * settings_.c() * settings_.bfield()), kfphi0);
0148   std::pair<float, float> helixrz(kfz0, kft);
0149
0150   // KF HLS uses HT mbin (which is binned q/Pt) to allow for scattering. So estimate it from tracklet.
0151   double chargeOverPt = helixrphi.first;
0152   int mBin = std::floor(TMTTsettings.houghNbinsPt() / 2) +
0153              std::floor((TMTTsettings.houghNbinsPt() / 2) * chargeOverPt / (1. / TMTTsettings.houghMinPt()));
0154   mBin = max(min(mBin, int(TMTTsettings.houghNbinsPt() - 1)), 0);  // protect precision issues.
0155   std::pair<unsigned int, unsigned int> celllocation(mBin, 1);
0156
0157   // Get range in z of tracks covered by this sector at chosen radius from beam-line
0158   const vector<double> etaRegions = TMTTsettings.etaRegions();
0159   const float chosenRofZ = TMTTsettings.chosenRofZ();
0160
0161   float kfzRef = kfz0 + chosenRofZ * kft;
0162
0163   unsigned int kf_eta_reg = 0;
0164   for (unsigned int iEtaSec = 1; iEtaSec < etaRegions.size() - 1; iEtaSec++) {  // Doesn't apply eta < 2.4 cut.
0165     const float etaMax = etaRegions[iEtaSec];
0166     const float zRefMax = chosenRofZ / tan(2. * atan(exp(-etaMax)));
0167     if (kfzRef > zRefMax)
0168       kf_eta_reg = iEtaSec;
0169   }
0170
0171   tmtt::L1track3D l1track3d(
0172       &TMTTsettings, TMTTstubs, celllocation, helixrphi, helixrz, kfd0, kf_phi_sec, kf_eta_reg, 1, false);
0173   unsigned int seedType = tracklet->getISeed();
0174   unsigned int numPS = tracklet->PSseed();  // Function PSseed() is out of date!
0175   l1track3d.setSeedLayerType(seedType);
0176   l1track3d.setSeedPS(numPS);
0177
0178   if (globals_->tmttKFParamsComb() == nullptr) {
0179     if (settings_.printDebugKF())
0180       edm::LogVerbatim("L1track") << "Will make KFParamsComb for " << settings_.nHelixPar() << " param fit";
0181     globals_->tmttKFParamsComb() = make_unique<tmtt::KFParamsComb>(&TMTTsettings, settings_.nHelixPar(), "KFfitter");
0182   }
0183
0184   tmtt::KFParamsComb& fitterKF = *globals_->tmttKFParamsComb();
0185
0186   // Call Kalman fit
0187   tmtt::L1fittedTrack fittedTrk = fitterKF.fit(l1track3d);
0188
0189   if (fittedTrk.accepted()) {
0190     tmtt::KFTrackletTrack trk = fittedTrk.returnKFTrackletTrack();
0191
0192     if (settings_.printDebugKF())
0193       edm::LogVerbatim("L1track") << "Done with Kalman fit. Pars: pt = " << trk.pt()
0194                                   << ", 1/2R = " << settings_.bfield() * 3 * trk.qOverPt() / 2000
0195                                   << ", phi0 = " << trk.phi0() << ", eta = " << trk.eta() << ", z0 = " << trk.z0()
0196                                   << ", chi2 = " << trk.chi2() << ", accepted = " << trk.accepted();
0197
0198     double d0, chi2rphi, phi0, qoverpt = -999;
0199     if (trk.done_bcon()) {
0200       d0 = trk.d0_bcon();
0201       chi2rphi = trk.chi2rphi_bcon();
0202       phi0 = trk.phi0_bcon();
0203       qoverpt = trk.qOverPt_bcon();
0204     } else {
0205       d0 = trk.d0();
0206       chi2rphi = trk.chi2rphi();
0207       phi0 = trk.phi0();
0208       qoverpt = trk.qOverPt();
0209     }
0210
0211     // Tracklet wants phi0 with respect to lower edge of sector, not global phi0.
0212     double phi0fit = reco::reduceRange(phi0 - iSector_ * 2 * M_PI / N_SECTOR + 0.5 * settings_.dphisectorHG());
0213     double rinvfit = 0.01 * settings_.c() * settings_.bfield() * qoverpt;
0214
0215     int irinvfit = floor(rinvfit / settings_.krinvpars());
0216     int iphi0fit = floor(phi0fit / settings_.kphi0pars());
0217     int itanlfit = floor(trk.tanLambda() / settings_.ktpars());
0218     int iz0fit = floor(trk.z0() / settings_.kz0pars());
0219     int id0fit = floor(d0 / settings_.kd0pars());
0220     int ichi2rphifit = chi2rphi / 16;
0221     int ichi2rzfit = trk.chi2rz() / 16;
0222
0223     const vector<const tmtt::Stub*>& stubsFromFit = trk.stubs();
0224     vector<const L1TStub*> l1stubsFromFit;
0225     for (const tmtt::Stub* s : stubsFromFit) {
0226       unsigned int IDf = s->index();
0227       const L1TStub* l1s = L1StubIndices.at(IDf);
0228       l1stubsFromFit.push_back(l1s);
0229     }
0230
0231     tracklet->setFitPars(rinvfit,
0232                          phi0fit,
0233                          d0,
0234                          trk.tanLambda(),
0235                          trk.z0(),
0236                          chi2rphi,
0237                          trk.chi2rz(),
0238                          rinvfit,
0239                          phi0fit,
0240                          d0,
0241                          trk.tanLambda(),
0242                          trk.z0(),
0243                          chi2rphi,
0244                          trk.chi2rz(),
0245                          irinvfit,
0246                          iphi0fit,
0247                          id0fit,
0248                          itanlfit,
0249                          iz0fit,
0250                          ichi2rphifit,
0251                          ichi2rzfit,
0252                          trk.hitPattern(),
0253                          l1stubsFromFit);
0254   } else {
0255     if (settings_.printDebugKF()) {
0256       edm::LogVerbatim("L1track") << "FitTrack:KF rejected track";
0257     }
0258   }
0259
0260   for (const tmtt::Stub* s : TMTTstubs) {
0261     delete s;
0262   }
0263 }
0264 #endif