Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-13 01:32:13 

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