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::reducePhiRange(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     if (fittedTrk.consistentSector()) {
0191       tmtt::KFTrackletTrack trk = fittedTrk.returnKFTrackletTrack();
0192
0193       if (settings_.printDebugKF())
0194         edm::LogVerbatim("L1track") << "Done with Kalman fit. Pars: pt = " << trk.pt()
0195                                     << ", 1/2R = " << settings_.bfield() * 3 * trk.qOverPt() / 2000
0196                                     << ", phi0 = " << trk.phi0() << ", eta = " << trk.eta() << ", z0 = " << trk.z0()
0197                                     << ", chi2 = " << trk.chi2() << ", accepted = " << trk.accepted();
0198
0199       double d0, chi2rphi, phi0, qoverpt = -999;
0200       if (trk.done_bcon()) {
0201         d0 = trk.d0_bcon();
0202         chi2rphi = trk.chi2rphi_bcon();
0203         phi0 = trk.phi0_bcon();
0204         qoverpt = trk.qOverPt_bcon();
0205       } else {
0206         d0 = trk.d0();
0207         chi2rphi = trk.chi2rphi();
0208         phi0 = trk.phi0();
0209         qoverpt = trk.qOverPt();
0210       }
0211
0212       // Tracklet wants phi0 with respect to lower edge of sector, not global phi0.
0213       double phi0fit = reco::reducePhiRange(phi0 - iSector_ * 2 * M_PI / N_SECTOR + 0.5 * settings_.dphisectorHG());
0214       double rinvfit = 0.01 * settings_.c() * settings_.bfield() * qoverpt;
0215
0216       int irinvfit = floor(rinvfit / settings_.krinvpars());
0217       int iphi0fit = floor(phi0fit / settings_.kphi0pars());
0218       int itanlfit = floor(trk.tanLambda() / settings_.ktpars());
0219       int iz0fit = floor(trk.z0() / settings_.kz0pars());
0220       int id0fit = floor(d0 / settings_.kd0pars());
0221       int ichi2rphifit = chi2rphi / 16;
0222       int ichi2rzfit = trk.chi2rz() / 16;
0223
0224       const vector<const tmtt::Stub*>& stubsFromFit = trk.stubs();
0225       vector<const L1TStub*> l1stubsFromFit;
0226       for (const tmtt::Stub* s : stubsFromFit) {
0227         unsigned int IDf = s->index();
0228         const L1TStub* l1s = L1StubIndices.at(IDf);
0229         l1stubsFromFit.push_back(l1s);
0230       }
0231
0232       tracklet->setFitPars(rinvfit,
0233                            phi0fit,
0234                            d0,
0235                            trk.tanLambda(),
0236                            trk.z0(),
0237                            chi2rphi,
0238                            trk.chi2rz(),
0239                            rinvfit,
0240                            phi0fit,
0241                            d0,
0242                            trk.tanLambda(),
0243                            trk.z0(),
0244                            chi2rphi,
0245                            trk.chi2rz(),
0246                            irinvfit,
0247                            iphi0fit,
0248                            id0fit,
0249                            itanlfit,
0250                            iz0fit,
0251                            ichi2rphifit,
0252                            ichi2rzfit,
0253                            trk.hitPattern(),
0254                            l1stubsFromFit);
0255     } else {
0256       if (settings_.printDebugKF()) {
0257         edm::LogVerbatim("L1track") << "FitTrack:KF rejected track";
0258       }
0259     }
0260   } else {
0261     if (settings_.printDebugKF()) {
0262       edm::LogVerbatim("L1track") << "FitTrack:KF rejected track";
0263     }
0264   }
0265
0266   for (const tmtt::Stub* s : TMTTstubs) {
0267     delete s;
0268   }
0269 }
0270 #endif