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 {
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 }
0058 // ------------ method called to for each event  ------------
0059 void CSCTFEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
0060 {
0061   using namespace csctf_analysis;
0062   ///////////////////////////////
0063   ////// Event Set-Up ///////////
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         //Reference tracks
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             //std::cout<<"\nMinPtSim= "<<minPtSim<<"\n";
0089             //std::cout<<"Ref_Pt = "<<refTrackTemp.getPt()<<"\n";
0090             referenceTrack->push_back(refTrackTemp);
0091         }
0092         }
0093
0094     }
0095
0096
0097   //dataType key: 0 - L1CSCTrack
0098   //              1 - L1MuRegionalCand
0099   //              2 - L1MuGMTExtendedCand
0100   if (dataType==1)//If you want the original production of sim csctf tracks without modes and ghost modes information/plots
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     //GMT Tracks
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 //If you want the new production of sim csctf tracks with modes and ghost modes information/plots
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         //Skips track if mode is not found in cutOnModes;
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   //////// Track Matching //////////////
0175   //////////////////////////////////////
0176   // Loop over all Reference tracks for an Event
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   //Fill Histograms
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())); //this line links the two tracks
0237         refTrack->fillSimvTFHist(*refTrack,refTrack->getMatchedTrack());
0238         refTrack->fillMatchHist();
0239         // resolution
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   ////// Ghost Finding ///////////
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 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called once each job just after ending the event loop  ------------
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 //define this as a plug-in
0407 DEFINE_FWK_MODULE(CSCTFEfficiency);
0408