CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch2/src/SUSYBSMAnalysis/HSCP/src/HSCPValidator.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    HSCP
00004 // Class:      HSCPValidator
00005 // 
00013 //
00014 // Original Author:  Seth Cooper,27 1-024,+41227672342,
00015 //         Created:  Wed Apr 14 14:27:52 CEST 2010
00016 // $Id: HSCPValidator.cc,v 1.9 2011/10/11 21:14:33 jiechen Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 #include <vector>
00024 #include <string>
00025 #include <map>
00026 
00027 // user include files
00028 #include "FWCore/ServiceRegistry/interface/Service.h"
00029 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00030 #include "DataFormats/Common/interface/Handle.h"
00031 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00032 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00033 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
00034 #include "DataFormats/RPCDigi/interface/RPCDigi.h"
00035 #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
00036 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
00037 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00038 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00039 #include <DataFormats/GeometrySurface/interface/LocalError.h>
00040 #include <DataFormats/GeometryVector/interface/LocalPoint.h>
00041 #include "DataFormats/GeometrySurface/interface/Surface.h"
00042 #include "DataFormats/DetId/interface/DetId.h"
00043 #include "DataFormats/Candidate/interface/Candidate.h"
00044 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00045 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00046 #include <DataFormats/GeometrySurface/interface/LocalError.h>
00047 #include <DataFormats/GeometryVector/interface/LocalPoint.h>
00048 #include "DataFormats/GeometrySurface/interface/Surface.h"
00049 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00050 #include "FastSimulation/Tracking/test/FastTrackAnalyzer.h"
00051 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
00052 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00053 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
00054 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00055 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
00056 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
00057 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00058 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00059 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00060 #include "MagneticField/Engine/interface/MagneticField.h"
00061 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00062 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00063 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
00064 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00065 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
00066 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00067 #include "SimDataFormats/Track/interface/SimTrack.h"
00068 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00069 #include "SimDataFormats/Vertex/interface/SimVertex.h"
00070 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00071 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00072 #include "SUSYBSMAnalysis/HSCP/interface/HSCPValidator.h"
00073 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00074 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00075 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00076 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00077 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00078 
00079 #include "DataFormats/TrackReco/interface/DeDxData.h"
00080 #include "DataFormats/TrackReco/interface/TrackToTrackMap.h"
00081 
00082 #include "TH1.h"
00083 #include "TGraph.h"
00084 #include "TCanvas.h"
00085 
00086 //
00087 // constants, enums and typedefs
00088 //
00089 
00090 //
00091 // static data member definitions
00092 //
00093 edm::Service<TFileService> fileService;
00094 
00095 //
00096 // constructors and destructor
00097 //
00098 HSCPValidator::HSCPValidator(const edm::ParameterSet& iConfig) :
00099   doGenPlots_ (iConfig.getParameter<bool>("MakeGenPlots")),
00100   doHLTPlots_ (iConfig.getParameter<bool>("MakeHLTPlots")),
00101   doSimTrackPlots_ (iConfig.getParameter<bool>("MakeSimTrackPlots")),
00102   doSimDigiPlots_ (iConfig.getParameter<bool>("MakeSimDigiPlots")),
00103   doRecoPlots_ (iConfig.getParameter<bool>("MakeRecoPlots")),
00104   label_ (iConfig.getParameter<edm::InputTag>("generatorLabel")),
00105   particleIds_ (iConfig.getParameter< std::vector<int> >("particleIds")),
00106   particleStatus_ (iConfig.getUntrackedParameter<int>("particleStatus",1)),
00107   ebSimHitTag_ (iConfig.getParameter<edm::InputTag>("EBSimHitCollection")),
00108   eeSimHitTag_ (iConfig.getParameter<edm::InputTag>("EESimHitCollection")),
00109   simTrackTag_ (iConfig.getParameter<edm::InputTag>("SimTrackCollection")),
00110   EBDigiCollection_ (iConfig.getParameter<edm::InputTag>("EBDigiCollection")),
00111   EEDigiCollection_ (iConfig.getParameter<edm::InputTag>("EEDigiCollection")),
00112   RPCRecHitTag_ (iConfig.getParameter<edm::InputTag>("RPCRecHitTag"))
00113 {
00114   //now do what ever initialization is needed
00115   // GEN
00116   particleEtaHist_ = fileService->make<TH1F>("particleEta","Eta of gen particle",100,-5,5);
00117   particlePhiHist_ = fileService->make<TH1F>("particlePhi","Phi of gen particle",180,-3.15,3.15);
00118   particlePHist_ = fileService->make<TH1F>("particleP","Momentum of gen particle",500,0,2000);
00119   particlePtHist_ = fileService->make<TH1F>("particlePt","P_{T} of gen particle",500,0,2000);
00120   particleMassHist_ = fileService->make<TH1F>("particleMass","Mass of gen particle",1000,0,2000);
00121   particleStatusHist_ = fileService->make<TH1F>("particleStatus","Status of gen particle",10,0,10);
00122   particleBetaHist_ = fileService->make<TH1F>("particleBeta","Beta of gen particle",100,0,1);
00123   particleBetaInverseHist_ = fileService->make<TH1F>("particleBetaInverse","1/#beta of gen particle",100,0,5);
00124   
00125   h_genhscp_met = fileService->make<TH1F>( "hscp_met"  , "missing E_{T} hscp" , 100, 0., 1500. );
00126   h_genhscp_met_nohscp = fileService->make<TH1F>( "hscp_met_nohscp"  , "missing E_{T} w/o hscp" , 100, 0., 1500. );
00127   h_genhscp_scaloret =  fileService->make<TH1F>( "hscp_scaloret"  , "scalor E_{T} sum" , 100, 0., 1500. );
00128   h_genhscp_scaloret_nohscp =  fileService->make<TH1F>( "hscp_scaloret_nohscp"  , "scalor E_{T} sum w/o hscp" , 100, 0., 1500. );
00129 
00130 
00131 
00132 
00133   //SIM track Info
00134   simTrackParticleEtaHist_ = fileService->make<TH1F>("simTrackParticleEta","Eta of simTrackParticle",100,-5,5);
00135   simTrackParticlePhiHist_ = fileService->make<TH1F>("simTrackParticlePhi","Phi of simTrackParticle",180,-3.15,3.15);
00136   simTrackParticlePHist_ = fileService->make<TH1F>("simTrackParticleP","Momentum of simTrackParticle",500,0,2000);
00137   simTrackParticlePtHist_ = fileService->make<TH1F>("simTrackParticlePt","P_{T} of simTrackParticle",500,0,2000);
00138   simTrackParticleBetaHist_ = fileService->make<TH1F>("simTrackParticleBeta","Beta of simTrackParticle",100,0,1);
00139 //reco track Info
00140 
00141   RecoHSCPPtVsGenPt= fileService->make<TH2F>("Recovsgenpt","RecovsGen",100,0,1000,100,0,1000);
00142   dedxVsp = fileService->make<TH2F>("dedxvsp","dedxvsp",100,0,1000,100,0,10);
00143 //HLT Info
00144   hltmet = fileService->make<TH1F>("HLT_MET","MET",3,-1,2);
00145   hltjet = fileService->make<TH1F>("HLT_JET","JET",3,-1,2);
00146   hltmu = fileService->make<TH1F>("HLT_Mu","Mu",3,-1,2);
00147 
00148 
00149   // SIM-DIGI: ECAL
00150   simHitsEcalEnergyHistEB_ = fileService->make<TH1F>("ecalEnergyOfSimHitsEB","HSCP SimTrack-matching SimHit energy EB [GeV]",125,-1,4);
00151   simHitsEcalEnergyHistEE_ = fileService->make<TH1F>("ecalEnergyOfSimHitsEE","HSCP SimTrack-matching SimHit energy EE [GeV]",125,-1,4);
00152   simHitsEcalTimeHistEB_ = fileService->make<TH1F>("ecalTimingOfSimHitsEB","HSCP SimTrack-matching SimHit time EB [ns]",115,-15,100);
00153   simHitsEcalTimeHistEE_ = fileService->make<TH1F>("ecalTimingOfSimHitsEE","HSCP SimTrack-matching SimHit time EE [ns]",115,-15,100);
00154   simHitsEcalNumHistEB_ = fileService->make<TH1F>("ecalNumberOfSimHitsEB","Number of HSCP SimTrack-matching EB sim hits in event",100,0,200);
00155   simHitsEcalNumHistEE_ = fileService->make<TH1F>("ecalNumberOfSimHitsEE","Number of HSCP SimTrack-matching EE sim hits in event",100,0,200);
00156   simHitsEcalEnergyVsTimeHistEB_ = fileService->make<TH2F>("ecalEnergyVsTimeOfSimHitsEB","Energy vs. time of HSCP SimTrack-matching EB sim hits in event",115,-15,100,125,-1,4);
00157   simHitsEcalEnergyVsTimeHistEE_ = fileService->make<TH2F>("ecalEnergyVsTimeOfSimHitsEE","Energy vs. time of HSCP SimTrack-matching EE sim hits in event",115,-15,100,125,-1,4);
00158   simHitsEcalDigiMatchEnergyHistEB_ = fileService->make<TH1F>("ecalEnergyOfDigiMatSimHitsEB","HSCP digi-matching SimHit energy EB [GeV]",125,-1,4);
00159   simHitsEcalDigiMatchEnergyHistEE_ = fileService->make<TH1F>("ecalEnergyOfDigiMatSimHitsEE","HSCP digi-matching SimHit energy EE [GeV]",125,-1,4);
00160   simHitsEcalDigiMatchTimeHistEB_ = fileService->make<TH1F>("ecalTimingOfDigiMatSimHitsEB","HSCP digi-matching SimHit time EB [ns]",115,-15,100);
00161   simHitsEcalDigiMatchTimeHistEE_ = fileService->make<TH1F>("ecalTimingOfDigiMatSimHitsEE","HSCP digi-matching SimHit time EE [ns]",115,-15,100);
00162   simHitsEcalDigiMatchEnergyVsTimeHistEB_ = fileService->make<TH2F>("ecalEnergyVsTimeOfDigiMatSimHitsEB","HSCP digi-matching EB SimHit energy vs. time",115,-15,100,125,-1,4);
00163   simHitsEcalDigiMatchEnergyVsTimeHistEE_ = fileService->make<TH2F>("ecalEnergyVsTimeOfDigiMatSimHitsEE","HSCP digi-matching EE SimHit energy vs. time",115,-15,100,125,-1,4);
00164   simHitsEcalDigiMatchIEtaHist_ = fileService->make<TH1F>("ecalIEtaOfDigiMatchSimHits","iEta of digi-matching Ecal simHits (EB)",171,-85,86);
00165   simHitsEcalDigiMatchIPhiHist_ = fileService->make<TH1F>("ecalIPhiOfDigiMatchSimHits","iPhi of digi-matching Ecal simHits (EB)",360,1,361);
00166   digisEcalNumHistEB_ = fileService->make<TH1F>("ecalDigisNumberEB","Number of EB digis matching simhits in event",200,0,1000);
00167   digisEcalNumHistEE_ = fileService->make<TH1F>("ecalDigisNumberEE","Number of EE digis matching simhits in event",200,0,1000);
00168   digiOccupancyMapEB_ = fileService->make<TH2F>("ecalDigiOccupancyMapEB","Occupancy of simhit-matching digis EB;i#phi;i#eta",360,1,361,171,-85,86);
00169   digiOccupancyMapEEP_ = fileService->make<TH2F>("ecalDigiOccupancyMapEEM","Occupancy of simhit-matching digis EEM;ix;iy",100,1,100,100,1,100);
00170   digiOccupancyMapEEM_ = fileService->make<TH2F>("ecalDigiOccupancyMapEEP","Occupancy of simhit-matching digis EEP;ix;iy",100,1,100,100,1,100);
00171 
00172   // SIM-DIGI: RPC
00173   residualsRPCRecHitSimDigis_ = fileService->make<TH1F>("residualsRPCRecHitSimDigis","HSCP SimHit - Clossest RPC RecHit",100,-5,5);
00174   efficiencyRPCRecHitSimDigis_ = fileService->make<TH1F>("efficiencyRPCRecHitSimDigis","HSCP SimHits RecHits Efficiency",2,-0.5,1.5);
00175   cluSizeDistribution_ = fileService->make<TH1F>("RPCCluSizeDistro","RPC HSCP CluSize Distribution",11,-0.5,10.5);
00176   rpcTimeOfFlightBarrel_[0] = fileService->make<TH1F>("RPCToFLayer1","RPC HSCP Time Of Flight Layer 1",50,5,100);
00177   rpcTimeOfFlightBarrel_[1] = fileService->make<TH1F>("RPCToFLayer2","RPC HSCP Time Of Flight Layer 2",50,5,100);
00178   rpcTimeOfFlightBarrel_[2] = fileService->make<TH1F>("RPCToFLayer3","RPC HSCP Time Of Flight Layer 3",50,5,100);
00179   rpcTimeOfFlightBarrel_[3] = fileService->make<TH1F>("RPCToFLayer4","RPC HSCP Time Of Flight Layer 4",50,5,100);
00180   rpcTimeOfFlightBarrel_[4] = fileService->make<TH1F>("RPCToFLayer5","RPC HSCP Time Of Flight Layer 5",50,5,100);
00181   rpcTimeOfFlightBarrel_[5] = fileService->make<TH1F>("RPCToFLayer6","RPC HSCP Time Of Flight Layer 6",50,5,100);
00182   rpcBXBarrel_[0] = fileService->make<TH1F>("RPCBXLayer1","RPC HSCP BX Layer 1",5,-0.5,4.5);
00183   rpcBXBarrel_[1] = fileService->make<TH1F>("RPCBXLayer2","RPC HSCP BX Layer 2",5,-0.5,4.5);
00184   rpcBXBarrel_[2] = fileService->make<TH1F>("RPCBXLayer3","RPC HSCP BX Layer 3",5,-0.5,4.5);
00185   rpcBXBarrel_[3] = fileService->make<TH1F>("RPCBXLayer4","RPC HSCP BX Layer 4",5,-0.5,4.5);
00186   rpcBXBarrel_[4] = fileService->make<TH1F>("RPCBXLayer5","RPC HSCP BX Layer 5",5,-0.5,4.5);
00187   rpcBXBarrel_[5] = fileService->make<TH1F>("RPCBXLayer6","RPC HSCP BX Layer 6",5,-0.5,4.5);
00188   rpcTimeOfFlightEndCap_[0]= fileService->make<TH1F>("RPCToFDisk1","RPC HSCP Time Of Flight Disk 1",50,5,100);
00189   rpcTimeOfFlightEndCap_[1]= fileService->make<TH1F>("RPCToFDisk2","RPC HSCP Time Of Flight Disk 2",50,5,100);
00190   rpcTimeOfFlightEndCap_[2]= fileService->make<TH1F>("RPCToFDisk3","RPC HSCP Time Of Flight Disk 3",50,5,100);
00191   rpcBXEndCap_[0] = fileService->make<TH1F>("RPCBXDisk1","RPC HSCP BX Disk 1",5,-0.5,4.5);
00192   rpcBXEndCap_[1] = fileService->make<TH1F>("RPCBXDisk2","RPC HSCP BX Disk 2",5,-0.5,4.5);
00193   rpcBXEndCap_[2] = fileService->make<TH1F>("RPCBXDisk3","RPC HSCP BX Disk 3",5,-0.5,4.5);
00194 }
00195 
00196 
00197 HSCPValidator::~HSCPValidator()
00198 {
00199  
00200    // do anything here that needs to be done at desctruction time
00201    // (e.g. close files, deallocate resources etc.)
00202 //   particleEtaHist_ = fileService->make<TH1F>("particleEta","Eta of gen particle",100,-5,5);
00203 //   particlePhiHist_ = fileService->make<TH1F>("particlePhi","Phi of gen particle",180,-3.15,3.15);
00204 //   particlePHist_ = fileService->make<TH1F>("particleP","Momentum of gen particle",500,0,2000);
00205 //   particlePtHist_ = fileService->make<TH1F>("particlePt","P_{T} of gen particle",500,0,2000);
00206 //   particleMassHist_ = fileService->make<TH1F>("particleMass","Mass of gen particle",1000,0,2000);
00207 //   particleStatusHist_ = fileService->make<TH1F>("particleStatus","Status of gen particle",10,0,10);
00208 //   particleBetaHist_ = fileService->make<TH1F>("particleBeta","Beta of gen particle",100,0,1);
00209 //   particleBetaInverseHist_ = fileService->make<TH1F>("particleBetaInverse","1/#beta of gen particle",100,0,5);
00210 
00211 }
00212 
00213 
00214 //
00215 // member functions
00216 //
00217 
00218 // ------------ method called to for each event  ------------
00219 void
00220 HSCPValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00221 {
00222   using namespace edm;
00223   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00224 
00225 
00226   if(doGenPlots_)
00227     makeGenPlots(iEvent);
00228   if(doSimTrackPlots_)
00229     makeSimTrackPlots(iEvent);
00230   if(doSimDigiPlots_){
00231     makeSimDigiPlotsECAL(iEvent);
00232     makeSimDigiPlotsRPC(iEvent);
00233   }
00234   if(doHLTPlots_){
00235     makeHLTPlots(iEvent);
00236   }
00237   if(doRecoPlots_){
00238      makeRecoPlots(iEvent);
00239   }
00240 }
00241 
00242 
00243 // ------------ method called once each job just before starting event loop  ------------
00244 void 
00245 HSCPValidator::beginJob()
00246 {
00247 }
00248 
00249 // ------------ method called once each job just after ending the event loop  ------------
00250 void 
00251 HSCPValidator::endJob()
00252 {
00253   std::string frequencies = "";
00254   for(std::map<int,int>::const_iterator itr = particleIdsFoundMap_.begin();
00255       itr != particleIdsFoundMap_.end(); ++itr)
00256   {
00257       frequencies+="PDG ID: ";
00258       frequencies+=intToString(itr->first);
00259       frequencies+=" Frequency: ";
00260       frequencies+=intToString(itr->second);
00261       frequencies+="\n";
00262   }
00263   std::cout << "Found PDGIds: " << "\n\n" << frequencies << std::endl;
00264 
00265 }
00266 
00267 // ------------- Make gen plots ---------------------------------------------------------
00268 void HSCPValidator::makeGenPlots(const edm::Event& iEvent)
00269 {
00270   using namespace edm;
00271 
00272   double missingpx=0;
00273   double missingpy=0;
00274   double missingpx_nohscp=0;
00275   double missingpy_nohscp=0;
00276   double scalorEt=0;
00277   double scalorEt_nohscp=0;
00278 
00279 
00280   Handle<HepMCProduct> evt;
00281   iEvent.getByLabel(label_, evt);
00282 
00283   HepMC::GenEvent * myGenEvent = new  HepMC::GenEvent(*(evt->GetEvent()));
00284   for(HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin();
00285       p != myGenEvent->particles_end(); ++p )
00286   {
00287 
00288     if((*p)->status() != particleStatus_)
00289       continue;
00290     //calculate MET(neutrino as MET)
00291     if(abs((*p)->pdg_id())!=12 && abs((*p)->pdg_id())!=14 && abs((*p)->pdg_id())!=16){ //for non-neutrino particles. 
00292        missingpx-=(*p)->momentum().px();
00293        missingpy-=(*p)->momentum().py();
00294        scalorEt+=(*p)->momentum().perp();
00295     }
00296 
00297     // Check if the particleId is in our R-hadron list
00298     std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),(*p)->pdg_id());
00299     if(partIdItr==particleIds_.end()){
00300        
00301        //calculate MET(neutrino+ HSCP as MET)
00302        if(abs((*p)->pdg_id())!=12 && abs((*p)->pdg_id())!=14 && abs((*p)->pdg_id())!=16){ //for non-neutrino particles. 
00303           missingpx_nohscp-=(*p)->momentum().px();
00304           missingpy_nohscp-=(*p)->momentum().py();
00305           scalorEt_nohscp+=(*p)->momentum().perp();
00306         }
00307     }
00308     else{
00309 
00310        particleStatusHist_->Fill((*p)->status());
00311     
00312        std::pair<std::map<int,int>::iterator,bool> pair = particleIdsFoundMap_.insert(std::make_pair<int,int>((*p)->pdg_id(),1));
00313        if(!pair.second)
00314        {
00315           ++(pair.first->second);
00316        }
00317        
00318        double mag = sqrt(pow((*p)->momentum().px(),2) + pow((*p)->momentum().py(),2) + pow((*p)->momentum().pz(),2) );
00319        particleEtaHist_->Fill((*p)->momentum().eta());
00320        particlePhiHist_->Fill((*p)->momentum().phi());
00321        particlePHist_->Fill(mag);
00322        particlePtHist_->Fill((*p)->momentum().perp());
00323        particleMassHist_->Fill((*p)->generated_mass());
00324        float particleP = mag;
00325        float particleM = (*p)->generated_mass();
00326        particleBetaHist_->Fill(particleP/sqrt(particleP*particleP+particleM*particleM));
00327        particleBetaInverseHist_->Fill(sqrt(particleP*particleP+particleM*particleM)/particleP);
00328     }
00329        
00330   }
00331 
00332   h_genhscp_met->Fill(sqrt(missingpx*missingpx+missingpy*missingpy));
00333   h_genhscp_met_nohscp->Fill(sqrt(missingpx_nohscp*missingpx_nohscp+missingpy_nohscp*missingpy_nohscp));
00334   h_genhscp_scaloret->Fill(scalorEt);
00335   h_genhscp_scaloret_nohscp->Fill(scalorEt_nohscp);
00336 
00337 
00338   delete myGenEvent; 
00339 
00340 
00341 
00342 }
00343 
00344 // ------------- Make SimTrack plots ---------------------------------------------------------
00345 void HSCPValidator::makeSimTrackPlots(const edm::Event& iEvent)
00346 {  using namespace edm;
00347   //get sim track infos
00348   Handle<edm::SimTrackContainer> simTracksHandle;
00349   iEvent.getByLabel("g4SimHits",simTracksHandle);
00350   const SimTrackContainer simTracks = *(simTracksHandle.product());
00351 
00352   SimTrackContainer::const_iterator simTrack;
00353 
00354   for (simTrack = simTracks.begin(); simTrack != simTracks.end(); ++simTrack){
00355      // Check if the particleId is in our list
00356      std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),simTrack->type());
00357      if(partIdItr==particleIds_.end()) continue;
00358 
00359      simTrackParticleEtaHist_->Fill((*simTrack).momentum().eta());
00360      simTrackParticlePhiHist_->Fill((*simTrack).momentum().phi());
00361      simTrackParticlePHist_->Fill((*simTrack).momentum().P());
00362      
00363      simTrackParticlePtHist_->Fill((*simTrack).momentum().pt());
00364      
00365      simTrackParticleBetaHist_->Fill((*simTrack).momentum().P()/(*simTrack).momentum().e());  
00366      
00367      // std::cout<<"Particle:"<<simTrack->type()<<" Charge:"<<simTrack->charge()<<std::endl;
00368 
00369   }
00370 }
00371 // ------------- Make HLT plots ---------------------------------------------------------
00372 void HSCPValidator::makeHLTPlots(const edm::Event& iEvent)
00373 {
00374   using namespace edm;
00375   //get HLT infos
00376      
00377 
00378       edm::TriggerResultsByName tr = iEvent.triggerResultsByName("HLT");
00379 
00380           if(!tr.isValid()){
00381         std::cout<<"Tirgger Results not available"<<std::endl;
00382       }
00383  
00384    edm::Handle< trigger::TriggerEvent > trEvHandle;
00385    iEvent.getByLabel("hltTriggerSummaryAOD", trEvHandle);
00386    trigger::TriggerEvent trEv = *trEvHandle;
00387 
00388 
00389    unsigned int TrIndex_Unknown     = tr.size();
00390 
00391 
00392    // HLT TRIGGER BASED ON 1 MUON!
00393    if(TrIndex_Unknown != tr.triggerIndex("HLT_Mu40_v1")){
00394       if(tr.accept(tr.triggerIndex("HLT_Mu40_v1"))) hltmu->Fill(1);
00395       else {hltmu->Fill(0);}
00396    }
00397    else{
00398       if(TrIndex_Unknown != tr.triggerIndex("HLT_Mu30_v1")){
00399          if(IncreasedTreshold(trEv, InputTag("hltSingleMu30L3Filtered30","","HLT"), 40,2.1, 1, false)) hltmu->Fill(1);
00400          else hltmu->Fill(0);
00401       }else{
00402          printf("BUG with HLT_Mu\n");
00403          std::cout<<"trigger names are : ";
00404          for(unsigned int i=0;i<tr.size();i++){
00405             std::cout<<" "<<tr.triggerName(i);
00406          }
00407          std::cout<<std::endl;       
00408       }
00409    }
00410 
00411 
00412  // HLT TRIGGER BASED ON MET!
00413    if(TrIndex_Unknown != tr.triggerIndex("HLT_PFMHT150_v3")){
00414       if(tr.accept(tr.triggerIndex("HLT_PFMHT150_v3")))hltmet->Fill(1);
00415       else hltmet->Fill(0);
00416    }else{
00417       if(TrIndex_Unknown != tr.triggerIndex("HLT_PFMHT150_v2")){
00418           if(tr.accept(tr.triggerIndex("HLT_PFMHT150_v2"))) hltmet->Fill(1);
00419           else hltmet->Fill(0);
00420       }
00421       else{
00422          printf("BUG with HLT_MET\n");
00423          
00424       }
00425    }
00426    
00427 
00428   // HLT TRIGGER BASED ON 1 JET!
00429    if(TrIndex_Unknown != tr.triggerIndex("HLT_Jet370_v1")){
00430        if(tr.accept(tr.triggerIndex("HLT_Jet370_v1")))hltjet->Fill(1);
00431        else   hltjet->Fill(0);
00432    }else{ 
00433       if(TrIndex_Unknown != tr.triggerIndex("HLT_Jet100U")){
00434          if(IncreasedTreshold(trEv, InputTag("hlt1jet100U","","HLT"), 140, 5.,1, false))hltjet->Fill(1);
00435          else   hltjet->Fill(0);
00436       }else{
00437          if(TrIndex_Unknown != tr.triggerIndex("HLT_Jet70U")){   
00438             if(IncreasedTreshold(trEv, InputTag("hlt1jet70U","","HLT"), 140, 5.,1, false))hltjet->Fill(1);
00439             else   hltjet->Fill(0);
00440          }else{
00441             if(TrIndex_Unknown != tr.triggerIndex("HLT_Jet50U")){
00442                if(IncreasedTreshold(trEv, InputTag("hlt1jet50U","","HLT"), 140,2.5, 1, false))hltjet->Fill(1);
00443                else   hltjet->Fill(0); 
00444             }else{
00445                printf("BUG with HLT_Jet\n");
00446                
00447             }
00448          }
00449       }
00450    }
00451    
00452 
00453 
00454 
00455   
00456 }
00457 
00458 // ------------- Make simDigi plots ECAL ------------------------------------------------
00459 void HSCPValidator::makeSimDigiPlotsECAL(const edm::Event& iEvent)
00460 {
00461   using namespace edm;
00462   // EB SimHits
00463   Handle<PCaloHitContainer> ebSimHits;
00464   iEvent.getByLabel(ebSimHitTag_, ebSimHits);
00465   if(!ebSimHits.isValid())
00466   {
00467     std::cout << "Cannot get EBSimHits from event!" << std::endl;
00468     return;
00469   }
00470   // EE SimHits
00471   Handle<PCaloHitContainer> eeSimHits;
00472   iEvent.getByLabel(eeSimHitTag_, eeSimHits);
00473   if(!eeSimHits.isValid())
00474   {
00475     std::cout << "Cannot get EESimHits from event!" << std::endl;
00476     return;
00477   }
00478   // SimTracks
00479   Handle<SimTrackContainer> simTracks;
00480   iEvent.getByLabel(simTrackTag_,simTracks);
00481   if(!simTracks.isValid())
00482   {
00483     std::cout << "Cannot get SimTracks from event!" << std::endl;
00484     return;
00485   }
00486   // EB Digis
00487   Handle<EBDigiCollection> ebDigis;
00488   iEvent.getByLabel(EBDigiCollection_,ebDigis);
00489   if(!ebDigis.isValid())
00490   {
00491     std::cout << "Cannot get EBDigis from event!" << std::endl;
00492     return;
00493   }
00494   // EE Digis
00495   Handle<EEDigiCollection> eeDigis;
00496   iEvent.getByLabel(EEDigiCollection_,eeDigis);
00497   if(!eeDigis.isValid())
00498   {
00499     std::cout << "Cannot get EEDigis from event!" << std::endl;
00500     return;
00501   }
00502 
00503   // EB first
00504   // 1) Look at SimTracks, getting only the HSCP tracks
00505   // 2) Match to PCaloHits
00506   // 3) Match to digis
00507   int numMatchedSimHitsEventEB = 0;
00508   int numMatchedDigisEventEB = 0;
00509   const PCaloHitContainer* phitsEB=0;
00510   phitsEB = ebSimHits.product();
00511   for(SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
00512   {
00513     // Check if the particleId is in our list
00514     std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),simTrack->type());
00515     if(partIdItr==particleIds_.end())
00516       continue;
00517 
00518     PCaloHitContainer mySimHitsEB;
00519     std::vector<EBDataFrame> myDigisEB;
00520 
00521     //int particleId = simTrack->type();
00522     int trackId = simTrack->trackId();
00523     PCaloHitContainer::const_iterator simHitItr = phitsEB->begin();
00524     while(simHitItr != phitsEB->end())
00525     {
00526       if(simHitItr->geantTrackId()==trackId)
00527         mySimHitsEB.push_back(*simHitItr);
00528       ++simHitItr;
00529     }
00530     if(mySimHitsEB.size()==0)
00531     {
00532       std::cout << "Could not find matching EB PCaloHits for SimTrack id: " << trackId << ".  Skipping this SimTrack" << std::endl;
00533       continue;
00534     }
00535 
00536     // Loop over matching PCaloHits
00537     for(simHitItr = mySimHitsEB.begin(); simHitItr != mySimHitsEB.end(); ++simHitItr)
00538     {
00539       simHitsEcalEnergyHistEB_->Fill(simHitItr->energy());
00540       simHitsEcalTimeHistEB_->Fill(simHitItr->time());
00541       simHitsEcalEnergyVsTimeHistEB_->Fill(simHitItr->time(),simHitItr->energy());
00542       EBDetId simHitId = EBDetId(simHitItr->id());
00543       std::cout << "SimHit DetId found: " << simHitId << " for PDGid: " << simTrack->type() << std::endl;
00544       //std::cout << "SimHit hashedIndex: " << simHitId.hashedIndex() << std::endl;
00545       std::cout << "SimHit energy: " << simHitItr->energy() << " time: " << simHitItr->time() << std::endl;
00546       ++numMatchedSimHitsEventEB;
00547 
00548       EBDigiCollection::const_iterator digiItr = ebDigis->begin();
00549       while(digiItr != ebDigis->end() && (digiItr->id() != simHitId))
00550         ++digiItr;
00551       if(digiItr==ebDigis->end())
00552       {
00553         // Commented out for debugging ease, Aug 3 2009
00554         std::cout << "Could not find simHit detId: " << simHitId << "in EBDigiCollection!" << std::endl;
00555         continue;
00556       }
00557       std::vector<EBDataFrame>::const_iterator myDigiItr = myDigisEB.begin();
00558       while(myDigiItr != myDigisEB.end() && (digiItr->id() != myDigiItr->id()))
00559         ++myDigiItr;
00560       if(myDigiItr!=myDigisEB.end())
00561         continue; // if this digi is already in the list, skip it
00562 
00563       ++numMatchedDigisEventEB;
00564       EBDataFrame df = *digiItr;
00565       myDigisEB.push_back(df);
00566       std::cout << "SAMPLE ADCs: " << "\t";
00567       for(int i=0; i<10;++i)
00568         std::cout << i << "\t";
00569       std::cout << std::endl << "\t\t";
00570       for(int i=0; i < df.size(); ++i)
00571       {
00572         std::cout << df.sample(i).adc() << "\t";
00573       }
00574       std::cout << std::endl << std::endl;
00575 
00576       simHitsEcalDigiMatchEnergyHistEB_->Fill(simHitItr->energy());
00577       simHitsEcalDigiMatchTimeHistEB_->Fill(simHitItr->time());
00578       simHitsEcalDigiMatchEnergyVsTimeHistEB_->Fill(simHitItr->time(),simHitItr->energy());
00579       simHitsEcalDigiMatchIEtaHist_->Fill(((EBDetId)digiItr->id()).ieta());
00580       simHitsEcalDigiMatchIPhiHist_->Fill(((EBDetId)digiItr->id()).iphi());
00581       digiOccupancyMapEB_->Fill(((EBDetId)digiItr->id()).iphi(),((EBDetId)digiItr->id()).ieta());
00582     }
00583   }
00584   simHitsEcalNumHistEB_->Fill(numMatchedSimHitsEventEB);
00585   digisEcalNumHistEB_->Fill(numMatchedDigisEventEB);
00586 
00587   // EE next
00588   int numMatchedSimHitsEventEE = 0;
00589   int numMatchedDigisEventEE = 0;
00590   const PCaloHitContainer* phitsEE=0;
00591   phitsEE = eeSimHits.product();
00592   for(SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
00593   {
00594     // Check if the particleId is in our list
00595     std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),simTrack->type());
00596     if(partIdItr==particleIds_.end())
00597       continue;
00598 
00599     PCaloHitContainer mySimHitsEE;
00600     std::vector<EEDataFrame> myDigisEE;
00601 
00602     //int particleId = simTrack->type();
00603     int trackId = simTrack->trackId();
00604     PCaloHitContainer::const_iterator simHitItr = phitsEE->begin();
00605     while(simHitItr != phitsEE->end())
00606     {
00607       if(simHitItr->geantTrackId()==trackId)
00608         mySimHitsEE.push_back(*simHitItr);
00609       ++simHitItr;
00610     }
00611     if(mySimHitsEE.size()==0)
00612     {
00613       std::cout << "Could not find matching EE PCaloHits for SimTrack id: " << trackId << ".  Skipping this SimTrack" << std::endl;
00614       continue;
00615     }
00616 
00617     // Loop over matching PCaloHits
00618     for(simHitItr = mySimHitsEE.begin(); simHitItr != mySimHitsEE.end(); ++simHitItr)
00619     {
00620       simHitsEcalEnergyHistEE_->Fill(simHitItr->energy());
00621       simHitsEcalTimeHistEE_->Fill(simHitItr->time());
00622       simHitsEcalEnergyVsTimeHistEE_->Fill(simHitItr->time(),simHitItr->energy());
00623       EEDetId simHitId = EEDetId(simHitItr->id());
00624       std::cout << "SimHit DetId found: " << simHitId << " for PDGid: " << simTrack->type() << std::endl;
00625       //std::cout << "SimHit hashedIndex: " << simHitId.hashedIndex() << std::endl;
00626       std::cout << "SimHit energy: " << simHitItr->energy() << " time: " << simHitItr->time() << std::endl;
00627       ++numMatchedSimHitsEventEE;
00628 
00629       EEDigiCollection::const_iterator digiItr = eeDigis->begin();
00630       while(digiItr != eeDigis->end() && (digiItr->id() != simHitId))
00631         ++digiItr;
00632       if(digiItr==eeDigis->end())
00633       {
00634         // Commented out for debugging ease, Aug 3 2009
00635         std::cout << "Could not find simHit detId: " << simHitId << "in EEDigiCollection!" << std::endl;
00636         continue;
00637       }
00638       std::vector<EEDataFrame>::const_iterator myDigiItr = myDigisEE.begin();
00639       while(myDigiItr != myDigisEE.end() && (digiItr->id() != myDigiItr->id()))
00640         ++myDigiItr;
00641       if(myDigiItr!=myDigisEE.end())
00642         continue; // if this digi is already in the list, skip it
00643 
00644       ++numMatchedDigisEventEE;
00645       EEDataFrame df = *digiItr;
00646       myDigisEE.push_back(df);
00647       std::cout << "SAMPLE ADCs: " << "\t";
00648       for(int i=0; i<10;++i)
00649         std::cout << i << "\t";
00650       std::cout << std::endl << "\t\t";
00651       for(int i=0; i < df.size(); ++i)
00652       {
00653         std::cout << df.sample(i).adc() << "\t";
00654       }
00655       std::cout << std::endl << std::endl;
00656 
00657       simHitsEcalDigiMatchEnergyHistEE_->Fill(simHitItr->energy());
00658       simHitsEcalDigiMatchTimeHistEE_->Fill(simHitItr->time());
00659       simHitsEcalDigiMatchEnergyVsTimeHistEE_->Fill(simHitItr->time(),simHitItr->energy());
00660       if(((EEDetId)digiItr->id()).zside() > 0)
00661         digiOccupancyMapEEP_->Fill(((EEDetId)digiItr->id()).ix(),((EEDetId)digiItr->id()).iy());
00662       else if(((EEDetId)digiItr->id()).zside() < 0)
00663         digiOccupancyMapEEM_->Fill(((EEDetId)digiItr->id()).ix(),((EEDetId)digiItr->id()).iy());
00664     }
00665   }
00666   simHitsEcalNumHistEE_->Fill(numMatchedSimHitsEventEE);
00667   digisEcalNumHistEE_->Fill(numMatchedDigisEventEE);
00668 
00669 }
00670 // ------------- Make Reco plots ---------------------------------------------------------
00671 void HSCPValidator::makeRecoPlots(const edm::Event& iEvent)
00672 {
00673   using namespace edm;
00674    using namespace reco;
00675 
00676   Handle<HepMCProduct> evt;
00677   iEvent.getByLabel(label_, evt);
00678 
00679   Handle<TrackCollection> tkTracks;
00680   iEvent.getByLabel("generalTracks",tkTracks);
00681   const reco::TrackCollection tkTC = *(tkTracks.product());
00682   
00683   Handle<ValueMap<DeDxData> >          dEdxTrackHandle;
00684   iEvent.getByLabel("dedxHarmonic2", dEdxTrackHandle);
00685   const ValueMap<DeDxData> dEdxTrack = *dEdxTrackHandle.product();
00686   
00687   for(size_t i=0; i<tkTracks->size(); i++){
00688      
00689      reco::TrackRef trkRef = reco::TrackRef(tkTracks, i);
00690         
00691      if(trkRef->pt()<5 || trkRef->normalizedChi2()>10) continue;
00692  
00693      double minR= 999;
00694      double hscpgenPt =-1;
00695      
00696      HepMC::GenEvent * myGenEvent = new  HepMC::GenEvent(*(evt->GetEvent()));
00697      for(HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin();
00698          p != myGenEvent->particles_end(); ++p )
00699      {
00700         
00701         if((*p)->status() != particleStatus_)
00702            continue;
00703         // Check if the particleId is in our R-hadron list
00704         std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),(*p)->pdg_id());
00705         if(partIdItr!=particleIds_.end()){
00706            
00707            //calculate DeltaR
00708            double distance =pow((*p)->momentum().eta()-trkRef->eta(),2)+pow((*p)->momentum().phi()-trkRef->phi(),2);
00709            distance =sqrt(distance);
00710            if(distance <minR ){
00711               minR = distance;
00712               hscpgenPt= (*p)->momentum().perp();
00713            }
00714         }
00715      }
00716      RecoHSCPPtVsGenPt->Fill(trkRef->pt(),hscpgenPt);
00717  
00718      delete myGenEvent;     
00719      double dedx = dEdxTrack[trkRef].dEdx();
00720      dedxVsp->Fill( trkRef->p(),dedx);
00721        
00722   }  
00723 
00724 }
00725 
00726 // ------------- Make simDigi plots RPC -------------------------------------------------
00727 void HSCPValidator::makeSimDigiPlotsRPC(const edm::Event& iEvent)
00728 {
00729   using namespace edm;
00730 
00731   //std::cout << " Getting the SimHits " <<std::endl;
00732   std::vector<Handle<edm::PSimHitContainer> > theSimHitContainers;
00733   iEvent.getManyByType(theSimHitContainers);
00734   //std::cout << " The Number of sim Hits is  " << theSimHitContainers.size() <<std::endl;
00735 
00736   Handle<RPCRecHitCollection> rpcRecHits;
00737   iEvent.getByLabel("rpcRecHits","",rpcRecHits);
00738 
00739   //SimTrack Stuff
00740   std::vector<PSimHit> theSimHits;
00741 
00742   for (int i = 0; i < int(theSimHitContainers.size()); i++){
00743     theSimHits.insert(theSimHits.end(),theSimHitContainers.at(i)->begin(),theSimHitContainers.at(i)->end());
00744   } 
00745 
00746 
00747   for (std::vector<PSimHit>::const_iterator iHit = theSimHits.begin(); iHit != theSimHits.end(); iHit++){
00748 
00749     std::vector<int>::const_iterator partIdItr = find(particleIds_.begin(),particleIds_.end(),(*iHit).particleType());
00750     if(partIdItr==particleIds_.end())
00751       continue;
00752 
00753     DetId theDetUnitId((*iHit).detUnitId());
00754 
00755     DetId simdetid= DetId((*iHit).detUnitId());
00756 
00757     if(simdetid.det()==DetId::Muon &&  simdetid.subdetId()== MuonSubdetId::RPC){//Only RPCs
00758 
00759       RPCDetId rollId(theDetUnitId);
00760       RPCGeomServ rpcsrv(rollId);
00761 
00762       //std::cout << " Reading the Roll"<<std::endl;
00763       const RPCRoll* rollasociated = rpcGeo->roll(rollId);
00764 
00765       //std::cout << " Getting the Surface"<<std::endl;
00766       const BoundPlane & RPCSurface = rollasociated->surface(); 
00767 
00768       GlobalPoint SimHitInGlobal = RPCSurface.toGlobal((*iHit).localPosition());
00769 
00770       std::cout<<"\t\t We have an RPC Sim Hit! in t="<<(*iHit).timeOfFlight()<<"ns "<<rpcsrv.name()<<" Global postition="<<SimHitInGlobal<<std::endl;
00771 
00772       int layer = 0;
00773 
00774       if(rollId.station()==1&&rollId.layer()==1) layer = 1;
00775       else if(rollId.station()==1&&rollId.layer()==2) layer = 2;
00776       else if(rollId.station()==2&&rollId.layer()==1) layer = 3;
00777       else if(rollId.station()==2&&rollId.layer()==2)  layer = 4;
00778       else if(rollId.station()==3) layer = 5;
00779       else if(rollId.station()==4) layer = 6;
00780 
00781       if(rollId.region()==0){
00782         rpcTimeOfFlightBarrel_[layer-1]->Fill((*iHit).timeOfFlight());
00783       }else{
00784         rpcTimeOfFlightEndCap_[rollId.station()-1]->Fill((*iHit).timeOfFlight());
00785       }
00786 
00787       std::cout<<"\t\t r="<<SimHitInGlobal.mag()<<" phi="<<SimHitInGlobal.phi()<<" eta="<<SimHitInGlobal.eta()<<std::endl;
00788 
00789       int cluSize = 0;
00790       int bx = 100;
00791       float minres = 3000.;
00792 
00793       std::cout<<"\t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00794 
00795       typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00796       rangeRecHits recHitCollection =  rpcRecHits->get(rollId);
00797       RPCRecHitCollection::const_iterator recHit;
00798 
00799       efficiencyRPCRecHitSimDigis_->Fill(0);
00800 
00801       for(recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++){
00802         LocalPoint recHitPos=recHit->localPosition();
00803         float res=(*iHit).localPosition().x()- recHitPos.x();
00804         if(fabs(res)<fabs(minres)){
00805           minres=res;
00806           cluSize = recHit->clusterSize();
00807           bx=recHit->BunchX();
00808           std::cout<<"\t New Min Res "<<res<<"cm."<<std::endl;
00809         }
00810       }
00811 
00812       if(minres<3000.){
00813         residualsRPCRecHitSimDigis_->Fill(minres);
00814         efficiencyRPCRecHitSimDigis_->Fill(1);
00815         cluSizeDistribution_->Fill(cluSize);
00816         if(rollId.region()==0) rpcBXBarrel_[layer-1]->Fill(bx);
00817         else rpcBXEndCap_[rollId.station()-1]->Fill(bx);
00818       }
00819     }
00820   }
00821 }
00822 
00823 // ------------- Convert int to string for printing -------------------------------------
00824 std::string HSCPValidator::intToString(int num)
00825 {
00826   using namespace std;
00827   ostringstream myStream;
00828   myStream << num << flush;
00829   return(myStream.str()); //returns the string form of the stringstream object
00830 }
00831 
00832 
00833 
00834 //------Increase trigger thresold----
00835 
00836 
00837 bool HSCPValidator::IncreasedTreshold(const trigger::TriggerEvent& trEv, const edm::InputTag& InputPath, double NewThreshold, double etaCut, int NObjectAboveThreshold, bool averageThreshold)
00838 {
00839    unsigned int filterIndex = trEv.filterIndex(InputPath);
00840    //if(filterIndex<trEv.sizeFilters())printf("SELECTED INDEX =%i --> %s    XXX   %s\n",filterIndex,trEv.filterTag(filterIndex).label().c_str(), trEv.filterTag(filterIndex).process().c_str());
00841 
00842    if (filterIndex<trEv.sizeFilters()){
00843       const trigger::Vids& VIDS(trEv.filterIds(filterIndex));
00844       const trigger::Keys& KEYS(trEv.filterKeys(filterIndex));
00845       const int nI(VIDS.size());
00846       const int nK(KEYS.size());
00847       assert(nI==nK);
00848       const int n(std::max(nI,nK));
00849       const trigger::TriggerObjectCollection& TOC(trEv.getObjects());
00850 
00851 
00852       if(!averageThreshold){
00853          int NObjectAboveThresholdObserved = 0;
00854          for (int i=0; i!=n; ++i) {
00855            if(TOC[KEYS[i]].pt()> NewThreshold && fabs(TOC[KEYS[i]].eta())<etaCut) NObjectAboveThresholdObserved++;
00856             //cout << "   " << i << " " << VIDS[i] << "/" << KEYS[i] << ": "<< TOC[KEYS[i]].id() << " " << TOC[KEYS[i]].pt() << " " << TOC[KEYS[i]].eta() << " " << TOC[KEYS[i]].phi() << " " << TOC[KEYS[i]].mass()<< endl;
00857          }
00858          if(NObjectAboveThresholdObserved>=NObjectAboveThreshold)return true;
00859 
00860       }else{
00861          std::vector<double> ObjPt;
00862 
00863          for (int i=0; i!=n; ++i) {
00864             ObjPt.push_back(TOC[KEYS[i]].pt());
00865             //cout << "   " << i << " " << VIDS[i] << "/" << KEYS[i] << ": "<< TOC[KEYS[i]].id() << " " << TOC[KEYS[i]].pt() << " " << TOC[KEYS[i]].eta() << " " << TOC[KEYS[i]].phi() << " " << TOC[KEYS[i]].mass()<< endl;
00866          }
00867          if((int)(ObjPt.size())<NObjectAboveThreshold)return false;
00868          std::sort(ObjPt.begin(), ObjPt.end());
00869 
00870          double Average = 0;
00871          for(int i=0; i<NObjectAboveThreshold;i++){
00872             Average+= ObjPt[ObjPt.size()-1-i];
00873          }Average/=NObjectAboveThreshold;
00874          //cout << "AVERAGE = " << Average << endl;
00875 
00876          if(Average>NewThreshold)return true;
00877       }
00878    }
00879    return false;
00880 }
00881 
00882 
00883 //define this as a plug-in
00884 DEFINE_FWK_MODULE(HSCPValidator);
Generated for CMSSW Reference Manual by doxygen 1.7.3