Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:17:16

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackProbabilityXMLtoDB
0004 // Class:      TrackProbabilityXMLtoDB
0005 //
0006 /**\class TrackProbabilityXMLtoDB TrackProbabilityXMLtoDB.cc RecoBTag/TrackProbabilityXMLtoDB/src/TrackProbabilityXMLtoDB.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Andrea Rizzi
0015 //         Created:  Wed Apr 12 11:12:49 CEST 2006
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <iostream>
0023 using namespace std;
0024 
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0035 #include "DataFormats/Math/interface/Vector3D.h"
0036 #include "DataFormats/Common/interface/Ref.h"
0037 #include "DataFormats/JetReco/interface/Jet.h"
0038 #include "DataFormats/JetReco/interface/CaloJet.h"
0039 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 
0042 //#include "RecoBTag/TrackProbability/interface/TrackClassFilterCategory.h"
0043 
0044 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0045 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0046 //#include "TrackProbabilityCalibratedHistogram.h"
0047 
0048 #include "RecoBTag/BTagTools/interface/SignedTransverseImpactParameter.h"
0049 #include "RecoBTag/BTagTools/interface/SignedImpactParameter3D.h"
0050 #include "RecoBTag/BTagTools/interface/SignedDecayLength3D.h"
0051 
0052 //CondFormats
0053 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
0054 
0055 #include "FWCore/Framework/interface/IOVSyncValue.h"
0056 #include "FWCore/ServiceRegistry/interface/Service.h"
0057 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0058 // Math
0059 #include "Math/GenVector/VectorUtil.h"
0060 #include "Math/GenVector/PxPyPzE4D.h"
0061 
0062 #include "RecoBTag/XMLCalibration/interface/AlgorithmCalibration.h"
0063 #include "RecoBTag/XMLCalibration/interface/CalibratedHistogramXML.h"
0064 #include "RecoBTag/TrackProbability/interface/TrackClassFilterCategory.h"
0065 #include "TrackingTools/IPTools/interface/IPTools.h"
0066 
0067 //#include "TH1F.h"
0068 //#include "TFile.h"
0069 
0070 #include <fstream>
0071 #include <iostream>
0072 
0073 using namespace reco;
0074 
0075 //
0076 // class decleration
0077 //
0078 
0079 class CalibrationSkeleton : public edm::one::EDAnalyzer<> {
0080 public:
0081   explicit CalibrationSkeleton(const edm::ParameterSet&);
0082 
0083   virtual void beginJob() {
0084     bool resetData = true;
0085     bool newBinning = false;
0086     edm::FileInPath f2d("RecoBTag/TrackProbability/data/2DHisto.xml");
0087     edm::FileInPath f3d("RecoBTag/TrackProbability/data/3DHisto.xml");
0088     calibrationNew =
0089         new AlgorithmCalibration<TrackClassFilterCategory, CalibratedHistogramXML>((f3d.fullPath()).c_str());
0090     calibration2dNew =
0091         new AlgorithmCalibration<TrackClassFilterCategory, CalibratedHistogramXML>((f2d.fullPath()).c_str());
0092     if (resetData) {
0093       vector<pair<TrackClassFilterCategory, CalibratedHistogramXML> > data = calibrationNew->categoriesWithData();
0094       vector<pair<TrackClassFilterCategory, CalibratedHistogramXML> > data2d = calibration2dNew->categoriesWithData();
0095       for (unsigned int i = 0; i < data.size(); i++) {
0096         if (newBinning)
0097           data[i].second = CalibratedHistogram(1000, 0, 50);
0098         else
0099           data[i].second = CalibratedHistogram();
0100       }
0101       for (unsigned int i = 0; i < data2d.size(); i++) {
0102         if (newBinning)
0103           data2d[i].second = CalibratedHistogram(1000, 0, 50);
0104         else
0105           data2d[i].second = CalibratedHistogram();
0106       }
0107     }
0108     //    calibrationNew->startCalibration();
0109     //  calibration2dNew->startCalibration();
0110   }
0111 
0112   virtual void endJob() {
0113     edm::Service<cond::service::PoolDBOutputService> mydbservice;
0114     if (!mydbservice.isAvailable())
0115       return;
0116 
0117     vector<pair<TrackClassFilterCategory, CalibratedHistogramXML> > data = calibrationNew->categoriesWithData();
0118     vector<pair<TrackClassFilterCategory, CalibratedHistogramXML> > data2d = calibration2dNew->categoriesWithData();
0119     TrackProbabilityCalibration calibration;
0120     TrackProbabilityCalibration calibration2d;
0121 
0122     for (unsigned int i = 0; i < data.size(); i++) {
0123       TrackProbabilityCalibration::Entry entry;
0124       entry.category = data[i].first.categoryData();
0125       entry.histogram = data[i].second;
0126       calibration.data.push_back(entry);
0127     }
0128     for (unsigned int i = 0; i < data2d.size(); i++) {
0129       TrackProbabilityCalibration::Entry entry;
0130       entry.category = data2d[i].first.categoryData();
0131       entry.histogram = data2d[i].second;
0132       calibration2d.data.push_back(entry);
0133     }
0134 
0135     mydbservice->createOneIOV(calibration, mydbservice->endOfTime(), "BTagTrackProbability3DRcd");
0136 
0137     mydbservice->createOneIOV(calibration2d, mydbservice->endOfTime(), "BTagTrackProbability2DRcd");
0138   }
0139 
0140   ~CalibrationSkeleton() {}
0141 
0142   virtual void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
0143 
0144 private:
0145   AlgorithmCalibration<TrackClassFilterCategory, CalibratedHistogramXML>* calibrationNew;
0146   AlgorithmCalibration<TrackClassFilterCategory, CalibratedHistogramXML>* calibration2dNew;
0147 
0148   int count;
0149   int ntracks;
0150   int m_cutPixelHits;
0151   int m_cutTotalHits;
0152   double m_cutMaxTIP;
0153   double m_cutMinPt;
0154   double m_cutMaxDecayLen;
0155   double m_cutMaxChiSquared;
0156   double m_cutMaxLIP;
0157   double m_cutMaxDistToAxis;
0158   double m_cutMinProb;
0159 
0160   edm::InputTag m_assoc;
0161   edm::InputTag m_jets;
0162   edm::InputTag m_primaryVertexProducer;
0163 };
0164 
0165 //
0166 // constructors and destructor
0167 //
0168 CalibrationSkeleton::CalibrationSkeleton(const edm::ParameterSet& parameters) {}
0169 
0170 void CalibrationSkeleton::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0171   using namespace edm;
0172   using namespace reco;
0173   using namespace std;
0174 
0175   Handle<reco::VertexCollection> primaryVertex;
0176   iEvent.getByLabel("offlinePrimaryVertices", primaryVertex);
0177 
0178   //***********************************************************************************
0179   //look at reco vertices
0180   const reco::Vertex* pv;
0181   bool newvertex = false;
0182 
0183   bool pvFound = (primaryVertex->size() != 0);
0184   if (pvFound) {
0185     pv = &(*primaryVertex->begin());
0186   } else {
0187     reco::Vertex::Error e;
0188     e(0, 0) = 0.0015 * 0.0015;
0189     e(1, 1) = 0.0015 * 0.0015;
0190     e(2, 2) = 15. * 15.;
0191     reco::Vertex::Point p(0, 0, 0);
0192     pv = new reco::Vertex(p, e, 1, 1, 1);
0193     newvertex = true;
0194   }
0195 
0196   //*************************************************************************
0197   //look at JetTracks
0198   edm::Handle<JetTracksAssociationCollection> associationHandle;
0199   iEvent.getByLabel("ak5JetTracksAssociatorAtVertex", associationHandle);
0200   reco::JetTracksAssociationCollection::const_iterator it = associationHandle->begin();
0201 
0202   //***********************************************************************************
0203   //mandatory for ip significance reco
0204   const TransientTrackBuilder* m_transientTrackBuilder_producer;
0205   edm::ESHandle<TransientTrackBuilder> builder;
0206   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
0207   m_transientTrackBuilder_producer = builder.product();
0208   //***********************************************************************************
0209 
0210   int i = 0;
0211   for (; it != associationHandle->end(); it++, i++) {
0212     const JetTracksAssociationRef& jetTracks = Ref<JetTracksAssociationCollection>(associationHandle, i);
0213 
0214     //GlobalVector direction(jetTracks->key->px(),jetTracks->key->py(),jetTracks->key->pz());
0215 
0216     //     if(jetTracks->second.size() <2 ) continue;
0217 
0218     bool directionWithTracks = false;
0219     TrackRefVector tracks = it->second;
0220     math::XYZVector jetMomentum = it->first->momentum() / 2.;
0221     if (directionWithTracks) {
0222       for (TrackRefVector::const_iterator itTrack = tracks.begin(); itTrack != tracks.end(); ++itTrack) {
0223         if ((**itTrack).numberOfValidHits() >= m_cutTotalHits)  //minimal quality cuts
0224           jetMomentum += (**itTrack).momentum();
0225       }
0226     } else {
0227       jetMomentum = it->first->momentum();
0228     }
0229     GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
0230 
0231     for (RefVector<reco::TrackCollection>::const_iterator itt = tracks.begin(); itt != tracks.end(); itt++) {
0232       //loop on jets and tracks
0233 
0234       TrackClassFilterInput input(**itt, *(it->first), *pv);
0235 
0236       const TransientTrack& transientTrack = builder->build(&(**itt));
0237       const_cast<CalibratedHistogramXML*>(calibrationNew->getCalibData(input))
0238           ->fill(IPTools::signedImpactParameter3D(transientTrack, direction, *pv).second.significance());
0239       //  calibration2dNew->getCalibData(input)->Fill(significance2D);
0240 
0241       //       calibrationNew->updateCalibration(input);
0242       //       calibration2dNew->updateCalibration(input);
0243     }
0244   }
0245 }
0246 //define this as a plug-in
0247 DEFINE_FWK_MODULE(CalibrationSkeleton);