File indexing completed on 2024-04-06 12:19:33
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include "L1Trigger/CSCTrackFinder/test/analysis/CSCTFEfficiency.h"
0022 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0023 #include "DataFormats/L1CSCTrackFinder/interface/L1CSCTrackCollection.h"
0024 #include "DataFormats/Scalers/interface/LumiScalers.h"
0025 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
0026 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
0027 #include <algorithm>
0028 #include "TStyle.h"
0029 using namespace std;
0030 int counterTFT;
0031 int counterRLimit;
0032 CSCTFEfficiency::CSCTFEfficiency(const edm::ParameterSet& iConfig) {
0033 usesResource(TFileService::kSharedResource);
0034
0035 inputTag = iConfig.getUntrackedParameter<edm::InputTag>("inputTag");
0036 dataType = iConfig.getUntrackedParameter<int>("type_of_data");
0037 minPtSim = iConfig.getUntrackedParameter<double>("MinPtSim");
0038 maxPtSim = iConfig.getUntrackedParameter<double>("MaxPtSim");
0039 minEtaSim = iConfig.getUntrackedParameter<double>("MinEtaSim");
0040 maxEtaSim = iConfig.getUntrackedParameter<double>("MaxEtaSim");
0041 minPtTF = iConfig.getUntrackedParameter<double>("MinPtTF");
0042 minQualityTF = iConfig.getUntrackedParameter<double>("MinQualityTF");
0043 ghostLoseParam = iConfig.getUntrackedParameter<std::string>("GhostLoseParam");
0044 inputData = iConfig.getUntrackedParameter<bool>("InputData");
0045 minMatchRParam = iConfig.getUntrackedParameter<double>("MinMatchR");
0046 statsFilename = iConfig.getUntrackedParameter<std::string>("StatsFilename");
0047 saveHistImages = iConfig.getUntrackedParameter<bool>("SaveHistImages");
0048 singleMuSample = iConfig.getUntrackedParameter<bool>("SingleMuSample");
0049 noRefTracks = iConfig.getUntrackedParameter<bool>("NoRefTracks");
0050 cutOnModes = iConfig.getUntrackedParameter<std::vector<unsigned> >("CutOnModes");
0051 nEvents = 0;
0052 configuration = &iConfig;
0053 }
0054 CSCTFEfficiency::~CSCTFEfficiency() {
0055
0056 }
0057
0058 void CSCTFEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0059 using namespace csctf_analysis;
0060
0061
0062
0063 std::vector<double>* Rlist = new std::vector<double>;
0064 std::vector<RefTrack>* referenceTrack = new std::vector<RefTrack>;
0065 std::vector<TFTrack>* trackFinderTrack = new std::vector<TFTrack>;
0066 std::vector<RefTrack>::iterator refTrack;
0067 std::vector<TFTrack>::iterator tfTrack;
0068
0069 if (!noRefTracks) {
0070
0071 edm::Handle<edm::SimTrackContainer> BaseSimTracks;
0072 iEvent.getByLabel("g4SimHits", BaseSimTracks);
0073 edm::SimTrackContainer::const_iterator BaseSimTrk;
0074 for (BaseSimTrk = BaseSimTracks->begin(); BaseSimTrk != BaseSimTracks->end(); BaseSimTrk++) {
0075 RefTrack refTrackTemp = RefTrack(*BaseSimTrk);
0076 double fabsEta = fabs(refTrackTemp.getEta());
0077 refTrackTemp.setHistList(refHistList);
0078 bool limits = refTrackTemp.getPt() < maxPtSim && refTrackTemp.getPt() > minPtSim && (fabsEta < maxEtaSim) &&
0079 (fabsEta > minEtaSim);
0080 bool isMuon = fabs(refTrackTemp.getType()) == 13;
0081 if (isMuon && limits) {
0082
0083
0084 referenceTrack->push_back(refTrackTemp);
0085 }
0086 }
0087 }
0088
0089
0090
0091
0092 if (dataType ==
0093 1)
0094 {
0095 edm::Handle<std::vector<L1MuRegionalCand> > trackFinderTracks;
0096 iEvent.getByLabel(inputTag, trackFinderTracks);
0097 std::vector<L1MuRegionalCand>::const_iterator BaseTFTrk;
0098 for (BaseTFTrk = trackFinderTracks->begin(); BaseTFTrk != trackFinderTracks->end(); BaseTFTrk++) {
0099 TFTrack tfTrackTemp = TFTrack(*BaseTFTrk);
0100 if (tfTrackTemp.getQuality() >= minQualityTF && tfTrackTemp.getPt() >= minPtTF) {
0101 tfTrackTemp.setHistList(tfHistList);
0102 trackFinderTrack->push_back(tfTrackTemp);
0103 }
0104 }
0105 } else if (dataType == 2) {
0106
0107 edm::Handle<L1MuGMTReadoutCollection> gmtreadoutCollectionH;
0108 iEvent.getByLabel(inputTag, gmtreadoutCollectionH);
0109 L1MuGMTReadoutCollection const* GMTrc = gmtreadoutCollectionH.product();
0110 vector<L1MuGMTReadoutRecord> gmt_records = GMTrc->getRecords();
0111 vector<L1MuGMTReadoutRecord>::const_iterator gmt_records_iter;
0112 for (gmt_records_iter = gmt_records.begin(); gmt_records_iter != gmt_records.end(); gmt_records_iter++) {
0113 vector<L1MuGMTExtendedCand>::const_iterator ExtCand_iter;
0114 vector<L1MuGMTExtendedCand> extendedCands = gmt_records_iter->getGMTCands();
0115 for (ExtCand_iter = extendedCands.begin(); ExtCand_iter != extendedCands.end(); ExtCand_iter++) {
0116 L1MuGMTExtendedCand gmtec = *ExtCand_iter;
0117 TFTrack tfTrackTemp = TFTrack(gmtec);
0118
0119 if (tfTrackTemp.getPt() >= minPtTF && fabs(tfTrackTemp.getEta()) <= maxEtaSim &&
0120 fabs(tfTrackTemp.getEta()) >= minEtaSim) {
0121 tfTrackTemp.setHistList(tfHistList);
0122 trackFinderTrack->push_back(tfTrackTemp);
0123 }
0124 }
0125 }
0126 } else
0127 {
0128 edm::Handle<L1CSCTrackCollection> trackFinderTracks;
0129 iEvent.getByLabel(inputTag, trackFinderTracks);
0130 L1CSCTrackCollection::const_iterator BaseTFTrk;
0131
0132 for (BaseTFTrk = trackFinderTracks->begin(); BaseTFTrk != trackFinderTracks->end(); BaseTFTrk++) {
0133 TFTrack tfTrackTemp = TFTrack(*BaseTFTrk, iSetup);
0134
0135
0136 if (cutOnModes.size() > 0) {
0137 std::vector<unsigned>::iterator found;
0138 found = std::find(cutOnModes.begin(), cutOnModes.end(), tfTrackTemp.getMode());
0139 if (found == cutOnModes.end())
0140 continue;
0141 }
0142
0143 if (tfTrackTemp.getQuality() >= minQualityTF && tfTrackTemp.getPt() >= minPtTF) {
0144 tfTrackTemp.setHistList(tfHistList);
0145 trackFinderTrack->push_back(tfTrackTemp);
0146 }
0147 }
0148 }
0149 multHistList->FillMultiplicityHist(trackFinderTrack);
0150 if (trackFinderTrack->size() >= 2) {
0151 rHistogram->fillR(trackFinderTrack->at(0), trackFinderTrack->at(1));
0152 }
0153
0154
0155
0156
0157
0158 unsigned int iRefTrack = 0;
0159 for (refTrack = referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++) {
0160 bool tftracksExist = trackFinderTrack->size() > 0;
0161 if (tftracksExist) {
0162 for (tfTrack = trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++) {
0163 Rlist->push_back(refTrack->distanceTo(&(*tfTrack)));
0164 }
0165 unsigned int iSmallR = minIndex(Rlist);
0166 double smallR = Rlist->at(iSmallR);
0167 bool bestMatch;
0168 bool oldMatch;
0169 if (trackFinderTrack->at(iSmallR).getMatched()) {
0170 bestMatch = smallR < trackFinderTrack->at(iSmallR).getR();
0171 oldMatch = true;
0172 } else if (smallR < minMatchRParam) {
0173 bestMatch = true;
0174 oldMatch = false;
0175 } else {
0176 bestMatch = false;
0177 oldMatch = false;
0178 }
0179
0180 if (bestMatch) {
0181 if (oldMatch) {
0182 int oldRefMatchi = trackFinderTrack->at(iSmallR).getMatchedIndex();
0183 referenceTrack->at(oldRefMatchi).unMatch();
0184 }
0185 int tfQ = trackFinderTrack->at(iSmallR).getQuality();
0186 double tfPt = trackFinderTrack->at(iSmallR).getPt();
0187 refTrack->matchedTo(iSmallR, smallR, tfQ, tfPt);
0188 refTrack->setQuality(tfQ);
0189 refTrack->setTFPt(tfPt);
0190 trackFinderTrack->at(iSmallR).matchedTo(iRefTrack, smallR);
0191 }
0192 }
0193 Rlist->clear();
0194 }
0195
0196
0197 int iHighest = 0;
0198 int i = 0;
0199 double highestPt = -100.0;
0200 for (refTrack = referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++) {
0201 refTrack->fillHist();
0202 if (refTrack->getPt() > highestPt) {
0203 iHighest = i;
0204 highestPt = refTrack->getPt();
0205 }
0206 if (refTrack->getMatched()) {
0207 refTrack->setMatch(trackFinderTrack->at(refTrack->getMatchedIndex()));
0208 refTrack->fillSimvTFHist(*refTrack, refTrack->getMatchedTrack());
0209 refTrack->fillMatchHist();
0210
0211 int TFTIndex = refTrack->getMatchedIndex();
0212 TFTrack tfTrk = trackFinderTrack->at(TFTIndex);
0213 resHistList->FillResolutionHist(*refTrack, tfTrk);
0214 }
0215 i++;
0216 }
0217 if (!referenceTrack->empty()) {
0218 referenceTrack->at(iHighest).fillRateHist();
0219 }
0220
0221 iHighest = 0;
0222 i = 0;
0223 highestPt = -100.0;
0224 for (tfTrack = trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++) {
0225 tfTrack->fillHist();
0226 if (tfTrack->getPt() > highestPt) {
0227 iHighest = i;
0228 highestPt = tfTrack->getPt();
0229 }
0230 if (tfTrack->getMatched()) {
0231 tfTrack->fillMatchHist();
0232 }
0233 }
0234 if (!trackFinderTrack->empty()) {
0235 trackFinderTrack->at(iHighest).fillRateHist();
0236 }
0237
0238
0239
0240
0241 if (singleMuSample) {
0242 if (trackFinderTrack->size() > 1) {
0243 std::vector<double>* Qlist = new std::vector<double>;
0244 for (tfTrack = trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++) {
0245 if (referenceTrack->size() > 0) {
0246 Rlist->push_back(referenceTrack->begin()->distanceTo(&(*tfTrack)));
0247 }
0248 Qlist->push_back(tfTrack->getQuality());
0249 }
0250 unsigned int iSmallR = minIndex(Rlist);
0251 unsigned int iSmallQ = minIndex(Qlist);
0252 if (referenceTrack->size() > 0) {
0253 if (ghostLoseParam == "Q") {
0254 trackFinderTrack->at(iSmallQ).fillGhostHist();
0255 } else if (ghostLoseParam == "R") {
0256 trackFinderTrack->at(iSmallR).fillGhostHist();
0257 } else {
0258 std::cout << "Warning: ghostLoseParam Not set to R or Q!" << std::endl;
0259 }
0260 referenceTrack->begin()->fillGhostHist();
0261 } else {
0262 std::cout << "Warning: Multiple TFTracks found, RefTrack NOT found!" << std::endl;
0263 trackFinderTrack->at(iSmallQ).fillGhostHist();
0264 }
0265 delete Qlist;
0266 }
0267 } else {
0268 unsigned int iTFTrack = 0;
0269 for (tfTrack = trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++) {
0270 bool reftracksExist = referenceTrack->size() > 0;
0271 if (reftracksExist) {
0272 for (refTrack = referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++) {
0273 Rlist->push_back(refTrack->distanceTo(&(*tfTrack)));
0274 }
0275 unsigned int iSmallR = minIndex(Rlist);
0276 double smallR = Rlist->at(iSmallR);
0277 RefTrack* matchedRefTrack = &referenceTrack->at(iSmallR);
0278 matchedRefTrack->ghostMatchedTo(*tfTrack, iTFTrack, smallR);
0279 Rlist->clear();
0280 iTFTrack++;
0281 }
0282 }
0283 for (refTrack = referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++) {
0284 if (refTrack->getGhost()) {
0285 refTrack->loseBestGhostCand(ghostLoseParam);
0286 std::vector<unsigned int>* ghosts;
0287 ghosts = refTrack->ghostMatchedIndecies();
0288 std::vector<unsigned int>::const_iterator iGhost;
0289 for (iGhost = ghosts->begin(); iGhost != ghosts->end(); iGhost++) {
0290 TFTrack* tempTFTrack = &trackFinderTrack->at(*iGhost);
0291 int tfQ = tempTFTrack->getQuality();
0292 refTrack->setQuality(tfQ);
0293 refTrack->fillGhostHist();
0294 tempTFTrack->fillGhostHist();
0295 }
0296 }
0297 }
0298 }
0299 delete Rlist;
0300 delete referenceTrack;
0301 delete trackFinderTrack;
0302 nEvents++;
0303 }
0304
0305 void CSCTFEfficiency::beginJob() {
0306 using namespace csctf_analysis;
0307
0308 gStyle->SetOptStat(0);
0309 tfHistList = new TrackHistogramList("TrackFinder", configuration);
0310 refHistList = new TrackHistogramList("Reference", configuration);
0311 effhistlist = new EffHistogramList("Efficiency", configuration);
0312 resHistList = new ResolutionHistogramList("Resolution", configuration);
0313 multHistList = new MultiplicityHistogramList();
0314 statFile = new StatisticsFile(statsFilename);
0315 rHistogram = new RHistogram("RHistograms");
0316
0317
0318 configuration = nullptr;
0319 }
0320
0321 void CSCTFEfficiency::endJob() {
0322 effhistlist->ComputeEff(refHistList);
0323 statFile->WriteStatistics(*tfHistList, *refHistList);
0324 if (saveHistImages) {
0325 effhistlist->Print();
0326 resHistList->Print();
0327 }
0328 delete tfHistList;
0329 delete refHistList;
0330 delete effhistlist;
0331 delete resHistList;
0332 delete statFile;
0333 }
0334 namespace csctf_analysis {
0335 unsigned int minIndex(const std::vector<int>* list) {
0336 unsigned int minI = 0;
0337 for (unsigned int i = 0; i < list->size(); i++) {
0338 if (list[i] < list[minI]) {
0339 minI = i;
0340 }
0341 }
0342 return minI;
0343 }
0344 unsigned int minIndex(const std::vector<double>* list) {
0345 unsigned int minI = 0;
0346 for (unsigned int i = 0; i < list->size(); i++) {
0347 if (list->at(i) < list->at(minI)) {
0348 minI = i;
0349 }
0350 }
0351 return minI;
0352 }
0353
0354 }
0355
0356 DEFINE_FWK_MODULE(CSCTFEfficiency);