00001 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionSeedFinder.h"
00002 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
00003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
00004
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006
00007 #include "MagneticField/Engine/interface/MagneticField.h"
00008
00009 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00010
00011 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00013
00014 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00015 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00016 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00017 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
00018
00019 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
00020 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00021
00022
00023 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00024 #include "CLHEP/Geometry/Point3D.h"
00025 #include "CLHEP/Geometry/Vector3D.h"
00026 #include "CLHEP/Geometry/Transform3D.h"
00027 #include <cfloat>
00028
00029
00030 OutInConversionSeedFinder::OutInConversionSeedFinder( const edm::ParameterSet& conf ): ConversionSeedFinder( conf ), conf_(conf)
00031 {
00032
00033 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR " << "\n";
00034
00035
00036 the2ndHitdphi_ = 0.03;
00037 the2ndHitdzConst_ = 5.;
00038 the2ndHitdznSigma_ = 2.;
00039 maxNumberOfOutInSeedsPerBC_=50;
00040
00041
00042 }
00043
00044
00045
00046
00047 OutInConversionSeedFinder::~OutInConversionSeedFinder() {
00048 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR " << "\n";
00049
00050 }
00051
00052
00053
00054 void OutInConversionSeedFinder::makeSeeds( const edm::Handle<edm::View<reco::CaloCluster> > & allBC ) const {
00055
00056
00057 theSeeds_.clear();
00058
00059
00060
00061
00062
00063
00064
00065
00066 findLayers();
00067
00068
00069
00070
00071 float theSCPhi=theSCPosition_.phi();
00072 float theSCEta=theSCPosition_.eta();
00073
00074
00075
00076
00077 reco::CaloClusterCollection::const_iterator bcItr;
00078 LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() All BC in the event " << "\n";
00079 for (unsigned i = 0; i < allBC->size(); ++i ) {
00080
00081
00082 nSeedsPerBC_=0;
00083
00084 theBCEnergy_=allBC->ptrAt(i)->energy();
00085 theBCPosition_ = GlobalPoint(allBC->ptrAt(i)->position().x(), allBC->ptrAt(i)->position().y(), allBC->ptrAt(i)->position().z() ) ;
00086 float theBcEta= theBCPosition_.eta();
00087 float theBcPhi= theBCPosition_.phi();
00088
00089
00090
00091
00092 if ( theBCEnergy_ < 1.5 ) continue;
00093
00094 LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() Passing the >=1.5 GeV cut BC eta " << theBcEta << " phi " << theBcPhi << " BC energy " << theBCEnergy_ << "\n";
00095
00096 if ( fabs(theBcEta-theSCEta) < 0.015 && fabs(theBcPhi-theSCPhi) < 0.3 ) {
00097 LogDebug("OutInConversionSeedFinder") << " OutInConversionSeedFinder::makeSeeds() in et and phi range passed to the analysis " << "\n";
00098 fillClusterSeeds( allBC->ptrAt(i) );
00099 }
00100
00101
00102 }
00103
00104
00105
00106
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143 }
00144
00145
00146
00147 void OutInConversionSeedFinder::makeSeeds( const reco::CaloClusterPtr& aBC ) const {
00148
00149 theSeeds_.clear();
00150
00151 findLayers();
00152
00153 float theSCPhi=theSCPosition_.phi();
00154 float theSCEta=theSCPosition_.eta();
00155
00156 nSeedsPerBC_=0;
00157
00158 theBCEnergy_=aBC->energy();
00159 theBCPosition_ = GlobalPoint(aBC->position().x(), aBC->position().y(), aBC->position().z() ) ;
00160 float theBcEta= theBCPosition_.eta();
00161 float theBcPhi= theBCPosition_.phi();
00162
00163
00164 if ( theBCEnergy_ < 1.5 ) return;
00165
00166 if ( fabs(theBcEta-theSCEta) < 0.015 && fabs(theBcPhi-theSCPhi) < 0.25 ) {
00167 fillClusterSeeds( aBC);
00168 }
00169
00170 }
00171
00172
00173
00174
00175
00176 void OutInConversionSeedFinder::fillClusterSeeds(const reco::CaloClusterPtr& bc) const {
00177
00178
00179 theFirstMeasurements_.clear();
00180 FreeTrajectoryState fts;
00181
00183 if ( makeTrackState(-1).second ) {
00184 fts = makeTrackState(-1).first;
00185 startSeed(fts);
00186 }
00188
00189 if ( makeTrackState(1).second ) {
00190 fts = makeTrackState(1).first;
00191 startSeed(fts);
00192 }
00193
00194 }
00195
00196
00197
00198 std::pair<FreeTrajectoryState,bool> OutInConversionSeedFinder::makeTrackState(int charge) const {
00199
00200 std::pair<FreeTrajectoryState,bool> result;
00201 result.second=false;
00202
00203
00204
00205
00206
00207
00208
00209
00210 GlobalPoint gpOrigine(theBeamSpot_.position().x(),theBeamSpot_.position().y(),theBeamSpot_.position().z());
00211 GlobalVector gvBcRadius = theBCPosition_ - gpOrigine ;
00212 HepGeom::Point3D<double> radiusBc(gvBcRadius.x(),gvBcRadius.y(),gvBcRadius.z()) ;
00213 HepGeom::Point3D<double> momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_ ;
00214
00215
00216 double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp() ;
00217 curvature /= 100. ;
00218
00219 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x() << " " << gpOrigine.y() << " " << gpOrigine.z() << " momentumWithoutCurvature " << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
00220
00221
00222 float R = theBCPosition_.perp();
00223 float r = gpOrigine.perp();
00224 float rho = 1./curvature;
00225
00226
00227 float d = sqrt(r*r+rho*rho);
00228 float u = rho + rho/d/d*(R*R-rho*rho) - r/d/d*sqrt((R*R-r*r+2*rho*R)*(R*R-r*r+2*rho*R));
00229
00230 if ( u <=R ) result.second=true;
00231
00232 double sinAlpha = 0.5*u/R;
00233 if ( sinAlpha>(1.-10*DBL_EPSILON) ) sinAlpha = 1.-10*DBL_EPSILON;
00234 else if ( sinAlpha<-(1.-10*DBL_EPSILON) ) sinAlpha = -(1.-10*DBL_EPSILON);
00235
00236 double newdphi = charge * asin( sinAlpha) ;
00237
00238 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R << " u/R " << u/R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
00239
00240 HepGeom::Transform3D rotation = HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double> (0., 0. ,1.));
00241
00242
00243 HepGeom::Point3D<double> momentumInTracker = momentumWithoutCurvature.transform(rotation) ;
00244 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState R " << R << " r " << r << " rho " << rho << " d " << d << " u " << u << " newdphi " << newdphi << " momentumInTracker " << momentumInTracker << "\n";
00245
00246 HepGeom::Point3D<double> hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z()) ;
00247
00248 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint " << hepStartingPoint << "\n";
00249
00250 hepStartingPoint.transform(rotation);
00251
00252 GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
00253
00254 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint << " calo position " << theBCPosition_ << "\n";
00255 GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
00256 GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_) );
00257
00258 AlgebraicSymMatrix55 m = AlgebraicMatrixID();
00259 m(0,0) = 0.1; m(1,1) = 0.1 ; m(2,2) = 0.1 ;
00260 m(3,3) = 0.1 ; m(4,4) = 0.1;
00261
00262
00263
00264 result.first= FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) ;
00265 return result;
00266
00267 }
00268
00269
00270 void OutInConversionSeedFinder::startSeed(const FreeTrajectoryState & fts) const {
00271
00272
00273
00274
00275
00276 std::vector<const DetLayer*> myLayers=layerList();
00277 if ( myLayers.size() > 3 ) {
00278
00279 for(unsigned int ilayer = myLayers.size()-1; ilayer >= myLayers.size()-2; --ilayer) {
00280 const DetLayer * layer = myLayers[ilayer];
00281
00282
00283
00284
00285
00286 float dphi = 0.030;
00287
00288 MeasurementEstimator * newEstimator = makeEstimator(layer, dphi);
00289
00290
00291
00292
00293 TSOS tsos(fts, layer->surface() );
00294
00295 LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed after TSOS tsos(fts, layer->surface() ) " << "\n";
00296
00297 LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
00298 theFirstMeasurements_ = theLayerMeasurements_.measurements( *layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
00299
00300
00301
00302 if(theFirstMeasurements_.size() > 1)
00303 LogDebug("OutInConversionSeedFinder") << " Found " << theFirstMeasurements_.size()-1 << " 1st hits in seed" << "\n";
00304
00305 delete newEstimator;
00306
00307 LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
00308
00309 for(unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
00310 TrajectoryMeasurement m1 = theFirstMeasurements_[i];
00311 if(m1.recHit()->isValid()) {
00312
00313
00314
00315 GlobalPoint hitPoint = m1.recHit()->globalPosition();
00316 LogDebug("OutInConversionSeedFinder") << " Valid hit at R " << m1.recHit()->globalPosition().perp() << " Z " << m1.recHit()->globalPosition().z() << " eta " << m1.recHit()->globalPosition().eta() << " phi " << m1.recHit()->globalPosition().phi() << " xyz " << m1.recHit()->globalPosition() << "\n";
00317
00318
00319 FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
00320 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed newfts " << newfts << "\n";
00321 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed propagationDirection after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
00322
00323
00324
00325 completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-1);
00326
00327 if(ilayer == myLayers.size()-1) {
00328 completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-2);
00329 }
00330 }
00331 }
00332
00333 }
00334 }
00335
00336
00337
00338
00339 }
00340
00341
00342
00343 MeasurementEstimator * OutInConversionSeedFinder::makeEstimator(const DetLayer * layer, float dphi) const {
00344
00345
00346
00347 MeasurementEstimator * newEstimator=0;
00348
00349 if (layer->location() == GeomDetEnumerators::barrel ) {
00350
00351 const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
00352 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Barrel r = " << barrelLayer->specificSurface().radius() << " " << "\n";
00353 float r = barrelLayer->specificSurface().radius();
00354 float zrange = 15.* (1.-r/theBCPosition_.perp());
00355 newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
00356 }
00357
00358
00359
00360 if (layer->location() == GeomDetEnumerators::endcap ) {
00361
00362 const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
00363 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius() << " R " << forwardLayer->specificSurface().outerRadius() << " Z " << forwardLayer->specificSurface().position().z() << "\n";
00364
00365 float zc = fabs(theBCPosition_.z());
00366 float z = fabs(forwardLayer->surface().position().z());
00367
00368 float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc*zc - 15.*zc);
00369 newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
00370 }
00371
00372
00373
00374
00375 return newEstimator;
00376 }
00377
00378
00379
00380
00381 void OutInConversionSeedFinder::completeSeed(const TrajectoryMeasurement & m1,
00382 FreeTrajectoryState & fts,
00383 const Propagator* propagator, int ilayer) const {
00384
00385
00386
00387 MeasurementEstimator * newEstimator=0;
00388 const DetLayer * layer = theLayerList_[ilayer];
00389
00390
00391 if ( layer->location() == GeomDetEnumerators::barrel ) {
00392
00393 LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
00394 float dz = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1.recHit()->globalPositionError().czz()
00395 + the2ndHitdzConst_*the2ndHitdzConst_);
00396 newEstimator =
00397 new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
00398 }
00399 else {
00400 LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
00401
00402
00403 float m1dr = sqrt(m1.recHit()->localPositionError().yy());
00404 float dr = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1dr*m1dr
00405 + the2ndHitdzConst_*the2ndHitdznSigma_);
00406
00407 newEstimator =
00408 new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
00409 }
00410
00411 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
00412
00413
00414 TSOS tsos(fts, layer->surface() );
00415 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection " << int(propagator->propagationDirection() ) << "\n";
00416 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
00417
00418 LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
00419 std::vector<TrajectoryMeasurement> measurements = theLayerMeasurements_.measurements( *layer, tsos, *propagator, *newEstimator);
00420
00421 delete newEstimator;
00422
00423 for(unsigned int i = 0; i < measurements.size(); ++i) {
00424 if( measurements[i].recHit()->isValid() ) {
00425 createSeed(m1, measurements[i]);
00426 }
00427 }
00428
00429
00430
00431
00432 }
00433
00434
00435
00436 void OutInConversionSeedFinder::createSeed(const TrajectoryMeasurement & m1,
00437 const TrajectoryMeasurement & m2) const {
00438
00439
00440
00441
00442 FreeTrajectoryState fts = createSeedFTS(m1, m2);
00443
00444
00445 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
00446 LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts <<"\n";
00447
00448
00449
00450
00451 TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts, m1.recHit()->det()->surface());
00452
00453
00454
00455
00456
00457
00458 if ( state1.isValid() ) {
00459 TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit() );
00460
00461 if ( updatedState1.isValid() ) {
00462 TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(*updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
00463
00464 if ( state2.isValid() ) {
00465
00466 TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit() );
00467 TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit() , m1.estimate(), m1.layer());
00468 TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit() , m2.estimate(), m2.layer());
00469
00470 edm::OwnVector<TrackingRecHit> myHits;
00471 myHits.push_back(meas1.recHit()->hit()->clone());
00472 myHits.push_back(meas2.recHit()->hit()->clone());
00473
00474 if ( nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_ ) return;
00475
00476 TrajectoryStateTransform tsTransform;
00477 PTrajectoryStateOnDet* ptsod= tsTransform.persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId() );
00478
00479 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed from state " << state2.globalPosition() << "\n";
00480 LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed ptsod " << ptsod->parameters().position() << " R " << ptsod->parameters().position().perp() << " phi " << ptsod->parameters().position().phi() << " eta " << ptsod->parameters().position().eta() << "\n" ;
00481
00482
00483
00484 theSeeds_.push_back(TrajectorySeed( *ptsod, myHits, oppositeToMomentum ));
00485 nSeedsPerBC_++;
00486
00487 delete ptsod;
00488
00489
00490 }
00491 }
00492 }
00493
00494
00495
00496
00497
00498 }
00499
00500
00501
00502
00503
00504 FreeTrajectoryState OutInConversionSeedFinder::createSeedFTS(const TrajectoryMeasurement & m1,
00505 const TrajectoryMeasurement & m2) const {
00506
00507
00508
00509 GlobalPoint xmeas = fixPointRadius(m1);
00510 GlobalPoint xvert = fixPointRadius(m2);
00511
00512
00513 float pt = theSCenergy_ * sin(theSCPosition_.theta());
00514 float pz = theSCenergy_ * cos(theSCPosition_.theta());
00515
00516
00517
00518
00519
00520
00521 int charge = m2.predictedState().charge();
00522
00523 double BInTesla = theMF_->inTesla(xmeas).z();
00524 GlobalVector xdiff = xmeas -xvert;
00525
00526 double phi= xdiff.phi();
00527 double pxOld = pt*cos(phi);
00528 double pyOld = pt*sin(phi);
00529 double RadCurv = 100*pt/(BInTesla*0.29979);
00530 double alpha = asin(0.5*xdiff.perp()/RadCurv);
00531 float ca = cos(charge*alpha);
00532 float sa = sin(charge*alpha);
00533 double pxNew = ca*pxOld + sa*pyOld;
00534 double pyNew = -sa*pxOld + ca*pyOld;
00535 GlobalVector pNew(pxNew, pyNew, pz);
00536
00537 GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_) );
00538
00539 AlgebraicSymMatrix55 m = AlgebraicMatrixID();
00540 m(0,0) = 0.05; m(1,1) = 0.02 ; m(2,2) = 0.007 ;
00541 m(3,3) = 10. ; m(4,4) = 10. ;
00542 return FreeTrajectoryState(gp, CurvilinearTrajectoryError(m));
00543
00544
00545 }
00546
00547
00548
00549
00550 GlobalPoint OutInConversionSeedFinder::fixPointRadius(const TrajectoryMeasurement& m1) const {
00551 GlobalPoint p1 = m1.recHit()->globalPosition();
00552 GlobalPoint p2;
00553 if(m1.layer()->location() == GeomDetEnumerators::barrel) {
00554 p2 = p1;
00555 } else {
00556 float z = p1.z();
00557 float phi = p1.phi();
00558 float theta = theSCPosition_.theta();
00559 float r = p1.z() * tan(theta);
00560 p2 = GlobalPoint(r*cos(phi), r*sin(phi), z);
00561
00562 }
00563 return p2;
00564 }
00565
00566
00567
00568
00569