| 1 |
#include <iostream>
|
| 2 |
|
| 3 |
#include "TLorentzVector.h"
|
| 4 |
|
| 5 |
#include "UserCode/DTDPGAnalysis/interface/HFDimuons.h"
|
| 6 |
|
| 7 |
//#include "AnalysisDataFormats/HeavyFlavorObjects/rootio/TAna01Event.hh"
|
| 8 |
//#include "AnalysisDataFormats/HeavyFlavorObjects/rootio/TAnaTrack.hh"
|
| 9 |
//#include "AnalysisDataFormats/HeavyFlavorObjects/rootio/TAnaCand.hh"
|
| 10 |
//#include "AnalysisDataFormats/HeavyFlavorObjects/rootio/TGenCand.hh"
|
| 11 |
|
| 12 |
#include "DataFormats/VertexReco/interface/VertexFwd.h"
|
| 13 |
#include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
|
| 14 |
#include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
|
| 15 |
#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
|
| 16 |
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
|
| 17 |
#include "CommonTools/Statistics/interface/ChiSquared.h"
|
| 18 |
|
| 19 |
#include "RecoVertex/KinematicFitPrimitives/interface/ParticleMass.h"
|
| 20 |
#include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
|
| 21 |
#include <RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h>
|
| 22 |
#include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
|
| 23 |
#include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
|
| 24 |
#include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
|
| 25 |
#include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
|
| 26 |
#include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
|
| 27 |
|
| 28 |
#include "DataFormats/Common/interface/Handle.h"
|
| 29 |
#include "DataFormats/Common/interface/Wrapper.h"
|
| 30 |
#include "DataFormats/Candidate/interface/Candidate.h"
|
| 31 |
#include "DataFormats/MuonReco/interface/MuonFwd.h"
|
| 32 |
#include "DataFormats/MuonReco/interface/Muon.h"
|
| 33 |
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
|
| 34 |
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
| 35 |
#include "DataFormats/TrackReco/interface/Track.h"
|
| 36 |
|
| 37 |
//extern TAna01Event *gHFEvent;
|
| 38 |
|
| 39 |
using namespace std;
|
| 40 |
using namespace reco;
|
| 41 |
using namespace edm;
|
| 42 |
|
| 43 |
#define MMUON 0.10566
|
| 44 |
|
| 45 |
// ----------------------------------------------------------------------
|
| 46 |
HFDimuons::HFDimuons(const edm::ParameterSet& iConfig) :
|
| 47 |
fVerbose(iConfig.getUntrackedParameter<int>("verbose", 0)),
|
| 48 |
fTracksLabel(iConfig.getUntrackedParameter<InputTag>("tracksLabel", InputTag("goodTracks"))),
|
| 49 |
fPrimaryVertexLabel(iConfig.getUntrackedParameter<InputTag>("PrimaryVertexLabel", InputTag("offlinePrimaryVertices"))),
|
| 50 |
fMuonsLabel(iConfig.getUntrackedParameter<InputTag>("muonsLabel")),
|
| 51 |
fMuonPt(iConfig.getUntrackedParameter<double>("muonPt", 1.0)),
|
| 52 |
fMassLow(iConfig.getUntrackedParameter<double>("massLow", 0.0)),
|
| 53 |
fMassHigh(iConfig.getUntrackedParameter<double>("massHigh", 14000.)){
|
| 54 |
using namespace std;
|
| 55 |
cout << "----------------------------------------------------------------------" << endl;
|
| 56 |
cout << "--- HFDimuons constructor" << endl;
|
| 57 |
cout << "--- tracksLabel: " << fTracksLabel << endl;
|
| 58 |
cout << "--- muonsLabel: " << fMuonsLabel << endl;
|
| 59 |
cout << "--- muonPt: " << fMuonPt << endl;
|
| 60 |
cout << "--- massLow: " << fMassLow << endl;
|
| 61 |
cout << "--- massHigh: " << fMassHigh << endl;
|
| 62 |
cout << "----------------------------------------------------------------------" << endl;
|
| 63 |
|
| 64 |
}
|
| 65 |
|
| 66 |
|
| 67 |
// ----------------------------------------------------------------------
|
| 68 |
HFDimuons::~HFDimuons() {
|
| 69 |
|
| 70 |
}
|
| 71 |
|
| 72 |
|
| 73 |
// ----------------------------------------------------------------------
|
| 74 |
void HFDimuons::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
|
| 75 |
|
| 76 |
// -- get the primary vertex
|
| 77 |
edm::Handle<VertexCollection> recoPrimaryVertexCollection;
|
| 78 |
iEvent.getByLabel(fPrimaryVertexLabel, recoPrimaryVertexCollection);
|
| 79 |
if(!recoPrimaryVertexCollection.isValid()) {
|
| 80 |
cout << "==>HFDimuons> No primary vertex collection found, skipping" << endl;
|
| 81 |
return;
|
| 82 |
}
|
| 83 |
const reco::VertexCollection vertices = *(recoPrimaryVertexCollection.product());
|
| 84 |
if (vertices.size() == 0) {
|
| 85 |
cout << "==>HFDimuons> No primary vertex found, skipping" << endl;
|
| 86 |
return;
|
| 87 |
}
|
| 88 |
//fPV = vertices[gHFEvent->fEventTag];
|
| 89 |
//if (fVerbose > 0) {
|
| 90 |
// cout << "HFDimuons: Taking vertex " << gHFEvent->fEventTag << " with ntracks = " << fPV.tracksSize() << endl;
|
| 91 |
//}
|
| 92 |
|
| 93 |
// -- get the collection of muons
|
| 94 |
Handle<MuonCollection> hMuons;
|
| 95 |
iEvent.getByLabel(fMuonsLabel, hMuons);
|
| 96 |
if (!hMuons.isValid()) {
|
| 97 |
cout << "==>HFDimuons> No valid MuonCollection with label "<< fMuonsLabel <<" found, skipping" << endl;
|
| 98 |
return;
|
| 99 |
}
|
| 100 |
|
| 101 |
// -- get the collection of tracks
|
| 102 |
edm::Handle<edm::View<reco::Track> > hTracks;
|
| 103 |
iEvent.getByLabel(fTracksLabel, hTracks);
|
| 104 |
if(!hTracks.isValid()) {
|
| 105 |
cout << "==>HFDimuons> No valid TrackCollection with label "<<fTracksLabel <<" found, skipping" << endl;
|
| 106 |
return;
|
| 107 |
}
|
| 108 |
|
| 109 |
// -- get the collection of muons and store their corresponding track indices
|
| 110 |
vector<int> muonIndices;
|
| 111 |
for (MuonCollection::const_iterator muon = hMuons->begin(); muon != hMuons->end(); ++muon) {
|
| 112 |
int im = muon->track().index();
|
| 113 |
if (im > 0) muonIndices.push_back(im);
|
| 114 |
}
|
| 115 |
if (fVerbose > 0) {
|
| 116 |
cout << "==>HFDimuons> nMuons = " << hMuons->size() << endl;
|
| 117 |
cout << "==>HFDimuons> nMuonIndices = " << muonIndices.size() << endl;
|
| 118 |
}
|
| 119 |
if (muonIndices.size() < 2) return;
|
| 120 |
|
| 121 |
// -- Transient tracks for vertexing
|
| 122 |
iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", fTTB);
|
| 123 |
if (!fTTB.isValid()) {
|
| 124 |
cout << " -->HFDimuons: Error: no TransientTrackBuilder found."<<endl;
|
| 125 |
return;
|
| 126 |
}
|
| 127 |
|
| 128 |
std::vector<reco::Track> fitTracks;
|
| 129 |
TLorentzVector dimuon, m1, m2;
|
| 130 |
//int iMuon1(-1), iMuon2(-1);
|
| 131 |
for (unsigned int imuon1 = 0; imuon1 < muonIndices.size()-1; ++imuon1) {
|
| 132 |
TrackBaseRef mu1TrackView(hTracks, muonIndices[imuon1]);
|
| 133 |
Track tMuon1(*mu1TrackView);
|
| 134 |
//iMuon1 = muonIndices[imuon1];
|
| 135 |
m1.SetPtEtaPhiM(tMuon1.pt(), tMuon1.eta(), tMuon1.phi(), MMUON);
|
| 136 |
if (tMuon1.pt() < fMuonPt) continue;
|
| 137 |
|
| 138 |
for (unsigned int imuon2 = imuon1 + 1; imuon2 < muonIndices.size(); ++imuon2) {
|
| 139 |
TrackBaseRef mu2TrackView(hTracks, muonIndices[imuon2]);
|
| 140 |
Track tMuon2(*mu2TrackView);
|
| 141 |
//iMuon2 = muonIndices[imuon2];
|
| 142 |
if (tMuon2.pt() < fMuonPt) continue;
|
| 143 |
|
| 144 |
m2.SetPtEtaPhiM(tMuon2.pt(), tMuon2.eta(), tMuon2.phi(), MMUON);
|
| 145 |
dimuon = m1 + m2;
|
| 146 |
|
| 147 |
if (dimuon.M() < fMassLow || dimuon.M() > fMassHigh) {
|
| 148 |
if (fVerbose > 0) {
|
| 149 |
cout << "==>HFDimuons> dimuon mass = " << dimuon.M() << ", skipping" << endl;
|
| 150 |
}
|
| 151 |
continue;
|
| 152 |
}
|
| 153 |
|
| 154 |
// -- Vertex the two muons only
|
| 155 |
//TAnaCand *pCand = gHFEvent->addCand();
|
| 156 |
//fitTracks.clear();
|
| 157 |
//fitTracks.push_back(tMuon1);
|
| 158 |
//fitTracks.push_back(tMuon2);
|
| 159 |
//doVertexFit(fitTracks, iMuon1, iMuon2, pCand);
|
| 160 |
|
| 161 |
}
|
| 162 |
}
|
| 163 |
}
|
| 164 |
|
| 165 |
/*
|
| 166 |
// ----------------------------------------------------------------------
|
| 167 |
void HFDimuons::doVertexFit(std::vector<reco::Track> &Tracks, int iMuon1, int iMuon2, TAnaCand *pCand){
|
| 168 |
|
| 169 |
Track tMuon1 = Tracks[0];
|
| 170 |
Track tMuon2 = Tracks[1];
|
| 171 |
|
| 172 |
std::vector<reco::TransientTrack> RecoTransientTrack;
|
| 173 |
RecoTransientTrack.clear();
|
| 174 |
RecoTransientTrack.push_back(fTTB->build(Tracks[0]));
|
| 175 |
RecoTransientTrack.push_back(fTTB->build(Tracks[1]));
|
| 176 |
|
| 177 |
// -- Do the vertexing
|
| 178 |
KalmanVertexFitter theFitter(true);
|
| 179 |
TransientVertex TransSecVtx = theFitter.vertex(RecoTransientTrack);
|
| 180 |
if (TransSecVtx.isValid() ) {
|
| 181 |
if (isnan(TransSecVtx.position().x())
|
| 182 |
|| isnan(TransSecVtx.position().y())
|
| 183 |
|| isnan(TransSecVtx.position().z()) ) {
|
| 184 |
cout << "==>HFDimuons> Something went wrong! SecVtx nan - continue ... " << endl;
|
| 185 |
pCand->fType = -1;
|
| 186 |
return;
|
| 187 |
}
|
| 188 |
} else {
|
| 189 |
cout << "==>HFDimuons> KVF failed! continue ..." << endl;
|
| 190 |
pCand->fType = -1;
|
| 191 |
return;
|
| 192 |
}
|
| 193 |
|
| 194 |
// -- Get refitted tracks
|
| 195 |
std::vector<reco::TransientTrack> refTT = TransSecVtx.refittedTracks();
|
| 196 |
std::vector<reco::Track> refT; refT.clear();
|
| 197 |
for(vector<reco::TransientTrack>::const_iterator i = refTT.begin(); i != refTT.end(); i++) {
|
| 198 |
const Track &ftt = i->track();
|
| 199 |
refT.push_back(ftt);
|
| 200 |
}
|
| 201 |
|
| 202 |
// -- Build composite
|
| 203 |
TLorentzVector comp, M1, M2;
|
| 204 |
M1.SetXYZM(refT[0].px(), refT[0].py(), refT[0].pz(), MMUON);
|
| 205 |
M2.SetXYZM(refT[1].px(), refT[1].py(), refT[1].pz(), MMUON);
|
| 206 |
comp = M1 + M2;
|
| 207 |
|
| 208 |
|
| 209 |
// -- Build vertex for ntuple
|
| 210 |
TAnaVertex anaVtx;
|
| 211 |
ChiSquared chi(TransSecVtx.totalChiSquared(), TransSecVtx.degreesOfFreedom());
|
| 212 |
anaVtx.setInfo(chi.value(), chi.degreesOfFreedom(), chi.probability(), 0, 0);
|
| 213 |
anaVtx.fPoint.SetXYZ(TransSecVtx.position().x(),
|
| 214 |
TransSecVtx.position().y(),
|
| 215 |
TransSecVtx.position().z());
|
| 216 |
|
| 217 |
anaVtx.addTrack(iMuon1);
|
| 218 |
anaVtx.addTrack(iMuon2);
|
| 219 |
|
| 220 |
VertexDistanceXY axy;
|
| 221 |
anaVtx.fDxy = axy.distance(fPV, TransSecVtx).value();
|
| 222 |
anaVtx.fDxyE = axy.distance(fPV, TransSecVtx).error();
|
| 223 |
VertexDistance3D a3d;
|
| 224 |
anaVtx.fD3d = a3d.distance(fPV, TransSecVtx).value();
|
| 225 |
anaVtx.fD3dE = a3d.distance(fPV, TransSecVtx).error();
|
| 226 |
|
| 227 |
|
| 228 |
// -- fill candidate
|
| 229 |
pCand->fPlab = comp.Vect();
|
| 230 |
pCand->fMass = comp.M();
|
| 231 |
pCand->fVtx = anaVtx;
|
| 232 |
pCand->fType = fType;
|
| 233 |
pCand->fDau1 = -1;
|
| 234 |
pCand->fDau2 = -1;
|
| 235 |
pCand->fSig1 = gHFEvent->nSigTracks();
|
| 236 |
pCand->fSig2 = pCand->fSig1 + 1;
|
| 237 |
|
| 238 |
// -- fill refitted sig tracks
|
| 239 |
TAnaTrack *pTrack = gHFEvent->addSigTrack();
|
| 240 |
pTrack->fMCID = tMuon1.charge()*13;
|
| 241 |
pTrack->fGenIndex = -1;
|
| 242 |
pTrack->fQ = tMuon1.charge();
|
| 243 |
pTrack->fPlab.SetXYZ(refT[0].px(),
|
| 244 |
refT[0].py(),
|
| 245 |
refT[0].pz()
|
| 246 |
);
|
| 247 |
pTrack->fIndex = iMuon1;
|
| 248 |
|
| 249 |
pTrack = gHFEvent->addSigTrack();
|
| 250 |
pTrack->fMCID = tMuon2.charge()*13;
|
| 251 |
pTrack->fGenIndex = -1;
|
| 252 |
pTrack->fQ = tMuon2.charge();
|
| 253 |
pTrack->fPlab.SetXYZ(refT[1].px(),
|
| 254 |
refT[1].py(),
|
| 255 |
refT[1].pz()
|
| 256 |
);
|
| 257 |
pTrack->fIndex = iMuon2;
|
| 258 |
}
|
| 259 |
*/
|
| 260 |
|
| 261 |
// ------------ method called once each job just before starting event loop ------------
|
| 262 |
void HFDimuons::beginJob(const edm::EventSetup& setup) {
|
| 263 |
}
|
| 264 |
|
| 265 |
// ------------ method called once each job just after ending the event loop ------------
|
| 266 |
void HFDimuons::endJob() {
|
| 267 |
}
|
| 268 |
|
| 269 |
//define this as a plug-in
|
| 270 |
//DEFINE_FWK_MODULE(HFDimuons);
|
