File indexing completed on 2021-02-14 14:21:47
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 {
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
0059 void CSCTFEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
0060 {
0061 using namespace csctf_analysis;
0062
0063
0064
0065 std::vector<double>* Rlist= new std::vector<double>;
0066 std::vector<RefTrack>* referenceTrack=new std::vector<RefTrack>;
0067 std::vector<TFTrack>* trackFinderTrack=new std::vector<TFTrack>;
0068 std::vector<RefTrack>::iterator refTrack;
0069 std::vector<TFTrack>::iterator tfTrack;
0070
0071
0072
0073 if(!noRefTracks)
0074 {
0075
0076 edm::Handle<edm::SimTrackContainer> BaseSimTracks;
0077 iEvent.getByLabel("g4SimHits",BaseSimTracks);
0078 edm::SimTrackContainer::const_iterator BaseSimTrk;
0079 for(BaseSimTrk=BaseSimTracks->begin(); BaseSimTrk != BaseSimTracks->end(); BaseSimTrk++)
0080 {
0081 RefTrack refTrackTemp= RefTrack(*BaseSimTrk);
0082 double fabsEta = fabs(refTrackTemp.getEta());
0083 refTrackTemp.setHistList(refHistList);
0084 bool limits = refTrackTemp.getPt() < maxPtSim && refTrackTemp.getPt() > minPtSim && (fabsEta < maxEtaSim) && (fabsEta > minEtaSim);
0085 bool isMuon = fabs(refTrackTemp.getType())==13;
0086 if(isMuon && limits)
0087 {
0088
0089
0090 referenceTrack->push_back(refTrackTemp);
0091 }
0092 }
0093
0094 }
0095
0096
0097
0098
0099
0100 if (dataType==1)
0101 {
0102 edm::Handle<std::vector<L1MuRegionalCand> > trackFinderTracks;
0103 iEvent.getByLabel(inputTag, trackFinderTracks);
0104 std::vector<L1MuRegionalCand>::const_iterator BaseTFTrk;
0105 for(BaseTFTrk=trackFinderTracks->begin();BaseTFTrk != trackFinderTracks->end(); BaseTFTrk++)
0106 {
0107 TFTrack tfTrackTemp= TFTrack(*BaseTFTrk);
0108 if(tfTrackTemp.getQuality() >= minQualityTF && tfTrackTemp.getPt() >= minPtTF )
0109 {
0110 tfTrackTemp.setHistList(tfHistList);
0111 trackFinderTrack->push_back(tfTrackTemp);
0112 }
0113 }
0114 }
0115 else if(dataType==2)
0116 {
0117
0118
0119 edm::Handle<L1MuGMTReadoutCollection> gmtreadoutCollectionH;
0120 iEvent.getByLabel(inputTag,gmtreadoutCollectionH);
0121 L1MuGMTReadoutCollection const* GMTrc = gmtreadoutCollectionH.product();
0122 vector<L1MuGMTReadoutRecord> gmt_records = GMTrc->getRecords();
0123 vector<L1MuGMTReadoutRecord>::const_iterator gmt_records_iter;
0124 for(gmt_records_iter=gmt_records.begin(); gmt_records_iter!=gmt_records.end(); gmt_records_iter++){
0125 vector<L1MuGMTExtendedCand>::const_iterator ExtCand_iter;
0126 vector<L1MuGMTExtendedCand> extendedCands = gmt_records_iter->getGMTCands();
0127 for(ExtCand_iter=extendedCands.begin();ExtCand_iter != extendedCands.end(); ExtCand_iter++){
0128
0129 L1MuGMTExtendedCand gmtec = *ExtCand_iter;
0130 TFTrack tfTrackTemp= TFTrack(gmtec);
0131
0132
0133 if(tfTrackTemp.getPt() >= minPtTF && fabs(tfTrackTemp.getEta())<=maxEtaSim && fabs(tfTrackTemp.getEta())>=minEtaSim){
0134 tfTrackTemp.setHistList(tfHistList);
0135 trackFinderTrack->push_back(tfTrackTemp);
0136 }
0137 }
0138 }
0139 }
0140 else
0141 {
0142 edm::Handle<L1CSCTrackCollection> trackFinderTracks;
0143 iEvent.getByLabel(inputTag,trackFinderTracks);
0144 L1CSCTrackCollection::const_iterator BaseTFTrk;
0145
0146 for(BaseTFTrk=trackFinderTracks->begin();BaseTFTrk != trackFinderTracks->end(); BaseTFTrk++)
0147 {
0148 TFTrack tfTrackTemp= TFTrack(*BaseTFTrk,iSetup);
0149
0150
0151 if(cutOnModes.size() > 0)
0152 {
0153 std::vector<unsigned>::iterator found;
0154 found = std::find(cutOnModes.begin(),cutOnModes.end(), tfTrackTemp.getMode());
0155 if(found == cutOnModes.end())
0156 continue;
0157 }
0158
0159
0160 if(tfTrackTemp.getQuality() >= minQualityTF && tfTrackTemp.getPt() >= minPtTF )
0161 {
0162 tfTrackTemp.setHistList(tfHistList);
0163 trackFinderTrack->push_back(tfTrackTemp);
0164 }
0165 }
0166 }
0167 multHistList->FillMultiplicityHist(trackFinderTrack);
0168 if(trackFinderTrack->size()>=2)
0169 {
0170 rHistogram->fillR(trackFinderTrack->at(0),trackFinderTrack->at(1));
0171 }
0172
0173
0174
0175
0176
0177 unsigned int iRefTrack=0;
0178 for(refTrack=referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++)
0179 {
0180 bool tftracksExist = trackFinderTrack->size() > 0;
0181 if (tftracksExist)
0182 {
0183 for(tfTrack=trackFinderTrack->begin();
0184 tfTrack != trackFinderTrack->end(); tfTrack++)
0185 {
0186 Rlist->push_back(refTrack->distanceTo(&(*tfTrack)));
0187 }
0188 unsigned int iSmallR = minIndex(Rlist);
0189 double smallR = Rlist->at(iSmallR);
0190 bool bestMatch;
0191 bool oldMatch;
0192 if (trackFinderTrack->at(iSmallR).getMatched())
0193 {
0194 bestMatch=smallR < trackFinderTrack->at(iSmallR).getR();
0195 oldMatch=true;
0196 }
0197 else if ( smallR < minMatchRParam )
0198 {
0199 bestMatch=true; oldMatch=false;
0200 }
0201 else
0202 {
0203 bestMatch = false; oldMatch = false;
0204 }
0205
0206 if (bestMatch)
0207 {
0208 if (oldMatch)
0209 {
0210 int oldRefMatchi = trackFinderTrack->at(iSmallR).getMatchedIndex();
0211 referenceTrack->at(oldRefMatchi).unMatch();
0212 }
0213 int tfQ=trackFinderTrack->at(iSmallR).getQuality();
0214 double tfPt=trackFinderTrack->at(iSmallR).getPt();
0215 refTrack->matchedTo(iSmallR, smallR,tfQ, tfPt);
0216 refTrack->setQuality(tfQ);
0217 refTrack->setTFPt(tfPt);
0218 trackFinderTrack->at(iSmallR).matchedTo(iRefTrack, smallR);
0219 }
0220 }
0221 Rlist->clear();
0222 }
0223
0224
0225
0226
0227 int iHighest=0; int i=0; double highestPt=-100.0;
0228 for(refTrack=referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++)
0229 {
0230
0231 refTrack->fillHist();
0232 if(refTrack->getPt()>highestPt){iHighest=i; highestPt=refTrack->getPt();}
0233 if (refTrack->getMatched())
0234 {
0235
0236 refTrack->setMatch(trackFinderTrack->at(refTrack->getMatchedIndex()));
0237 refTrack->fillSimvTFHist(*refTrack,refTrack->getMatchedTrack());
0238 refTrack->fillMatchHist();
0239
0240 int TFTIndex = refTrack->getMatchedIndex();
0241 TFTrack tfTrk = trackFinderTrack->at( TFTIndex );
0242 resHistList->FillResolutionHist( *refTrack, tfTrk );
0243 }
0244 i++;
0245
0246
0247 }
0248 if(!referenceTrack->empty()){referenceTrack->at(iHighest).fillRateHist();}
0249
0250 iHighest=0; i=0; highestPt=-100.0;
0251 for(tfTrack=trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++)
0252 {
0253
0254 tfTrack->fillHist();
0255 if(tfTrack->getPt()>highestPt){iHighest=i; highestPt=tfTrack->getPt();}
0256 if(tfTrack->getMatched())
0257 {
0258 tfTrack->fillMatchHist();
0259 }
0260
0261 }
0262 if(!trackFinderTrack->empty()){trackFinderTrack->at(iHighest).fillRateHist();}
0263
0264
0265
0266
0267 if(singleMuSample)
0268 {
0269 if(trackFinderTrack->size()>1)
0270 {
0271 std::vector<double>* Qlist= new std::vector<double>;
0272 for(tfTrack=trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++)
0273 {
0274 if(referenceTrack->size()>0)
0275 {
0276 Rlist->push_back(referenceTrack->begin()->distanceTo(&(*tfTrack)));
0277 }
0278 Qlist->push_back(tfTrack->getQuality());
0279 }
0280 unsigned int iSmallR = minIndex(Rlist);
0281 unsigned int iSmallQ = minIndex(Qlist);
0282 if(referenceTrack->size()>0)
0283 {
0284 if(ghostLoseParam=="Q")
0285 {
0286 trackFinderTrack->at(iSmallQ).fillGhostHist();
0287 }
0288 else if(ghostLoseParam=="R")
0289 {
0290 trackFinderTrack->at(iSmallR).fillGhostHist();
0291 }
0292 else
0293 {
0294 std::cout << "Warning: ghostLoseParam Not set to R or Q!" << std::endl;
0295 }
0296 referenceTrack->begin()->fillGhostHist();
0297 }
0298 else
0299 {
0300 std::cout << "Warning: Multiple TFTracks found, RefTrack NOT found!" << std::endl;
0301 trackFinderTrack->at(iSmallQ).fillGhostHist();
0302 }
0303 delete Qlist;
0304 }
0305 }
0306 else
0307 {
0308 unsigned int iTFTrack=0;
0309 for(tfTrack=trackFinderTrack->begin(); tfTrack != trackFinderTrack->end(); tfTrack++)
0310 {
0311 bool reftracksExist = referenceTrack->size() > 0;
0312 if(reftracksExist)
0313 {
0314 for(refTrack=referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++)
0315 {
0316 Rlist->push_back(refTrack->distanceTo(&(*tfTrack)));
0317 }
0318 unsigned int iSmallR = minIndex(Rlist);
0319 double smallR = Rlist->at(iSmallR);
0320 RefTrack* matchedRefTrack = &referenceTrack->at(iSmallR);
0321 matchedRefTrack->ghostMatchedTo(*tfTrack,iTFTrack,smallR);
0322 Rlist->clear();
0323 iTFTrack++;
0324 }
0325 }
0326 for(refTrack=referenceTrack->begin(); refTrack != referenceTrack->end(); refTrack++)
0327 {
0328 if(refTrack->getGhost())
0329 {
0330 refTrack->loseBestGhostCand(ghostLoseParam);
0331 std::vector<unsigned int>* ghosts;
0332 ghosts = refTrack->ghostMatchedIndecies();
0333 std::vector<unsigned int>::const_iterator iGhost;
0334 for(iGhost = ghosts->begin();iGhost != ghosts->end(); iGhost++)
0335 {
0336 TFTrack* tempTFTrack = &trackFinderTrack->at(*iGhost);
0337 int tfQ = tempTFTrack->getQuality();
0338 refTrack->setQuality(tfQ);
0339 refTrack->fillGhostHist();
0340 tempTFTrack->fillGhostHist();
0341 }
0342 }
0343 }
0344 }
0345 delete Rlist;
0346 delete referenceTrack;
0347 delete trackFinderTrack;
0348 nEvents++;
0349 }
0350
0351 void CSCTFEfficiency::beginJob()
0352 {
0353 using namespace csctf_analysis;
0354
0355 gStyle->SetOptStat(0);
0356 tfHistList = new TrackHistogramList("TrackFinder",configuration);
0357 refHistList = new TrackHistogramList("Reference",configuration);
0358 effhistlist = new EffHistogramList("Efficiency",configuration);
0359 resHistList = new ResolutionHistogramList("Resolution",configuration);
0360 multHistList = new MultiplicityHistogramList();
0361 statFile = new StatisticsFile( statsFilename );
0362 rHistogram = new RHistogram("RHistograms");
0363 }
0364
0365 void CSCTFEfficiency::endJob()
0366 {
0367 effhistlist->ComputeEff(refHistList);
0368 statFile->WriteStatistics( *tfHistList, *refHistList);
0369 if(saveHistImages)
0370 {
0371 effhistlist->Print();
0372 resHistList->Print();
0373 }
0374 delete tfHistList;
0375 delete refHistList;
0376 delete effhistlist;
0377 delete resHistList;
0378 delete statFile;
0379 }
0380 namespace csctf_analysis
0381 {
0382 unsigned int minIndex(const std::vector<int>* list){
0383 unsigned int minI=0;
0384 for(unsigned int i=0;i<list->size();i++)
0385 {
0386 if(list[i]<list[minI])
0387 {
0388 minI=i;
0389 }
0390 }
0391 return minI;
0392 }
0393 unsigned int minIndex(const std::vector<double>* list){
0394 unsigned int minI=0;
0395 for(unsigned int i=0;i<list->size();i++)
0396 {
0397 if(list->at(i)<list->at(minI))
0398 {
0399 minI=i;
0400 }
0401 }
0402 return minI;
0403 }
0404
0405 }
0406
0407 DEFINE_FWK_MODULE(CSCTFEfficiency);
0408