test/analysis/CSCTFEfficiency.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CSCTFEfficiency
0004 // Class:      CSCTFEfficiency
0005 //
0006 /**\class CSCTFEfficiency CSCTFEfficiency.cc jhugon/CSCTFEfficiency/src/CSCTFEfficiency.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Justin Hugon,Ivan's graduate student,jhugon@phys.ufl.edu
0015 //         Created:  Thu Jun 10 10:40:10 EDT 2010
0016 // $Id: CSCTFEfficiency.cc,v 1.2 2012/02/10 10:54:57 jhugon Exp $
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;     // Added By Daniel 07/02
0031 int counterRLimit;  // Added By Daniel 07/02
0032 CSCTFEfficiency::CSCTFEfficiency(const edm::ParameterSet& iConfig) {
0033   usesResource(TFileService::kSharedResource);
0034   //now do what ever initialization is needed
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 // ------------ method called to for each event  ------------
0058 void CSCTFEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0059   using namespace csctf_analysis;
0060   ///////////////////////////////
0061   ////// Event Set-Up ///////////
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     //Reference tracks
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         //std::cout<<"\nMinPtSim= "<<minPtSim<<"\n";
0083         //std::cout<<"Ref_Pt = "<<refTrackTemp.getPt()<<"\n";
0084         referenceTrack->push_back(refTrackTemp);
0085       }
0086     }
0087   }
0088
0089   //dataType key: 0 - L1CSCTrack
0090   //              1 - L1MuRegionalCand
0091   //              2 - L1MuGMTExtendedCand
0092   if (dataType ==
0093       1)  //If you want the original production of sim csctf tracks without modes and ghost modes information/plots
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     //GMT Tracks
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  //If you want the new production of sim csctf tracks with modes and ghost modes information/plots
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       //Skips track if mode is not found in cutOnModes;
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   //////// Track Matching //////////////
0156   //////////////////////////////////////
0157   // Loop over all Reference tracks for an Event
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   //Fill Histograms
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()));  //this line links the two tracks
0208       refTrack->fillSimvTFHist(*refTrack, refTrack->getMatchedTrack());
0209       refTrack->fillMatchHist();
0210       // resolution
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   ////// Ghost Finding ///////////
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 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called once each job just after ending the event loop  ------------
0318 void CSCTFEfficiency::endJob() {
0319   effhistlist->ComputeEff(refHistList);
0320   statFile->WriteStatistics(*tfHistList, *refHistList);
0321   if (saveHistImages) {
0322     effhistlist->Print();
0323     resHistList->Print();
0324   }
0325   delete tfHistList;
0326   delete refHistList;
0327   delete effhistlist;
0328   delete resHistList;
0329   delete statFile;
0330 }
0331 namespace csctf_analysis {
0332   unsigned int minIndex(const std::vector<int>* list) {
0333     unsigned int minI = 0;
0334     for (unsigned int i = 0; i < list->size(); i++) {
0335       if (list[i] < list[minI]) {
0336         minI = i;
0337       }
0338     }
0339     return minI;
0340   }
0341   unsigned int minIndex(const std::vector<double>* list) {
0342     unsigned int minI = 0;
0343     for (unsigned int i = 0; i < list->size(); i++) {
0344       if (list->at(i) < list->at(minI)) {
0345         minI = i;
0346       }
0347     }
0348     return minI;
0349   }
0350
0351 }  // namespace csctf_analysis
0352 //define this as a plug-in
0353 DEFINE_FWK_MODULE(CSCTFEfficiency);