| 1 |
bendavid |
1.1 |
// $Id: ProducerConversionsKinematic.cc,v 1.23 2010/01/18 14:41:33 bendavid Exp $
|
| 2 |
|
|
|
| 3 |
|
|
#include "MitEdm/Producers/interface/ProducerConversionsKinematic.h"
|
| 4 |
|
|
#include "DataFormats/Common/interface/Handle.h"
|
| 5 |
|
|
#include "DataFormats/TrackReco/interface/Track.h"
|
| 6 |
|
|
#include "DataFormats/TrackReco/interface/TrackFwd.h"
|
| 7 |
|
|
#include "DataFormats/VertexReco/interface/Vertex.h"
|
| 8 |
|
|
#include "DataFormats/VertexReco/interface/VertexFwd.h"
|
| 9 |
|
|
#include "MagneticField/Engine/interface/MagneticField.h"
|
| 10 |
|
|
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
|
| 11 |
|
|
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
|
| 12 |
|
|
#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
|
| 13 |
|
|
#include "MitEdm/Producers/interface/HitDropperRecord.h"
|
| 14 |
|
|
#include "MitEdm/Producers/interface/HitDropper.h"
|
| 15 |
|
|
#include "MitEdm/DataFormats/interface/Types.h"
|
| 16 |
|
|
#include "MitEdm/DataFormats/interface/Collections.h"
|
| 17 |
|
|
#include "MitEdm/DataFormats/interface/DecayPart.h"
|
| 18 |
|
|
#include "MitEdm/DataFormats/interface/StablePart.h"
|
| 19 |
|
|
#include "MitEdm/DataFormats/interface/StableData.h"
|
| 20 |
|
|
#include "MitEdm/VertexFitInterface/interface/MvfInterface.h"
|
| 21 |
|
|
#include "MitEdm/VertexFitInterface/interface/TrackParameters.h"
|
| 22 |
|
|
#include <TMath.h>
|
| 23 |
|
|
|
| 24 |
|
|
//#include "DataFormats/GeometrySurface/interface/ReferenceCounted.h"
|
| 25 |
|
|
//#include "RecoVertex/VertexPrimitives/interface/VertexState.h"
|
| 26 |
|
|
//#include "DataFormats/VertexReco/interface/Vertex.h"
|
| 27 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/KinematicVertex.h"
|
| 28 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticle.h"
|
| 29 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/KinematicTree.h"
|
| 30 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/TransientTrackKinematicParticle.h"
|
| 31 |
|
|
#include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
|
| 32 |
|
|
|
| 33 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/ParticleMass.h"
|
| 34 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
|
| 35 |
|
|
#include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h"
|
| 36 |
|
|
#include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
|
| 37 |
|
|
#include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
|
| 38 |
|
|
#include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
|
| 39 |
|
|
#include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
|
| 40 |
|
|
#include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
|
| 41 |
|
|
#include "RecoVertex/KinematicFit/interface/ColinearityKinematicConstraint.h"
|
| 42 |
|
|
#include "RecoVertex/KinematicFit/interface/MultiTrackMassKinematicConstraint.h"
|
| 43 |
|
|
#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
|
| 44 |
|
|
#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
|
| 45 |
|
|
#include "DataFormats/Math/interface/deltaPhi.h"
|
| 46 |
|
|
|
| 47 |
|
|
|
| 48 |
|
|
|
| 49 |
|
|
using namespace std;
|
| 50 |
|
|
using namespace edm;
|
| 51 |
|
|
using namespace reco;
|
| 52 |
|
|
using namespace mitedm;
|
| 53 |
|
|
using namespace mithep;
|
| 54 |
|
|
|
| 55 |
|
|
//--------------------------------------------------------------------------------------------------
|
| 56 |
|
|
ProducerConversionsKinematic::ProducerConversionsKinematic(const ParameterSet& cfg) :
|
| 57 |
|
|
BaseCandProducer (cfg),
|
| 58 |
|
|
iStables1_ (cfg.getUntrackedParameter<string>("iStables1","")),
|
| 59 |
|
|
iStables2_ (cfg.getUntrackedParameter<string>("iStables2","")),
|
| 60 |
|
|
iPVertexes_ (cfg.getUntrackedParameter<string>("iPVertexes","offlinePrimaryVerticesWithBS")),
|
| 61 |
|
|
usePVertex_ (cfg.getUntrackedParameter<bool> ("usePVertex",true)),
|
| 62 |
|
|
convConstraint_ (cfg.getUntrackedParameter<bool> ("convConstraint",false)),
|
| 63 |
|
|
convConstraint3D_(cfg.getUntrackedParameter<bool> ("convConstraint3D",true)),
|
| 64 |
|
|
rhoMin_ (cfg.getUntrackedParameter<double>("rhoMin",0.0)),
|
| 65 |
|
|
useRhoMin_ (cfg.getUntrackedParameter<bool> ("useRhoMin",true)),
|
| 66 |
|
|
useHitDropper_ (cfg.getUntrackedParameter<bool> ("useHitDropper",true)),
|
| 67 |
|
|
applyChargeConstraint_(cfg.getUntrackedParameter<bool> ("applyChargeConstraint",false))
|
| 68 |
|
|
{
|
| 69 |
|
|
// Constructor.
|
| 70 |
|
|
|
| 71 |
|
|
produces<DecayPartCol>();
|
| 72 |
|
|
}
|
| 73 |
|
|
|
| 74 |
|
|
//--------------------------------------------------------------------------------------------------
|
| 75 |
|
|
ProducerConversionsKinematic::~ProducerConversionsKinematic()
|
| 76 |
|
|
{
|
| 77 |
|
|
// Destructor.
|
| 78 |
|
|
}
|
| 79 |
|
|
|
| 80 |
|
|
//--------------------------------------------------------------------------------------------------
|
| 81 |
|
|
void ProducerConversionsKinematic::produce(Event &evt, const EventSetup &setup)
|
| 82 |
|
|
{
|
| 83 |
|
|
// Produce our DecayPartCol.
|
| 84 |
|
|
|
| 85 |
|
|
// First input collection
|
| 86 |
|
|
Handle<StablePartCol> hStables1;
|
| 87 |
|
|
if (!GetProduct(iStables1_, hStables1, evt)) {
|
| 88 |
|
|
printf("Stable collection 1 not found in ProducerConversionsKinematic\n");
|
| 89 |
|
|
return;
|
| 90 |
|
|
}
|
| 91 |
|
|
const StablePartCol *pS1 = hStables1.product();
|
| 92 |
|
|
// Second input collection
|
| 93 |
|
|
Handle<StablePartCol> hStables2;
|
| 94 |
|
|
if (!GetProduct(iStables2_, hStables2, evt)) {
|
| 95 |
|
|
printf("Stable collection 2 not found in ProducerConversionsKinematic\n");
|
| 96 |
|
|
return;
|
| 97 |
|
|
}
|
| 98 |
|
|
const StablePartCol *pS2 = hStables2.product();
|
| 99 |
|
|
|
| 100 |
|
|
// Get hit dropper
|
| 101 |
|
|
ESHandle<HitDropper> hDropper;
|
| 102 |
|
|
setup.get<HitDropperRecord>().get("HitDropper",hDropper);
|
| 103 |
|
|
const HitDropper *dropper = hDropper.product();
|
| 104 |
|
|
|
| 105 |
|
|
//Get Magnetic Field from event setup, taking value at (0,0,0)
|
| 106 |
|
|
edm::ESHandle<MagneticField> magneticField;
|
| 107 |
|
|
setup.get<IdealMagneticFieldRecord>().get(magneticField);
|
| 108 |
|
|
const double bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
|
| 109 |
|
|
|
| 110 |
|
|
edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
|
| 111 |
|
|
setup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
|
| 112 |
|
|
const TransientTrackBuilder *transientTrackBuilder = hTransientTrackBuilder.product();
|
| 113 |
|
|
|
| 114 |
|
|
//construct intermediate collection of TrackParameters in mvf format for vertex fit
|
| 115 |
|
|
std::vector<reco::TransientTrack> ttrks1;
|
| 116 |
|
|
for (UInt_t i = 0; i<pS1->size(); ++i) {
|
| 117 |
|
|
const reco::Track *t = pS1->at(i).track();
|
| 118 |
|
|
const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
|
| 119 |
|
|
ttrks1.push_back(ttrk);
|
| 120 |
|
|
}
|
| 121 |
|
|
|
| 122 |
|
|
std::vector<reco::TransientTrack> ttrks2;
|
| 123 |
|
|
if (iStables1_ == iStables2_) {
|
| 124 |
|
|
ttrks2 = ttrks1;
|
| 125 |
|
|
}
|
| 126 |
|
|
else for (UInt_t i = 0; i<pS2->size(); ++i) {
|
| 127 |
|
|
const reco::Track *t = pS2->at(i).track();
|
| 128 |
|
|
const reco::TransientTrack ttrk = transientTrackBuilder->build(t);
|
| 129 |
|
|
ttrks2.push_back(ttrk);
|
| 130 |
|
|
}
|
| 131 |
|
|
|
| 132 |
|
|
|
| 133 |
|
|
// Create the output collection
|
| 134 |
|
|
auto_ptr<DecayPartCol> pD(new DecayPartCol());
|
| 135 |
|
|
|
| 136 |
|
|
ClosestApproachInRPhi helixIntersector;
|
| 137 |
|
|
|
| 138 |
|
|
int nFits = 0;
|
| 139 |
|
|
|
| 140 |
|
|
//printf("S1 size = %i\n", pS1->size());
|
| 141 |
|
|
|
| 142 |
|
|
// Simple double loop
|
| 143 |
|
|
for (UInt_t i = 0; i<pS1->size(); ++i) {
|
| 144 |
|
|
const StablePart &s1 = pS1->at(i);
|
| 145 |
|
|
|
| 146 |
|
|
const reco::Track * t1 = s1.track();
|
| 147 |
|
|
// const TrackParameters &trkPar1 = trkPars1.at(i);
|
| 148 |
|
|
const reco::TransientTrack &ttrk1 = ttrks1.at(i);
|
| 149 |
|
|
|
| 150 |
|
|
UInt_t j;
|
| 151 |
|
|
if (iStables1_ == iStables2_)
|
| 152 |
|
|
j = i+1;
|
| 153 |
|
|
else
|
| 154 |
|
|
j = 0;
|
| 155 |
|
|
|
| 156 |
|
|
TrajectoryStateTransform tsTransform;
|
| 157 |
|
|
|
| 158 |
|
|
FreeTrajectoryState initialState1 = tsTransform.initialFreeState(*s1.track(),&*magneticField);
|
| 159 |
|
|
|
| 160 |
|
|
for (; j<pS2->size(); ++j) {
|
| 161 |
|
|
const StablePart &s2 = pS2->at(j);
|
| 162 |
|
|
|
| 163 |
|
|
//Do fast helix fit to check if there's any hope
|
| 164 |
|
|
const reco::Track * t2 = s2.track();
|
| 165 |
|
|
// const TrackParameters &trkPar2 = trkPars2.at(j);
|
| 166 |
|
|
const reco::TransientTrack &ttrk2 = ttrks2.at(j);
|
| 167 |
|
|
|
| 168 |
|
|
int trackCharge = t1->charge() + t2->charge();
|
| 169 |
|
|
|
| 170 |
|
|
double dR0 = 0.0;
|
| 171 |
|
|
double dZ0 = 1e6;
|
| 172 |
|
|
double angle1 = 0.0;
|
| 173 |
|
|
double angle2 = 0.0;
|
| 174 |
|
|
double angle1simple = 0.0;
|
| 175 |
|
|
double angle2simple = 0.0;
|
| 176 |
|
|
double phi1helix = 0.0;
|
| 177 |
|
|
double phi2helix = 0.0;
|
| 178 |
|
|
double eta1helix = 0.0;
|
| 179 |
|
|
double eta2helix = 0.0;
|
| 180 |
|
|
if (!applyChargeConstraint_ || trackCharge==0) {
|
| 181 |
|
|
FreeTrajectoryState initialState2 = tsTransform.initialFreeState(*s2.track(),&*magneticField);
|
| 182 |
|
|
helixIntersector.calculate(initialState1, initialState2);
|
| 183 |
|
|
if (helixIntersector.status()) {
|
| 184 |
|
|
dR0 = helixIntersector.crossingPoint().perp();
|
| 185 |
|
|
//std::pair<GlobalPoint, GlobalPoint> intPoints = helixIntersector.points();
|
| 186 |
|
|
//dZ0 = intPoints.first.z() - intPoints.second.z();
|
| 187 |
|
|
|
| 188 |
|
|
std::pair <GlobalTrajectoryParameters, GlobalTrajectoryParameters > trajpars = helixIntersector.trajectoryParameters();
|
| 189 |
|
|
dZ0 = trajpars.first.position().z() - trajpars.second.position().z();
|
| 190 |
|
|
angle1simple = trajpars.first.momentum().phi() - t1->phi();
|
| 191 |
|
|
angle2simple = trajpars.second.momentum().phi() - t2->phi();
|
| 192 |
|
|
phi1helix = trajpars.first.momentum().phi();
|
| 193 |
|
|
phi2helix = trajpars.second.momentum().phi();
|
| 194 |
|
|
eta1helix = trajpars.first.momentum().eta();
|
| 195 |
|
|
eta2helix = trajpars.second.momentum().eta();
|
| 196 |
|
|
|
| 197 |
|
|
TwoTrackMinimumDistance trackmin;
|
| 198 |
|
|
trackmin.calculate(initialState1,initialState2);
|
| 199 |
|
|
if (trackmin.status()) {
|
| 200 |
|
|
angle1 = trackmin.firstAngle();
|
| 201 |
|
|
angle2 = trackmin.secondAngle();
|
| 202 |
|
|
}
|
| 203 |
|
|
//initialFitPoint = helixIntersector.crossingPoint();
|
| 204 |
|
|
}
|
| 205 |
|
|
}
|
| 206 |
|
|
|
| 207 |
|
|
double prob = -99.0;
|
| 208 |
|
|
int fitStatus = 0;
|
| 209 |
|
|
RefCountedKinematicTree myTree;
|
| 210 |
|
|
RefCountedKinematicVertex gamma_dec_vertex;
|
| 211 |
|
|
RefCountedKinematicParticle the_photon;
|
| 212 |
|
|
if ( (!applyChargeConstraint_ || trackCharge==0) && (!useRhoMin_ || dR0 > rhoMin_) ) {
|
| 213 |
|
|
|
| 214 |
|
|
// Vertex fit now, possibly with conversion constraint
|
| 215 |
|
|
nFits++;
|
| 216 |
|
|
|
| 217 |
|
|
|
| 218 |
|
|
|
| 219 |
|
|
|
| 220 |
|
|
|
| 221 |
|
|
|
| 222 |
|
|
|
| 223 |
|
|
bool cmsFitStatus = false;
|
| 224 |
|
|
|
| 225 |
|
|
float sigma = 0.00000000001;
|
| 226 |
|
|
float chi = 0.;
|
| 227 |
|
|
float ndf = 0.;
|
| 228 |
|
|
|
| 229 |
|
|
edm::ParameterSet pSet;
|
| 230 |
|
|
pSet.addParameter<double>("maxDistance", 1e-3);//0.001
|
| 231 |
|
|
pSet.addParameter<double>("maxOfInitialValue",1.4) ;//1.4
|
| 232 |
|
|
pSet.addParameter<int>("maxNbrOfIterations", 40);//40
|
| 233 |
|
|
|
| 234 |
|
|
KinematicParticleFactoryFromTransientTrack pFactory;
|
| 235 |
|
|
|
| 236 |
|
|
vector<RefCountedKinematicParticle> particles;
|
| 237 |
|
|
|
| 238 |
|
|
|
| 239 |
|
|
|
| 240 |
|
|
particles.push_back(pFactory.particle (ttrk1,s1.mass(),chi,ndf,sigma));
|
| 241 |
|
|
particles.push_back(pFactory.particle (ttrk2,s2.mass(),chi,ndf,sigma));
|
| 242 |
|
|
|
| 243 |
|
|
MultiTrackKinematicConstraint *constr = new ColinearityKinematicConstraint(ColinearityKinematicConstraint::PhiTheta);
|
| 244 |
|
|
//MultiTrackKinematicConstraint *constr = new ColinearityKinematicConstraint(ColinearityKinematicConstraint::Phi);
|
| 245 |
|
|
|
| 246 |
|
|
// ParticleMass pmass(kmass);
|
| 247 |
|
|
//MultiTrackKinematicConstraint *constr = new TwoTrackMassKinematicConstraint(pmass);
|
| 248 |
|
|
|
| 249 |
|
|
KinematicConstrainedVertexFitter kcvFitter;
|
| 250 |
|
|
kcvFitter.setParameters(pSet);
|
| 251 |
|
|
//RefCountedKinematicTree myTree = kcvFitter.fit(particles, constr);
|
| 252 |
|
|
//RefCountedKinematicTree myTree = kcvFitter.fit(particles, constr, &initialFitPoint);
|
| 253 |
|
|
//RefCountedKinematicTree myTree = kcvFitter.fit(particles);
|
| 254 |
|
|
|
| 255 |
|
|
|
| 256 |
|
|
|
| 257 |
|
|
|
| 258 |
|
|
//if (cdfProb>1e-6) {
|
| 259 |
|
|
myTree = kcvFitter.fit(particles, constr);
|
| 260 |
|
|
//}
|
| 261 |
|
|
|
| 262 |
|
|
|
| 263 |
|
|
|
| 264 |
|
|
|
| 265 |
|
|
|
| 266 |
|
|
if( myTree->isValid() ) {
|
| 267 |
|
|
myTree->movePointerToTheTop();
|
| 268 |
|
|
the_photon = myTree->currentParticle();
|
| 269 |
|
|
if (the_photon->currentState().isValid()){
|
| 270 |
|
|
//const ParticleMass photon_mass = the_photon->currentState().mass();
|
| 271 |
|
|
gamma_dec_vertex = myTree->currentDecayVertex();
|
| 272 |
|
|
if( gamma_dec_vertex->vertexIsValid() ){
|
| 273 |
|
|
cmsFitStatus = true;
|
| 274 |
|
|
//const float chi2Prob = ChiSquaredProbability(gamma_dec_vertex->chiSquared(), gamma_dec_vertex->degreesOfFreedom());
|
| 275 |
|
|
double cmsChi2 = gamma_dec_vertex->chiSquared();
|
| 276 |
|
|
double cmsNdof = gamma_dec_vertex->degreesOfFreedom();
|
| 277 |
|
|
//cmsR = gamma_dec_vertex->position().perp();
|
| 278 |
|
|
//cmsPhi = gamma_dec_vertex->position().phi();
|
| 279 |
|
|
prob = TMath::Prob(cmsChi2,cmsNdof);
|
| 280 |
|
|
}
|
| 281 |
|
|
}
|
| 282 |
|
|
|
| 283 |
|
|
|
| 284 |
|
|
|
| 285 |
|
|
}
|
| 286 |
|
|
|
| 287 |
|
|
|
| 288 |
|
|
|
| 289 |
|
|
|
| 290 |
|
|
}
|
| 291 |
|
|
|
| 292 |
|
|
if (prob>1e-10) {
|
| 293 |
|
|
DecayPart *d = new DecayPart(oPid_,DecayPart::Fast);
|
| 294 |
|
|
|
| 295 |
|
|
// Update temporarily some of the quantities (prob, chi2, nDoF, mass, lxy, pt, fourMomentum)
|
| 296 |
|
|
d->setProb(prob);
|
| 297 |
|
|
d->setChi2(gamma_dec_vertex->chiSquared());
|
| 298 |
|
|
d->setNdof(gamma_dec_vertex->degreesOfFreedom());
|
| 299 |
|
|
|
| 300 |
|
|
GlobalVector momphoton = the_photon->currentState().globalMomentum();
|
| 301 |
|
|
ThreeVector gammaMom(momphoton);
|
| 302 |
|
|
double gammamass = the_photon->currentState().mass();
|
| 303 |
|
|
double photonenergy = sqrt(gammaMom.R()*gammaMom.R() + gammamass*gammamass);
|
| 304 |
|
|
|
| 305 |
|
|
FourVector p4Fitted(gammaMom.x(),gammaMom.y(),gammaMom.z(),photonenergy);
|
| 306 |
|
|
//p4Fitted += fit.getTrackP4(1);
|
| 307 |
|
|
//p4Fitted += fit.getTrackP4(2);
|
| 308 |
|
|
d->setFourMomentum(p4Fitted);
|
| 309 |
|
|
ThreeVector vtxPos(gamma_dec_vertex->position());
|
| 310 |
|
|
d->setPosition(vtxPos);
|
| 311 |
|
|
//d->setError (fit.getErrorMatrix(MultiVertexFitterD::VERTEX_1));
|
| 312 |
|
|
//float mass, massErr;
|
| 313 |
|
|
//const int trksIds[2] = { 1, 2 };
|
| 314 |
|
|
//mass = fit.getMass(2,trksIds,massErr);
|
| 315 |
|
|
|
| 316 |
|
|
|
| 317 |
|
|
ThreeVector p3Fitted(p4Fitted.px(), p4Fitted.py(), p4Fitted.pz());
|
| 318 |
|
|
|
| 319 |
|
|
// Get decay length in xy plane
|
| 320 |
|
|
// float dl, dlErr;
|
| 321 |
|
|
// dl = fit.getDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
|
| 322 |
|
|
// p3Fitted, dlErr);
|
| 323 |
|
|
//
|
| 324 |
|
|
// // Get Z decay length
|
| 325 |
|
|
// float dlz, dlzErr;
|
| 326 |
|
|
// dlz = fit.getZDecayLength(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
|
| 327 |
|
|
// p3Fitted, dlzErr);
|
| 328 |
|
|
//
|
| 329 |
|
|
// // Get impact parameter
|
| 330 |
|
|
// float dxy, dxyErr;
|
| 331 |
|
|
// dxy = fit.getImpactPar(MultiVertexFitterD::PRIMARY_VERTEX, MultiVertexFitterD::VERTEX_1,
|
| 332 |
|
|
// p3Fitted, dxyErr);
|
| 333 |
|
|
|
| 334 |
|
|
double dlz = vtxPos.z()*TMath::Abs(p4Fitted.Pz())/p4Fitted.Pz();
|
| 335 |
|
|
|
| 336 |
|
|
ThreeVector momPerp(p4Fitted.Px(),p4Fitted.Py(),0);
|
| 337 |
|
|
ThreeVector posPerp(vtxPos.x(),vtxPos.y(),0);
|
| 338 |
|
|
double dl = momPerp.Dot(posPerp)/momPerp.R();
|
| 339 |
|
|
|
| 340 |
|
|
|
| 341 |
|
|
double dxy = momPerp.Cross(vtxPos).Z()/p4Fitted.Pt();
|
| 342 |
|
|
|
| 343 |
|
|
|
| 344 |
|
|
|
| 345 |
|
|
BasePartPtr ptr1(hStables1,i);
|
| 346 |
|
|
BasePartPtr ptr2(hStables2,j);
|
| 347 |
|
|
|
| 348 |
|
|
std::vector< RefCountedKinematicParticle > bs_children = myTree->finalStateParticles();
|
| 349 |
|
|
|
| 350 |
|
|
GlobalVector mom1 = bs_children.at(0)->currentState().globalMomentum();
|
| 351 |
|
|
GlobalVector mom2 = bs_children.at(1)->currentState().globalMomentum();
|
| 352 |
|
|
|
| 353 |
|
|
const ThreeVector trkMom1(mom1.x(),mom1.y(),mom1.z());
|
| 354 |
|
|
const ThreeVector trkMom2(mom2.x(),mom2.y(),mom2.z());
|
| 355 |
|
|
|
| 356 |
|
|
double deltaphi = reco::deltaPhi(trkMom1.Phi(),trkMom2.Phi());
|
| 357 |
|
|
double etaratio = trkMom1.Eta()/trkMom2.Eta();
|
| 358 |
|
|
|
| 359 |
|
|
//if (prob>1e-6 && TMath::Abs(p4Fitted.Eta()) <0.01) {
|
| 360 |
|
|
if (prob>1e-6 && deltaphi > 0.1) {
|
| 361 |
|
|
printf("kin track 1 : qoverp = %5f, pt = %5f, theta = %5f, eta = %5f, phi = %5f, dxy = %5f, dsz = %5f, dz = %5f\n",t1->qoverp(),t1->pt(),t1->theta(),t1->eta(),t1->phi(), t1->dxy(), t1->dsz(), t1->dz());
|
| 362 |
|
|
printf("kin track 1 : qoverperr = %5f, pterr = %5f, thetaerr = %5f, etaerr = %5f, phierr = %5f, dxyerr = %5f, dszerr = %5f, dzerr = %5f\n",t1->qoverpError(),t1->ptError(),t1->thetaError(),t1->etaError(),t1->phiError(), t1->dxyError(), t1->dszError(), t1->dzError());
|
| 363 |
|
|
printf("kin track 2 : qoverp = %5f, pt = %5f, theta = %5f, eta = %5f, phi = %5f, dxy = %5f, dsz = %5f, dz = %5f\n",t2->qoverp(),t2->pt(),t2->theta(),t2->eta(),t2->phi(), t2->dxy(), t2->dsz(), t2->dz());
|
| 364 |
|
|
printf("kin track 2 : qoverperr = %5f, pterr = %5f, thetaerr = %5f, etaerr = %5f, phierr = %5f, dxyerr = %5f, dszerr = %5f, dzerr = %5f\n",t2->qoverpError(),t2->ptError(),t2->thetaError(),t2->etaError(),t2->phiError(), t2->dxyError(), t2->dszError(), t2->dzError());
|
| 365 |
|
|
printf("kin conversion pt = %5f, eta = %5f, theta = %5f, phi = %5f\n",p4Fitted.Pt(),p4Fitted.Eta(),p4Fitted.Theta(),p4Fitted.Phi());
|
| 366 |
|
|
printf("kin trkMom1 pt = %5f, eta = %5f, theta = %5f, phi = %5f\n",trkMom1.Rho(),trkMom1.Eta(),trkMom1.Theta(),trkMom1.Phi());
|
| 367 |
|
|
printf("kin trkMom2 pt = %5f, eta = %5f, theta = %5f, phi = %5f\n",trkMom2.Rho(),trkMom2.Eta(),trkMom2.Theta(),trkMom2.Phi());
|
| 368 |
|
|
printf("kin helix1 eta = %5f, phi = %5f\n",phi1helix,eta1helix);
|
| 369 |
|
|
printf("kin helix2 eta = %5f, phi = %5f\n",phi2helix,eta2helix);
|
| 370 |
|
|
printf("kin dR0 = %5f, dZ0 = %5f, firstAngle = %5f, secondAngle = %5f, angle1simple = %5f, angle2simple = %5f, deltaphi = %5f, etaratio = %5f\n",dR0, dZ0,angle1,angle2,angle1simple,angle2simple,deltaphi,etaratio);
|
| 371 |
|
|
}
|
| 372 |
|
|
|
| 373 |
|
|
|
| 374 |
|
|
StableData c1(mom1.x(),mom1.y(), mom1.z(), ptr1);
|
| 375 |
|
|
StableData c2(mom2.x(),mom2.y(), mom2.z(), ptr2);
|
| 376 |
|
|
|
| 377 |
|
|
|
| 378 |
|
|
// Build corrected HitPattern for StableData, removing hits before the fit vertex
|
| 379 |
|
|
if (useHitDropper_) {
|
| 380 |
|
|
std::pair<reco::HitPattern,uint> hits1 = dropper->CorrectedHitsAOD(s1.track(), vtxPos, trkMom1, 1.0, 1.0);
|
| 381 |
|
|
std::pair<reco::HitPattern,uint> hits2 = dropper->CorrectedHitsAOD(s2.track(), vtxPos, trkMom2, 1.0, 1.0);
|
| 382 |
|
|
|
| 383 |
|
|
c1.SetHits(hits1.first);
|
| 384 |
|
|
c2.SetHits(hits2.first);
|
| 385 |
|
|
c1.SetHitsFilled();
|
| 386 |
|
|
c2.SetHitsFilled();
|
| 387 |
|
|
c1.SetNWrongHits(hits1.second);
|
| 388 |
|
|
c2.SetNWrongHits(hits2.second);
|
| 389 |
|
|
|
| 390 |
|
|
reco::HitPattern sharedHits = dropper->SharedHits(s1.track(),s2.track());
|
| 391 |
|
|
d->setSharedHits(sharedHits);
|
| 392 |
|
|
}
|
| 393 |
|
|
|
| 394 |
|
|
d->addStableChild(c1);
|
| 395 |
|
|
d->addStableChild(c2);
|
| 396 |
|
|
|
| 397 |
|
|
|
| 398 |
|
|
// d->setFittedMass (mass);
|
| 399 |
|
|
// d->setFittedMassError(massErr);
|
| 400 |
|
|
|
| 401 |
|
|
d->setLxy(dl);
|
| 402 |
|
|
// d->setLxyError(dlErr);
|
| 403 |
|
|
d->setLxyToPv(dl);
|
| 404 |
|
|
// d->setLxyToPvError(dlErr);
|
| 405 |
|
|
|
| 406 |
|
|
d->setLz(dlz);
|
| 407 |
|
|
// d->setLzError(dlzErr);
|
| 408 |
|
|
d->setLzToPv(dlz);
|
| 409 |
|
|
// d->setLzToPvError(dlzErr);
|
| 410 |
|
|
|
| 411 |
|
|
d->setDxy(dxy);
|
| 412 |
|
|
// d->setDxyError(dxyErr);
|
| 413 |
|
|
d->setDxyToPv(dxy);
|
| 414 |
|
|
// d->setDxyToPvError(dxyErr);
|
| 415 |
|
|
|
| 416 |
|
|
// if (usePVertex_)
|
| 417 |
|
|
// d->setPrimaryVertex(vPtr);
|
| 418 |
|
|
|
| 419 |
|
|
// Put the result into our collection
|
| 420 |
|
|
pD->push_back(*d);
|
| 421 |
|
|
|
| 422 |
|
|
delete d;
|
| 423 |
|
|
}
|
| 424 |
|
|
}
|
| 425 |
|
|
}
|
| 426 |
|
|
|
| 427 |
|
|
//printf("nConversionFits = %i\n",nFits);
|
| 428 |
|
|
|
| 429 |
|
|
// Write the collection even if it is empty
|
| 430 |
|
|
if (0) {
|
| 431 |
|
|
cout << " ProducerConversionsKinematic::produce - " << pD->size() << " entries collection created -"
|
| 432 |
|
|
<< " (Pid: " << oPid_ << ")\n";
|
| 433 |
|
|
}
|
| 434 |
|
|
evt.put(pD);
|
| 435 |
|
|
}
|
| 436 |
|
|
|
| 437 |
|
|
//define this as a plug-in
|
| 438 |
|
|
DEFINE_FWK_MODULE(ProducerConversionsKinematic);
|
