Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:58 

0001 #include "L1Trigger/L1TMuonEndCap/interface/PtAssignmentEngine2017.h"
0002 #include "L1Trigger/L1TMuonEndCap/interface/PtAssignmentEngineAux2017.h"
0003 
0004 #include <iostream>
0005 #include <sstream>
0006 
0007 const PtAssignmentEngineAux2017& PtAssignmentEngine2017::aux() const {
0008   static const PtAssignmentEngineAux2017
0009       instance;  // KK: arguable design solution, but const qualifier makes it thread-safe anyway
0010   return instance;
0011 }
0012 
0013 float PtAssignmentEngine2017::scale_pt(const float pt, const int mode) const {
0014   emtf_assert(ptLUTVersion_ >= 6);
0015 
0016   float pt_xml = -99;
0017   float pt_scale = -99;
0018 
0019   // Scaling to achieve 90% efficency at any given L1 pT threshold
0020   // For now, a linear scaling based on SingleMu-quality (modes 11, 13, 14, 15), CSC+RPC tracks
0021   // Should maybe scale each mode differently in the future - AWB 31.05.17
0022 
0023   // TRG       = (1.2 + 0.015*TRG) * XML
0024   // TRG       = 1.2*XML / (1 - 0.015*XML)
0025   // TRG / XML = 1.2 / (1 - 0.015*XML)
0026   if (ptLUTVersion_ >= 9) {  // LUTs with lower pT scale for 2023, deployed in May 2023
0027     pt_xml = fmin(20., pt);  // Maximum scale set by muons with XML pT = 20 GeV (scaled pT ~31 GeV)
0028     pt_scale = 1.07 / (1 - 0.015 * pt_xml);
0029   } else if (ptLUTVersion_ == 8) {  // First "physics" LUTs for 2022, will be deployed in June 2022
0030     pt_xml = fmin(20., pt);         // Maximum scale set by muons with XML pT = 20 GeV (scaled pT ~32 GeV)
0031     pt_scale = 1.13 / (1 - 0.015 * pt_xml);
0032   } else if (ptLUTVersion_ >= 6) {  // First "physics" LUTs for 2017, deployed June 7
0033     pt_xml = fmin(20., pt);         // Maximum scale set by muons with XML pT = 20 GeV (scaled pT ~35 GeV)
0034     pt_scale = 1.2 / (1 - 0.015 * pt_xml);
0035   }
0036 
0037   return pt_scale;
0038 }
0039 
0040 float PtAssignmentEngine2017::unscale_pt(const float pt, const int mode) const {
0041   emtf_assert(ptLUTVersion_ >= 6);
0042 
0043   float pt_unscale = -99;
0044 
0045   if (ptLUTVersion_ >= 9) {  // LUTs with lower pT scale for 2023, deployed in May 2023
0046     pt_unscale = 1 / (1.07 + 0.015 * pt);
0047     pt_unscale = fmax(pt_unscale, (1 - 0.015 * 20) / 1.07);
0048   } else if (ptLUTVersion_ == 8) {  // First "physics" LUTs for 2022, will be deployed in June 2022
0049     pt_unscale = 1 / (1.13 + 0.015 * pt);
0050     pt_unscale = fmax(pt_unscale, (1 - 0.015 * 20) / 1.13);
0051   } else if (ptLUTVersion_ >= 6) {  // First "physics" LUTs for 2017, deployed June 7
0052     pt_unscale = 1 / (1.2 + 0.015 * pt);
0053     pt_unscale = fmax(pt_unscale, (1 - 0.015 * 20) / 1.2);
0054   }
0055 
0056   return pt_unscale;
0057 }
0058 
0059 PtAssignmentEngine::address_t PtAssignmentEngine2017::calculate_address(const EMTFTrack& track) const {
0060   address_t address = 0;
0061 
0062   EMTFPtLUT data = track.PtLUT();
0063 
0064   int mode = track.Mode();
0065   int theta = track.Theta_fp();
0066   int endcap = track.Endcap();
0067   int nHits = (mode / 8) + ((mode % 8) / 4) + ((mode % 4) / 2) + ((mode % 2) / 1);
0068   emtf_assert(nHits > 1 && nHits < 5);
0069 
0070   // 'A' is first station in the track, 'B' the second, etc.
0071   int mode_ID = -1;
0072   int iA = -1, iB = -1, iC = -1, iD = -1;
0073   int iAB, iAC, iAD, iBC, iCD;
0074 
0075   switch (mode) {  // Indices for quantities by station or station pair
0076     case 15:
0077       mode_ID = 0b1;
0078       iA = 0;
0079       iB = 1;
0080       iC = 2;
0081       iD = 3;
0082       break;
0083     case 14:
0084       mode_ID = 0b11;
0085       iA = 0;
0086       iB = 1;
0087       iC = 2;
0088       break;
0089     case 13:
0090       mode_ID = 0b10;
0091       iA = 0;
0092       iB = 1;
0093       iC = 3;
0094       break;
0095     case 11:
0096       mode_ID = 0b01;
0097       iA = 0;
0098       iB = 2;
0099       iC = 3;
0100       break;
0101     case 7:
0102       mode_ID = 0b1;
0103       iA = 1;
0104       iB = 2;
0105       iC = 3;
0106       break;
0107     case 12:
0108       mode_ID = 0b111;
0109       iA = 0;
0110       iB = 1;
0111       break;
0112     case 10:
0113       mode_ID = 0b110;
0114       iA = 0;
0115       iB = 2;
0116       break;
0117     case 9:
0118       mode_ID = 0b101;
0119       iA = 0;
0120       iB = 3;
0121       break;
0122     case 6:
0123       mode_ID = 0b100;
0124       iA = 1;
0125       iB = 2;
0126       break;
0127     case 5:
0128       mode_ID = 0b011;
0129       iA = 1;
0130       iB = 3;
0131       break;
0132     case 3:
0133       mode_ID = 0b010;
0134       iA = 2;
0135       iB = 3;
0136       break;
0137     default:
0138       break;
0139   }
0140   iAB = (iA >= 0 && iB >= 0) ? iA + iB - (iA == 0) : -1;
0141   iAC = (iA >= 0 && iC >= 0) ? iA + iC - (iA == 0) : -1;
0142   iAD = (iA >= 0 && iD >= 0) ? 2 : -1;
0143   iBC = (iB >= 0 && iC >= 0) ? iB + iC : -1;
0144   iCD = (iC >= 0 && iD >= 0) ? 5 : -1;
0145 
0146   // Fill variable info from pT LUT data
0147   int st1_ring2 = data.st1_ring2;
0148   if (nHits == 4) {
0149     int dPhiAB = data.delta_ph[iAB];
0150     int dPhiBC = data.delta_ph[iBC];
0151     int dPhiCD = data.delta_ph[iCD];
0152     int sPhiAB = data.sign_ph[iAB];
0153     int sPhiBC = (data.sign_ph[iBC] == sPhiAB);
0154     int sPhiCD = (data.sign_ph[iCD] == sPhiAB);
0155     int dTheta = data.delta_th[iAD] * (data.sign_th[iAD] ? 1 : -1);
0156     int frA = data.fr[iA];
0157     int clctA = data.cpattern[iA];
0158     int clctB = data.cpattern[iB];
0159     int clctC = data.cpattern[iC];
0160     int clctD = data.cpattern[iD];
0161 
0162     // Convert variables to words for pT LUT address
0163     dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
0164     dPhiBC = aux().getNLBdPhiBin(dPhiBC, 5, 256);
0165     dPhiCD = aux().getNLBdPhiBin(dPhiCD, 4, 256);
0166     dTheta = aux().getdTheta(dTheta, 2);
0167     // Combines track theta, stations with RPC hits, and station 1 bend information
0168     int mode15_8b = aux().get8bMode15(theta, st1_ring2, endcap, (sPhiAB == 1 ? 1 : -1), clctA, clctB, clctC, clctD);
0169 
0170     // Form the pT LUT address
0171     address |= (dPhiAB & ((1 << 7) - 1)) << (0);
0172     address |= (dPhiBC & ((1 << 5) - 1)) << (0 + 7);
0173     address |= (dPhiCD & ((1 << 4) - 1)) << (0 + 7 + 5);
0174     address |= (sPhiBC & ((1 << 1) - 1)) << (0 + 7 + 5 + 4);
0175     address |= (sPhiCD & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1);
0176     address |= (dTheta & ((1 << 2) - 1)) << (0 + 7 + 5 + 4 + 1 + 1);
0177     address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2);
0178     address |= (mode15_8b & ((1 << 8) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1);
0179     address |= (mode_ID & ((1 << 1) - 1)) << (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1 + 8);
0180     emtf_assert(address < pow(2, 30) && address >= pow(2, 29));
0181   } else if (nHits == 3) {
0182     int dPhiAB = data.delta_ph[iAB];
0183     int dPhiBC = data.delta_ph[iBC];
0184     int sPhiAB = data.sign_ph[iAB];
0185     int sPhiBC = (data.sign_ph[iBC] == sPhiAB);
0186     int dTheta = data.delta_th[iAC] * (data.sign_th[iAC] ? 1 : -1);
0187     int frA = data.fr[iA];
0188     int frB = data.fr[iB];
0189     int clctA = data.cpattern[iA];
0190     int clctB = data.cpattern[iB];
0191     int clctC = data.cpattern[iC];
0192 
0193     // Convert variables to words for pT LUT address
0194     dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
0195     dPhiBC = aux().getNLBdPhiBin(dPhiBC, 5, 256);
0196     dTheta = aux().getdTheta(dTheta, 3);
0197     // Identifies which stations have RPC hits in 3-station tracks
0198     int rpc_2b = aux().get2bRPC(clctA, clctB, clctC);  // Have to use un-compressed CLCT words
0199     clctA = aux().getCLCT(clctA, endcap, (sPhiAB == 1 ? 1 : -1), 2);
0200     theta = aux().getTheta(theta, st1_ring2, 5);
0201 
0202     // Form the pT LUT address
0203     address |= (dPhiAB & ((1 << 7) - 1)) << (0);
0204     address |= (dPhiBC & ((1 << 5) - 1)) << (0 + 7);
0205     address |= (sPhiBC & ((1 << 1) - 1)) << (0 + 7 + 5);
0206     address |= (dTheta & ((1 << 3) - 1)) << (0 + 7 + 5 + 1);
0207     address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3);
0208     int bit = 0;
0209     if (mode != 7) {
0210       address |= (frB & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3 + 1);
0211       bit = 1;
0212     }
0213     address |= (clctA & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit);
0214     address |= (rpc_2b & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2);
0215     address |= (theta & ((1 << 5) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2);
0216     if (mode != 7) {
0217       address |= (mode_ID & ((1 << 2) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5);
0218       emtf_assert(address < pow(2, 29) && address >= pow(2, 27));
0219     } else {
0220       address |= (mode_ID & ((1 << 1) - 1)) << (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5);
0221       emtf_assert(address < pow(2, 27) && address >= pow(2, 26));
0222     }
0223 
0224   } else if (nHits == 2) {
0225     int dPhiAB = data.delta_ph[iAB];
0226     int sPhiAB = data.sign_ph[iAB];
0227     int dTheta = data.delta_th[iAB] * (data.sign_th[iAB] ? 1 : -1);
0228     int frA = data.fr[iA];
0229     int frB = data.fr[iB];
0230     int clctA = data.cpattern[iA];
0231     int clctB = data.cpattern[iB];
0232 
0233     // Convert variables to words for pT LUT address
0234     dPhiAB = aux().getNLBdPhiBin(dPhiAB, 7, 512);
0235     dTheta = aux().getdTheta(dTheta, 3);
0236     clctA = aux().getCLCT(clctA, endcap, (sPhiAB == 1 ? 1 : -1), 3);
0237     clctB = aux().getCLCT(clctB, endcap, (sPhiAB == 1 ? 1 : -1), 3);
0238     theta = aux().getTheta(theta, st1_ring2, 5);
0239 
0240     // Form the pT LUT address
0241     address |= (dPhiAB & ((1 << 7) - 1)) << (0);
0242     address |= (dTheta & ((1 << 3) - 1)) << (0 + 7);
0243     address |= (frA & ((1 << 1) - 1)) << (0 + 7 + 3);
0244     address |= (frB & ((1 << 1) - 1)) << (0 + 7 + 3 + 1);
0245     address |= (clctA & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1);
0246     address |= (clctB & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1 + 3);
0247     address |= (theta & ((1 << 5) - 1)) << (0 + 7 + 3 + 1 + 1 + 3 + 3);
0248     address |= (mode_ID & ((1 << 3) - 1)) << (0 + 7 + 3 + 1 + 1 + 3 + 3 + 5);
0249     emtf_assert(address < pow(2, 26) && address >= pow(2, 24));
0250   }
0251   return address;
0252 }  // End function: PtAssignmentEngine2017::calculate_address()
0253 
0254 // Calculate XML pT from address
0255 float PtAssignmentEngine2017::calculate_pt_xml(const address_t& address) const {
0256   // std::cout << "Inside calculate_pt_xml, examining address: ";
0257   // for (int j = 0; j < 30; j++) {
0258   //   std::cout << ((address >> (29 - j)) & 0x1);
0259   //   if ((j % 5) == 4) std::cout << "  ";
0260   // }
0261   // std::cout << std::endl;
0262 
0263   float pt_xml = 0.;
0264   int nHits = -1, mode = -1;
0265 
0266   emtf_assert(address < pow(2, 30));
0267   if (address >= pow(2, 29)) {
0268     nHits = 4;
0269     mode = 15;
0270   } else if (address >= pow(2, 27)) {
0271     nHits = 3;
0272   } else if (address >= pow(2, 26)) {
0273     nHits = 3;
0274     mode = 7;
0275   } else if (address >= pow(2, 24)) {
0276     nHits = 2;
0277   } else
0278     return pt_xml;
0279 
0280   // Variables to unpack from the pT LUT address
0281   int mode_ID, theta, dTheta;
0282   int dPhiAB, dPhiBC = -1, dPhiCD = -1;
0283   int sPhiAB, sPhiBC = -1, sPhiCD = -1;
0284   int frA, frB = -1;
0285   int clctA, clctB = -1;
0286   int rpcA, rpcB, rpcC, rpcD;
0287   int endcap = 1;
0288   sPhiAB = 1;          // Assume positive endcap and dPhiAB for unpacking CLCT bend
0289   int mode15_8b = -1;  // Combines track theta, stations with RPC hits, and station 1 bend information
0290   int rpc_2b = -1;     // Identifies which stations have RPC hits in 3-station tracks
0291 
0292   // Unpack variable words from the pT LUT address
0293   if (nHits == 4) {
0294     dPhiAB = (address >> (0) & ((1 << 7) - 1));
0295     dPhiBC = (address >> (0 + 7) & ((1 << 5) - 1));
0296     dPhiCD = (address >> (0 + 7 + 5) & ((1 << 4) - 1));
0297     sPhiBC = (address >> (0 + 7 + 5 + 4) & ((1 << 1) - 1));
0298     sPhiCD = (address >> (0 + 7 + 5 + 4 + 1) & ((1 << 1) - 1));
0299     dTheta = (address >> (0 + 7 + 5 + 4 + 1 + 1) & ((1 << 2) - 1));
0300     frA = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2) & ((1 << 1) - 1));
0301     mode15_8b = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1) & ((1 << 8) - 1));
0302     mode_ID = (address >> (0 + 7 + 5 + 4 + 1 + 1 + 2 + 1 + 8) & ((1 << 1) - 1));
0303     emtf_assert(address < pow(2, 30));
0304   } else if (nHits == 3) {
0305     dPhiAB = (address >> (0) & ((1 << 7) - 1));
0306     dPhiBC = (address >> (0 + 7) & ((1 << 5) - 1));
0307     sPhiBC = (address >> (0 + 7 + 5) & ((1 << 1) - 1));
0308     dTheta = (address >> (0 + 7 + 5 + 1) & ((1 << 3) - 1));
0309     frA = (address >> (0 + 7 + 5 + 1 + 3) & ((1 << 1) - 1));
0310     int bit = 0;
0311     if (mode != 7) {
0312       frB = (address >> (0 + 7 + 5 + 1 + 3 + 1) & ((1 << 1) - 1));
0313       bit = 1;
0314     }
0315     clctA = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit) & ((1 << 2) - 1));
0316     rpc_2b = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2) & ((1 << 2) - 1));
0317     theta = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2) & ((1 << 5) - 1));
0318     if (mode != 7) {
0319       mode_ID = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5) & ((1 << 2) - 1));
0320       emtf_assert(address < pow(2, 29));
0321     } else {
0322       mode_ID = (address >> (0 + 7 + 5 + 1 + 3 + 1 + bit + 2 + 2 + 5) & ((1 << 1) - 1));
0323       emtf_assert(address < pow(2, 27));
0324     }
0325   } else if (nHits == 2) {
0326     dPhiAB = (address >> (0) & ((1 << 7) - 1));
0327     dTheta = (address >> (0 + 7) & ((1 << 3) - 1));
0328     frA = (address >> (0 + 7 + 3) & ((1 << 1) - 1));
0329     frB = (address >> (0 + 7 + 3 + 1) & ((1 << 1) - 1));
0330     clctA = (address >> (0 + 7 + 3 + 1 + 1) & ((1 << 3) - 1));
0331     clctB = (address >> (0 + 7 + 3 + 1 + 1 + 3) & ((1 << 3) - 1));
0332     theta = (address >> (0 + 7 + 3 + 1 + 1 + 3 + 3) & ((1 << 5) - 1));
0333     mode_ID = (address >> (0 + 7 + 3 + 1 + 1 + 3 + 3 + 5) & ((1 << 3) - 1));
0334     emtf_assert(address < pow(2, 26));
0335   }
0336 
0337   // Infer track mode (and stations with hits) from mode_ID
0338   if (nHits == 3 && mode != 7) {
0339     switch (mode_ID) {
0340       case 0b11:
0341         mode = 14;
0342         break;
0343       case 0b10:
0344         mode = 13;
0345         break;
0346       case 0b01:
0347         mode = 11;
0348         break;
0349       default:
0350         break;
0351     }
0352   } else if (nHits == 2) {
0353     switch (mode_ID) {
0354       case 0b111:
0355         mode = 12;
0356         break;
0357       case 0b110:
0358         mode = 10;
0359         break;
0360       case 0b101:
0361         mode = 9;
0362         break;
0363       case 0b100:
0364         mode = 6;
0365         break;
0366       case 0b011:
0367         mode = 5;
0368         break;
0369       case 0b010:
0370         mode = 3;
0371         break;
0372       default:
0373         break;
0374     }
0375   }
0376 
0377   emtf_assert(mode > 0);
0378 
0379   // Un-compress words from address
0380   // For most variables (e.g. theta, dTheta, CLCT) don't need to unpack, since compressed version was used in training
0381   int St1_ring2 = -1;
0382   if (nHits == 4) {
0383     dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
0384     dPhiBC = aux().getdPhiFromBin(dPhiBC, 5, 256) * (sPhiBC == 1 ? 1 : -1);
0385     dPhiCD = aux().getdPhiFromBin(dPhiCD, 4, 256) * (sPhiCD == 1 ? 1 : -1);
0386     aux().unpack8bMode15(mode15_8b, theta, St1_ring2, endcap, (sPhiAB == 1 ? 1 : -1), clctA, rpcA, rpcB, rpcC, rpcD);
0387 
0388     // // Check bit-wise compression / de-compression
0389     // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 2), 2) );
0390   } else if (nHits == 3) {
0391     dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
0392     dPhiBC = aux().getdPhiFromBin(dPhiBC, 5, 256) * (sPhiBC == 1 ? 1 : -1);
0393     St1_ring2 = aux().unpackSt1Ring2(theta, 5);
0394     aux().unpack2bRPC(rpc_2b, rpcA, rpcB, rpcC);
0395 
0396     // // Check bit-wise compression / de-compression
0397     // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 3), 3) );
0398     // emtf_assert( clctA  == aux().getCLCT( aux().unpackCLCT( clctA, endcap, (sPhiAB == 1 ? 1 : -1), 2),
0399     //                               endcap, (sPhiAB == 1 ? 1 : -1), 2) );
0400     // int theta_unp = theta;
0401     // aux().unpackTheta( theta_unp, St1_ring2, 5 );
0402     // emtf_assert( theta == aux().getTheta(theta_unp, St1_ring2, 5) );
0403   } else if (nHits == 2) {
0404     dPhiAB = aux().getdPhiFromBin(dPhiAB, 7, 512);
0405     St1_ring2 = aux().unpackSt1Ring2(theta, 5);
0406 
0407     // // Check bit-wise compression / de-compression
0408     // emtf_assert( dTheta == aux().getdTheta( aux().unpackdTheta( dTheta, 3), 3) );
0409     // emtf_assert( clctA  == aux().getCLCT( aux().unpackCLCT( clctA, endcap, (sPhiAB == 1 ? 1 : -1), 3),
0410     //                               endcap, (sPhiAB == 1 ? 1 : -1), 3) );
0411     // emtf_assert( clctB  == aux().getCLCT( aux().unpackCLCT( clctB, endcap, (sPhiAB == 1 ? 1 : -1), 3),
0412     //                               endcap, (sPhiAB == 1 ? 1 : -1), 3) );
0413     // int theta_unp = theta;
0414     // aux().unpackTheta( theta_unp, St1_ring2, 5 );
0415     // emtf_assert( theta == aux().getTheta(theta_unp, St1_ring2, 5) );
0416   }
0417 
0418   // Fill vectors of variables for XMLs
0419   // KK: sequence of variables here should exaclty match <Variables> block produced by TMVA
0420   std::vector<int> predictors;
0421 
0422   // Variables for input to XMLs
0423   int dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi;
0424 
0425   // Convert words into variables for XMLs
0426   if (nHits == 4) {
0427     predictors = {
0428         theta, St1_ring2, dPhiAB, dPhiBC, dPhiCD, dPhiAB + dPhiBC, dPhiAB + dPhiBC + dPhiCD, dPhiBC + dPhiCD, frA, clctA};
0429 
0430     aux().calcDeltaPhiSums(dPhiSum4,
0431                            dPhiSum4A,
0432                            dPhiSum3,
0433                            dPhiSum3A,
0434                            outStPhi,
0435                            dPhiAB,
0436                            dPhiAB + dPhiBC,
0437                            dPhiAB + dPhiBC + dPhiCD,
0438                            dPhiBC,
0439                            dPhiBC + dPhiCD,
0440                            dPhiCD);
0441 
0442     int tmp[10] = {dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi, dTheta, rpcA, rpcB, rpcC, rpcD};
0443     predictors.insert(predictors.end(), tmp, tmp + 10);
0444   } else if (nHits == 3) {
0445     if (mode == 14)
0446       predictors = {theta, St1_ring2, dPhiAB, dPhiBC, dPhiAB + dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
0447     else if (mode == 13)
0448       predictors = {theta, St1_ring2, dPhiAB, dPhiAB + dPhiBC, dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
0449     else if (mode == 11)
0450       predictors = {theta, St1_ring2, dPhiBC, dPhiAB, dPhiAB + dPhiBC, frA, frB, clctA, dTheta, rpcA, rpcB, rpcC};
0451     else if (mode == 7)
0452       predictors = {theta, dPhiAB, dPhiBC, dPhiAB + dPhiBC, frA, clctA, dTheta, rpcA, rpcB, rpcC};
0453   } else if (nHits == 2 && mode >= 8) {
0454     predictors = {theta, St1_ring2, dPhiAB, frA, frB, clctA, clctB, dTheta, (clctA == 0), (clctB == 0)};
0455   } else if (nHits == 2 && mode < 8) {
0456     predictors = {theta, dPhiAB, frA, frB, clctA, clctB, dTheta, (clctA == 0), (clctB == 0)};
0457   } else {
0458     emtf_assert(false && "Incorrect nHits or mode");
0459   }
0460 
0461   // Retreive pT from XMLs
0462   std::vector<double> tree_data(predictors.cbegin(), predictors.cend());
0463 
0464   auto tree_event = std::make_unique<emtf::Event>();
0465   tree_event->predictedValue = 0;
0466   tree_event->data = tree_data;
0467 
0468   // forests_.at(mode).predictEvent(tree_event.get(), 400);
0469   emtf::Forest& forest = const_cast<emtf::Forest&>(forests_.at(mode));
0470   forest.predictEvent(tree_event.get(), 400);
0471 
0472   // // Adjust this for different XMLs
0473   if (ptLUTVersion_ >= 8) {  // Run 3 2022 BDT uses log2(pT) target
0474     float log2_pt = tree_event->predictedValue;
0475     pt_xml = pow(2, fmax(0.0, log2_pt));  // Protect against negative values
0476   } else if (ptLUTVersion_ >= 6) {        // Run 2 2017/2018 BDTs use 1/pT target
0477     float inv_pt = tree_event->predictedValue;
0478     pt_xml = 1.0 / fmax(0.001, inv_pt);  // Protect against negative values
0479   }
0480 
0481   return pt_xml;
0482 
0483 }  // End function: float PtAssignmentEngine2017::calculate_pt_xml(const address_t& address)
0484 
0485 // Calculate XML pT directly from track quantities, without forming an address
0486 float PtAssignmentEngine2017::calculate_pt_xml(const EMTFTrack& track) const {
0487   float pt_xml = 0.;
0488 
0489   EMTFPtLUT data = track.PtLUT();
0490 
0491   auto contain = [](const std::vector<int>& vec, int elem) {
0492     return (std::find(vec.begin(), vec.end(), elem) != vec.end());
0493   };
0494 
0495   int endcap = track.Endcap();
0496   int mode = track.Mode();
0497   int theta = track.Theta_fp();
0498   int phi = track.Phi_fp();
0499   if (!contain(allowedModes_, mode))
0500     return pt_xml;
0501 
0502   // Which stations have hits
0503   int st1 = (mode >= 8);
0504   int st2 = ((mode % 8) >= 4);
0505   int st3 = ((mode % 4) >= 2);
0506   int st4 = ((mode % 2) == 1);
0507 
0508   // Variables for input to XMLs
0509   int dPhi_12, dPhi_13, dPhi_14, dPhi_23, dPhi_24, dPhi_34, dPhiSign;
0510   int dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi;
0511   int dTh_12, dTh_13, dTh_14, dTh_23, dTh_24, dTh_34;
0512   int FR_1, FR_2, FR_3, FR_4;
0513   int bend_1, bend_2, bend_3, bend_4;
0514   int RPC_1, RPC_2, RPC_3, RPC_4;
0515   int St1_ring2 = data.st1_ring2;
0516 
0517   int ph1 = -99, ph2 = -99, ph3 = -99, ph4 = -99;
0518   int th1 = -99, th2 = -99, th3 = -99, th4 = -99;
0519   int pat1 = -99, pat2 = -99, pat3 = -99, pat4 = -99;
0520 
0521   // Compute the original phi and theta coordinates
0522   if (st2) {
0523     ph2 = phi;    // Track phi is from station 2 (if it exists), otherwise 3 or 4
0524     th2 = theta;  // Likewise for track theta
0525     if (st1)
0526       ph1 = ph2 - data.delta_ph[0] * (data.sign_ph[0] ? 1 : -1);
0527     if (st1)
0528       th1 = th2 - data.delta_th[0] * (data.sign_th[0] ? 1 : -1);
0529     if (st3)
0530       ph3 = ph2 + data.delta_ph[3] * (data.sign_ph[3] ? 1 : -1);
0531     // Important that phi be from adjacent station pairs (see note below)
0532     if (st3 && st4)
0533       ph4 = ph3 + data.delta_ph[5] * (data.sign_ph[5] ? 1 : -1);
0534     else if (st4)
0535       ph4 = ph2 + data.delta_ph[4] * (data.sign_ph[4] ? 1 : -1);
0536     // Important that theta be from first-last station pair, not adjacent pairs: delta_th values are "best" for each pair, but
0537     // thanks to duplicated CSC LCTs, are not necessarily consistent (or physical) between pairs or between delta_th and delta_ph.
0538     // This is an artifact of the firmware implementation of deltas: see src/AngleCalculation.cc.
0539     if (st1 && st3)
0540       th3 = th1 + data.delta_th[1] * (data.sign_th[1] ? 1 : -1);
0541     else if (st3)
0542       th3 = th2 + data.delta_th[3] * (data.sign_th[3] ? 1 : -1);
0543     if (st1 && st4)
0544       th4 = th1 + data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
0545     else if (st4)
0546       th4 = th2 + data.delta_th[4] * (data.sign_th[4] ? 1 : -1);
0547   } else if (st3) {
0548     ph3 = phi;
0549     th3 = theta;
0550     if (st1)
0551       ph1 = ph3 - data.delta_ph[1] * (data.sign_ph[1] ? 1 : -1);
0552     if (st1)
0553       th1 = th3 - data.delta_th[1] * (data.sign_th[1] ? 1 : -1);
0554     if (st4)
0555       ph4 = ph3 + data.delta_ph[5] * (data.sign_ph[5] ? 1 : -1);
0556     if (st1 && st4)
0557       th4 = th1 + data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
0558     else if (st4)
0559       th4 = th3 + data.delta_th[5] * (data.sign_th[5] ? 1 : -1);
0560   } else if (st4) {
0561     ph4 = phi;
0562     th4 = theta;
0563     if (st1)
0564       ph1 = ph4 - data.delta_ph[2] * (data.sign_ph[2] ? 1 : -1);
0565     if (st1)
0566       th1 = th4 - data.delta_th[2] * (data.sign_th[2] ? 1 : -1);
0567   }
0568 
0569   if (st1)
0570     pat1 = data.cpattern[0];
0571   if (st2)
0572     pat2 = data.cpattern[1];
0573   if (st3)
0574     pat3 = data.cpattern[2];
0575   if (st4)
0576     pat4 = data.cpattern[3];
0577 
0578   // BEGIN: Identical (almost) to BDT training code in EMTFPtAssign2017/PtRegression_Apr_2017.C
0579 
0580   theta = aux().calcTrackTheta(th1, th2, th3, th4, St1_ring2, mode, true);
0581 
0582   aux().calcDeltaPhis(dPhi_12,
0583                       dPhi_13,
0584                       dPhi_14,
0585                       dPhi_23,
0586                       dPhi_24,
0587                       dPhi_34,
0588                       dPhiSign,
0589                       dPhiSum4,
0590                       dPhiSum4A,
0591                       dPhiSum3,
0592                       dPhiSum3A,
0593                       outStPhi,
0594                       ph1,
0595                       ph2,
0596                       ph3,
0597                       ph4,
0598                       mode,
0599                       true);
0600 
0601   aux().calcDeltaThetas(dTh_12, dTh_13, dTh_14, dTh_23, dTh_24, dTh_34, th1, th2, th3, th4, mode, true);
0602 
0603   FR_1 = (st1 ? data.fr[0] : -99);
0604   FR_2 = (st2 ? data.fr[1] : -99);
0605   FR_3 = (st3 ? data.fr[2] : -99);
0606   FR_4 = (st4 ? data.fr[3] : -99);
0607 
0608   aux().calcBends(bend_1, bend_2, bend_3, bend_4, pat1, pat2, pat3, pat4, dPhiSign, endcap, mode, true);
0609 
0610   RPC_1 = (st1 ? (pat1 == 0) : -99);
0611   RPC_2 = (st2 ? (pat2 == 0) : -99);
0612   RPC_3 = (st3 ? (pat3 == 0) : -99);
0613   RPC_4 = (st4 ? (pat4 == 0) : -99);
0614 
0615   aux().calcRPCs(RPC_1, RPC_2, RPC_3, RPC_4, mode, St1_ring2, theta, true);
0616 
0617   // END: Identical (almost) to BDT training code in EMTFPtAssign2017/PtRegression_Apr_2017.C
0618 
0619   // Fill vectors of variables for XMLs
0620   // KK: sequence of variables here should exaclty match <Variables> block produced by TMVA
0621   std::vector<int> predictors;
0622   switch (mode) {
0623     case 15:  // 1-2-3-4
0624       predictors = {theta,    St1_ring2, dPhi_12,  dPhi_23,   dPhi_34,  dPhi_13, dPhi_14, dPhi_24, FR_1,  bend_1,
0625                     dPhiSum4, dPhiSum4A, dPhiSum3, dPhiSum3A, outStPhi, dTh_14,  RPC_1,   RPC_2,   RPC_3, RPC_4};
0626       break;
0627     case 14:  // 1-2-3
0628       predictors = {theta, St1_ring2, dPhi_12, dPhi_23, dPhi_13, FR_1, FR_2, bend_1, dTh_13, RPC_1, RPC_2, RPC_3};
0629       break;
0630     case 13:  // 1-2-4
0631       predictors = {theta, St1_ring2, dPhi_12, dPhi_14, dPhi_24, FR_1, FR_2, bend_1, dTh_14, RPC_1, RPC_2, RPC_4};
0632       break;
0633     case 11:  // 1-3-4
0634       predictors = {theta, St1_ring2, dPhi_34, dPhi_13, dPhi_14, FR_1, FR_3, bend_1, dTh_14, RPC_1, RPC_3, RPC_4};
0635       break;
0636     case 7:  // 2-3-4
0637       predictors = {theta, dPhi_23, dPhi_34, dPhi_24, FR_2, bend_2, dTh_24, RPC_2, RPC_3, RPC_4};
0638       break;
0639     case 12:  // 1-2
0640       predictors = {theta, St1_ring2, dPhi_12, FR_1, FR_2, bend_1, bend_2, dTh_12, RPC_1, RPC_2};
0641       break;
0642     case 10:  // 1-3
0643       predictors = {theta, St1_ring2, dPhi_13, FR_1, FR_3, bend_1, bend_3, dTh_13, RPC_1, RPC_3};
0644       break;
0645     case 9:  // 1-4
0646       predictors = {theta, St1_ring2, dPhi_14, FR_1, FR_4, bend_1, bend_4, dTh_14, RPC_1, RPC_4};
0647       break;
0648     case 6:  // 2-3
0649       predictors = {theta, dPhi_23, FR_2, FR_3, bend_2, bend_3, dTh_23, RPC_2, RPC_3};
0650       break;
0651     case 5:  // 2-4
0652       predictors = {theta, dPhi_24, FR_2, FR_4, bend_2, bend_4, dTh_24, RPC_2, RPC_4};
0653       break;
0654     case 3:  // 3-4
0655       predictors = {theta, dPhi_34, FR_3, FR_4, bend_3, bend_4, dTh_34, RPC_3, RPC_4};
0656       break;
0657   }
0658 
0659   std::vector<double> tree_data(predictors.cbegin(), predictors.cend());
0660 
0661   auto tree_event = std::make_unique<emtf::Event>();
0662   tree_event->predictedValue = 0;
0663   tree_event->data = tree_data;
0664 
0665   // forests_.at(mode).predictEvent(tree_event.get(), 400);
0666   emtf::Forest& forest = const_cast<emtf::Forest&>(forests_.at(mode));
0667   forest.predictEvent(tree_event.get(), 400);
0668 
0669   // // Adjust this for different XMLs
0670   if (ptLUTVersion_ >= 8) {  // Run 3 2022 BDT uses log2(pT) target
0671     float log2_pt = tree_event->predictedValue;
0672     pt_xml = pow(2, fmax(0.0, log2_pt));  // Protect against negative values
0673   } else if (ptLUTVersion_ >= 6) {        // Run 2 2017/2018 BDTs use 1/pT target
0674     float inv_pt = tree_event->predictedValue;
0675     pt_xml = 1.0 / fmax(0.001, inv_pt);  // Protect against negative values
0676   }
0677 
0678   return pt_xml;
0679 
0680 }  // End function: float PtAssignmentEngine2017::calculate_pt_xml(const EMTFTrack& track)
[ 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!