test/deprecated/CSCAnodeLCTAnalyzer.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/CSCLayer.h"
0005 #include "L1Trigger/CSCTriggerPrimitives/plugins/CSCAnodeLCTAnalyzer.h"
0006 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0007
0008 using namespace std;
0009
0010 //-----------------
0011 // Static variables
0012 //-----------------
0013
0014 bool CSCAnodeLCTAnalyzer::debug = true;
0015 bool CSCAnodeLCTAnalyzer::doME1A = true;
0016
0017 vector<CSCAnodeLayerInfo> CSCAnodeLCTAnalyzer::getSimInfo(const CSCALCTDigi& alct,
0018                                                           const CSCDetId& alctId,
0019                                                           const CSCWireDigiCollection* wiredc,
0020                                                           const edm::PSimHitContainer* allSimHits) {
0021   // Fills vector of CSCAnodeLayerInfo objects.  There can be up to 6 such
0022   // objects (one per layer); they contain the list of wire digis used to
0023   // build a given ALCT, and the list of associated (closest) SimHits.
0024   // Filling is done in two steps: first, we construct the list of wire
0025   // digis; next, find associated SimHits.
0026   vector<CSCAnodeLayerInfo> alctInfo = lctDigis(alct, alctId, wiredc);
0027
0028   // Sanity checks.
0029   if (alctInfo.size() > CSCConstants::NUM_LAYERS) {
0030     throw cms::Exception("CSCAnodeLCTAnalyzer") << "+++ Number of CSCAnodeLayerInfo objects, " << alctInfo.size()
0031                                                 << ", exceeds max expected, " << CSCConstants::NUM_LAYERS << " +++\n";
0032   }
0033   // not a good check for high PU
0034   //if (alctInfo.size() != (unsigned)alct.getQuality()+3) {
0035   //  edm::LogWarning("L1CSCTPEmulatorWrongValues")
0036   //    << "+++ Warning: mismatch between ALCT quality, " << alct.getQuality()
0037   //    << ", and the number of layers with digis, " << alctInfo.size()
0038   //    << ", in alctInfo! +++\n";
0039   //}
0040
0041   // Find the closest SimHit to each Digi.
0042   vector<CSCAnodeLayerInfo>::iterator pali;
0043   for (pali = alctInfo.begin(); pali != alctInfo.end(); pali++) {
0044     digiSimHitAssociator(*pali, allSimHits);
0045   }
0046
0047   return alctInfo;
0048 }
0049
0050 vector<CSCAnodeLayerInfo> CSCAnodeLCTAnalyzer::lctDigis(const CSCALCTDigi& alct,
0051                                                         const CSCDetId& alctId,
0052                                                         const CSCWireDigiCollection* wiredc) {
0053   // Function to find a list of WireDigis used to create an LCT.
0054   // The list of WireDigis is stored in a class called CSCLayerInfo which
0055   // contains the layerId's of the stored WireDigis as well as the actual digis
0056   // themselves.
0057   CSCAnodeLayerInfo tempInfo;
0058   vector<CSCAnodeLayerInfo> vectInfo;
0059
0060   // Inquire the alct for its pattern and key wiregroup.
0061   int alct_pattern = 0;
0062   if (!alct.getAccelerator())
0063     alct_pattern = alct.getCollisionB() + 1;
0064   int alct_keywire = alct.getKeyWG();
0065   int alct_bx = alct.getBX();
0066
0067   // Choose pattern envelope.
0068   int MESelection = (alctId.station() < 3) ? 0 : 1;
0069
0070   if (debug) {
0071     LogDebug("lctDigis") << "\nlctDigis: ALCT keywire = " << alct_keywire << "; alctId:"
0072                          << " endcap " << alctId.endcap() << ", station " << alctId.station() << ", ring "
0073                          << alctId.ring() << ", chamber " << alctId.chamber();
0074   }
0075
0076   for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; i_layer++) {
0077     map<int, CSCWireDigi> digiMap;
0078     // Clear tempInfo values every iteration before using.
0079     tempInfo.clear();
0080
0081     // ALCTs belong to a chamber, so their layerId is 0.  Wire digis belong
0082     // to layers, so in order to access them we need to add layer number to
0083     // ALCT's endcap, chamber, etc.
0084     CSCDetId layerId(alctId.endcap(), alctId.station(), alctId.ring(), alctId.chamber(), i_layer + 1);
0085     // Preselection of Digis: right layer and bx.
0086     preselectDigis(alct_bx, layerId, wiredc, digiMap);
0087
0088     // In case of ME1/1, one can also look for digis in ME1/A.
0089     // Keep "on" by defailt since the resolution should not be different
0090     // from that in ME1/B.
0091     if (doME1A) {
0092       if (alctId.station() == 1 && alctId.ring() == 1) {
0093         CSCDetId layerId_me1a(alctId.endcap(), alctId.station(), 4, alctId.chamber(), i_layer + 1);
0094         preselectDigis(alct_bx, layerId_me1a, wiredc, digiMap);
0095       }
0096     }
0097
0098     // Loop over all the wires in a pattern.
0099     for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; ++i_layer) {
0100       for (int i_wire = 0; i_wire < CSCConstants::ALCT_PATTERN_WIDTH; ++i_wire) {
0101         // check the mask
0102         if (CSCPatternBank::alct_pattern_legacy_[alct_pattern][i_layer][i_wire]) {
0103           int wire = alct_keywire + CSCPatternBank::alct_keywire_offset_[MESelection][i_wire];
0104           if (wire >= 0 && wire < CSCConstants::MAX_NUM_WIREGROUPS) {
0105             // Check if there is a "good" Digi on this wire.
0106             if (digiMap.count(wire) > 0) {
0107               tempInfo.setId(layerId);               // store the layer of this object
0108               tempInfo.addComponent(digiMap[wire]);  // and the RecDigi
0109               if (debug)
0110                 LogDebug("lctDigis") << " Digi on ALCT: wire group " << digiMap[wire].getWireGroup();
0111             }
0112           }
0113         }
0114       }
0115     }
0116
0117     // Save results for each non-empty layer.
0118     if (tempInfo.getId().layer() != 0) {
0119       vectInfo.push_back(tempInfo);
0120     }
0121   }
0122   return vectInfo;
0123 }
0124
0125 void CSCAnodeLCTAnalyzer::preselectDigis(const int alct_bx,
0126                                          const CSCDetId& layerId,
0127                                          const CSCWireDigiCollection* wiredc,
0128                                          map<int, CSCWireDigi>& digiMap) {
0129   // Preselection of Digis: right layer and bx.
0130
0131   // Parameters defining time window for accepting hits; should come from
0132   // configuration file eventually.
0133   const int fifo_tbins = 16;
0134   const int drift_delay = 2;
0135   const int hit_persist = 6;  // not a config. parameter, just const
0136
0137   const CSCWireDigiCollection::Range rwired = wiredc->get(layerId);
0138   for (CSCWireDigiCollection::const_iterator digiIt = rwired.first; digiIt != rwired.second; ++digiIt) {
0139     if (debug)
0140       LogDebug("lctDigis") << "Wire digi: layer " << layerId.layer() - 1 << (*digiIt);
0141     int bx_time = (*digiIt).getTimeBin();
0142     if (bx_time >= 0 && bx_time < fifo_tbins) {
0143       // Do not use digis which could not have contributed to a given ALCT.
0144       int latch_bx = alct_bx + drift_delay;
0145       if (bx_time <= latch_bx - hit_persist || bx_time > latch_bx) {
0146         if (debug)
0147           LogDebug("lctDigis") << "Late wire digi: layer " << layerId.layer() - 1 << " " << (*digiIt) << " skipping...";
0148         continue;
0149       }
0150
0151       int i_wire = (*digiIt).getWireGroup() - 1;
0152
0153       // If there is more than one digi on the same wire, pick the one
0154       // which occurred earlier.
0155       if (digiMap.count(i_wire) > 0) {
0156         if (digiMap[i_wire].getTimeBin() > bx_time) {
0157           if (debug) {
0158             LogDebug("lctDigis") << " Replacing good wire digi on wire " << i_wire;
0159           }
0160           digiMap.erase(i_wire);
0161         }
0162       }
0163
0164       digiMap[i_wire] = *digiIt;
0165       if (debug) {
0166         LogDebug("lctDigis") << " Good wire digi: wire group " << i_wire;
0167       }
0168     }
0169   }
0170 }
0171
0172 void CSCAnodeLCTAnalyzer::digiSimHitAssociator(CSCAnodeLayerInfo& info, const edm::PSimHitContainer* allSimHits) {
0173   // This routine matches up the closest simHit to every digi on a given layer.
0174   // Iit is possible to have up to 3 digis contribute to an LCT on a given
0175   // layer.  In a primitive algorithm used now more than one digi on a layer
0176   // can be associated with the same simHit.
0177   vector<PSimHit> simHits;
0178
0179   vector<CSCWireDigi> thisLayerDigis = info.getRecDigis();
0180   if (!thisLayerDigis.empty()) {
0181     CSCDetId layerId = info.getId();
0182     bool me11 = (layerId.station() == 1) && (layerId.ring() == 1);
0183
0184     // Get simHits in this layer.
0185     for (edm::PSimHitContainer::const_iterator simHitIt = allSimHits->begin(); simHitIt != allSimHits->end();
0186          simHitIt++) {
0187       // Find detId where simHit is located.
0188       CSCDetId hitId = (CSCDetId)(*simHitIt).detUnitId();
0189       if (hitId == layerId)
0190         simHits.push_back(*simHitIt);
0191       if (me11) {
0192         CSCDetId layerId_me1a(layerId.endcap(), layerId.station(), 4, layerId.chamber(), layerId.layer());
0193         if (hitId == layerId_me1a)
0194           simHits.push_back(*simHitIt);
0195       }
0196     }
0197
0198     if (!simHits.empty()) {
0199       ostringstream strstrm;
0200       if (debug) {
0201         strstrm << "\nLayer " << layerId.layer() << " has " << simHits.size() << " SimHit(s); eta value(s) = ";
0202       }
0203
0204       // Get the wire number for every digi and convert to eta.
0205       for (vector<CSCWireDigi>::const_iterator prd = thisLayerDigis.begin(); prd != thisLayerDigis.end(); prd++) {
0206         double deltaEtaMin = 999.;
0207         double bestHitEta = 999.;
0208         PSimHit* bestHit = nullptr;
0209
0210         int wiregroup = prd->getWireGroup();  // counted from 1
0211         double digiEta = getWGEta(layerId, wiregroup - 1);
0212
0213         const CSCLayer* csclayer = geom_->layer(layerId);
0214         for (vector<PSimHit>::iterator psh = simHits.begin(); psh != simHits.end(); psh++) {
0215           // Get the local eta for the simHit.
0216           LocalPoint hitLP = psh->localPosition();
0217           GlobalPoint hitGP = csclayer->toGlobal(hitLP);
0218           double hitEta = hitGP.eta();
0219           if (debug)
0220             strstrm << hitEta << " ";
0221           // Find the lowest deltaEta and store the associated simHit.
0222           double deltaEta = fabs(hitEta - digiEta);
0223           if (deltaEta < deltaEtaMin) {
0224             deltaEtaMin = deltaEta;
0225             bestHit = &(*psh);
0226             bestHitEta = hitEta;
0227           }
0228         }
0229         if (debug) {
0230           strstrm << "\nDigi eta: " << digiEta << ", closest SimHit eta: " << bestHitEta
0231                   << ", particle type: " << bestHit->particleType();
0232           //strstrm << "\nlocal position:" << bestHit->localPosition();
0233         }
0234         info.addComponent(*bestHit);
0235       }
0236       if (debug) {
0237         LogDebug("digiSimHitAssociator") << strstrm.str();
0238       }
0239     }
0240   }
0241 }
0242
0243 int CSCAnodeLCTAnalyzer::nearestWG(const vector<CSCAnodeLayerInfo>& allLayerInfo,
0244                                    double& closestPhi,
0245                                    double& closestEta) {
0246   // Function to set the simulated values for comparison to the reconstructed.
0247   // It first tries to look for the SimHit in the key layer.  If it is
0248   // unsuccessful, it loops over all layers and looks for an associated
0249   // hit in any one of the layers.  First instance of a hit gives a calculation
0250   // for eta.
0251   int nearestWG = -999;
0252   PSimHit matchedHit;
0253   bool hit_found = false;
0254   CSCDetId layerId;
0255
0256   vector<CSCAnodeLayerInfo>::const_iterator pli;
0257   for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0258     // For ALCT search, the key layer is the 3rd one, counting from 1.
0259     if (pli->getId().layer() == CSCConstants::KEY_ALCT_LAYER) {
0260       vector<PSimHit> thisLayerHits = pli->getSimHits();
0261       if (!thisLayerHits.empty()) {
0262         // There can be only one RecDigi (and therefore only one SimHit)
0263         // in a key layer.
0264         if (thisLayerHits.size() != 1) {
0265           edm::LogWarning("L1CSCTPEmulatorWrongValues")
0266               << "+++ Warning: " << thisLayerHits.size() << " SimHits in key layer " << CSCConstants::KEY_ALCT_LAYER
0267               << "! +++ \n";
0268           for (unsigned i = 0; i < thisLayerHits.size(); i++) {
0269             edm::LogWarning("L1CSCTPEmulatorWrongValues") << " SimHit # " << i << thisLayerHits[i] << "\n";
0270           }
0271         }
0272         matchedHit = thisLayerHits[0];
0273         layerId = pli->getId();
0274         hit_found = true;
0275         break;
0276       }
0277     }
0278   }
0279
0280   if (!hit_found) {
0281     for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0282       // if there is any occurrence of simHit size greater that zero, use this.
0283       if (!(pli->getRecDigis()).empty() && !(pli->getSimHits()).empty()) {
0284         // Always use the first SimHit for now.
0285         vector<PSimHit> thisLayerHits = pli->getSimHits();
0286         matchedHit = thisLayerHits[0];
0287         layerId = pli->getId();
0288         hit_found = true;
0289         break;
0290       }
0291     }
0292   }
0293
0294   // Set the eta if there were any hits found.
0295   if (hit_found) {
0296     const CSCLayer* csclayer = geom_->layer(layerId);
0297     const CSCLayerGeometry* layerGeom = csclayer->geometry();
0298     int nearestW = layerGeom->nearestWire(matchedHit.localPosition());
0299     nearestWG = layerGeom->wireGroup(nearestW);
0300     // Wire groups in ALCTs are counted starting from 0, whereas they
0301     // are counted from 1 in MC-related info.
0302     nearestWG -= 1;
0303     if (nearestWG < 0 || nearestWG >= CSCConstants::MAX_NUM_WIREGROUPS) {
0304       edm::LogWarning("L1CSCTPEmulatorWrongInput")
0305           << "+++ Warning: nearest wire group, " << nearestWG << ", is not in [0-" << CSCConstants::MAX_NUM_WIREGROUPS
0306           << ") interval +++\n";
0307     }
0308
0309     GlobalPoint thisPoint = csclayer->toGlobal(matchedHit.localPosition());
0310     closestPhi = thisPoint.phi();
0311     closestEta = thisPoint.eta();
0312     ostringstream strstrm;
0313     if (debug)
0314       strstrm << "Matched anode phi: " << closestPhi;
0315     if (closestPhi < 0.) {
0316       closestPhi += 2. * M_PI;
0317       if (debug)
0318         strstrm << " (" << closestPhi << ")";
0319     }
0320     if (debug) {
0321       strstrm << " eta: " << closestEta << " on a layer " << layerId.layer() << " (1-6);"
0322               << " nearest wire group: " << nearestWG;
0323       LogDebug("nearestWG") << strstrm.str();
0324     }
0325   }
0326
0327   return nearestWG;
0328 }
0329
0330 void CSCAnodeLCTAnalyzer::setGeometry(const CSCGeometry* geom) { geom_ = geom; }
0331
0332 double CSCAnodeLCTAnalyzer::getWGEta(const CSCDetId& layerId, const int wiregroup) {
0333   // Returns eta position of a given wiregroup.
0334   if (wiregroup < 0 || wiregroup >= CSCConstants::MAX_NUM_WIREGROUPS) {
0335     edm::LogWarning("L1CSCTPEmulatorWrongInput") << "+++ Warning: wire group, " << wiregroup << ", is not in [0-"
0336                                                  << CSCConstants::MAX_NUM_WIREGROUPS << ") interval +++\n";
0337   }
0338
0339   const CSCLayer* csclayer = geom_->layer(layerId);
0340   const CSCLayerGeometry* layerGeom = csclayer->geometry();
0341   LocalPoint digiLP = layerGeom->localCenterOfWireGroup(wiregroup + 1);
0342   //int wirePerWG = layerGeom->numberOfWiresPerGroup(wiregroup+1);
0343   //float middleW = layerGeom->middleWireOfGroup(wiregroup+1);
0344   //float ywire = layerGeom->yOfWire(middleW, 0.);
0345   //digiLP = LocalPoint(0., ywire);
0346   GlobalPoint digiGP = csclayer->toGlobal(digiLP);
0347   double eta = digiGP.eta();
0348
0349   return eta;
0350 }