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