Reference Manual

Go to the documentation of this file.
00001 #include "DQMOffline/JetMET/interface/JPTJetAnalyzer.h"
00002 #include "DQMServices/Core/interface/MonitorElement.h"
00003 #include "DQMServices/Core/interface/DQMStore.h"
00004 #include "DataFormats/TrackReco/interface/Track.h"
00005 #include "DataFormats/JetReco/interface/CaloJet.h"
00006 #include "FWCore/Framework/interface/Event.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "DataFormats/Common/interface/Handle.h"
00011 #include "FWCore/Framework/interface/ESHandle.h"
00012 #include "DataFormats/Math/interface/deltaR.h"
00013 #include "RecoJets/JetAssociationAlgorithms/interface/JetTracksAssociationDRCalo.h"
00014 #include "DataFormats/Math/interface/Point3D.h"
00015 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00016 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00017 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00018 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00019 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00020 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00021 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
00022 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
00023 #include "DataFormats/DetId/interface/DetId.h"
00024 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00025 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
00026 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
00027 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
00028 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
00029 #include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
00030 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
00031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00032 #include "RecoJets/JetProducers/interface/JetIDHelper.h"
00033 #include <cmath>
00034 #include <string>
00035 #include <memory>
00036 #include <vector>
00037 
00038 namespace jptJetAnalysis {
00039   
00040   // Helpper class to propagate tracks to the calo surface using the same implementation as the JetTrackAssociator
00041   class TrackPropagatorToCalo
00042   {
00043    public:
00044     TrackPropagatorToCalo();
00045     void update(const edm::EventSetup& eventSetup);
00046     math::XYZPoint impactPoint(const reco::Track& track) const;
00047    private:
00048     const MagneticField* magneticField_;
00049     const Propagator* propagator_;
00050     uint32_t magneticFieldCacheId_;
00051     uint32_t propagatorCacheId_;
00052   };
00053   
00054   // Helper class to calculate strip signal to noise and manage the necessary ES objects
00055   class StripSignalOverNoiseCalculator
00056   {
00057    public:
00058     StripSignalOverNoiseCalculator(const std::string& theQualityLabel = std::string(""));
00059     void update(const edm::EventSetup& eventSetup);
00060     double signalOverNoise(const SiStripCluster& cluster,
00061                            const uint32_t& id) const;
00062     double operator () (const SiStripCluster& cluster,
00063                         const uint32_t& id) const
00064     { return signalOverNoise(cluster,id); }
00065    private:
00066     const std::string qualityLabel_;
00067     const SiStripQuality* quality_;
00068     const SiStripNoises* noise_;
00069     const SiStripGain* gain_;
00070     uint32_t qualityCacheId_;
00071     uint32_t noiseCacheId_;
00072     uint32_t gainCacheId_;
00073   };
00074   
00075 }
00076 
00077 const char* JPTJetAnalyzer::messageLoggerCatregory = "JetPlusTrackDQM";
00078 
00079 JPTJetAnalyzer::JPTJetAnalyzer(const edm::ParameterSet& config)
00080   : histogramPath_(config.getParameter<std::string>("HistogramPath")),
00081     verbose_(config.getUntrackedParameter<bool>("PrintDebugMessages",false)),
00082     writeDQMStore_(config.getUntrackedParameter<bool>("WriteDQMStore")),
00083     dqmStoreFileName_(),
00084     n90HitsMin_(config.getParameter<int>("n90HitsMin")),
00085     fHPDMax_(config.getParameter<double>("fHPDMax")),
00086     resEMFMin_(config.getParameter<double>("resEMFMin")),
00087     correctedPtMin_(config.getParameter<double>("correctedPtThreshold")),
00088     trackPropagator_(new jptJetAnalysis::TrackPropagatorToCalo),
00089     sOverNCalculator_(new jptJetAnalysis::StripSignalOverNoiseCalculator),
00090     jetID_(new reco::helper::JetIDHelper(config.getParameter<edm::ParameterSet>("JetIDParams")))
00091 {
00092   //print config to debug log
00093   std::ostringstream debugStream;
00094   if (verbose_) {
00095     debugStream << "Configuration for JPTJetAnalyzer: " << std::endl
00096                 << "\tHistogramPath: " << histogramPath_ << std::endl
00097                 << "\tPrintDebugMessages? " << (verbose_ ? "yes" : "no") << std::endl;
00098   }
00099   if (writeDQMStore_) {
00100     dqmStoreFileName_ = config.getUntrackedParameter<std::string>("DQMStoreFileName");
00101   }
00102   
00103   //don't generate debug mesages if debug is disabled
00104   std::ostringstream* pDebugStream = (verbose_ ? &debugStream : NULL);
00105   
00106   //get histogram configuration
00107   getConfigForHistogram("E",config,pDebugStream);
00108   getConfigForHistogram("Et",config,pDebugStream);
00109   getConfigForHistogram("P",config,pDebugStream);
00110   /*  getConfigForHistogram("Mass",config,pDebugStream);  */
00111   getConfigForHistogram("Pt",config,pDebugStream);
00112   getConfigForHistogram("Pt1",config,pDebugStream);
00113   getConfigForHistogram("Pt2",config,pDebugStream);
00114   getConfigForHistogram("Pt3",config,pDebugStream);
00115   getConfigForHistogram("Px",config,pDebugStream);
00116   getConfigForHistogram("Py",config,pDebugStream);
00117   getConfigForHistogram("Pz",config,pDebugStream);
00118   getConfigForHistogram("Eta",config,pDebugStream);
00119   getConfigForHistogram("Phi",config,pDebugStream);
00120   /*  getConfigForHistogram("deltaEta",config,pDebugStream);
00121       getConfigForHistogram("deltaPhi",config,pDebugStream);  */
00122   getConfigForHistogram("PhiVsEta",config,pDebugStream);
00123   /*  getConfigForHistogram("N90Hits",config,pDebugStream);  
00124   getConfigForHistogram("fHPD",config,pDebugStream);
00125   getConfigForHistogram("ResEMF",config,pDebugStream);   
00126   getConfigForHistogram("fRBX",config,pDebugStream); 
00127   getConfigForHistogram("TrackSiStripHitStoN",config,pDebugStream); */
00128   getConfigForHistogram("InCaloTrackDirectionJetDR",config,pDebugStream);
00129   getConfigForHistogram("OutCaloTrackDirectionJetDR",config,pDebugStream);
00130   getConfigForHistogram("InVertexTrackImpactPointJetDR",config,pDebugStream);
00131   getConfigForHistogram("OutVertexTrackImpactPointJetDR",config,pDebugStream);
00132   /*  getConfigForHistogram("PtFractionInCone",config,pDebugStream);
00133   getConfigForHistogram("PtFractionInConeVsJetRawEt",config,pDebugStream);
00134   getConfigForHistogram("PtFractionInConeVsJetEta",config,pDebugStream);  */
00135   getConfigForHistogram("nTracks",config,pDebugStream);
00136   getConfigForHistogram("nTracksVsJetEt",config,pDebugStream);
00137   getConfigForHistogram("nTracksVsJetEta",config,pDebugStream);
00138   /*  getConfigForHistogram("CorrFactor",config,pDebugStream);
00139   getConfigForHistogram("CorrFactorVsJetEt",config,pDebugStream);
00140   getConfigForHistogram("CorrFactorVsJetEta",config,pDebugStream);  
00141   getConfigForHistogram("ZSPCorrFactor",config,pDebugStream);
00142   getConfigForHistogram("ZSPCorrFactorVsJetEt",config,pDebugStream);
00143   getConfigForHistogram("ZSPCorrFactorVsJetEta",config,pDebugStream);
00144   getConfigForHistogram("JPTCorrFactor",config,pDebugStream);
00145   getConfigForHistogram("JPTCorrFactorVsJetEt",config,pDebugStream);
00146   getConfigForHistogram("JPTCorrFactorVsJetEta",config,pDebugStream);  */
00147   getConfigForTrackHistograms("AllPions",config,pDebugStream);
00148   getConfigForTrackHistograms("InCaloInVertexPions",config,pDebugStream);
00149   getConfigForTrackHistograms("InCaloOutVertexPions",config,pDebugStream);
00150   getConfigForTrackHistograms("OutCaloInVertexPions",config,pDebugStream);
00151   getConfigForTrackHistograms("AllMuons",config,pDebugStream);
00152   getConfigForTrackHistograms("InCaloInVertexMuons",config,pDebugStream);
00153   getConfigForTrackHistograms("InCaloOutVertexMuons",config,pDebugStream);
00154   getConfigForTrackHistograms("OutCaloInVertexMuons",config,pDebugStream); 
00155   getConfigForTrackHistograms("AllElectrons",config,pDebugStream);
00156   getConfigForTrackHistograms("InCaloInVertexElectrons",config,pDebugStream);
00157   getConfigForTrackHistograms("InCaloOutVertexElectrons",config,pDebugStream);
00158   getConfigForTrackHistograms("OutCaloInVertexElectrons",config,pDebugStream);
00159   
00160   if (verbose_) LogTrace(messageLoggerCatregory) << debugStream.str();
00161 }
00162 
00163 JPTJetAnalyzer::~JPTJetAnalyzer()
00164 {}
00165 
00166 void JPTJetAnalyzer::beginJob(DQMStore* dqmStore)
00167 {
00168   dqm_ = dqmStore;
00169   //book histograms
00170   dqmStore->setCurrentFolder(histogramPath_);
00171   bookHistograms(dqmStore);
00172 }
00173 
00174 
00175 void JPTJetAnalyzer::analyze(const edm::Event&             event,
00176                              const edm::EventSetup&        eventSetup,
00177                              const reco::JPTJetCollection& jptJets,
00178                              const int                     numPV)
00179 {
00180   double pt1 = 0;
00181   double pt2 = 0;
00182   double pt3 = 0;
00183 
00184   for (reco::JPTJetCollection::const_iterator iJet = jptJets.begin(); iJet != jptJets.end(); ++iJet) {
00185     analyze(event, eventSetup, *iJet, pt1, pt2, pt3, numPV);
00186   }
00187   if (pt1) fillHistogram(JetPt1_,pt1);
00188   if (pt2) fillHistogram(JetPt2_,pt2);
00189   if (pt3) fillHistogram(JetPt3_,pt3);
00190 }
00191 
00192 void JPTJetAnalyzer::analyze(const edm::Event&      event,
00193                              const edm::EventSetup& eventSetup,
00194                              const reco::JPTJet&    jptJet,
00195                              double&                pt1,
00196                              double&                pt2,
00197                              double&                pt3,
00198                              const int              numPV)
00199 {
00200   //update the track propagator and strip noise calculator
00201   trackPropagator_->update(eventSetup);
00202   sOverNCalculator_->update(eventSetup);
00203   
00204   //make corrected jets 
00205   /* */
00206   //  const double factorZSP = jptJet.getZSPCor();  
00207   const reco::Jet& rawJet = *jptJet.getCaloJetRef();
00208   /* */ 
00209   //  const double factorZSPJPT = jptJet.energy()/rawJet.energy(); 
00210   
00211   //check jet is correctable by JPT
00212   if ( fabs(rawJet.eta()) > 2.1) return;
00213   
00214   //get consitiuents of jet
00215   const reco::TrackRefVector& pionsInVertexInCalo = jptJet.getPionsInVertexInCalo();
00216   const reco::TrackRefVector& pionsInVertexOutCalo = jptJet.getPionsInVertexOutCalo();
00217   const reco::TrackRefVector& pionsOutVertexInCalo = jptJet.getPionsOutVertexInCalo();
00218   const reco::TrackRefVector& muonsInVertexInCalo = jptJet.getMuonsInVertexInCalo();
00219   const reco::TrackRefVector& muonsInVertexOutCalo = jptJet.getMuonsInVertexOutCalo();
00220   const reco::TrackRefVector& muonsOutVertexInCalo = jptJet.getMuonsOutVertexInCalo();
00221   const reco::TrackRefVector& electronsInVertexInCalo = jptJet.getElecsInVertexInCalo();
00222   const reco::TrackRefVector& electronsInVertexOutCalo = jptJet.getElecsInVertexOutCalo();
00223   const reco::TrackRefVector& electronsOutVertexInCalo = jptJet.getElecsOutVertexInCalo();
00224   
00225   //check pt against highest values
00226   const double pt = jptJet.pt();
00227   if (pt > pt1) {
00228     pt3 = pt2;
00229     pt2 = pt1;
00230     pt1 = pt;
00231   } else if (pt > pt2) {
00232     pt3 = pt2;
00233     pt2 = pt;
00234   } else if (pt > pt3) {
00235     pt3 = pt;
00236   }
00237   
00238   //fill histograms
00239   try {
00240     const reco::CaloJet& rawCaloJet = dynamic_cast<const reco::CaloJet&>(rawJet);
00241     jetID_->calculate(event,rawCaloJet);
00242   } catch (const std::bad_cast&) {
00243     edm::LogError(messageLoggerCatregory) << "Failed to cast raw jet to CaloJet. JPT Jet does not appear to have been built from a CaloJet. "
00244                                           << "Histograms not filled. ";
00245     return;
00246   }
00247   const bool idPassed = ( (jetID_->n90Hits() >= n90HitsMin_) && (jetID_->fHPD() < fHPDMax_) && (jetID_->restrictedEMF() >= resEMFMin_) );
00248   /*  fillHistogram(JetN90Hits_,jetID_->n90Hits());  
00249   fillHistogram(JetfHPD_,jetID_->fHPD());
00250   fillHistogram(JetResEMF_,jetID_->restrictedEMF());  
00251   fillHistogram(JetfRBX_,jetID_->fRBX());  */
00252   if (idPassed) {
00253     //    const double deltaEta = jptJet.eta() - rawJet.eta();
00254     //    const double deltaPhi = jptJet.phi() - rawJet.phi();
00255     if (jptJet.pt() > correctedPtMin_) {
00256       fillHistogram(JetE_,jptJet.energy());
00257       fillHistogram(JetEt_,jptJet.et());
00258       fillHistogram(JetP_,jptJet.p());
00259       /*      fillHistogram(JetMass_,jptJet.mass());  */
00260       fillHistogram(JetPt_,jptJet.pt());
00261       fillHistogram(JetPx_,jptJet.px());
00262       fillHistogram(JetPy_,jptJet.py());
00263       fillHistogram(JetPz_,jptJet.pz());
00264 
00265       fillHistogram(JetEta_,jptJet.eta());
00266       fillHistogram(JetPhi_,jptJet.phi());
00267       /*      fillHistogram(JetDeltaEta_,deltaEta);
00268               fillHistogram(JetDeltaPhi_,deltaPhi);    */
00269       fillHistogram(JetPhiVsEta_,jptJet.phi(),jptJet.eta());
00270 
00271 
00272       const uint16_t totalTracks = jptJet.chargedMultiplicity();
00273       fillHistogram(NTracksPerJetHisto_,totalTracks);
00274       fillHistogram(NTracksPerJetVsJetEtHisto_,rawJet.et(),totalTracks);
00275       fillHistogram(NTracksPerJetVsJetEtaHisto_,rawJet.eta(),totalTracks);
00276       /*      fillHistogram(CorrFactorHisto_,factorZSPJPT);
00277       fillHistogram(CorrFactorVsJetEtHisto_,rawJet.et(),factorZSPJPT);
00278       fillHistogram(CorrFactorVsJetEtaHisto_,rawJet.eta(),factorZSPJPT);
00279       fillHistogram(ZSPCorrFactorHisto_,factorZSP);
00280       fillHistogram(ZSPCorrFactorVsJetEtHisto_,rawJet.et(),factorZSP);
00281       fillHistogram(ZSPCorrFactorVsJetEtaHisto_,rawJet.eta(),factorZSP); */
00282       /*     const double factorJPT = factorZSPJPT / factorZSP;   */
00283       /*      fillHistogram(JPTCorrFactorHisto_,factorJPT);
00284       fillHistogram(JPTCorrFactorVsJetEtHisto_,rawJet.et(),factorJPT);
00285       fillHistogram(JPTCorrFactorVsJetEtaHisto_,rawJet.eta(),factorJPT);  */
00286       //      const double ptFractionInCone = findPtFractionInCone(pionsInVertexInCalo,pionsInVertexOutCalo);
00287       /*      fillHistogram(PtFractionInConeHisto_,ptFractionInCone);
00288       fillHistogram(PtFractionInConeVsJetRawEtHisto_,rawJet.et(),ptFractionInCone);
00289       fillHistogram(PtFractionInConeVsJetEtaHisto_,rawJet.eta(),ptFractionInCone); */
00290       //fill track level histograms
00291       fillTrackHistograms(allPionHistograms_,inCaloInVertexPionHistograms_,inCaloOutVertexPionHistograms_,outCaloInVertexPionHistograms_,
00292                           pionsInVertexInCalo,pionsOutVertexInCalo,pionsInVertexOutCalo,rawJet);
00293       fillTrackHistograms(allMuonHistograms_,inCaloInVertexMuonHistograms_,inCaloOutVertexMuonHistograms_,outCaloInVertexMuonHistograms_,
00294                           muonsInVertexInCalo,muonsOutVertexInCalo,muonsInVertexOutCalo,rawJet);
00295       fillTrackHistograms(allElectronHistograms_,inCaloInVertexElectronHistograms_,inCaloOutVertexElectronHistograms_,outCaloInVertexElectronHistograms_,
00296                           electronsInVertexInCalo,electronsOutVertexInCalo,electronsInVertexOutCalo,rawJet);
00297     }
00298   }
00299 }
00300 
00301 void JPTJetAnalyzer::endJob()
00302 {
00303   if(writeDQMStore_) dqm_->save(dqmStoreFileName_);
00304 }
00305 
00306 void JPTJetAnalyzer::getConfigForHistogram(const std::string& configName, const edm::ParameterSet& psetContainingConfigPSet,
00307                                            std::ostringstream* pDebugStream)
00308 {
00309   const std::string psetName = configName+std::string("HistogramConfig");
00310   if (!psetContainingConfigPSet.exists(psetName)) {
00311     edm::LogWarning(messageLoggerCatregory) << "Histogram " << configName << " config not found" << std::endl;
00312     histogramConfig_[configName] = HistogramConfig();
00313   } else {
00314     if (pDebugStream) {
00315       (*pDebugStream) << "Histogram " << configName << " config found and loaded" << std::endl;
00316     }
00317     const edm::ParameterSet& pset = psetContainingConfigPSet.getParameter<edm::ParameterSet>(psetName);
00318     const bool enabled = (pset.exists("Enabled") ? pset.getParameter<bool>("Enabled") : true);
00319     if (!enabled) {
00320       histogramConfig_[configName] = HistogramConfig();
00321       if (pDebugStream) {
00322         (*pDebugStream) << "\tHistogram: " << configName << " Disabled" << std::endl;
00323       }
00324     } else {
00325       const unsigned int nBins = (pset.exists("NBins") ? pset.getParameter<unsigned int>("NBins") : 0);
00326       const double min = (pset.exists("Min") ? pset.getParameter<double>("Min") : 0);
00327       const double max = (pset.exists("Max") ? pset.getParameter<double>("Max") : 0);
00328       const unsigned int nBinsY = (pset.exists("NBinsY") ? pset.getParameter<unsigned int>("NBinsY") : 0);
00329       const double minY = (pset.exists("MinY") ? pset.getParameter<double>("MinY") : 0);
00330       const double maxY = (pset.exists("MaxY") ? pset.getParameter<double>("MaxY") : 0);
00331       if (nBins) {
00332         if (pDebugStream) {
00333           (*pDebugStream) << "\tHistogram: " << configName << "\tEnabled"
00334                           << "\tNBins: " << nBins << "\tMin: " << min << "\tMax: " << max;
00335           if (nBinsY) (*pDebugStream) << "\tNBinsY: " << nBinsY << "\tMinY: " << minY << "\tMaxY: " << maxY;
00336           (*pDebugStream) << std::endl;
00337         }
00338         if (nBinsY) {
00339           histogramConfig_[configName] = HistogramConfig(nBins,min,max,nBinsY,minY,maxY);
00340         } else {
00341           histogramConfig_[configName] = HistogramConfig(nBins,min,max);
00342         }
00343       }
00344     }
00345   }
00346 }
00347 
00348 void JPTJetAnalyzer::getConfigForTrackHistograms(const std::string& tag, const edm::ParameterSet& psetContainingConfigPSet,
00349                                                   std::ostringstream* pDebugStream)
00350 {
00351   getConfigForHistogram("n"+tag+"TracksPerJet",psetContainingConfigPSet,pDebugStream);
00352   getConfigForHistogram(tag+"TrackPt",psetContainingConfigPSet,pDebugStream);
00353   getConfigForHistogram(tag+"TrackPhi",psetContainingConfigPSet,pDebugStream);
00354   getConfigForHistogram(tag+"TrackEta",psetContainingConfigPSet,pDebugStream);
00355   /*  getConfigForHistogram(tag+"TrackNHits",psetContainingConfigPSet,pDebugStream);
00356       getConfigForHistogram(tag+"TrackNLayers",psetContainingConfigPSet,pDebugStream);
00357       getConfigForHistogram(tag+"TrackNLayers",psetContainingConfigPSet,pDebugStream);     */
00358   getConfigForHistogram(tag+"TrackPtVsEta",psetContainingConfigPSet,pDebugStream);
00359   /*  getConfigForHistogram(tag+"TrackDz",psetContainingConfigPSet,pDebugStream);
00360   getConfigForHistogram(tag+"TrackDxy",psetContainingConfigPSet,pDebugStream); */
00361 }
00362 
00363 MonitorElement* JPTJetAnalyzer::bookHistogram(const std::string& name, const std::string& title, const std::string& xAxisTitle, DQMStore* dqm)
00364 {
00365   std::map<std::string,HistogramConfig>::const_iterator configIterator = histogramConfig_.find(name);
00366   if (configIterator == histogramConfig_.end()) {
00367     edm::LogWarning(messageLoggerCatregory) << "Trying to book histogram with name " << name << " when no config was not retrieved from ParameterSet";
00368     return NULL;
00369   }
00370   const HistogramConfig& histoConfig = (*configIterator).second;
00371   if (histoConfig.enabled) {
00372     MonitorElement* histo = dqm->book1D(name,title,histoConfig.nBins,histoConfig.min,histoConfig.max);
00373     histo->setAxisTitle(xAxisTitle,1);
00374     return histo;
00375   } else {
00376     return NULL;
00377   }
00378 }
00379 
00380 MonitorElement* JPTJetAnalyzer::book2DHistogram(const std::string& name, const std::string& title,
00381                                                 const std::string& xAxisTitle, const std::string& yAxisTitle, DQMStore* dqm)
00382 {
00383   std::map<std::string,HistogramConfig>::const_iterator configIterator = histogramConfig_.find(name);
00384   if (configIterator == histogramConfig_.end()) {
00385     edm::LogWarning(messageLoggerCatregory) << "Trying to book histogram with name " << name << " when no config was not retrieved from ParameterSet";
00386     return NULL;
00387   }
00388   const HistogramConfig& histoConfig = (*configIterator).second;
00389   if (histoConfig.enabled) {
00390     MonitorElement* histo = dqm->book2D(name,title,histoConfig.nBins,histoConfig.min,histoConfig.max,histoConfig.nBinsY,histoConfig.minY,histoConfig.maxY);
00391     histo->setAxisTitle(xAxisTitle,1);
00392     histo->setAxisTitle(yAxisTitle,2);
00393     return histo;
00394   } else {
00395     return NULL;
00396   }
00397 }
00398 
00399 MonitorElement* JPTJetAnalyzer::bookProfile(const std::string& name, const std::string& title,
00400                                             const std::string& xAxisTitle, const std::string& yAxisTitle, DQMStore* dqm)
00401 {
00402   std::map<std::string,HistogramConfig>::const_iterator configIterator = histogramConfig_.find(name);
00403   if (configIterator == histogramConfig_.end()) {
00404     edm::LogWarning(messageLoggerCatregory) << "Trying to book histogram with name " << name << " when no config was not retrieved from ParameterSet";
00405     return NULL;
00406   }
00407   const HistogramConfig& histoConfig = (*configIterator).second;
00408   if (histoConfig.enabled) {
00409     TProfile* underlyingRootObject = new TProfile(name.c_str(),title.c_str(),histoConfig.nBins,histoConfig.min,histoConfig.max);
00410     MonitorElement* histo = dqm->bookProfile(name,underlyingRootObject);
00411     histo->setAxisTitle(xAxisTitle,1);
00412     histo->setAxisTitle(yAxisTitle,2);
00413     return histo;
00414   } else {
00415     return NULL;
00416   }
00417 }
00418 
00419 void JPTJetAnalyzer::bookHistograms(DQMStore* dqm)
00420 {
00421   JetE_        = bookHistogram("E","Corrected Jet Energy","E /GeV",dqm);
00422   JetEt_       = bookHistogram("Et","Corrected Jet E_{T}","E_{T} /GeV",dqm);
00423   JetP_        = bookHistogram("P","Corrected Jet Momentum","p /GeV/c",dqm);
00424   /*  JetMass_     = bookHistogram("Mass","Jet mass","Mass /GeV/c^{2}",dqm);  */
00425   JetPt_       = bookHistogram("Pt","Jet p_{T}","p_{T} /GeV/c",dqm);
00426   JetPt1_      = bookHistogram("Pt1","1st Jet p_{T}","p_{T} /GeV/c",dqm);
00427   JetPt2_      = bookHistogram("Pt2","2nd Jet p_{T}","p_{T} /GeV/c",dqm);
00428   JetPt3_      = bookHistogram("Pt3","3rd Jet p_{T}","p_{T} /GeV/c",dqm);
00429   JetPx_       = bookHistogram("Px","Jet p_{X}","p_{X} /GeV/c",dqm);
00430   JetPy_       = bookHistogram("Py","Jet p_{Y}","p_{Y} /GeV/c",dqm);
00431   JetPz_       = bookHistogram("Pz","Jet p_{Z}","p_{Z} /GeV/c",dqm);
00432   JetEta_      = bookHistogram("Eta","Jet #eta","#eta",dqm);
00433   JetPhi_      = bookHistogram("Phi","Jet #phi","#phi",dqm);
00434   /*  JetDeltaEta_ = bookHistogram("deltaEta","Change in #eta","#Delta #eta",dqm);
00435       JetDeltaPhi_ = bookHistogram("deltaPhi","Change in #phi","#Delta #phi",dqm);  */
00436   JetPhiVsEta_ = book2DHistogram("PhiVsEta","Corrected jet #phi vs #eta","jet #phi","jet #eta",dqm);
00437   /*  JetN90Hits_  = bookHistogram("N90Hits","Jet N90Hits","N90 Hits",dqm);  
00438   JetfHPD_     = bookHistogram("fHPD","Jet fHPD","fHPD",dqm);
00439   JetResEMF_   = bookHistogram("ResEMF","Jet restricted EM fraction","restricted EMF",dqm); 
00440   JetfRBX_     = bookHistogram("fRBX","Jet fRBX","fRBX",dqm); 
00441 
00442 
00443   TrackSiStripHitStoNHisto_            = bookHistogram("TrackSiStripHitStoN","Signal to noise of track SiStrip hits","S/N",dqm);  */
00444   InCaloTrackDirectionJetDRHisto_      = bookHistogram("InCaloTrackDirectionJetDR",
00445                                                        "#Delta R between track direrction at vertex and jet axis (track in cone at calo)","#Delta R",dqm);
00446   OutCaloTrackDirectionJetDRHisto_     = bookHistogram("OutCaloTrackDirectionJetDR",
00447                                                        "#Delta R between track direrction at vertex and jet axis (track out of cone at calo)","#Delta R",dqm);
00448   InVertexTrackImpactPointJetDRHisto_  = bookHistogram("InVertexTrackImpactPointJetDR",
00449                                                        "#Delta R between track impact point on calo and jet axis (track in cone at vertex)","#Delta R",dqm);
00450   OutVertexTrackImpactPointJetDRHisto_ = bookHistogram("OutVertexTrackImpactPointJetDR",
00451                                                        "#Delta R between track impact point on calo and jet axis (track out of cone at vertex)","#Delta R",dqm);
00452   
00453   NTracksPerJetHisto_         = bookHistogram("nTracks","Number of tracks for correction per jet","n tracks",dqm);
00454   NTracksPerJetVsJetEtHisto_  = bookProfile("nTracksVsJetEt","Number of tracks for correction per jet vs jet raw E_{T}","Jet raw E_{T} /GeV","n Tracks",dqm);
00455   NTracksPerJetVsJetEtaHisto_ = bookProfile("nTracksVsJetEta","Number of tracks for correction per jet vs jet #eta","Jet #eta","n Tracks",dqm);
00456   
00457   /*  PtFractionInConeHisto_           = bookHistogram("PtFractionInCone","#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}}",
00458                                                  "#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}}",dqm);   
00459   PtFractionInConeVsJetRawEtHisto_ = bookProfile("PtFractionInConeVsJetRawEt","#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}} vs jet raw E_{T}",
00460                                                  "Jet raw E_{T} / GeV","#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}}",dqm);
00461   PtFractionInConeVsJetEtaHisto_   = bookProfile("PtFractionInConeVsJetEta","#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}} vs jet #eta",
00462                                                  "Jet #eta","#frac{p_{T}^{in-cone}}{p_{T}^{in-cone}+p_{T}^{out-of-cone}}",dqm);
00463   
00464   CorrFactorHisto_            = bookHistogram("CorrFactor","Correction factor","Correction factor",dqm);
00465   CorrFactorVsJetEtHisto_     = bookProfile("CorrFactorVsJetEt","Correction factor vs jet raw E_{T}","Jet raw E_{T}","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);
00466   CorrFactorVsJetEtaHisto_    = bookProfile("CorrFactorVsJetEta","Correction factor vs jet #eta","Jet #eta","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);
00467   ZSPCorrFactorHisto_         = bookHistogram("ZSPCorrFactor","Correction factor from ZSP step","Correction factor",dqm);
00468   ZSPCorrFactorVsJetEtHisto_  = bookProfile("ZSPCorrFactorVsJetEt","Correction factor from ZSP step vs jet raw E_{T}",
00469                                             "Jet raw E_{T}","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);
00470   ZSPCorrFactorVsJetEtaHisto_ = bookProfile("ZSPCorrFactorVsJetEta","Correction factor from ZSP step vs jet #eta",
00471                                             "Jet #eta","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);
00472   JPTCorrFactorHisto_         = bookHistogram("JPTCorrFactor","Correction factor from JPT step","Correction factor",dqm);
00473   JPTCorrFactorVsJetEtHisto_  = bookProfile("JPTCorrFactorVsJetEt","Correction factor from JPT step vs jet raw E_{T}",
00474                                             "Jet raw E_{T}","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);
00475   JPTCorrFactorVsJetEtaHisto_ = bookProfile("JPTCorrFactorVsJetEta","Correction factor from JPT step vs jet #eta",
00476   "Jet #eta","#frac{E_{T}^{corr}}{E_{T}^{raw}}",dqm);   */
00477 
00478 
00479   
00480   bookTrackHistograms(&allPionHistograms_,"AllPions","pion",NULL,NULL,dqm);
00481   bookTrackHistograms(&inCaloInVertexPionHistograms_,"InCaloInVertexPions","pions in cone at calo and vertex",
00482                       InCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00483   bookTrackHistograms(&inCaloOutVertexPionHistograms_,"InCaloOutVertexPions","pions in cone at calo but out at vertex",
00484                       InCaloTrackDirectionJetDRHisto_,OutVertexTrackImpactPointJetDRHisto_,dqm);
00485   bookTrackHistograms(&outCaloInVertexPionHistograms_,"OutCaloInVertexPions","pions out of cone at calo but in at vertex",
00486                       OutCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00487   
00488   bookTrackHistograms(&allMuonHistograms_,"AllMuons","muon",NULL,NULL,dqm);
00489   bookTrackHistograms(&inCaloInVertexMuonHistograms_,"InCaloInVertexMuons","muons in cone at calo and vertex",
00490                       InCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00491   bookTrackHistograms(&inCaloOutVertexMuonHistograms_,"InCaloOutVertexMuons","muons in cone at calo but out at vertex",
00492                       InCaloTrackDirectionJetDRHisto_,OutVertexTrackImpactPointJetDRHisto_,dqm);
00493   bookTrackHistograms(&outCaloInVertexMuonHistograms_,"OutCaloInVertexMuons","muons out of cone at calo but in at vertex",
00494                       OutCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00495   
00496   bookTrackHistograms(&allElectronHistograms_,"AllElectrons","electron",NULL,NULL,dqm);
00497   bookTrackHistograms(&inCaloInVertexElectronHistograms_,"InCaloInVertexElectrons","electrons in cone at calo and vertex",
00498                       InCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00499   bookTrackHistograms(&inCaloOutVertexElectronHistograms_,"InCaloOutVertexElectrons","electrons in cone at calo but out at vertex",
00500                       InCaloTrackDirectionJetDRHisto_,OutVertexTrackImpactPointJetDRHisto_,dqm);
00501   bookTrackHistograms(&outCaloInVertexElectronHistograms_,"OutCaloInVertexElectrons","electrons out of cone at calo but in at vertex",
00502                       OutCaloTrackDirectionJetDRHisto_,InVertexTrackImpactPointJetDRHisto_,dqm);
00503 }
00504 
00505 void JPTJetAnalyzer::bookTrackHistograms(TrackHistograms* histos, const std::string& tag, const std::string& titleTag,
00506                                          MonitorElement* trackDirectionJetDRHisto, MonitorElement* trackImpactPointJetDRHisto, DQMStore* dqm)
00507 {
00508   histos->nTracksHisto = bookHistogram("n"+tag+"TracksPerJet","Number of "+titleTag+" tracks per jet","n Tracks",dqm);
00509   histos->ptHisto = bookHistogram(tag+"TrackPt",titleTag+" p_{T}","p_{T} /GeV/c",dqm);
00510   histos->phiHisto = bookHistogram(tag+"TrackPhi",titleTag+" track #phi","#phi",dqm);
00511   histos->etaHisto = bookHistogram(tag+"TrackEta",titleTag+" track #eta","#eta",dqm);
00512   /*  histos->nHitsHisto = bookHistogram(tag+"TrackNHits",titleTag+" track N hits","N hits",dqm);
00513       histos->nLayersHisto = bookHistogram(tag+"TrackNLayers",titleTag+" track N layers with hits","N layers",dqm);    */
00514   histos->ptVsEtaHisto = bookProfile(tag+"TrackPtVsEta",titleTag+" track p_{T} vs #eta","#eta","p_{T} /GeV/c",dqm);
00515   /*  histos->dzHisto = bookHistogram(tag+"TrackDz",titleTag+" track dz","dz /cm",dqm);
00516   histos->dxyHisto = bookHistogram(tag+"TrackDxy",titleTag+" track dxy","dxy /cm",dqm); */
00517   histos->trackDirectionJetDRHisto = trackDirectionJetDRHisto;
00518   histos->trackImpactPointJetDRHisto = trackImpactPointJetDRHisto;
00519 }
00520 
00521 void JPTJetAnalyzer::fillTrackHistograms(TrackHistograms& allTracksHistos, TrackHistograms& inCaloInVertexHistos,
00522                                          TrackHistograms& inCaloOutVertexHistos, TrackHistograms& outCaloInVertexHistos,
00523                                          const reco::TrackRefVector& inVertexInCalo,
00524                                          const reco::TrackRefVector& outVertexInCalo,
00525                                          const reco::TrackRefVector& inVertexOutCalo,
00526                                          const reco::Jet& rawJet)
00527 {
00528   fillTrackHistograms(allTracksHistos,inVertexInCalo,rawJet);
00529   fillTrackHistograms(allTracksHistos,outVertexInCalo,rawJet);
00530   fillTrackHistograms(allTracksHistos,inVertexOutCalo,rawJet);
00531   fillHistogram(allTracksHistos.nTracksHisto,inVertexInCalo.size()+outVertexInCalo.size()+inVertexOutCalo.size());
00532   fillSiStripSoNForTracks(inVertexInCalo);
00533   fillSiStripSoNForTracks(outVertexInCalo);
00534   fillSiStripSoNForTracks(inVertexOutCalo);
00535   fillTrackHistograms(inCaloInVertexHistos,inVertexInCalo,rawJet);
00536   fillHistogram(inCaloInVertexHistos.nTracksHisto,inVertexInCalo.size());
00537   fillTrackHistograms(inCaloOutVertexHistos,outVertexInCalo,rawJet);
00538   fillHistogram(inCaloOutVertexHistos.nTracksHisto,outVertexInCalo.size());
00539   fillTrackHistograms(outCaloInVertexHistos,inVertexOutCalo,rawJet);
00540   fillHistogram(outCaloInVertexHistos.nTracksHisto,inVertexOutCalo.size());
00541 }
00542 
00543 void JPTJetAnalyzer::fillTrackHistograms(TrackHistograms& histos, const reco::TrackRefVector& tracks, const reco::Jet& rawJet)
00544 {
00545   const reco::TrackRefVector::const_iterator tracksEnd = tracks.end();
00546   for (reco::TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracksEnd; ++iTrack) {
00547     const reco::Track& track = **iTrack;
00548     const double pt = track.pt();
00549     const double phi = track.phi();
00550     const double eta = track.eta();
00551     //    const unsigned int nHits = track.found();
00552     //    const unsigned int nLayers = track.hitPattern().trackerLayersWithMeasurement();
00553     //    const double dz = track.dz();
00554     //    const double dxy = track.dxy();  
00555     fillHistogram(histos.ptHisto,pt);
00556     fillHistogram(histos.phiHisto,phi);
00557     fillHistogram(histos.etaHisto,eta);
00558     /*    fillHistogram(histos.nHitsHisto,nHits);
00559           fillHistogram(histos.nLayersHisto,nLayers);   */
00560     fillHistogram(histos.ptVsEtaHisto,eta,pt);
00561     /*    fillHistogram(histos.dzHisto,dz);
00562     fillHistogram(histos.dxyHisto,dxy);   */
00563     const double trackDirectionJetDR = deltaR(rawJet,track);
00564     fillHistogram(histos.trackDirectionJetDRHisto,trackDirectionJetDR);
00565     const double impactPointJetDR = deltaR(rawJet,trackPropagator_->impactPoint(track));
00566     fillHistogram(histos.trackImpactPointJetDRHisto,impactPointJetDR);
00567   }
00568 }
00569 
00570 void JPTJetAnalyzer::fillSiStripSoNForTracks(const reco::TrackRefVector& tracks)
00571 {
00572   const reco::TrackRefVector::const_iterator tracksEnd = tracks.end();
00573   for (reco::TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracksEnd; ++iTrack) {
00574     const trackingRecHit_iterator trackRecHitsEnd = (*iTrack)->recHitsEnd();
00575     for (trackingRecHit_iterator iHit = (*iTrack)->recHitsBegin(); iHit != trackRecHitsEnd; ++iHit) {
00576       fillSiStripHitSoN(**iHit);
00577     }
00578   }
00579 }
00580 
00581 void JPTJetAnalyzer::fillSiStripHitSoN(const TrackingRecHit& hit)
00582 {
00583   //check it is an SiStrip hit
00584   DetId detId(hit.geographicalId());
00585   if (!( (detId.det() == DetId::Tracker) &&
00586         ( (detId.subdetId() == SiStripDetId::TIB) ||
00587           (detId.subdetId() == SiStripDetId::TID) ||
00588           (detId.subdetId() == SiStripDetId::TOB) ||
00589           (detId.subdetId() == SiStripDetId::TEC)
00590         )
00591      )) return;
00592 
00593 
00594   // J.Piedra, 2012/09/24
00595   //try to determine the type of the hit
00596   //  const TrackingRecHit* pHit = &hit;
00597   //  const SiStripRecHit2D* pRecHit2D = dynamic_cast<const SiStripRecHit2D*>(pHit);
00598   //  const SiStripMatchedRecHit2D* pMatchedRecHit2D = dynamic_cast<const SiStripMatchedRecHit2D*>(pHit);
00599   //  const ProjectedSiStripRecHit2D* pProjctedRecHit2D = dynamic_cast<const ProjectedSiStripRecHit2D*>(pHit);
00600   //  const InvalidTrackingRecHit* pInvalidHit = dynamic_cast<const InvalidTrackingRecHit*>(pHit);
00601 
00602   //  //fill signal to noise for appropriate hit
00603   //  if (pMatchedRecHit2D) {
00604   //    fillSiStripHitSoNForSingleHit(pMatchedRecHit2D->monoHit());
00605   //    fillSiStripHitSoNForSingleHit(pMatchedRecHit2D->stereoHit());
00606   //  } else if (pProjctedRecHit2D) {
00607   //    fillSiStripHitSoNForSingleHit(pProjctedRecHit2D->originalHit());
00608   //  } else if (pRecHit2D) {
00609   //    fillSiStripHitSoNForSingleHit(*pRecHit2D);
00610   //  } else if (pInvalidHit) {
00611   //    return;
00612   //  } else {
00613   //    edm::LogInfo(messageLoggerCatregory) << "Hit on det ID " << hit.geographicalId().rawId() << " cannot be converted to a strip hit";
00614   //  }
00615 }
00616 
00617 
00618 // J.Piedra, 2012/09/24
00619 // void JPTJetAnalyzer::fillSiStripHitSoNForSingleHit(const SiStripRecHit2D& hit)
00620 // {
00621 //   //get the cluster
00622 //   const SiStripCluster* cluster = NULL;
00623 //   const SiStripRecHit2D::ClusterRegionalRef& regionalClusterRef = hit.cluster_regional();
00624 //   const SiStripRecHit2D::ClusterRef& normalClusterRef = hit.cluster();
00625 //   if (regionalClusterRef.isNonnull()) {
00626 //     cluster = &*regionalClusterRef;
00627 //   } else if (normalClusterRef.isNonnull()) {
00628 //     cluster = &*normalClusterRef;
00629 //   } else {
00630 //     edm::LogError(messageLoggerCatregory) << "Unable to get cluster from SiStripRecHit2D with det ID " << hit.geographicalId().rawId();
00631 //     return;
00632 //   }
00633 //   //calculate signal to noise for cluster
00634 //   //  const double sOverN = (*sOverNCalculator_)(*cluster,hit.geographicalId());
00635 //   //fill histogram
00636 //   /*  fillHistogram(TrackSiStripHitStoNHisto_,sOverN);  */
00637 // }
00638 
00639 
00640 double JPTJetAnalyzer::findPtFractionInCone(const reco::TrackRefVector& inConeTracks, const reco::TrackRefVector& outOfConeTracks)
00641 {
00642   double totalPt = 0;
00643   double inConePt = 0;
00644   const reco::TrackRefVector::const_iterator inConeTracksEnd = inConeTracks.end();
00645   for (reco::TrackRefVector::const_iterator iInConeTrack = inConeTracks.begin(); iInConeTrack != inConeTracksEnd; ++iInConeTrack) {
00646     const double pt = (*iInConeTrack)->pt();
00647     totalPt += pt;
00648     inConePt += pt;
00649   }
00650   const reco::TrackRefVector::const_iterator outOfConeTracksEnd = outOfConeTracks.end();
00651   for (reco::TrackRefVector::const_iterator iOutOfConeTrack = outOfConeTracks.begin(); iOutOfConeTrack != outOfConeTracksEnd; ++iOutOfConeTrack) {
00652     const double pt = (*iOutOfConeTrack)->pt();
00653     totalPt += pt;
00654   }
00655   if (totalPt) return inConePt/totalPt;
00656   //return 0 if there are no tracks at all
00657   else return 0;
00658 }
00659 
00660 
00661 
00662 JPTJetAnalyzer::HistogramConfig::HistogramConfig()
00663   : enabled(false),
00664     nBins(0),
00665     min(0),
00666     max(0),
00667     nBinsY(0),
00668     minY(0),
00669     maxY(0)
00670 {}
00671 
00672 JPTJetAnalyzer::HistogramConfig::HistogramConfig(const unsigned int theNBins, const double theMin, const double theMax)
00673   : enabled(true),
00674     nBins(theNBins),
00675     min(theMin),
00676     max(theMax),
00677     nBinsY(0),
00678     minY(0),
00679     maxY(0)
00680 {}
00681 
00682 JPTJetAnalyzer::HistogramConfig::HistogramConfig(const unsigned int theNBinsX, const double theMinX, const double theMaxX,
00683                                                  const unsigned int theNBinsY, const double theMinY, const double theMaxY)
00684   : enabled(true),
00685     nBins(theNBinsX),
00686     min(theMinX),
00687     max(theMaxX),
00688     nBinsY(theNBinsY),
00689     minY(theMinY),
00690     maxY(theMaxY)
00691 {}
00692 
00693 JPTJetAnalyzer::TrackHistograms::TrackHistograms()
00694   : ptHisto(NULL),
00695     phiHisto(NULL),
00696     etaHisto(NULL),
00697     /*    nHitsHisto(NULL),
00698           nLayersHisto(NULL), */
00699     ptVsEtaHisto(NULL),
00700     /*    dzHisto(NULL),
00701     dxyHisto(NULL),  */
00702     trackDirectionJetDRHisto(NULL),
00703     trackImpactPointJetDRHisto(NULL)
00704 {}
00705 
00706 JPTJetAnalyzer::TrackHistograms::TrackHistograms(MonitorElement* theNTracksHisto, 
00707 MonitorElement* thePtHisto, MonitorElement* thePhiHisto, MonitorElement* theEtaHisto, 
00708 /*  MonitorElement* theNHitsHisto, MonitorElement* theNLayersHisto, */ 
00709 MonitorElement* thePtVsEtaHisto, 
00710 /* MonitorElement* theDzHisto, MonitorElement* theDxyHisto, */ 
00711 MonitorElement* theTrackDirectionJetDRHisto, MonitorElement* theTrackImpactPointJetDRHisto)
00712   : nTracksHisto(theNTracksHisto),
00713     ptHisto(thePtHisto),
00714     phiHisto(thePhiHisto),
00715     etaHisto(theEtaHisto),
00716     /*    nHitsHisto(theNHitsHisto),
00717           nLayersHisto(theNLayersHisto),  */
00718     ptVsEtaHisto(thePtVsEtaHisto),
00719     /*    dzHisto(theDzHisto),
00720     dxyHisto(theDxyHisto),  */
00721     trackDirectionJetDRHisto(theTrackDirectionJetDRHisto),
00722     trackImpactPointJetDRHisto(theTrackImpactPointJetDRHisto)
00723 {}
00724 
00725 namespace jptJetAnalysis {
00726   
00727   TrackPropagatorToCalo::TrackPropagatorToCalo()
00728     : magneticField_(NULL),
00729       propagator_(NULL),
00730       magneticFieldCacheId_(0),
00731       propagatorCacheId_(0)
00732     {}
00733   
00734   void TrackPropagatorToCalo::update(const edm::EventSetup& eventSetup)
00735   {
00736     //update magnetic filed if necessary
00737     const IdealMagneticFieldRecord& magneticFieldRecord = eventSetup.get<IdealMagneticFieldRecord>();
00738     const uint32_t newMagneticFieldCacheId = magneticFieldRecord.cacheIdentifier();
00739     if ((newMagneticFieldCacheId != magneticFieldCacheId_) || !magneticField_) {
00740       edm::ESHandle<MagneticField> magneticFieldHandle;
00741       magneticFieldRecord.get(magneticFieldHandle);
00742       magneticField_ = magneticFieldHandle.product();
00743       magneticFieldCacheId_ = newMagneticFieldCacheId;
00744     }
00745     //update propagator if necessary
00746     const TrackingComponentsRecord& trackingComponentsRecord = eventSetup.get<TrackingComponentsRecord>();
00747     const uint32_t newPropagatorCacheId = trackingComponentsRecord.cacheIdentifier();
00748     if ((propagatorCacheId_ != newPropagatorCacheId) || !propagator_) {
00749       edm::ESHandle<Propagator> propagatorHandle;
00750       trackingComponentsRecord.get("SteppingHelixPropagatorAlong",propagatorHandle);
00751       propagator_ = propagatorHandle.product();
00752       propagatorCacheId_ = newPropagatorCacheId;
00753     }
00754   }
00755   
00756   inline math::XYZPoint TrackPropagatorToCalo::impactPoint(const reco::Track& track) const
00757   {
00758     return JetTracksAssociationDRCalo::propagateTrackToCalorimeter(track,*magneticField_,*propagator_);
00759   }
00760   
00761   StripSignalOverNoiseCalculator::StripSignalOverNoiseCalculator(const std::string& theQualityLabel)
00762     : qualityLabel_(theQualityLabel),
00763       quality_(NULL),
00764       noise_(NULL),
00765       gain_(NULL),
00766       qualityCacheId_(0),
00767       noiseCacheId_(0),
00768       gainCacheId_(0)
00769     {}
00770   
00771   void StripSignalOverNoiseCalculator::update(const edm::EventSetup& eventSetup)
00772   {
00773     //update the quality if necessary
00774     const SiStripQualityRcd& qualityRecord = eventSetup.get<SiStripQualityRcd>();
00775     const uint32_t newQualityCacheId = qualityRecord.cacheIdentifier();
00776     if ((newQualityCacheId != qualityCacheId_) || !quality_) {
00777       edm::ESHandle<SiStripQuality> qualityHandle;
00778       qualityRecord.get(qualityLabel_,qualityHandle);
00779       quality_ = qualityHandle.product();
00780       qualityCacheId_ = newQualityCacheId;
00781     }
00782     //update the noise if necessary
00783     const SiStripNoisesRcd& noiseRecord = eventSetup.get<SiStripNoisesRcd>();
00784     const uint32_t newNoiseCacheId = noiseRecord.cacheIdentifier();
00785     if ((newNoiseCacheId != noiseCacheId_) || !noise_) {
00786       edm::ESHandle<SiStripNoises> noiseHandle;
00787       noiseRecord.get(noiseHandle);
00788       noise_ = noiseHandle.product();
00789       noiseCacheId_ = newNoiseCacheId;
00790     }
00791     //update the gain if necessary
00792     const SiStripGainRcd& gainRecord = eventSetup.get<SiStripGainRcd>();
00793     const uint32_t newGainCacheId = gainRecord.cacheIdentifier();
00794     if ((newGainCacheId != gainCacheId_) || !gain_) {
00795       edm::ESHandle<SiStripGain> gainHandle;
00796       gainRecord.get(gainHandle);
00797       gain_ = gainHandle.product();
00798       gainCacheId_ = newGainCacheId;
00799     }
00800   }
00801   
00802   double StripSignalOverNoiseCalculator::signalOverNoise(const SiStripCluster& cluster,
00803                                                          const uint32_t& detId) const
00804   {
00805     //const uint32_t detId = cluster.geographicalId();
00806     
00807     const uint16_t firstStrip = cluster.firstStrip();
00808     const SiStripQuality::Range& qualityRange = quality_->getRange(detId);
00809     const SiStripNoises::Range& noiseRange = noise_->getRange(detId);
00810     const SiStripApvGain::Range& gainRange = gain_->getRange(detId);
00811     double signal = 0;
00812     double noise2 = 0;
00813     unsigned int nNonZeroStrips = 0;
00814     const std::vector<uint8_t>& clusterAmplitudes = cluster.amplitudes();
00815     const std::vector<uint8_t>::const_iterator clusterAmplitudesEnd = clusterAmplitudes.end();
00816     const std::vector<uint8_t>::const_iterator clusterAmplitudesBegin = clusterAmplitudes.begin();
00817     for (std::vector<uint8_t>::const_iterator iAmp = clusterAmplitudesBegin; iAmp != clusterAmplitudesEnd; ++iAmp) {
00818       const uint8_t adc = *iAmp;
00819       const uint16_t strip = iAmp-clusterAmplitudesBegin+firstStrip;
00820       const bool stripBad = quality_->IsStripBad(qualityRange,strip);
00821       const double noise = noise_->getNoise(strip,noiseRange);
00822       const double gain = gain_->getStripGain(strip,gainRange);
00823       signal += adc;
00824       if (adc) ++nNonZeroStrips;
00825       const double noiseContrib = (stripBad ? 0 : noise/gain);
00826       noise2 += noiseContrib*noiseContrib;
00827     }
00828     const double noise = sqrt(noise2/nNonZeroStrips);
00829     if (noise) return signal/noise;
00830     else return 0;
00831   }
00832   
00833 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/DQMOffline/JetMET/src/JPTJetAnalyzer.cc
