Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    IsolatedTracksCone
00004 // Class:      IsolatedTracksCone
00005 // 
00006 
00015 //
00016 // Original Author: Jim Hirschauer (adaptation of Seema Sharma's
00017 // IsolatedTracksNew)
00018 //         Created:  Thu Nov  6 15:30:40 CST 2008
00019 //
00020 //
00021 #define _DEBUG_QUIET
00022 
00023 #include "Calibration/IsolatedParticles/plugins/IsolatedTracksCone.h"
00024 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
00025 
00026 #include "DataFormats/HLTReco/interface/TriggerObject.h"
00027 #include "FWCore/Common/interface/TriggerNames.h"
00028 #include "FWCore/Framework/interface/LuminosityBlock.h"
00029 #include "DataFormats/Common/interface/TriggerResults.h"
00030 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
00031 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
00032 
00033 IsolatedTracksCone::IsolatedTracksCone(const edm::ParameterSet& iConfig) {
00034 
00035   //now do what ever initialization is needed
00036   doMC            = iConfig.getUntrackedParameter<bool>  ("DoMC", false); 
00037   myverbose_      = 
00038     iConfig.getUntrackedParameter<int>( "Verbosity", 5 );
00039   useJetTrigger_  = 
00040     iConfig.getUntrackedParameter<bool>( "useJetTrigger", false);
00041   drLeadJetVeto_  = 
00042     iConfig.getUntrackedParameter<double>( "drLeadJetVeto",  1.2 );
00043   ptMinLeadJet_   = 
00044     iConfig.getUntrackedParameter<double>( "ptMinLeadJet",  15.0 );
00045 
00046   debugTrks_          = 
00047     iConfig.getUntrackedParameter<int>("DebugTracks");
00048   printTrkHitPattern_ = 
00049     iConfig.getUntrackedParameter<bool>("PrintTrkHitPattern");
00050   
00051   minTrackP_     = 
00052     iConfig.getUntrackedParameter<double>( "minTrackP", 10.0);
00053   maxTrackEta_    = 
00054     iConfig.getUntrackedParameter<double>( "maxTrackEta", 5.0);
00055   maxNearTrackP_  = 
00056     iConfig.getUntrackedParameter<double>( "maxNearTrackP", 1.0);
00057 
00058   debugEcalSimInfo_   = 
00059     iConfig.getUntrackedParameter<int>("DebugEcalSimInfo");
00060 
00061   applyEcalIsolation_ = 
00062     iConfig.getUntrackedParameter<bool>("ApplyEcalIsolation");
00063 
00064   _L1extraTauJetSource = 
00065     iConfig.getParameter<edm::InputTag>("L1extraTauJetSource");
00066   _L1extraCenJetSource = 
00067     iConfig.getParameter<edm::InputTag>("L1extraCenJetSource");
00068   _L1extraFwdJetSource = 
00069     iConfig.getParameter<edm::InputTag>("L1extraFwdJetSource");
00070 
00071   edm::ParameterSet parameters = 
00072     iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00073   parameters_.loadParameters( parameters );
00074   trackAssociator_ =  new TrackDetectorAssociator();
00075   trackAssociator_->useDefaultPropagator();
00076 
00077   if(myverbose_>=0) {
00078     std::cout <<"Parameters read from config file \n" 
00079               << "myverbose_ "          << myverbose_     << "\n"     
00080               << "useJetTrigger_ "      << useJetTrigger_ << "\n"
00081               << "drLeadJetVeto_ "      << drLeadJetVeto_ << "\n"
00082               << "minTrackP_ "         << minTrackP_    << "\n"
00083               << "maxTrackEta_ "        << maxTrackEta_   << "\n"
00084               << "maxNearTrackP_ "      << maxNearTrackP_ 
00085               << std::endl;
00086   }
00087 }
00088 
00089 IsolatedTracksCone::~IsolatedTracksCone() {
00090   delete  trackAssociator_;
00091 }
00092 
00093 void IsolatedTracksCone::analyze(const edm::Event& iEvent, 
00094                                  const edm::EventSetup& iSetup) {
00095   
00096   unsigned int irun = (unsigned int)iEvent.id().run();
00097   unsigned int ilum = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
00098   unsigned int ievt = (unsigned int)iEvent.id().event();
00099 
00100   
00101   
00102   clearTrackVectors();
00103   
00104   // check the L1 objects
00105   bool   L1Pass = false;
00106   leadL1JetPT=-999, leadL1JetEta=-999,  leadL1JetPhi=-999;
00107   if( !useJetTrigger_ ) {
00108     L1Pass = true;
00109   } else {
00110     edm::Handle<l1extra::L1JetParticleCollection> l1TauHandle;
00111     iEvent.getByLabel(_L1extraTauJetSource,l1TauHandle);
00112     l1extra::L1JetParticleCollection::const_iterator itr;
00113     for(itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr ) 
00114     {
00115       if( itr->pt()>leadL1JetPT ) {
00116         leadL1JetPT  = itr->pt();
00117         leadL1JetEta = itr->eta();
00118         leadL1JetPhi = itr->phi();
00119       }
00120     }
00121     edm::Handle<l1extra::L1JetParticleCollection> l1CenJetHandle;
00122     iEvent.getByLabel(_L1extraCenJetSource,l1CenJetHandle);
00123     for( itr = l1CenJetHandle->begin();  itr != l1CenJetHandle->end(); ++itr ) 
00124     {
00125       if( itr->pt()>leadL1JetPT ) {
00126         leadL1JetPT  = itr->pt();
00127         leadL1JetEta = itr->eta();
00128         leadL1JetPhi = itr->phi();
00129       }
00130     }
00131     edm::Handle<l1extra::L1JetParticleCollection> l1FwdJetHandle;
00132     iEvent.getByLabel(_L1extraFwdJetSource,l1FwdJetHandle);
00133     for( itr = l1FwdJetHandle->begin();  itr != l1FwdJetHandle->end(); ++itr ) 
00134     {
00135       if( itr->pt()>leadL1JetPT ) {
00136         leadL1JetPT  = itr->pt();
00137         leadL1JetEta = itr->eta();
00138         leadL1JetPhi = itr->phi();
00139       }
00140     }
00141     if(leadL1JetPT>ptMinLeadJet_) 
00142     {
00143       L1Pass = true;
00144     }
00145     
00146   }
00147   
00148 
00150   // Break now if L1Pass is false
00152   //   if (!L1Pass) {
00153   //     nEVT_failL1++;
00154   //     //    std::cout << "L1Pass is false : " << L1Pass << std::endl;
00155   //     return;  
00156   //   }
00157   
00159   // Get the collection handles
00161   
00162 
00163   edm::ESHandle<CaloGeometry> pG;
00164   iSetup.get<CaloGeometryRecord>().get(pG);
00165   const CaloGeometry* geo = pG.product();
00166   const CaloSubdetectorGeometry* gEB = 
00167     geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00168   const CaloSubdetectorGeometry* gEE = 
00169     geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00170   const CaloSubdetectorGeometry* gHB = 
00171     geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00172   const CaloSubdetectorGeometry* gHE = 
00173     geo->getSubdetectorGeometry(DetId::Hcal,HcalEndcap);
00174   
00175   edm::ESHandle<CaloTopology> theCaloTopology;
00176   iSetup.get<CaloTopologyRecord>().get(theCaloTopology); 
00177   const CaloTopology *caloTopology = theCaloTopology.product();
00178   //  const CaloSubdetectorTopology* theEBTopology   = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel);
00179   //  const CaloSubdetectorTopology* theEETopology   = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalEndcap);
00180   
00181   edm::ESHandle<HcalTopology> htopo;
00182   iSetup.get<IdealGeometryRecord>().get(htopo);
00183   const HcalTopology* theHBHETopology = htopo.product();
00184   
00185   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
00186   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
00187   iEvent.getByLabel("ecalRecHit","EcalRecHitsEB",barrelRecHitsHandle);
00188   iEvent.getByLabel("ecalRecHit","EcalRecHitsEE",endcapRecHitsHandle);
00189   
00190   edm::Handle<HBHERecHitCollection> hbhe;
00191   iEvent.getByLabel("hbhereco",hbhe);
00192   const HBHERecHitCollection Hithbhe = *(hbhe.product());
00193 
00194   edm::Handle<reco::TrackCollection> trkCollection;
00195   iEvent.getByLabel("generalTracks", trkCollection);
00196   reco::TrackCollection::const_iterator trkItr;
00197   if(debugTrks_>1){
00198     std::cout << "Track Collection: " << std::endl;
00199     std::cout << "Number of Tracks " << trkCollection->size() << std::endl;
00200   }
00201   std::string theTrackQuality = "highPurity";
00202   reco::TrackBase::TrackQuality trackQuality_=
00203     reco::TrackBase::qualityByName(theTrackQuality);
00204   
00205   //get Handles to SimTracks and SimHits
00206   edm::Handle<edm::SimTrackContainer> SimTk;
00207   if (doMC) iEvent.getByLabel("g4SimHits",SimTk);
00208   edm::SimTrackContainer::const_iterator simTrkItr;
00209 
00210   edm::Handle<edm::SimVertexContainer> SimVtx;
00211   if (doMC) iEvent.getByLabel("g4SimHits",SimVtx);
00212   edm::SimVertexContainer::const_iterator vtxItr = SimVtx->begin();
00213 
00214   //get Handles to PCaloHitContainers of eb/ee/hbhe
00215   edm::Handle<edm::PCaloHitContainer> pcaloeb;
00216   if (doMC) iEvent.getByLabel("g4SimHits", "EcalHitsEB", pcaloeb);
00217 
00218   edm::Handle<edm::PCaloHitContainer> pcaloee;
00219   if (doMC) iEvent.getByLabel("g4SimHits", "EcalHitsEE", pcaloee);
00220 
00221   edm::Handle<edm::PCaloHitContainer> pcalohh;
00222   if (doMC) iEvent.getByLabel("g4SimHits", "HcalHits", pcalohh);
00223   
00224   
00225   
00227   // Get HLT_IsoTrackHB/HE Information
00229     
00230   edm::InputTag theTriggerResultsLabel ("TriggerResults","","HLT");
00231   edm::Handle<edm::TriggerResults> triggerResults;
00232   iEvent.getByLabel( theTriggerResultsLabel, triggerResults);
00233 
00234   
00235 
00236   std::vector<int> v_hlTriggers;
00237   int hltHB(-99);
00238   int hltHE(-99);
00239   int hltL1Jet15                        (-99);
00240   int hltJet30                  (-99);
00241   int hltJet50                  (-99);
00242   int hltJet80                  (-99);
00243   int hltJet110                 (-99);
00244   int hltJet140                 (-99);
00245   int hltJet180                 (-99);
00246   int hltL1SingleEG5            (-99);
00247   int hltZeroBias               (-99);
00248   int hltMinBiasHcal            (-99);
00249   int hltMinBiasEcal            (-99);
00250   int hltMinBiasPixel           (-99);
00251   int hltSingleIsoTau30_Trk5    (-99);
00252   int hltDoubleLooseIsoTau15_Trk5(-99);
00253 
00254   if (triggerResults.isValid()) {
00255     
00256     const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
00257     // TriggerNames class  triggerNames.init(*triggerResults);
00258     
00259 
00260     for (unsigned int i=0; i<triggerResults->size(); i++){
00261       //      std::cout << "triggerNames.triggerName(" << i << ") = " << triggerNames.triggerName(i) << std::endl;
00262       if (triggerNames.triggerName(i) == "HLT_IsoTrackHE_1E31") hltHE = triggerResults->accept(i);
00263       if (triggerNames.triggerName(i) == "HLT_IsoTrackHB_1E31") hltHB = triggerResults->accept(i);
00264       if (triggerNames.triggerName(i) == "HLT_L1Jet15"                  ) hltL1Jet15                     = triggerResults->accept(i);
00265       if (triggerNames.triggerName(i) == "HLT_Jet30"                    ) hltJet30                       = triggerResults->accept(i);
00266       if (triggerNames.triggerName(i) == "HLT_Jet50"                    ) hltJet50                       = triggerResults->accept(i);
00267       if (triggerNames.triggerName(i) == "HLT_Jet80"                    ) hltJet80                       = triggerResults->accept(i);
00268       if (triggerNames.triggerName(i) == "HLT_Jet110"                   ) hltJet110                      = triggerResults->accept(i);
00269       if (triggerNames.triggerName(i) == "HLT_Jet140"                   ) hltJet140                      = triggerResults->accept(i);
00270       if (triggerNames.triggerName(i) == "HLT_Jet180"                   ) hltJet180                      = triggerResults->accept(i);
00271       if (triggerNames.triggerName(i) == "HLT_L1SingleEG5"              ) hltL1SingleEG5                 = triggerResults->accept(i);
00272       if (triggerNames.triggerName(i) == "HLT_ZeroBias"                 ) hltZeroBias                    = triggerResults->accept(i);
00273       if (triggerNames.triggerName(i) == "HLT_MinBiasHcal"              ) hltMinBiasHcal                 = triggerResults->accept(i);
00274       if (triggerNames.triggerName(i) == "HLT_MinBiasEcal"              ) hltMinBiasEcal                 = triggerResults->accept(i);
00275       if (triggerNames.triggerName(i) == "HLT_MinBiasPixel"             ) hltMinBiasPixel                = triggerResults->accept(i);
00276       if (triggerNames.triggerName(i) == "HLT_SingleIsoTau30_Trk5"      ) hltSingleIsoTau30_Trk5         = triggerResults->accept(i);
00277       if (triggerNames.triggerName(i) == "HLT_DoubleLooseIsoTau15_Trk5" ) hltDoubleLooseIsoTau15_Trk5    = triggerResults->accept(i);
00278     }
00279   }
00280   
00281     
00282 
00283   
00285   // Primary loop over tracks
00287   TrackerHitAssociator* associate=0;
00288   if (doMC) associate = new TrackerHitAssociator(iEvent);
00289   
00290 
00291   nTRK      = 0;
00292   nRawTRK   = 0;
00293   nFailPt   = 0;
00294   nFailEta  = 0;
00295   nFailHighPurityQaul = 0;
00296   nMissEcal = 0;
00297   nMissHcal = 0;
00298 
00299   for( trkItr = trkCollection->begin(); 
00300        trkItr != trkCollection->end(); ++trkItr)
00301   {
00302 
00303     nRawTRK++;
00304 
00305     const reco::Track* pTrack = &(*trkItr);
00306 
00308     // Check for min Pt and max Eta P
00310 
00311     bool trkQual  = pTrack->quality(trackQuality_);
00312     bool goodPt   = pTrack->p()>minTrackP_;
00313     bool goodEta  = std::abs(pTrack->momentum().eta())<maxTrackEta_;
00314 
00315     double eta1       = pTrack->momentum().eta();
00316     double phi1       = pTrack->momentum().phi();
00317     double pt1        = pTrack->pt();
00318     double p1         = pTrack->p();
00319 
00320 
00321     if (!goodEta){
00322       nFailEta++;
00323     }
00324     if (!goodPt){
00325       nFailPt++;
00326     }
00327     if (!trkQual){
00328       nFailHighPurityQaul++;
00329     }
00330     
00331     hRawPt ->Fill(pt1 );
00332     hRawP  ->Fill(p1  );
00333     hRawEta->Fill(eta1);
00334     hRawPhi->Fill(phi1);
00335       
00336     if( !goodEta || !goodPt || !trkQual ) continue; // Skip to next track
00337     
00339     // Find track trajectory
00341     
00342       
00343     const FreeTrajectoryState fts1 = 
00344       trackAssociator_->getFreeTrajectoryState(iSetup, *pTrack);
00345     
00346       
00347     TrackDetMatchInfo info1 = 
00348       trackAssociator_->associate(iEvent, iSetup, fts1, parameters_);
00349     
00350       
00351 
00353     // First confirm track makes it to Hcal
00355 
00356     if (info1.trkGlobPosAtHcal.x()==0 && 
00357         info1.trkGlobPosAtHcal.y()==0 && 
00358         info1.trkGlobPosAtHcal.z()==0) 
00359     {
00360       nMissHcal++;
00361       continue;      
00362     }
00363     
00364     const GlobalPoint hpoint1(info1.trkGlobPosAtHcal.x(),
00365                               info1.trkGlobPosAtHcal.y(),
00366                               info1.trkGlobPosAtHcal.z());
00367 
00368       
00369 
00371     // Get basic quantities
00373 
00374     const reco::HitPattern& hitp = pTrack->hitPattern();
00375     int nLayersCrossed = hitp.trackerLayersWithMeasurement();        
00376     int nOuterHits     = hitp.stripTOBLayersWithMeasurement()
00377       +hitp.stripTECLayersWithMeasurement() ;
00378 
00379     
00380     double simP = 0;
00381     if (doMC) {
00382       edm::SimTrackContainer::const_iterator matchedSimTrk = 
00383         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
00384       simP = matchedSimTrk->momentum().P();
00385     }
00387     // Get Ecal Point
00389 
00390     const GlobalPoint point1(info1.trkGlobPosAtEcal.x(),
00391                              info1.trkGlobPosAtEcal.y(),
00392                              info1.trkGlobPosAtEcal.z());
00393     
00394 
00395     // Sanity check that track hits Ecal
00396 
00397     if (info1.trkGlobPosAtEcal.x()==0 && 
00398         info1.trkGlobPosAtEcal.y()==0 && 
00399         info1.trkGlobPosAtEcal.z()==0) 
00400     {
00401       std::cout << "Track doesn't reach Ecal." << std::endl;
00402       nMissEcal++;
00403       continue;
00404     }
00405 
00406     // Get Track Momentum - make sure you have latest version of TrackDetMatchInfo
00407     
00408     GlobalVector trackMomAtEcal = info1.trkMomAtEcal;
00409     GlobalVector trackMomAtHcal = info1.trkMomAtHcal;
00410 
00412     // If using Jet trigger, get distance from leading jet
00414 
00415     double drFromLeadJet = 999.0;
00416     if( useJetTrigger_ ) {
00417       double dphi = DeltaPhi(phi1, leadL1JetPhi);
00418       double deta = eta1 - leadL1JetEta;
00419       drFromLeadJet =  sqrt(dphi*dphi + deta*deta);
00420     }
00421     
00422 
00424     // Define Arrays for sizes of Charge, Neut Iso Radii and
00425     // Clustering Cone Radius.
00427 
00428     const int a_size = 7;
00429     double a_coneR[a_size];
00430     double a_charIsoR[a_size];
00431     double a_neutIsoR[a_size];
00432 
00433     a_coneR[0] = 17.49; // = area of 2x2
00434     a_coneR[1] = 26.23; // = area of 3x3
00435     a_coneR[2] = 30.61;
00436     a_coneR[3] = 34.98; // = area of 4x4
00437     a_coneR[4] = 39.35;
00438     a_coneR[5] = 43.72; // = area of 5x5
00439     a_coneR[6] = 52.46; // = area of 6x6
00440     
00441     for (int i=0; i<a_size; i++){
00442       a_charIsoR[i] = a_coneR[i]+28.9; // 28.9 gives 55.1 for 3x3 benchmark 
00443       a_neutIsoR[i] = a_charIsoR[i]*0.726; // Ecal radius = 0.726*Hcal radius
00444     }
00445     
00447     // Do Neutral Iso in radius on Ecal surface.
00449   
00450     // NxN cluster
00451     double e3x3=-999.0;
00452     double trkEcalEne =-999.0;
00453     if(std::abs(point1.eta())<1.479) {
00454       const DetId isoCell = gEB->getClosestCell(point1);
00455       e3x3   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, geo, caloTopology, 1,1);  
00456       trkEcalEne   = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
00457     } else {
00458       const DetId isoCell = gEE->getClosestCell(point1);
00459       e3x3   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, geo, caloTopology, 1,1);
00460       trkEcalEne   = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
00461     }
00462 
00463     // Cone cluster
00464 
00465     // Set up array of cone sizes for MIP cut
00466     const int a_mip_size = 5;
00467     double a_mipR[a_mip_size];
00468     a_mipR[0] = 3.84; // = area of 3x3 ecal
00469     a_mipR[1] = 14.0;
00470     a_mipR[2] = 19.0;
00471     a_mipR[3] = 24.0;
00472     a_mipR[4] = 9.0;  // = from standard analyzer
00473 
00474     std::vector<double> v_eDR;
00475     for (int i = 0 ; i < a_size ; i++){
00476       int nRH_eDR = 0;
00477 
00478       // Cone in ecal
00479       double eDR = spr::eCone_ecal(geo, 
00480                                    barrelRecHitsHandle, 
00481                                    endcapRecHitsHandle, 
00482                                    hpoint1, point1, 
00483                                    a_neutIsoR[i],  
00484                                    trackMomAtEcal, nRH_eDR);
00485       v_eDR.push_back(eDR);
00486       
00487     }
00488 
00489     std::vector<double> v_eMipDR;
00490     for (int i = 0 ; i < a_mip_size ; i++){
00491       int nRH_eMipDR = 0;
00492       double eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle, 
00493                                       endcapRecHitsHandle, 
00494                                       hpoint1, point1, 
00495                                       a_mipR[i], trackMomAtEcal, nRH_eMipDR);
00496       
00497       v_eMipDR.push_back(eMipDR);
00498     }
00499     
00500             
00502     // Do charge isolation in radius at Hcal surface for 5 different
00503     // radii defined above in a_charIso
00505     
00506     std::vector<double> v_hmaxNearP_goodTrk;
00507     std::vector<double> v_hmaxNearP        ;
00508     std::vector<int>    v_hnNearTRKs       ;
00509     std::vector<int>    v_hnLayers_maxNearP;
00510     std::vector<int>    v_htrkQual_maxNearP;
00511 
00512     std::vector<double> v_cone_hmaxNearP_goodTrk;
00513     std::vector<double> v_cone_hmaxNearP        ;
00514     std::vector<int>    v_cone_hnNearTRKs       ;
00515     std::vector<int>    v_cone_hnLayers_maxNearP;
00516     std::vector<int>    v_cone_htrkQual_maxNearP;
00517 
00518     for (int i = 0 ; i < a_size ; i++){
00519 
00520       double hmaxNearP         = -999.0;
00521       int    hnNearTRKs        =  0;
00522       int    hnLayers_maxNearP =  0;
00523       int    htrkQual_maxNearP = -1; 
00524       double hmaxNearP_goodTrk = -999.0;
00525 
00526       double conehmaxNearP         = -999.0;
00527       int    conehnNearTRKs        =  0;
00528       int    conehnLayers_maxNearP =  0;
00529       int    conehtrkQual_maxNearP = -1; 
00530       double conehmaxNearP_goodTrk = -999.0;
00531 
00532       conehmaxNearP = spr::coneChargeIsolation(iEvent, iSetup, 
00533                                                trkItr, trkCollection, 
00534                                                *trackAssociator_, parameters_, 
00535                                                theTrackQuality, conehnNearTRKs,
00536                                                conehnLayers_maxNearP,
00537                                                conehtrkQual_maxNearP, 
00538                                                conehmaxNearP_goodTrk, 
00539                                                hpoint1, trackMomAtHcal, 
00540                                                a_charIsoR[i]); 
00541 
00542       v_hmaxNearP_goodTrk.push_back(hmaxNearP_goodTrk);
00543       v_hmaxNearP        .push_back(hmaxNearP        );
00544       v_hnNearTRKs       .push_back(hnNearTRKs       );
00545       v_hnLayers_maxNearP.push_back(hnLayers_maxNearP);
00546       v_htrkQual_maxNearP.push_back(htrkQual_maxNearP);
00547 
00548       v_cone_hmaxNearP_goodTrk.push_back(conehmaxNearP_goodTrk);
00549       v_cone_hmaxNearP        .push_back(conehmaxNearP        );
00550       v_cone_hnNearTRKs       .push_back(conehnNearTRKs       );
00551       v_cone_hnLayers_maxNearP.push_back(conehnLayers_maxNearP);
00552       v_cone_htrkQual_maxNearP.push_back(conehtrkQual_maxNearP);
00553       
00554     }
00555     
00556     double h3x3=-999.0, h5x5=-999.0;
00557     double hsim3x3=-999.0, hsim5x5=-999.0, trkHcalEne=-999.0;
00558     std::map<std::string, double> hsimInfo3x3, hsimInfo5x5;
00559     double distFromHotCell_h3x3 = -99.;
00560     int ietaFromHotCell_h3x3    = -99;
00561     int iphiFromHotCell_h3x3    = -99;
00562     double distFromHotCell_h5x5 = -99.;
00563     int ietaFromHotCell_h5x5    = -99;
00564     int iphiFromHotCell_h5x5    = -99;
00565 
00566     GlobalPoint gPosHotCell_h3x3(0.,0.,0.);
00567     GlobalPoint gPosHotCell_h5x5(0.,0.,0.);
00568 
00569     int nRH_h3x3(0), nRH_h5x5(0);
00570 
00571     // Hcal Energy Clustering
00572     
00573     // Get closetcell for ietaFromHotCell and iphiFromHotCell
00574     DetId ClosestCell;
00575     if( std::abs(pTrack->eta())<1.392 ) {
00576       ClosestCell = gHB->getClosestCell(hpoint1);
00577     } else {
00578       ClosestCell = gHE->getClosestCell(hpoint1);
00579     }
00580     // Transform into HcalDetId so that I can get ieta, iphi later.
00581     HcalDetId ClosestCell_HcalDetId(ClosestCell.rawId());
00582 
00583     // Using NxN Cluster
00584     std::vector<int>    v_RH_h3x3_ieta;
00585     std::vector<int>    v_RH_h3x3_iphi;
00586     std::vector<double> v_RH_h3x3_ene;
00587     std::vector<int>    v_RH_h5x5_ieta;
00588     std::vector<int>    v_RH_h5x5_iphi;
00589     std::vector<double> v_RH_h5x5_ene;
00590     
00591     
00592     h3x3 = spr::eHCALmatrix(geo, theHBHETopology, ClosestCell, hbhe,1,1,
00593                             nRH_h3x3, v_RH_h3x3_ieta, v_RH_h3x3_iphi, v_RH_h3x3_ene, 
00594                             gPosHotCell_h3x3);
00595     distFromHotCell_h3x3 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h3x3);
00596     
00597     h5x5 = spr::eHCALmatrix(geo, theHBHETopology, ClosestCell, hbhe,2,2,
00598                             nRH_h5x5, v_RH_h5x5_ieta, v_RH_h5x5_iphi,  v_RH_h5x5_ene, 
00599                             gPosHotCell_h5x5);
00600     distFromHotCell_h5x5 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h5x5);
00601     
00602 
00603     //     double heta = (double)hpoint1.eta();
00604     //     double hphi = (double)hpoint1.phi();
00605     std::vector<int> multiplicity3x3;
00606     std::vector<int> multiplicity5x5;
00607     if (doMC) {
00608       hsim3x3   = spr::eHCALmatrix(theHBHETopology, ClosestCell, 
00609                                    pcalohh,1,1);
00610       hsim5x5   = spr::eHCALmatrix(theHBHETopology, ClosestCell,
00611                                    pcalohh,2,2);
00612 
00613       hsimInfo3x3 = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, multiplicity3x3);
00614       hsimInfo5x5 = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 2,2, multiplicity5x5);
00615     
00616       // Get energy from all simhits in hcal associated with iso track
00617       trkHcalEne   = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
00618     }
00619 
00620     // Finally for cones of varying radii.
00621     std::vector<double> v_hsimInfoConeMatched;
00622     std::vector<double> v_hsimInfoConeRest   ;
00623     std::vector<double> v_hsimInfoConePhoton  ;
00624     std::vector<double> v_hsimInfoConeNeutHad;
00625     std::vector<double> v_hsimInfoConeCharHad;
00626     std::vector<double> v_hsimInfoConePdgMatched;
00627     std::vector<double> v_hsimInfoConeTotal  ;
00628 
00629     std::vector<int> v_hsimInfoConeNMatched;
00630     std::vector<int> v_hsimInfoConeNTotal  ;
00631     std::vector<int> v_hsimInfoConeNNeutHad;
00632     std::vector<int> v_hsimInfoConeNCharHad;
00633     std::vector<int> v_hsimInfoConeNPhoton ;
00634     std::vector<int> v_hsimInfoConeNRest   ;
00635 
00636     std::vector<double> v_hsimCone           ;
00637     std::vector<double> v_hCone              ;
00638     
00639     std::vector<int>    v_nRecHitsCone       ;
00640     std::vector<int>    v_nSimHitsCone       ;
00641 
00642     std::vector<double> v_distFromHotCell;
00643     std::vector<int> v_ietaFromHotCell;
00644     std::vector<int> v_iphiFromHotCell;
00645     GlobalPoint gposHotCell(0.,0.,0.);
00646     
00647     
00648     std::vector<int>    v_RH_r26_ieta;
00649     std::vector<int>    v_RH_r26_iphi;
00650     std::vector<double> v_RH_r26_ene;
00651     std::vector<int>    v_RH_r44_ieta;
00652     std::vector<int>    v_RH_r44_iphi;
00653     std::vector<double> v_RH_r44_ene;
00654     
00655     
00656         
00657     for (int i = 0 ; i < a_size ; i++){
00658 
00659       
00660       std::map<std::string, double> hsimInfoCone;
00661       double hsimCone = -999.0, hCone = -999.0;
00662       double distFromHotCell = -99.0;
00663       int ietaFromHotCell = -99;
00664       int iphiFromHotCell = -99;
00665       int ietaHotCell = -99;
00666       int iphiHotCell = -99;
00667       int nRecHitsCone = -999;
00668       int nSimHitsCone = -999;
00669 
00670       std::vector<int> multiplicityCone;
00671       std::vector<DetId> coneRecHitDetIds;
00672       if (doMC) 
00673         hsimCone = spr::eCone_hcal(geo, pcalohh, hpoint1, point1, 
00674                                    a_coneR[i], trackMomAtHcal, nSimHitsCone);
00675       
00676       // If needed, get ieta and iphi of rechits for cones of 23.25
00677       // and for hitmap for debugging
00678       bool makeHitmaps = false;
00679       if (a_coneR[i] == 26.23 && makeHitmaps)
00680       {
00681         
00682         hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1, 
00683                                 a_coneR[i], trackMomAtHcal,nRecHitsCone,
00684                                 v_RH_r26_ieta, v_RH_r26_iphi,  v_RH_r26_ene, 
00685                                 coneRecHitDetIds, distFromHotCell, 
00686                                 ietaHotCell, iphiHotCell, gposHotCell);
00687       } 
00688       else if (a_coneR[i] == 43.72 && makeHitmaps)
00689       {
00690         
00691         hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1, 
00692                                 a_coneR[i], trackMomAtHcal,nRecHitsCone,
00693                                 v_RH_r44_ieta, v_RH_r44_iphi,  v_RH_r44_ene, 
00694                                 coneRecHitDetIds, distFromHotCell, 
00695                                 ietaHotCell, iphiHotCell, gposHotCell);
00696       } 
00697       else 
00698       {
00699         
00700         hCone = spr::eCone_hcal(geo, hbhe, hpoint1, point1, 
00701                                 a_coneR[i], trackMomAtHcal, nRecHitsCone, 
00702                                 coneRecHitDetIds, distFromHotCell, 
00703                                 ietaHotCell, iphiHotCell, gposHotCell);
00704       }
00705 
00706       
00707       
00708       if (ietaHotCell != 99){
00709         ietaFromHotCell = ietaHotCell-ClosestCell_HcalDetId.ieta();
00710         iphiFromHotCell = iphiHotCell-ClosestCell_HcalDetId.iphi();
00711       }
00712       
00713       // SimHits NOT matched to RecHits
00714       if (doMC) {
00715         hsimInfoCone = spr::eHCALSimInfoCone(iEvent,pcalohh, SimTk, SimVtx, pTrack, *associate, geo, hpoint1, point1, a_coneR[i], trackMomAtHcal, multiplicityCone);      
00716       
00717       
00718           
00719         // SimHits matched to RecHits
00720         //       hsimInfoCone = spr::eHCALSimInfoCone(iEvent,pcalohh, SimTk, SimVtx, 
00721         //                                                 pTrack, *associate, 
00722         //                                                 geo, hpoint1, point1, 
00723         //                                                 a_coneR[i], trackMomAtHcal, 
00724         //                                                 multiplicityCone,
00725         //                                                 coneRecHitDetIds);      
00726       
00727         v_hsimInfoConeMatched   .push_back(hsimInfoCone["eMatched"   ]);
00728         v_hsimInfoConeRest      .push_back(hsimInfoCone["eRest"      ]);
00729         v_hsimInfoConePhoton    .push_back(hsimInfoCone["eGamma"     ]);
00730         v_hsimInfoConeNeutHad   .push_back(hsimInfoCone["eNeutralHad"]);
00731         v_hsimInfoConeCharHad   .push_back(hsimInfoCone["eChargedHad"]);
00732         v_hsimInfoConePdgMatched.push_back(hsimInfoCone["pdgMatched" ]);
00733         v_hsimInfoConeTotal     .push_back(hsimInfoCone["eTotal"     ]);
00734 
00735         v_hsimInfoConeNMatched  .push_back(multiplicityCone.at(0));
00736 
00737         v_hsimInfoConeNTotal      .push_back(multiplicityCone.at(1));
00738         v_hsimInfoConeNNeutHad    .push_back(multiplicityCone.at(2));
00739         v_hsimInfoConeNCharHad    .push_back(multiplicityCone.at(3));
00740         v_hsimInfoConeNPhoton     .push_back(multiplicityCone.at(4));
00741         v_hsimInfoConeNRest       .push_back(multiplicityCone.at(5));
00742 
00743         v_hsimCone                .push_back(hsimCone                   );
00744         v_nSimHitsCone            .push_back(nSimHitsCone               );
00745       }
00746       v_hCone                   .push_back(hCone                      );
00747       v_nRecHitsCone            .push_back(nRecHitsCone               );
00748       
00749       v_distFromHotCell         .push_back(distFromHotCell            );
00750       v_ietaFromHotCell         .push_back(ietaFromHotCell            );
00751       v_iphiFromHotCell         .push_back(iphiFromHotCell            );
00752       
00753           
00754     }
00755     
00756  
00758     // Fill Vectors that go into root file
00760 
00761     t_v_hnNearTRKs        ->push_back(v_hnNearTRKs                 );
00762     t_v_hnLayers_maxNearP ->push_back(v_hnLayers_maxNearP          );
00763     t_v_htrkQual_maxNearP ->push_back(v_htrkQual_maxNearP          );
00764     t_v_hmaxNearP_goodTrk ->push_back(v_hmaxNearP_goodTrk          ); 
00765     t_v_hmaxNearP         ->push_back(v_hmaxNearP                  );    
00766 
00767     t_v_cone_hnNearTRKs        ->push_back(v_cone_hnNearTRKs       );
00768     t_v_cone_hnLayers_maxNearP ->push_back(v_cone_hnLayers_maxNearP);
00769     t_v_cone_htrkQual_maxNearP ->push_back(v_cone_htrkQual_maxNearP);
00770     t_v_cone_hmaxNearP_goodTrk ->push_back(v_cone_hmaxNearP_goodTrk); 
00771     t_v_cone_hmaxNearP         ->push_back(v_cone_hmaxNearP        );    
00772 
00773     //    t_hScale            ->push_back(hScale                     );
00774     t_trkNOuterHits     ->push_back(nOuterHits                 );
00775     t_trkNLayersCrossed ->push_back(nLayersCrossed             );
00776     t_dtFromLeadJet     ->push_back(drFromLeadJet              );
00777     t_trkP              ->push_back(p1                         );
00778     t_trkPt             ->push_back(pt1                        );
00779     t_trkEta            ->push_back(eta1                       );
00780     t_trkPhi            ->push_back(phi1                       );
00781 
00782     t_e3x3              ->push_back(e3x3                       );
00783     t_v_eDR             ->push_back(v_eDR                      );
00784     t_v_eMipDR          ->push_back(v_eMipDR                   );
00785 
00786     t_h3x3              ->push_back(h3x3                       );
00787     t_h5x5              ->push_back(h5x5                       );
00788     t_nRH_h3x3          ->push_back(nRH_h3x3                   );
00789     t_nRH_h5x5          ->push_back(nRH_h5x5                   );
00790     
00791     t_v_RH_h3x3_ieta    ->push_back(v_RH_h3x3_ieta);
00792     t_v_RH_h3x3_iphi    ->push_back(v_RH_h3x3_iphi);
00793     t_v_RH_h3x3_ene     ->push_back(v_RH_h3x3_ene);
00794     t_v_RH_h5x5_ieta    ->push_back(v_RH_h5x5_ieta);
00795     t_v_RH_h5x5_iphi    ->push_back(v_RH_h5x5_iphi);
00796     t_v_RH_h5x5_ene     ->push_back(v_RH_h5x5_ene);
00797 
00798     if (doMC) {
00799       t_simP              ->push_back(simP                       );
00800       t_hsim3x3           ->push_back(hsim3x3                    );
00801       t_hsim5x5           ->push_back(hsim5x5                    );
00802 
00803       t_hsim3x3Matched    ->push_back(hsimInfo3x3["eMatched"]    );
00804       t_hsim5x5Matched    ->push_back(hsimInfo5x5["eMatched"]    );
00805       t_hsim3x3Rest       ->push_back(hsimInfo3x3["eRest"]       );
00806       t_hsim5x5Rest       ->push_back(hsimInfo5x5["eRest"]       );
00807       t_hsim3x3Photon     ->push_back(hsimInfo3x3["eGamma"]      );
00808       t_hsim5x5Photon     ->push_back(hsimInfo5x5["eGamma"]      );
00809       t_hsim3x3NeutHad    ->push_back(hsimInfo3x3["eNeutralHad"] );
00810       t_hsim5x5NeutHad    ->push_back(hsimInfo5x5["eNeutralHad"] );
00811       t_hsim3x3CharHad    ->push_back(hsimInfo3x3["eChargedHad"] );
00812       t_hsim5x5CharHad    ->push_back(hsimInfo5x5["eChargedHad"] );
00813       t_hsim3x3Total      ->push_back(hsimInfo3x3["eTotal"]      );
00814       t_hsim5x5Total      ->push_back(hsimInfo5x5["eTotal"]      );
00815       t_hsim3x3PdgMatched ->push_back(hsimInfo3x3["pdgMatched"]  );
00816       t_hsim5x5PdgMatched ->push_back(hsimInfo5x5["pdgMatched"]  );
00817 
00818       t_hsim3x3NMatched   ->push_back(multiplicity3x3.at(0));
00819       t_hsim3x3NTotal     ->push_back(multiplicity3x3.at(1));
00820       t_hsim3x3NNeutHad   ->push_back(multiplicity3x3.at(2));
00821       t_hsim3x3NCharHad   ->push_back(multiplicity3x3.at(3));
00822       t_hsim3x3NPhoton    ->push_back(multiplicity3x3.at(4));
00823       t_hsim3x3NRest      ->push_back(multiplicity3x3.at(5));
00824 
00825       t_hsim5x5NMatched   ->push_back(multiplicity5x5.at(0));
00826       t_hsim5x5NTotal     ->push_back(multiplicity5x5.at(1));
00827       t_hsim5x5NNeutHad   ->push_back(multiplicity5x5.at(2));
00828       t_hsim5x5NCharHad   ->push_back(multiplicity5x5.at(3));
00829       t_hsim5x5NPhoton    ->push_back(multiplicity5x5.at(4));
00830       t_hsim5x5NRest      ->push_back(multiplicity5x5.at(5));
00831     }
00832     
00833     t_distFromHotCell_h3x3->push_back(distFromHotCell_h3x3);
00834     t_ietaFromHotCell_h3x3->push_back(ietaFromHotCell_h3x3);
00835     t_iphiFromHotCell_h3x3->push_back(iphiFromHotCell_h3x3);
00836     t_distFromHotCell_h5x5->push_back(distFromHotCell_h5x5);
00837     t_ietaFromHotCell_h5x5->push_back(ietaFromHotCell_h5x5);
00838     t_iphiFromHotCell_h5x5->push_back(iphiFromHotCell_h5x5);
00839     
00840     if (doMC) {
00841       t_trkHcalEne                ->push_back(trkHcalEne                 );
00842       t_trkEcalEne                ->push_back(trkEcalEne                 );
00843 
00844       t_v_hsimInfoConeMatched     ->push_back(v_hsimInfoConeMatched   );
00845       t_v_hsimInfoConeRest        ->push_back(v_hsimInfoConeRest      );
00846       t_v_hsimInfoConePhoton      ->push_back(v_hsimInfoConePhoton   );
00847       t_v_hsimInfoConeNeutHad     ->push_back(v_hsimInfoConeNeutHad   );
00848       t_v_hsimInfoConeCharHad     ->push_back(v_hsimInfoConeCharHad   );
00849       t_v_hsimInfoConePdgMatched  ->push_back(v_hsimInfoConePdgMatched);  
00850       t_v_hsimInfoConeTotal       ->push_back(v_hsimInfoConeTotal     );
00851 
00852       t_v_hsimInfoConeNMatched    ->push_back(v_hsimInfoConeNMatched  );
00853       t_v_hsimInfoConeNTotal      ->push_back(v_hsimInfoConeNTotal    );
00854       t_v_hsimInfoConeNNeutHad    ->push_back(v_hsimInfoConeNNeutHad  );
00855       t_v_hsimInfoConeNCharHad    ->push_back(v_hsimInfoConeNCharHad  );
00856       t_v_hsimInfoConeNPhoton     ->push_back(v_hsimInfoConeNPhoton    );
00857       t_v_hsimInfoConeNRest       ->push_back(v_hsimInfoConeNRest     );  
00858 
00859       t_v_hsimCone    ->push_back(v_hsimCone    );
00860       t_v_hCone       ->push_back(v_hCone       );
00861       t_v_nRecHitsCone->push_back(v_nRecHitsCone);
00862       t_v_nSimHitsCone->push_back(v_nSimHitsCone);
00863     }
00864 
00865     
00866     t_v_distFromHotCell->push_back(v_distFromHotCell);
00867     t_v_ietaFromHotCell->push_back(v_ietaFromHotCell);
00868     t_v_iphiFromHotCell->push_back(v_iphiFromHotCell);
00869     
00870     t_v_RH_r26_ieta    ->push_back(v_RH_r26_ieta);
00871     t_v_RH_r26_iphi    ->push_back(v_RH_r26_iphi);
00872     t_v_RH_r26_ene     ->push_back(v_RH_r26_ene);
00873     t_v_RH_r44_ieta    ->push_back(v_RH_r44_ieta);
00874     t_v_RH_r44_iphi    ->push_back(v_RH_r44_iphi);
00875     t_v_RH_r44_ene     ->push_back(v_RH_r44_ene);
00876 
00877     
00878 
00879     t_v_hlTriggers    ->push_back(v_hlTriggers);
00880     t_hltHB ->push_back(hltHB);
00881     t_hltHE ->push_back(hltHE);
00882     t_hltL1Jet15                 ->push_back(hltL1Jet15                 );
00883     t_hltJet30                   ->push_back(hltJet30                   );
00884     t_hltJet50                   ->push_back(hltJet50                   );
00885     t_hltJet80                   ->push_back(hltJet80                   );
00886     t_hltJet110                  ->push_back(hltJet110                  );
00887     t_hltJet140                  ->push_back(hltJet140                  );
00888     t_hltJet180                  ->push_back(hltJet180                  );
00889     t_hltL1SingleEG5             ->push_back(hltL1SingleEG5             );
00890     t_hltZeroBias                ->push_back(hltZeroBias                );
00891     t_hltMinBiasHcal             ->push_back(hltMinBiasHcal             );
00892     t_hltMinBiasEcal             ->push_back(hltMinBiasEcal             );
00893     t_hltMinBiasPixel            ->push_back(hltMinBiasPixel            );
00894     t_hltSingleIsoTau30_Trk5     ->push_back(hltSingleIsoTau30_Trk5     );
00895     t_hltDoubleLooseIsoTau15_Trk5->push_back(hltDoubleLooseIsoTau15_Trk5);
00896 
00897     t_irun->push_back(irun);
00898     t_ievt->push_back(ievt);
00899     t_ilum->push_back(ilum);
00900     
00901     nTRK++;
00902     
00903     
00904   } // Loop over track collection
00905   
00906     //  std::cout << "nEVT = " << nEVT << std::endl;
00907   
00908   ntp->Fill();
00909   nEVT++;
00910   
00911   
00912   delete associate;
00913 }
00914   
00915 
00916 
00917 void IsolatedTracksCone::beginJob(const edm::EventSetup&) {
00918 
00919   //   hbScale = 120.0;
00920   //   heScale = 135.0;
00921   nEVT=0;
00922   nEVT_failL1=0;
00923   nTRK=0;
00924   
00925   double tempgen_TH[22] = { 0.0,  1.0,  2.0,  3.0,  4.0,  
00926                             5.0,  6.0,  7.0,  8.0,  9.0, 
00927                             10.0, 12.0, 15.0, 20.0, 25.0, 
00928                             30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 100};
00929   for(int i=0; i<22; i++)  genPartPBins[i]  = tempgen_TH[i];
00930 
00931 
00932   double tempgen_Eta[5] = {0.0, 0.5, 1.1, 1.7, 2.0};
00933   for(int i=0; i<5; i++) genPartEtaBins[i] = tempgen_Eta[i];
00934 
00935   t_v_hnNearTRKs           = new std::vector<std::vector<int> >   (); 
00936   t_v_hnLayers_maxNearP    = new std::vector<std::vector<int> >   (); 
00937   t_v_htrkQual_maxNearP    = new std::vector<std::vector<int> >  (); 
00938   t_v_hmaxNearP_goodTrk    = new std::vector<std::vector<double> >();
00939   t_v_hmaxNearP            = new std::vector<std::vector<double> >();
00940                                                             
00941   t_v_cone_hnNearTRKs           = new std::vector<std::vector<int> >   (); 
00942   t_v_cone_hnLayers_maxNearP    = new std::vector<std::vector<int> >   (); 
00943   t_v_cone_htrkQual_maxNearP    = new std::vector<std::vector<int> >  (); 
00944   t_v_cone_hmaxNearP_goodTrk    = new std::vector<std::vector<double> >();
00945   t_v_cone_hmaxNearP            = new std::vector<std::vector<double> >();
00946                                                             
00947   //  t_hScale             = new std::vector<double>();
00948   t_trkNOuterHits      = new std::vector<double>();
00949   t_trkNLayersCrossed  = new std::vector<double>();
00950   t_dtFromLeadJet      = new std::vector<double>();
00951   t_trkP               = new std::vector<double>();
00952   t_trkPt              = new std::vector<double>();
00953   t_trkEta             = new std::vector<double>();
00954   t_trkPhi             = new std::vector<double>();
00955 
00956   t_e3x3               = new std::vector<double>();
00957   t_v_eDR              = new std::vector<std::vector<double> >();
00958   t_v_eMipDR           = new std::vector<std::vector<double> >();
00959 
00960   t_h3x3               = new std::vector<double>();
00961   t_h5x5               = new std::vector<double>();
00962 
00963   t_nRH_h3x3           = new std::vector<double>();
00964   t_nRH_h5x5           = new std::vector<double>();
00965 
00966   if (doMC) {
00967     t_simP               = new std::vector<double>();
00968     t_hsim3x3            = new std::vector<double>();
00969     t_hsim5x5            = new std::vector<double>();
00970 
00971     t_hsim3x3Matched     = new std::vector<double>();
00972     t_hsim5x5Matched     = new std::vector<double>();
00973     t_hsim3x3Rest        = new std::vector<double>();
00974     t_hsim5x5Rest        = new std::vector<double>();
00975     t_hsim3x3Photon      = new std::vector<double>();
00976     t_hsim5x5Photon      = new std::vector<double>();
00977     t_hsim3x3NeutHad     = new std::vector<double>();
00978     t_hsim5x5NeutHad     = new std::vector<double>();
00979     t_hsim3x3CharHad     = new std::vector<double>();
00980     t_hsim5x5CharHad     = new std::vector<double>();
00981     t_hsim3x3PdgMatched  = new std::vector<double>();
00982     t_hsim5x5PdgMatched  = new std::vector<double>();
00983     t_hsim3x3Total       = new std::vector<double>();
00984     t_hsim5x5Total       = new std::vector<double>();
00985 
00986     t_hsim3x3NMatched  = new std::vector<int>();
00987     t_hsim3x3NTotal    = new std::vector<int>();
00988     t_hsim3x3NNeutHad  = new std::vector<int>();
00989     t_hsim3x3NCharHad  = new std::vector<int>();
00990     t_hsim3x3NPhoton   = new std::vector<int>();
00991     t_hsim3x3NRest     = new std::vector<int>();
00992 
00993     t_hsim5x5NMatched  = new std::vector<int>();
00994     t_hsim5x5NTotal    = new std::vector<int>();
00995     t_hsim5x5NNeutHad  = new std::vector<int>();
00996     t_hsim5x5NCharHad  = new std::vector<int>();
00997     t_hsim5x5NPhoton   = new std::vector<int>();
00998     t_hsim5x5NRest     = new std::vector<int>();
00999 
01000     t_trkHcalEne         = new std::vector<double>();
01001     t_trkEcalEne         = new std::vector<double>();
01002   }
01003 
01004   t_distFromHotCell_h3x3  = new std::vector<double>();
01005   t_ietaFromHotCell_h3x3  = new std::vector<int>();
01006   t_iphiFromHotCell_h3x3  = new std::vector<int>();
01007   t_distFromHotCell_h5x5  = new std::vector<double>();
01008   t_ietaFromHotCell_h5x5  = new std::vector<int>();
01009   t_iphiFromHotCell_h5x5  = new std::vector<int>();
01010 
01011   if (doMC) {
01012     t_v_hsimInfoConeMatched   = new std::vector<std::vector<double> >();
01013     t_v_hsimInfoConeRest      = new std::vector<std::vector<double> >();
01014     t_v_hsimInfoConePhoton     = new std::vector<std::vector<double> >();
01015     t_v_hsimInfoConeNeutHad   = new std::vector<std::vector<double> >();
01016     t_v_hsimInfoConeCharHad   = new std::vector<std::vector<double> >();
01017     t_v_hsimInfoConePdgMatched= new std::vector<std::vector<double> >();
01018     t_v_hsimInfoConeTotal     = new std::vector<std::vector<double> >();
01019 
01020     t_v_hsimInfoConeNMatched  = new std::vector<std::vector<int> >();
01021     t_v_hsimInfoConeNTotal    = new std::vector<std::vector<int> >();
01022     t_v_hsimInfoConeNNeutHad  = new std::vector<std::vector<int> >();
01023     t_v_hsimInfoConeNCharHad  = new std::vector<std::vector<int> >();
01024     t_v_hsimInfoConeNPhoton    = new std::vector<std::vector<int> >();
01025     t_v_hsimInfoConeNRest     = new std::vector<std::vector<int> >();
01026 
01027     t_v_hsimCone    = new std::vector<std::vector<double> >();
01028   }
01029 
01030   t_v_hCone       = new std::vector<std::vector<double> >();
01031   t_v_nRecHitsCone= new std::vector<std::vector<int> >();
01032   t_v_nSimHitsCone= new std::vector<std::vector<int> >();
01033 
01034   t_v_distFromHotCell= new std::vector<std::vector<double> >();
01035   t_v_ietaFromHotCell= new std::vector<std::vector<int> >();
01036   t_v_iphiFromHotCell= new std::vector<std::vector<int> >();
01037 
01038   t_v_RH_h3x3_ieta = new std::vector<std::vector<int> >();
01039   t_v_RH_h3x3_iphi = new std::vector<std::vector<int> >();
01040   t_v_RH_h3x3_ene  = new std::vector<std::vector<double> >();
01041   t_v_RH_h5x5_ieta = new std::vector<std::vector<int> >();
01042   t_v_RH_h5x5_iphi = new std::vector<std::vector<int> >();
01043   t_v_RH_h5x5_ene  = new std::vector<std::vector<double> >();
01044   t_v_RH_r26_ieta  = new std::vector<std::vector<int> >();
01045   t_v_RH_r26_iphi  = new std::vector<std::vector<int> >();
01046   t_v_RH_r26_ene   = new std::vector<std::vector<double> >();
01047   t_v_RH_r44_ieta  = new std::vector<std::vector<int> >();
01048   t_v_RH_r44_iphi  = new std::vector<std::vector<int> >();
01049   t_v_RH_r44_ene   = new std::vector<std::vector<double> >();
01050 
01051 
01052   t_v_hlTriggers    = new std::vector<std::vector<int> >();
01053 
01054   t_hltHE                         = new std::vector<int>();
01055   t_hltHB                         = new std::vector<int>();
01056   t_hltL1Jet15                    = new std::vector<int>();
01057   t_hltJet30                      = new std::vector<int>();
01058   t_hltJet50                      = new std::vector<int>();
01059   t_hltJet80                      = new std::vector<int>();
01060   t_hltJet110                     = new std::vector<int>();
01061   t_hltJet140                     = new std::vector<int>();
01062   t_hltJet180                     = new std::vector<int>();
01063   t_hltL1SingleEG5                = new std::vector<int>();
01064   t_hltZeroBias                   = new std::vector<int>();
01065   t_hltMinBiasHcal                = new std::vector<int>();
01066   t_hltMinBiasEcal                = new std::vector<int>();
01067   t_hltMinBiasPixel               = new std::vector<int>();
01068   t_hltSingleIsoTau30_Trk5        = new std::vector<int>();
01069   t_hltDoubleLooseIsoTau15_Trk5   = new std::vector<int>();
01070   
01071 
01072   t_irun = new std::vector<unsigned int>();
01073   t_ievt = new std::vector<unsigned int>();
01074   t_ilum = new std::vector<unsigned int>();
01075 
01076   BuildTree();
01077 }
01078 
01079 void IsolatedTracksCone::clearTrackVectors() {
01080 
01081   t_v_hnNearTRKs          ->clear();   
01082   t_v_hnLayers_maxNearP   ->clear();   
01083   t_v_htrkQual_maxNearP   ->clear();   
01084   t_v_hmaxNearP_goodTrk   ->clear();
01085   t_v_hmaxNearP           ->clear();
01086 
01087   t_v_cone_hnNearTRKs          ->clear();   
01088   t_v_cone_hnLayers_maxNearP   ->clear();   
01089   t_v_cone_htrkQual_maxNearP   ->clear();   
01090   t_v_cone_hmaxNearP_goodTrk   ->clear();
01091   t_v_cone_hmaxNearP           ->clear();
01092 
01093   //  t_hScale             ->clear();
01094   t_trkNOuterHits      ->clear();
01095   t_trkNLayersCrossed  ->clear();
01096   t_dtFromLeadJet      ->clear();
01097   t_trkP               ->clear();
01098   t_trkPt              ->clear();
01099   t_trkEta             ->clear();
01100   t_trkPhi             ->clear();
01101   t_e3x3               ->clear();
01102   t_v_eDR              ->clear();
01103   t_v_eMipDR           ->clear();
01104   t_h3x3               ->clear();
01105   t_h5x5               ->clear();
01106   t_nRH_h3x3           ->clear();
01107   t_nRH_h5x5           ->clear();
01108 
01109   if (doMC) {
01110     t_simP               ->clear();
01111     t_hsim3x3            ->clear();
01112     t_hsim5x5            ->clear();
01113     t_hsim3x3Matched     ->clear();
01114     t_hsim5x5Matched     ->clear();
01115     t_hsim3x3Rest        ->clear();
01116     t_hsim5x5Rest        ->clear();
01117     t_hsim3x3Photon      ->clear();
01118     t_hsim5x5Photon      ->clear();
01119     t_hsim3x3NeutHad     ->clear();
01120     t_hsim5x5NeutHad     ->clear();
01121     t_hsim3x3CharHad     ->clear();
01122     t_hsim5x5CharHad     ->clear();
01123     t_hsim3x3PdgMatched  ->clear();
01124     t_hsim5x5PdgMatched  ->clear();
01125     t_hsim3x3Total       ->clear();
01126     t_hsim5x5Total       ->clear();
01127 
01128     t_hsim3x3NMatched    ->clear();
01129     t_hsim3x3NTotal      ->clear();
01130     t_hsim3x3NNeutHad    ->clear();
01131     t_hsim3x3NCharHad    ->clear();
01132     t_hsim3x3NPhoton      ->clear();
01133     t_hsim3x3NRest       ->clear();
01134 
01135     t_hsim5x5NMatched    ->clear();
01136     t_hsim5x5NTotal      ->clear();
01137     t_hsim5x5NNeutHad    ->clear();
01138     t_hsim5x5NCharHad    ->clear();
01139     t_hsim5x5NPhoton      ->clear();
01140     t_hsim5x5NRest       ->clear();
01141 
01142     t_trkHcalEne         ->clear();
01143     t_trkEcalEne         ->clear();
01144   }
01145 
01146   t_distFromHotCell_h3x3  ->clear();
01147   t_ietaFromHotCell_h3x3  ->clear();
01148   t_iphiFromHotCell_h3x3  ->clear();
01149   t_distFromHotCell_h5x5  ->clear();
01150   t_ietaFromHotCell_h5x5  ->clear();
01151   t_iphiFromHotCell_h5x5  ->clear();
01152 
01153   if (doMC) {
01154     t_v_hsimInfoConeMatched   ->clear();
01155     t_v_hsimInfoConeRest      ->clear();
01156     t_v_hsimInfoConePhoton     ->clear();
01157     t_v_hsimInfoConeNeutHad   ->clear();
01158     t_v_hsimInfoConeCharHad   ->clear();
01159     t_v_hsimInfoConePdgMatched->clear();  
01160     t_v_hsimInfoConeTotal     ->clear();
01161 
01162     t_v_hsimInfoConeNMatched   ->clear();
01163     t_v_hsimInfoConeNRest      ->clear();
01164     t_v_hsimInfoConeNPhoton     ->clear();
01165     t_v_hsimInfoConeNNeutHad   ->clear();
01166     t_v_hsimInfoConeNCharHad   ->clear();
01167     t_v_hsimInfoConeNTotal     ->clear();
01168 
01169     t_v_hsimCone    ->clear();
01170   }
01171 
01172   t_v_hCone       ->clear();
01173   t_v_nRecHitsCone->clear();
01174   t_v_nSimHitsCone->clear();
01175 
01176   t_v_distFromHotCell->clear();
01177   t_v_ietaFromHotCell->clear();
01178   t_v_iphiFromHotCell->clear();
01179   
01180   t_v_RH_h3x3_ieta ->clear();
01181   t_v_RH_h3x3_iphi ->clear();
01182   t_v_RH_h3x3_ene ->clear();
01183   t_v_RH_h5x5_ieta ->clear();
01184   t_v_RH_h5x5_iphi ->clear();
01185   t_v_RH_h5x5_ene ->clear();
01186   t_v_RH_r26_ieta  ->clear();
01187   t_v_RH_r26_iphi  ->clear();
01188   t_v_RH_r26_ene  ->clear();
01189   t_v_RH_r44_ieta  ->clear();
01190   t_v_RH_r44_iphi  ->clear();
01191   t_v_RH_r44_ene  ->clear();
01192 
01193   t_v_hlTriggers  ->clear();
01194   t_hltHB->clear();
01195   t_hltHE->clear();
01196   t_hltL1Jet15               ->clear();
01197   t_hltJet30                 ->clear();
01198   t_hltJet50                 ->clear();
01199   t_hltJet80                 ->clear();
01200   t_hltJet110                ->clear();
01201   t_hltJet140                ->clear();
01202   t_hltJet180                ->clear();
01203   t_hltL1SingleEG5           ->clear();
01204   t_hltZeroBias              ->clear();
01205   t_hltMinBiasHcal           ->clear();
01206   t_hltMinBiasEcal           ->clear();
01207   t_hltMinBiasPixel          ->clear();
01208   t_hltSingleIsoTau30_Trk5     ->clear();
01209   t_hltDoubleLooseIsoTau15_Trk5->clear();
01210   
01211   t_irun->clear();
01212   t_ievt->clear();
01213   t_ilum->clear();
01214 
01215 
01216 
01217 }
01218 
01219 // ----- method called once each job just after ending the event loop ----
01220 void IsolatedTracksCone::endJob() {
01221 
01222   std::cout << "Number of Events Processed " << nEVT << std::endl;
01223   
01224 }
01225 
01226 
01227 void IsolatedTracksCone::BuildTree(){
01228 
01229 
01230   hRawPt  = fs->make<TH1F>("hRawPt ", "hRawPt ",  100,  0.0, 100.0);
01231   hRawP   = fs->make<TH1F>("hRawP  ", "hRawP  ",  100,  0.0, 100.0);
01232   hRawEta = fs->make<TH1F>("hRawEta", "hRawEta",   15,  0.0,   3.0);
01233   hRawPhi = fs->make<TH1F>("hRawPhi", "hRawPhi",  100, -3.2,   3.2);
01234 
01235   ntp = fs->make<TTree>("ntp", "ntp");
01236 
01237   
01238   // Counters
01239   ntp->Branch("nEVT"        , &nEVT        , "nEVT/I"        );
01240   ntp->Branch("leadL1JetPT" , &leadL1JetPT , "leadL1JetPT/D" );
01241   ntp->Branch("leadL1JetEta", &leadL1JetEta, "leadL1JetEta/D");
01242   ntp->Branch("leadL1JetPhi", &leadL1JetPhi, "leadL1JetPhi/D");
01243   ntp->Branch("nTRK",         &nTRK,         "nTRK/I");
01244   ntp->Branch("nRawTRK"            , &nRawTRK            ,"nRawTRK/I"            );
01245   ntp->Branch("nFailHighPurityQaul", &nFailHighPurityQaul,"nFailHighPurityQaul/I");
01246   ntp->Branch("nFailPt"            , &nFailPt            ,"nFailPt/I"            );
01247   ntp->Branch("nFailEta"           , &nFailEta           ,"nFailEta/I"           );
01248   ntp->Branch("nMissEcal"          , &nMissEcal          ,"nMissEcal/I"          );
01249   ntp->Branch("nMissHcal"          , &nMissHcal          ,"nMissHcal/I"          );
01250 
01251   ntp->Branch("hnNearTRKs"          ,"vector<vector<int> >   ",&t_v_hnNearTRKs          );
01252   ntp->Branch("hnLayers_maxNearP"   ,"vector<vector<int> >   ",&t_v_hnLayers_maxNearP   );
01253   ntp->Branch("htrkQual_maxNearP"   ,"vector<vector<int> >   ",&t_v_htrkQual_maxNearP   );
01254   ntp->Branch("hmaxNearP_goodTrk"   ,"vector<vector<double> >",&t_v_hmaxNearP_goodTrk   );
01255   ntp->Branch("hmaxNearP"           ,"vector<vector<double> >",&t_v_hmaxNearP           );
01256 
01257   ntp->Branch("cone_hnNearTRKs"       ,"vector<vector<int> >   ",&t_v_cone_hnNearTRKs       );
01258   ntp->Branch("cone_hnLayers_maxNearP","vector<vector<int> >   ",&t_v_cone_hnLayers_maxNearP);
01259   ntp->Branch("cone_htrkQual_maxNearP","vector<vector<int> >   ",&t_v_cone_htrkQual_maxNearP);
01260   ntp->Branch("cone_hmaxNearP_goodTrk","vector<vector<double> >",&t_v_cone_hmaxNearP_goodTrk);
01261   ntp->Branch("cone_hmaxNearP"        ,"vector<vector<double> >",&t_v_cone_hmaxNearP        );
01262                                                                     
01263   //  ntp->Branch("hScale"           , "vector<double>", &t_hScale           );
01264   ntp->Branch("trkNOuterHits"    , "vector<double>", &t_trkNOuterHits    );
01265   ntp->Branch("trkNLayersCrossed", "vector<double>", &t_trkNLayersCrossed);
01266   ntp->Branch("drFromLeadJet"    , "vector<double>", &t_dtFromLeadJet    );
01267   ntp->Branch("trkP"             , "vector<double>", &t_trkP             );
01268   ntp->Branch("trkPt"            , "vector<double>", &t_trkPt            );
01269   ntp->Branch("trkEta"           , "vector<double>", &t_trkEta           );
01270   ntp->Branch("trkPhi"           , "vector<double>", &t_trkPhi           );
01271   ntp->Branch("e3x3"             , "vector<double>", &t_e3x3             );
01272 
01273   ntp->Branch("e3x3"          , "vector<double>"         , &t_e3x3 );
01274   ntp->Branch("v_eDR"         , "vector<vector<double> >", &t_v_eDR);
01275   ntp->Branch("v_eMipDR"      , "vector<vector<double> >", &t_v_eMipDR);
01276 
01277   ntp->Branch("h3x3"             , "vector<double>", &t_h3x3             );
01278   ntp->Branch("h5x5"             , "vector<double>", &t_h5x5             );
01279   ntp->Branch("nRH_h3x3"         , "vector<double>", &t_nRH_h3x3         );
01280   ntp->Branch("nRH_h5x5"         , "vector<double>", &t_nRH_h5x5         );
01281 
01282   if (doMC) {
01283     ntp->Branch("simP"             , "vector<double>", &t_simP             );
01284     ntp->Branch("hsim3x3"          , "vector<double>", &t_hsim3x3          );
01285     ntp->Branch("hsim5x5"          , "vector<double>", &t_hsim5x5          );
01286 
01287     ntp->Branch("hsim3x3Matched"   , "vector<double>", &t_hsim3x3Matched   );
01288     ntp->Branch("hsim5x5Matched"   , "vector<double>", &t_hsim5x5Matched   );
01289     ntp->Branch("hsim3x3Rest"      , "vector<double>", &t_hsim3x3Rest      );
01290     ntp->Branch("hsim5x5Rest"      , "vector<double>", &t_hsim5x5Rest      );
01291     ntp->Branch("hsim3x3Photon"    , "vector<double>", &t_hsim3x3Photon    );
01292     ntp->Branch("hsim5x5Photon"    , "vector<double>", &t_hsim5x5Photon    );
01293     ntp->Branch("hsim3x3NeutHad"   , "vector<double>", &t_hsim3x3NeutHad   );
01294     ntp->Branch("hsim5x5NeutHad"   , "vector<double>", &t_hsim5x5NeutHad   );
01295     ntp->Branch("hsim3x3CharHad"   , "vector<double>", &t_hsim3x3CharHad   );
01296     ntp->Branch("hsim5x5CharHad"   , "vector<double>", &t_hsim5x5CharHad   );
01297     ntp->Branch("hsim3x3PdgMatched", "vector<double>", &t_hsim3x3PdgMatched);
01298     ntp->Branch("hsim5x5PdgMatched", "vector<double>", &t_hsim5x5PdgMatched);
01299     ntp->Branch("hsim3x3Total"     , "vector<double>", &t_hsim3x3Total     );
01300     ntp->Branch("hsim5x5Total"     , "vector<double>", &t_hsim5x5Total     );
01301 
01302     ntp->Branch("hsim3x3NMatched"   , "vector<int>", &t_hsim3x3NMatched   );
01303     ntp->Branch("hsim3x3NRest"      , "vector<int>", &t_hsim3x3NRest      );
01304     ntp->Branch("hsim3x3NPhoton"    , "vector<int>", &t_hsim3x3NPhoton    );
01305     ntp->Branch("hsim3x3NNeutHad"   , "vector<int>", &t_hsim3x3NNeutHad   );
01306     ntp->Branch("hsim3x3NCharHad"   , "vector<int>", &t_hsim3x3NCharHad   );
01307     ntp->Branch("hsim3x3NTotal"     , "vector<int>", &t_hsim3x3NTotal     );
01308 
01309     ntp->Branch("hsim5x5NMatched"   , "vector<int>", &t_hsim5x5NMatched   );
01310     ntp->Branch("hsim5x5NRest"      , "vector<int>", &t_hsim5x5NRest      );
01311     ntp->Branch("hsim5x5NPhoton"    , "vector<int>", &t_hsim5x5NPhoton    );
01312     ntp->Branch("hsim5x5NNeutHad"   , "vector<int>", &t_hsim5x5NNeutHad   );
01313     ntp->Branch("hsim5x5NCharHad"   , "vector<int>", &t_hsim5x5NCharHad   );
01314     ntp->Branch("hsim5x5NTotal"     , "vector<int>", &t_hsim5x5NTotal     );
01315 
01316     ntp->Branch("trkHcalEne"       , "vector<double>", &t_trkHcalEne       );
01317     ntp->Branch("trkEcalEne"       , "vector<double>", &t_trkEcalEne       );
01318   }
01319 
01320   ntp->Branch("distFromHotCell_h3x3", "vector<double>", &t_distFromHotCell_h3x3);
01321   ntp->Branch("ietaFromHotCell_h3x3", "vector<int>", &t_ietaFromHotCell_h3x3);
01322   ntp->Branch("iphiFromHotCell_h3x3", "vector<int>", &t_iphiFromHotCell_h3x3);
01323   ntp->Branch("distFromHotCell_h5x5", "vector<double>", &t_distFromHotCell_h5x5);
01324   ntp->Branch("ietaFromHotCell_h5x5", "vector<int>", &t_ietaFromHotCell_h5x5);
01325   ntp->Branch("iphiFromHotCell_h5x5", "vector<int>", &t_iphiFromHotCell_h5x5);
01326 
01327   if (doMC) {
01328     ntp->Branch("v_hsimInfoConeMatched"   ,"vector<vector<double> >",&t_v_hsimInfoConeMatched   );
01329     ntp->Branch("v_hsimInfoConeRest"      ,"vector<vector<double> >",&t_v_hsimInfoConeRest      );
01330     ntp->Branch("v_hsimInfoConePhoton"     ,"vector<vector<double> >",&t_v_hsimInfoConePhoton     );
01331     ntp->Branch("v_hsimInfoConeNeutHad"   ,"vector<vector<double> >",&t_v_hsimInfoConeNeutHad   );
01332     ntp->Branch("v_hsimInfoConeCharHad"   ,"vector<vector<double> >",&t_v_hsimInfoConeCharHad   );
01333     ntp->Branch("v_hsimInfoConePdgMatched","vector<vector<double> >",&t_v_hsimInfoConePdgMatched);
01334     ntp->Branch("v_hsimInfoConeTotal"     ,"vector<vector<double> >",&t_v_hsimInfoConeTotal     );
01335 
01336     ntp->Branch("v_hsimInfoConeNMatched"  ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNMatched   );
01337     ntp->Branch("v_hsimInfoConeNRest"     ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNRest      );
01338     ntp->Branch("v_hsimInfoConeNPhoton"    ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNPhoton     );
01339     ntp->Branch("v_hsimInfoConeNNeutHad"  ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNNeutHad   );
01340     ntp->Branch("v_hsimInfoConeNCharHad"  ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNCharHad   );
01341     ntp->Branch("v_hsimInfoConeNTotal"    ,"vector<vector<int> >"   ,&t_v_hsimInfoConeNTotal     );
01342 
01343     ntp->Branch("v_hsimCone"     ,"vector<vector<double> >",&t_v_hsimCone    );
01344   }
01345 
01346   ntp->Branch("v_hCone"        ,"vector<vector<double> >",&t_v_hCone       );
01347   ntp->Branch("v_nRecHitsCone" ,"vector<vector<int> >"   ,&t_v_nRecHitsCone);
01348   ntp->Branch("v_nSimHitsCone" ,"vector<vector<int> >"   ,&t_v_nSimHitsCone);
01349 
01350   ntp->Branch("v_distFromHotCell"        ,"vector<vector<double> >",&t_v_distFromHotCell       );
01351   ntp->Branch("v_ietaFromHotCell"        ,"vector<vector<int> >",&t_v_ietaFromHotCell       );
01352   ntp->Branch("v_iphiFromHotCell"        ,"vector<vector<int> >",&t_v_iphiFromHotCell       );
01353 
01354   ntp->Branch("v_RH_h3x3_ieta" ,"vector<vector<int> >"   ,&t_v_RH_h3x3_ieta);
01355   ntp->Branch("v_RH_h3x3_iphi" ,"vector<vector<int> >"   ,&t_v_RH_h3x3_iphi);
01356   ntp->Branch("v_RH_h3x3_ene"  ,"vector<vector<double> >",&t_v_RH_h3x3_ene );
01357   ntp->Branch("v_RH_h5x5_ieta" ,"vector<vector<int> >"   ,&t_v_RH_h5x5_ieta);
01358   ntp->Branch("v_RH_h5x5_iphi" ,"vector<vector<int> >"   ,&t_v_RH_h5x5_iphi);
01359   ntp->Branch("v_RH_h5x5_ene"  ,"vector<vector<double> >",&t_v_RH_h5x5_ene );
01360   ntp->Branch("v_RH_r26_ieta"  ,"vector<vector<int> >"   ,&t_v_RH_r26_ieta );
01361   ntp->Branch("v_RH_r26_iphi"  ,"vector<vector<int> >"   ,&t_v_RH_r26_iphi );
01362   ntp->Branch("v_RH_r26_ene"   ,"vector<vector<double> >",&t_v_RH_r26_ene  );
01363   ntp->Branch("v_RH_r44_ieta"  ,"vector<vector<int> >"   ,&t_v_RH_r44_ieta );
01364   ntp->Branch("v_RH_r44_iphi"  ,"vector<vector<int> >"   ,&t_v_RH_r44_iphi );
01365   ntp->Branch("v_RH_r44_ene"   ,"vector<vector<double> >",&t_v_RH_r44_ene  );
01366 
01367   ntp->Branch("v_hlTriggers" ,"vector<vector<int> >",&t_v_hlTriggers);
01368   ntp->Branch("v_hltHB" ,"vector<int>",&t_hltHB);
01369   ntp->Branch("v_hltHE" ,"vector<int>",&t_hltHE);
01370   ntp->Branch("v_hltL1Jet15" ,"vector<int>",&t_hltL1Jet15                    );
01371   ntp->Branch("v_hltJet30" ,"vector<int>",&t_hltJet30                );
01372   ntp->Branch("v_hltJet50" ,"vector<int>",&t_hltJet50                );
01373   ntp->Branch("v_hltJet80" ,"vector<int>",&t_hltJet80                );
01374   ntp->Branch("v_hltJet110" ,"vector<int>",&t_hltJet110              );
01375   ntp->Branch("v_hltJet140" ,"vector<int>",&t_hltJet140              );
01376   ntp->Branch("v_hltJet180" ,"vector<int>",&t_hltJet180              );
01377   ntp->Branch("v_hltL1SingleEG5" ,"vector<int>",&t_hltL1SingleEG5            );
01378   ntp->Branch("v_hltZeroBias" ,"vector<int>",&t_hltZeroBias                  );
01379   ntp->Branch("v_hltMinBiasHcal" ,"vector<int>",&t_hltMinBiasHcal            );
01380   ntp->Branch("v_hltMinBiasEcal" ,"vector<int>",&t_hltMinBiasEcal            );
01381   ntp->Branch("v_hltMinBiasPixel" ,"vector<int>",&t_hltMinBiasPixel          );
01382   ntp->Branch("v_hltSingleIsoTau30_Trk5" ,"vector<int>",&t_hltSingleIsoTau30_Trk5     );
01383   ntp->Branch("v_hltDoubleLooseIsoTau15_Trk5" ,"vector<int>",&t_hltDoubleLooseIsoTau15_Trk5);
01384 
01385   ntp->Branch("irun" ,"vector<unsigned int>", &t_irun);
01386   ntp->Branch("ievt" ,"vector<unsigned int>", &t_ievt);
01387   ntp->Branch("ilum" ,"vector<unsigned int>", &t_ilum);
01388    
01389 }
01390 
01391 
01392 void IsolatedTracksCone::printTrack(const reco::Track* pTrack) {
01393   
01394   std::string theTrackQuality = "highPurity";
01395   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
01396 
01397   std::cout << " Reference Point " << pTrack->referencePoint() <<"\n"
01398             << " TrackMmentum " << pTrack->momentum()
01399             << " (pt,eta,phi)(" << pTrack->pt()<<","<<pTrack->eta()<<","<<pTrack->phi()<<")"
01400             << " p " << pTrack->p() << "\n"
01401             << " Normalized chi2 " << pTrack->normalizedChi2() <<"  charge " << pTrack->charge()
01402             << " qoverp() " << pTrack->qoverp() <<"+-" << pTrack->qoverpError()
01403             << " d0 " << pTrack->d0() << "\n"
01404             << " NValidHits " << pTrack->numberOfValidHits() << "  NLostHits " << pTrack->numberOfLostHits()
01405             << " TrackQuality " << pTrack->qualityName(trackQuality_) << " " << pTrack->quality(trackQuality_) 
01406             << std::endl;
01407   
01408   if( printTrkHitPattern_ ) {
01409     const reco::HitPattern& p = pTrack->hitPattern();
01410     
01411     for (int i=0; i<p.numberOfHits(); i++) {
01412       p.printHitPattern(i, std::cout);
01413     }
01414   }
01415 
01416 }
01417 
01418 double IsolatedTracksCone::DeltaPhi(double v1, double v2) {
01419   // Computes the correctly normalized phi difference
01420   // v1, v2 = phi of object 1 and 2
01421   
01422   double pi    = 3.141592654;
01423   double twopi = 6.283185307;
01424   
01425   double diff = std::abs(v2 - v1);
01426   double corr = twopi - diff;
01427   if (diff < pi){ return diff;} else { return corr;} 
01428 }
01429 
01430 
01431 double IsolatedTracksCone::DeltaR(double eta1, double phi1, 
01432                                   double eta2, double phi2) {
01433   double deta = eta1 - eta2;
01434   double dphi = DeltaPhi(phi1, phi2);
01435   return std::sqrt(deta*deta + dphi*dphi);
01436 }
01437 
01438 
01439 //define this as a plug-in
01440 DEFINE_FWK_MODULE(IsolatedTracksCone);
01441           
01442           
01443           
01444           
01445           
01446 
01447           
01448           
01449 
01450 
01451 
01452 
01453 
01454 
01455
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Calibration/IsolatedParticles/plugins/IsolatedTracksCone.cc
