test/deprecated/CSCCathodeLCTAnalyzer.cc

0001 #include "FWCore/Utilities/interface/Exception.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0004 #include "Geometry/CSCGeometry/interface/OffsetRadialStripTopology.h"
0005 #include "L1Trigger/CSCTriggerPrimitives/plugins/CSCCathodeLCTAnalyzer.h"
0006
0007 using namespace std;
0008
0009 //-----------------
0010 // Static variables
0011 //-----------------
0012
0013 bool CSCCathodeLCTAnalyzer::debug = true;
0014 bool CSCCathodeLCTAnalyzer::doME1A = false;
0015
0016 vector<CSCCathodeLayerInfo> CSCCathodeLCTAnalyzer::getSimInfo(const CSCCLCTDigi& clct,
0017                                                               const CSCDetId& clctId,
0018                                                               const CSCComparatorDigiCollection* compdc,
0019                                                               const edm::PSimHitContainer* allSimHits) {
0020   // Fills vector of CSCCathodeLayerInfo objects.  There can be up to 6 such
0021   // objects (one per layer); they contain the list of comparator digis used to
0022   // build a given CLCT, and the list of associated (closest) SimHits.
0023   // Filling is done in two steps: first, we construct the list of comparator
0024   // digis; next, find associated SimHits.
0025   vector<CSCCathodeLayerInfo> clctInfo = lctDigis(clct, clctId, compdc);
0026
0027   // Sanity checks.
0028   if (clctInfo.size() > CSCConstants::NUM_LAYERS) {
0029     throw cms::Exception("CSCCathodeLCTAnalyzer") << "+++ Number of CSCCathodeLayerInfo objects, " << clctInfo.size()
0030                                                   << ", exceeds max expected, " << CSCConstants::NUM_LAYERS << " +++\n";
0031   }
0032   //  not a good check for high PU
0033   //if (clctInfo.size() != (unsigned)clct.getQuality()) {
0034   //  edm::LogWarning("L1CSCTPEmulatorWrongValues")
0035   //    << "+++ Warning: mismatch between CLCT quality, " << clct.getQuality()
0036   //    << ", and the number of layers with digis, " << clctInfo.size()
0037   //    << ", in clctInfo! +++\n";
0038   //}
0039
0040   // Find the closest SimHit to each Digi.
0041   vector<CSCCathodeLayerInfo>::iterator pcli;
0042   for (pcli = clctInfo.begin(); pcli != clctInfo.end(); pcli++) {
0043     digiSimHitAssociator(*pcli, allSimHits);
0044   }
0045
0046   return clctInfo;
0047 }
0048
0049 vector<CSCCathodeLayerInfo> CSCCathodeLCTAnalyzer::lctDigis(const CSCCLCTDigi& clct,
0050                                                             const CSCDetId& clctId,
0051                                                             const CSCComparatorDigiCollection* compdc) {
0052   // Function to find a list of ComparatorDigis used to create an LCT.
0053   // The list of ComparatorDigis is stored in a class called CSCLayerInfo which
0054   // contains the layerId's of the stored ComparatorDigis as well as the actual
0055   // digis themselves.
0056   int hfstripDigis[CSCConstants::MAX_NUM_HALF_STRIPS_RUN1_TRIGGER];
0057   int time[CSCConstants::MAX_NUM_STRIPS_RUN1], comp[CSCConstants::MAX_NUM_STRIPS_RUN1];
0058   int digiNum[CSCConstants::MAX_NUM_STRIPS_RUN1];
0059   int digiId = -999;
0060   CSCCathodeLayerInfo tempInfo;
0061   vector<CSCCathodeLayerInfo> vectInfo;
0062
0063   // Inquire the clct for its key half-strip, strip type and pattern number.
0064   int clct_keystrip = clct.getKeyStrip();
0065   int clct_stripType = clct.getStripType();
0066   int clct_pattern = clct.getPattern();
0067   int clct_bx = clct.getBX();
0068
0069   // Re-scale the key di-strip number to the count used in the search.
0070   if (clct_stripType == 0)
0071     clct_keystrip /= 4;
0072
0073   if (debug) {
0074     char stype = (clct_stripType == 0) ? 'D' : 'H';
0075     LogDebug("lctDigis") << "\nlctDigis: CLCT keystrip = " << clct_keystrip << " (" << stype << ")"
0076                          << " pattern = " << clct_pattern << "; Id:"
0077                          << " endcap " << clctId.endcap() << ", station " << clctId.station() << ", ring "
0078                          << clctId.ring() << ", chamber " << clctId.chamber();
0079   }
0080
0081   // 'Staggering' for key layer.  Needed for TMB07 version.
0082   const int key_layer = 3;  // counting from 1.
0083   CSCDetId layerId(clctId.endcap(), clctId.station(), clctId.ring(), clctId.chamber(), key_layer);
0084   const CSCLayer* csclayer = geom_->layer(layerId);
0085   int key_stagger = (csclayer->geometry()->stagger() + 1) / 2;
0086
0087   for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; i_layer++) {
0088     // Clear tempInfo values every iteration before using.
0089     tempInfo.clear();
0090
0091     // @ Switch to maps eventually
0092     vector<CSCComparatorDigi> digiMap;
0093     int digi_num = 0;
0094     for (int i_hstrip = 0; i_hstrip < CSCConstants::MAX_NUM_HALF_STRIPS_RUN1_TRIGGER; i_hstrip++) {
0095       hfstripDigis[i_hstrip] = -999;
0096     }
0097     for (int i_strip = 0; i_strip < CSCConstants::MAX_NUM_STRIPS_RUN1; i_strip++) {
0098       time[i_strip] = -999;
0099       comp[i_strip] = 0;
0100       digiNum[i_strip] = -999;
0101     }
0102
0103     // CLCTs belong to a chamber, so their layerId is 0.  Comparator digis
0104     // belong to layers, so in order to access them we need to add layer
0105     // number to CLCT's endcap, chamber, etc.
0106     CSCDetId layerId(clctId.endcap(), clctId.station(), clctId.ring(), clctId.chamber(), i_layer + 1);
0107
0108     // Preselection of Digis: right layer and bx.
0109     digi_num += preselectDigis(clct_bx, layerId, compdc, digiMap, hfstripDigis, time, comp, digiNum);
0110
0111     // In case of ME1/1, one can also look for digis in ME1/A.
0112     // Skip them for now since the resolution of CLCTs in ME1/A is
0113     // terrible (strips are ganged; channel numbers translated to be
0114     // in CFEB=4).
0115     if (doME1A) {
0116       if (clctId.station() == 1 && clctId.ring() == 1) {
0117         CSCDetId layerId_me1a(clctId.endcap(), clctId.station(), 4, clctId.chamber(), i_layer + 1);
0118         digi_num += preselectDigis(clct_bx, layerId_me1a, compdc, digiMap, hfstripDigis, time, comp, digiNum);
0119       }
0120     }
0121
0122     // Loop over all the strips in a pattern.
0123     // FIXME !!!
0124     int strip;
0125     for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; i_layer++) {
0126       for (int i_strip = 0; i_strip < CSCConstants::CLCT_PATTERN_WIDTH; i_strip++) {
0127         strip = clct_keystrip + key_stagger + CSCPatternBank::clct_pattern_offset_[i_strip];
0128         if (strip >= 0 && strip < CSCConstants::MAX_NUM_HALF_STRIPS_RUN1_TRIGGER) {
0129           digiId = hfstripDigis[strip];
0130           // halfstripDigis contains the digi numbers
0131           // that were carried through the different transformations
0132           // to keystrip. 0 means there was no Digi.
0133           if (digiId >= 0) {
0134             tempInfo.setId(layerId);                 // store the layer of this object
0135             tempInfo.addComponent(digiMap[digiId]);  // and the RecDigi
0136             if (debug) {
0137               LogTrace("lctDigis") << " Digi on CLCT: strip/comp/time " << digiMap[digiId];
0138             }
0139           }
0140         }
0141       }
0142     }
0143
0144     // Save results for each non-empty layer.
0145     if (tempInfo.getId().layer() != 0) {
0146       vectInfo.push_back(tempInfo);
0147     }
0148   }
0149   return vectInfo;
0150 }
0151
0152 int CSCCathodeLCTAnalyzer::preselectDigis(const int clct_bx,
0153                                           const CSCDetId& layerId,
0154                                           const CSCComparatorDigiCollection* compdc,
0155                                           vector<CSCComparatorDigi>& digiMap,
0156                                           int hfstripDigis[CSCConstants::MAX_NUM_HALF_STRIPS_RUN1_TRIGGER],
0157                                           int time[CSCConstants::MAX_NUM_STRIPS_RUN1],
0158                                           int comp[CSCConstants::MAX_NUM_STRIPS_RUN1],
0159                                           int digiNum[CSCConstants::MAX_NUM_STRIPS_RUN1]) {
0160   // Preselection of Digis: right layer and bx.
0161   int digi_num = 0;
0162
0163   // Parameters defining time window for accepting hits; should come from
0164   // configuration file eventually.
0165   const int fifo_tbins = 12;
0166   const int hit_persist = 4;
0167   const int drift_delay = 2;
0168
0169   // 'Staggering' for this layer.
0170   const CSCLayer* csclayer = geom_->layer(layerId);
0171   int stagger = (csclayer->geometry()->stagger() + 1) / 2;
0172
0173   bool me1a = (layerId.station() == 1) && (layerId.ring() == 4);
0174
0175   const CSCComparatorDigiCollection::Range rcompd = compdc->get(layerId);
0176   for (CSCComparatorDigiCollection::const_iterator digiIt = rcompd.first; digiIt != rcompd.second; ++digiIt) {
0177     if (debug)
0178       LogDebug("lctDigis") << "Comparator digi: layer " << layerId.layer() - 1
0179                            << " strip/comparator/time =" << (*digiIt);
0180
0181     if ((*digiIt).getComparator() == 0 || (*digiIt).getComparator() == 1) {
0182       int bx_time = (*digiIt).getTimeBin();
0183       if (bx_time >= 0 && bx_time < fifo_tbins) {
0184         // Do not use digis which could not have contributed to a given CLCT.
0185         int latch_bx = clct_bx + drift_delay;
0186         if (bx_time <= latch_bx - hit_persist || bx_time > latch_bx) {
0187           if (debug)
0188             LogDebug("lctDigis") << "Late comparator digi: layer " << layerId.layer() - 1
0189                                  << " strip/comparator/time =" << (*digiIt) << " skipping...";
0190           continue;
0191         }
0192
0193         // If there is more than one digi on the same strip, pick the one
0194         // which occurred earlier.
0195         int i_strip = (*digiIt).getStrip() - 1;  // starting from 0
0196         if (me1a && i_strip < 16) {
0197           // Move ME1/A comparators from CFEB=0 to CFEB=4 if this has not
0198           // been done already.
0199           i_strip += 64;
0200         }
0201
0202         if (time[i_strip] <= 0 || time[i_strip] > bx_time) {
0203           // @ Switch to maps; check for match in time.
0204           int i_hfstrip = 2 * i_strip + (*digiIt).getComparator() + stagger;
0205           hfstripDigis[i_hfstrip] = digi_num;
0206
0207           // Arrays for distrip stagger
0208           comp[i_strip] = (*digiIt).getComparator();
0209           time[i_strip] = bx_time;
0210           digiNum[i_strip] = digi_num;
0211
0212           if (debug)
0213             LogDebug("lctDigis") << "digi_num = " << digi_num << " half-strip = " << i_hfstrip
0214                                  << " strip = " << i_strip;
0215         }
0216       }
0217     }
0218     digiMap.push_back(*digiIt);
0219     digi_num++;
0220   }
0221
0222   return digi_num;
0223 }
0224
0225 void CSCCathodeLCTAnalyzer::digiSimHitAssociator(CSCCathodeLayerInfo& info, const edm::PSimHitContainer* allSimHits) {
0226   // This routine matches up the closest simHit to every digi on a given layer.
0227   // Note that it is possible to have up to 2 digis contribute to an LCT
0228   // on a given layer. In a primitive algorithm used now more than one digi on
0229   // a layer can be associated with the same simHit.
0230   vector<PSimHit> simHits;
0231
0232   vector<CSCComparatorDigi> thisLayerDigis = info.getRecDigis();
0233   if (!thisLayerDigis.empty()) {
0234     CSCDetId layerId = info.getId();
0235     bool me11 = (layerId.station() == 1) && (layerId.ring() == 1);
0236
0237     // Get simHits in this layer.
0238     for (edm::PSimHitContainer::const_iterator simHitIt = allSimHits->begin(); simHitIt != allSimHits->end();
0239          simHitIt++) {
0240       // Find detId where simHit is located.
0241       CSCDetId hitId = (CSCDetId)(*simHitIt).detUnitId();
0242       if (hitId == layerId)
0243         simHits.push_back(*simHitIt);
0244       if (me11) {
0245         CSCDetId layerId_me1a(layerId.endcap(), layerId.station(), 4, layerId.chamber(), layerId.layer());
0246         if (hitId == layerId_me1a)
0247           simHits.push_back(*simHitIt);
0248       }
0249     }
0250
0251     if (!simHits.empty()) {
0252       ostringstream strstrm;
0253       if (debug) {
0254         strstrm << "\nLayer " << layerId.layer() << " has " << simHits.size() << " SimHit(s); phi value(s) = ";
0255       }
0256
0257       // Get the strip number for every digi and convert to phi.
0258       for (vector<CSCComparatorDigi>::iterator prd = thisLayerDigis.begin(); prd != thisLayerDigis.end(); prd++) {
0259         double deltaPhiMin = 999.;
0260         double bestHitPhi = 999.;
0261         PSimHit* bestHit = nullptr;
0262
0263         int strip = prd->getStrip();
0264         double digiPhi = getStripPhi(layerId, strip - 0.5);
0265
0266         const CSCLayer* csclayer = geom_->layer(layerId);
0267         for (vector<PSimHit>::iterator psh = simHits.begin(); psh != simHits.end(); psh++) {
0268           // Get the local phi for the simHit.
0269           LocalPoint hitLP = psh->localPosition();
0270           GlobalPoint hitGP = csclayer->toGlobal(hitLP);
0271           double hitPhi = hitGP.phi();
0272           if (debug)
0273             strstrm << hitPhi << " ";
0274           // Find the lowest deltaPhi and store the associated simHit.
0275           // Check to see if comparison is on pi border, and if so, make
0276           // corrections.
0277           double deltaPhi = 999.;
0278           if (fabs(hitPhi - digiPhi) < M_PI) {
0279             deltaPhi = fabs(hitPhi - digiPhi);
0280           } else {
0281             if (debug)
0282               LogDebug("digiSimHitAssociator") << "SimHit and Digi are close to edge!!!";
0283             deltaPhi = fabs(fabs(hitPhi - digiPhi) - 2. * M_PI);
0284           }
0285           if (deltaPhi < deltaPhiMin) {
0286             deltaPhiMin = deltaPhi;
0287             bestHit = &(*psh);
0288             bestHitPhi = hitPhi;
0289           }
0290         }
0291         if (debug) {
0292           strstrm << "\nDigi Phi: " << digiPhi << ", closest SimHit phi: " << bestHitPhi
0293                   << ", particle type: " << bestHit->particleType();
0294         }
0295         info.addComponent(*bestHit);
0296       }
0297       if (debug) {
0298         LogTrace("digiSimHitAssociator") << strstrm.str();
0299       }
0300     }
0301   }
0302 }
0303
0304 int CSCCathodeLCTAnalyzer::nearestHS(const vector<CSCCathodeLayerInfo>& allLayerInfo,
0305                                      double& closestPhi,
0306                                      double& closestEta) {
0307   // Function to set the simulated values for comparison to the reconstructed.
0308   // It first tries to look for the SimHit in the key layer.  If it is
0309   // unsuccessful, it loops over all layers and looks for an associated
0310   // hit in any one of the layers.  First instance of a hit gives a calculation
0311   // for phi.
0312   int nearestHS = -999;
0313   PSimHit matchedHit;
0314   bool hit_found = false;
0315   CSCDetId layerId;
0316   static const int key_layer = CSCConstants::KEY_CLCT_LAYER;
0317
0318   vector<CSCCathodeLayerInfo>::const_iterator pli;
0319   for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0320     if (pli->getId().layer() == key_layer) {
0321       vector<PSimHit> thisLayerHits = pli->getSimHits();
0322       if (!thisLayerHits.empty()) {
0323         // There can be one RecDigi (and therefore only one SimHit)
0324         // in a keylayer.
0325         if (thisLayerHits.size() != 1) {
0326           edm::LogWarning("L1CSCTPEmulatorWrongValues")
0327               << "+++ Warning: " << thisLayerHits.size() << " SimHits in key layer " << key_layer << "! +++ \n";
0328           for (unsigned i = 0; i < thisLayerHits.size(); i++) {
0329             edm::LogWarning("L1CSCTPEmulatorWrongValues") << " SimHit # " << i << thisLayerHits[i] << "\n";
0330           }
0331         }
0332         matchedHit = thisLayerHits[0];
0333         layerId = pli->getId();
0334         hit_found = true;
0335         break;
0336       }
0337     }
0338   }
0339
0340   if (hit_found == false) {
0341     for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0342       // if there is any occurrence of simHit size greater that zero, use this.
0343       if (!(pli->getRecDigis()).empty() && !(pli->getSimHits()).empty()) {
0344         // Always use the first SimHit for now.
0345         matchedHit = (pli->getSimHits())[0];
0346         layerId = pli->getId();
0347         hit_found = true;
0348         break;
0349       }
0350     }
0351   }
0352
0353   // Set phi and eta if there were any hits found.
0354   if (hit_found) {
0355     const CSCLayer* csclayer = geom_->layer(layerId);
0356     const CSCLayerGeometry* layerGeom = csclayer->geometry();
0357     //nearestStrip = layerGeom->nearestStrip(matchedHit.localPosition());
0358     // Float in units of the strip (angular) width.  From RadialStripTopology
0359     // comments: "Strip in which a given LocalPoint lies. This is a float
0360     // which represents the fractional strip position within the detector.
0361     // Returns zero if the LocalPoint falls at the extreme low edge of the
0362     // detector or BELOW, and float(nstrips) if it falls at the extreme high
0363     // edge or ABOVE."
0364     float strip = layerGeom->topology()->strip(matchedHit.localPosition());
0365
0366     // Should be in the interval [0-MAX_STRIPS).  I see (rarely) cases when
0367     // strip = nearestStrip = MAX_STRIPS; do not know how to handle them.
0368     int nearestStrip = static_cast<int>(strip);
0369     if (nearestStrip < 0 || nearestStrip >= CSCConstants::MAX_NUM_STRIPS_RUN1) {
0370       edm::LogWarning("L1CSCTPEmulatorWrongInput")
0371           << "+++ Warning: nearest strip, " << nearestStrip << ", is not in [0-" << CSCConstants::MAX_NUM_STRIPS_RUN1
0372           << ") interval; strip = " << strip << " +++\n";
0373     }
0374     // Left/right half of the strip.
0375     int comp = ((strip - nearestStrip) < 0.5) ? 0 : 1;
0376     nearestHS = 2 * nearestStrip + comp;
0377
0378     GlobalPoint thisPoint = csclayer->toGlobal(matchedHit.localPosition());
0379     closestPhi = thisPoint.phi();
0380     closestEta = thisPoint.eta();
0381     ostringstream strstrm;
0382     if (debug)
0383       strstrm << "Matched cathode phi: " << closestPhi;
0384     if (closestPhi < 0.) {
0385       closestPhi += 2. * M_PI;
0386       if (debug)
0387         strstrm << "(" << closestPhi << ")";
0388     }
0389     if (debug) {
0390       strstrm << " eta: " << closestEta << " on a layer " << layerId.layer() << " (1-6);"
0391               << " nearest strip: " << nearestStrip;
0392       LogDebug("nearestHS") << strstrm.str();
0393     }
0394   }
0395
0396   return nearestHS;
0397 }
0398
0399 void CSCCathodeLCTAnalyzer::setGeometry(const CSCGeometry* geom) { geom_ = geom; }
0400
0401 double CSCCathodeLCTAnalyzer::getStripPhi(const CSCDetId& layerId, const float strip) {
0402   // Returns phi position of a given strip.
0403   if (strip < 0. || strip >= CSCConstants::MAX_NUM_STRIPS_RUN1) {
0404     edm::LogWarning("L1CSCTPEmulatorWrongInput") << "+++ Warning: strip, " << strip << ", is not in [0-"
0405                                                  << CSCConstants::MAX_NUM_STRIPS_RUN1 << ") interval +++\n";
0406   }
0407
0408   const CSCLayer* csclayer = geom_->layer(layerId);
0409   const CSCLayerGeometry* layerGeom = csclayer->geometry();
0410
0411   // Position at the center of the strip.
0412   LocalPoint digiLP = layerGeom->topology()->localPosition(strip);
0413   // The alternative calculation gives exactly the same answer.
0414   // double ystrip = 0.0;
0415   // double xstrip = topology->xOfStrip(strip, ystrip);
0416   // LocalPoint  digiLP = LocalPoint(xstrip, ystrip);
0417   GlobalPoint digiGP = csclayer->toGlobal(digiLP);
0418   double phi = digiGP.phi();
0419
0420   return phi;
0421 }