|
|
|
#include <DiMuonSeedGeneratorHIC.h>
Definition at line 69 of file DiMuonSeedGeneratorHIC.h.
| typedef std::vector<DiMuonTrajectorySeed> cms::DiMuonSeedGeneratorHIC::SeedContainer |
Definition at line 72 of file DiMuonSeedGeneratorHIC.h.
| typedef SeedContainer::iterator cms::DiMuonSeedGeneratorHIC::SeedIterator |
Definition at line 73 of file DiMuonSeedGeneratorHIC.h.
| cms::DiMuonSeedGeneratorHIC::DiMuonSeedGeneratorHIC | ( | edm::InputTag | , |
| const MagneticField * | , | ||
| const GeometricSearchTracker * | , | ||
| const HICConst * | , | ||
| const std::string | , | ||
| int | aMult | ||
| ) |
Definition at line 16 of file DiMuonSeedGeneratorHIC.cc.
References autoMagneticFieldProducer_cfi::magfield, and oppositeToMomentum.
{
// initialization
theTracker = theTracker0;
theHICConst = hh;
theLowMult = aMult;
magfield = magfield0;
rphirecHitsTag = rphirecHitsTag0;
TTRHbuilder = 0;
builderName = bb;
thePropagator=new PropagatorWithMaterial(oppositeToMomentum,0.1057,&(*magfield) );
}
| virtual cms::DiMuonSeedGeneratorHIC::~DiMuonSeedGeneratorHIC | ( | ) | [inline, virtual] |
Definition at line 91 of file DiMuonSeedGeneratorHIC.h.
{}
| pair< TrajectoryMeasurement, bool > cms::DiMuonSeedGeneratorHIC::barrelUpdateSeed | ( | const FreeTrajectoryState & | FTSOLD, |
| const TrajectoryMeasurement & | tm | ||
| ) | const [private] |
Definition at line 165 of file DiMuonSeedGeneratorHIC.cc.
References GlobalTrajectoryParameters::charge(), funct::cos(), gather_cfg::cout, TrajectoryMeasurement::forwardPredictedState(), TrajectoryStateOnSurface::freeTrajectoryState(), TrajectoryMeasurement::layer(), m, autoMagneticFieldProducer_cfi::magfield, edm::max(), edm::min(), GlobalTrajectoryParameters::momentum(), FreeTrajectoryState::parameters(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, pi, GlobalTrajectoryParameters::position(), ptmin, TrajectoryMeasurement::recHit(), funct::sin(), mathSSE::sqrt(), GeometricSearchDet::surface(), funct::tan(), theta(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
{
bool good=false;
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::barrelUpdateSeed::BarrelSeed "<<std::endl;
#endif
const DetLayer* dl = tm.layer();
// std::cout<<" BarrelSeed 0"<<std::endl;
const TransientTrackingRecHit::ConstRecHitPointer rh = tm.recHit();
// std::cout<<" BarrelSeed 1"<<std::endl;
if(!(rh->isValid()))
{
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::barrelUpdateSeed::hit is not valid "<<std::endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::barrelUpdateSeed::hit is valid "<<std::endl;
#endif
// std::cout<<" BarrelSeed 2"<<std::endl;
FreeTrajectoryState FTS = *(tm.forwardPredictedState().freeTrajectoryState());
//
// Define local variables.
//
int imin = 0;
int imax = 0;
int imin1 = 0;
int imax1 = 0;
double phi = FTSOLD.parameters().position().phi();
double pt = FTS.parameters().momentum().perp();
double aCharge = FTS.parameters().charge();
//AlgebraicSymMatrix55 e = FTS.curvilinearError().matrix();
double dptup = 0.35*pt;
double dptdown = 0.7*pt;
double ptshift = 0.22*pt;
// std::cout<<" BarrelSeed 3 "<<std::endl;
//
// Calculate a bin for lower and upper boundary of PT interval available for track.
//
int imax0 = (int)((pt+ptshift+dptup-theHICConst->ptboun)/theHICConst->step) + 1;
int imin0 = (int)((pt+ptshift-dptdown-theHICConst->ptboun)/theHICConst->step) + 1;
if( imin0 < 1 ) imin0 = 1;
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::barrelUpdateSeed::imin0,imax0 "<<imin0<<" "<<imax0<<" pt,dpt "<<pt+ptshift<<" "<<dptup<<" "<<dptdown<<std::endl;
#endif
double dens,df,ptmax,ptmin;
GlobalPoint realhit = (*rh).globalPosition();
df = fabs(realhit.phi() - phi);
double pi=4.*atan(1.);
double twopi=8.*atan(1.);
if(df > pi) df = twopi-df;
if(df > 1.e-5)
{
dens = 1./df;
} //end if
else
{
dens = 100000.;
} // end else
// std::cout<<" Phi rh "<<realhit.phi()<<" phumu "<<phi<<" df "<<df<<" dens "<<dens<<std::endl;
//
// Calculate new imin, imax, pt (works till 20GeV/c)
// It is necessary to parametrized for different Pt value with some step (to be done)
//
ptmax = (dens-(double)(theHICConst->phias[26]))/(double)(theHICConst->phibs[26]) + theHICConst->ptbmax;
ptmin = (dens-(double)(theHICConst->phiai[26]))/(double)(theHICConst->phibi[26]) + theHICConst->ptbmax;
#ifdef DEBUG
std::cout<<" Phias,phibs,phiai,phibi "<<theHICConst->phias[26]<<" "<<theHICConst->phibs[26]<<" "<<
theHICConst->phiai[26]<<" "<<theHICConst->phibi[26]<<" "<<theHICConst->ptbmax<<std::endl;
std::cout<<" ptmin= "<<ptmin<<" ptmax "<<ptmax<<std::endl;
std::cout<<" ptboun "<<theHICConst->ptboun<<" "<<theHICConst->step<<std::endl;
#endif
imax = (int)((ptmax-theHICConst->ptboun)/theHICConst->step)+1;
imin = (int)((ptmin-theHICConst->ptboun)/theHICConst->step)+1;
if(imin > imax) {
#ifdef DEBUG
std::cout<<" imin>imax "<<imin<<" "<<imax<<std::endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);}
if(imax < 1) {
#ifdef DEBUG
std::cout<<"imax < 1 "<<imax<<std::endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);}
imin1 = max(imin,imin0);
imax1 = min(imax,imax0);
if(imin1 > imax1) {
#ifdef DEBUG
std::cout<<" imin,imax "<<imin<<" "<<imax<<std::endl;
std::cout<<" imin,imax "<<imin0<<" "<<imax0<<std::endl;
std::cout<<" imin1>imax1 "<<imin1<<" "<<imax1<<std::endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
//
// Define new trajectory.
//
double ptnew = theHICConst->ptboun + theHICConst->step * (imax1 + imin1)/2. - theHICConst->step/2.; // recalculated PT of track
//
// new theta angle of track
//
double dfmax = 1./((double)(theHICConst->phias[26])+(double)(theHICConst->phibs[26])*(ptnew-theHICConst->ptbmax));
double dfmin = 1./((double)(theHICConst->phiai[26])+(double)(theHICConst->phibi[26])*(ptnew-theHICConst->ptbmax));
double dfcalc = fabs(dfmax+dfmin)/2.;
double phinew = phi+aCharge*dfcalc;
//
// Recalculate phi, and Z.
//
double rad = 100.*ptnew/(0.3*4.);
double alf = 2.*asin(realhit.perp()/rad);
double alfnew = phinew - aCharge*alf;
//
// Fill GlobalPoint,GlobalVector
//
double delx = realhit.z()-theHICConst->zvert;
double delr = sqrt( realhit.y()*realhit.y()+realhit.x()*realhit.x() );
double theta = atan2(delr,delx);
// std::cout<<" Point 3 "<<std::endl;
//
// Each trajectory in tracker starts from real point
// GlobalPoint xnew0( realhit.perp()*cos(phinew), realhit.perp()*sin(phinew), realhit.z() );
GlobalPoint xnew0( realhit.x(), realhit.y(), realhit.z() );
GlobalVector pnew0(ptnew*cos(alfnew),ptnew*sin(alfnew),ptnew/tan(theta));
AlgebraicSymMatrix55 m;
m(0,0) = 0.5*ptnew; m(1,1) = theHICConst->phiro[12];
m(2,2) = theHICConst->tetro[12];
m(3,3) = theHICConst->phiro[12];
m(4,4) = theHICConst->tetro[12];
TrajectoryStateOnSurface updatedTsosOnDet=TrajectoryStateOnSurface
( GlobalTrajectoryParameters( xnew0, pnew0, (int)aCharge, &(*magfield) ),
CurvilinearTrajectoryError(m), dl->surface() );
float estimate = 1.;
TrajectoryMeasurement newtm(tm.forwardPredictedState(), updatedTsosOnDet, rh, estimate, dl );
good=true;
pair<TrajectoryMeasurement,bool> newtmr(newtm,good);
// std::cout<<" Barrel newtm estimate= "<<newtmr.first.estimate()<<" "<<newtmr.second<<std::endl;
return newtmr;
}
| pair< TrajectoryMeasurement, bool > cms::DiMuonSeedGeneratorHIC::forwardUpdateSeed | ( | const FreeTrajectoryState & | FTSOLD, |
| const TrajectoryMeasurement & | tm | ||
| ) | const [private] |
Definition at line 331 of file DiMuonSeedGeneratorHIC.cc.
References GlobalTrajectoryParameters::charge(), funct::cos(), gather_cfg::cout, PV3DBase< T, PVType, FrameType >::eta(), TrajectoryMeasurement::forwardPredictedState(), TrajectoryStateOnSurface::freeTrajectoryState(), TrajectoryMeasurement::layer(), m, autoMagneticFieldProducer_cfi::magfield, GlobalTrajectoryParameters::momentum(), FreeTrajectoryState::parameters(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, pi, GlobalTrajectoryParameters::position(), ptmin, TrajectoryMeasurement::recHit(), funct::sin(), mathSSE::sqrt(), GeometricSearchDet::surface(), funct::tan(), theta(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
{
bool good=false;
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::forwardUpdateSeed::EndcapSeed::start "<<std::endl;
#endif
const DetLayer* dl = tm.layer();
// std::cout<<" EndcapSeed 0"<<std::endl;
const TransientTrackingRecHit::ConstRecHitPointer rh = tm.recHit();
// std::cout<<" EndcapSeed 1"<<std::endl;
if(!(rh->isValid()))
{
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::forwardUpdateSeed::EndcapSeed::EndcapSeed::hit is not valid "<<std::endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
#ifdef DEBUG
std::cout<<" DiMuonSeedGeneratorHIC::forwardUpdateSeed::EndcapSeed::EndcapSeed::valid "<<std::endl;
#endif
FreeTrajectoryState FTS = *(tm.forwardPredictedState().freeTrajectoryState());
//
// Define local variables.
//
double phi = FTSOLD.parameters().position().phi();
double aCharge = FTS.parameters().charge();
//AlgebraicSymMatrix55 e = FTS.curvilinearError().matrix();
double pt = FTS.parameters().momentum().perp();
// double pz = FTS.parameters().momentum().z();
// double dpt = 0.6*pt;
// std::cout<<" Point 0 "<<std::endl;
double pi=4.*atan(1.);
double twopi=8.*atan(1.);
GlobalPoint realhit = rh->globalPosition();
// std::cout<<" Point 1 "<<std::endl;
double df = fabs(realhit.phi() - phi);
if(df > pi) df = twopi-df;
#ifdef DEBUG
cout<<" DiMuonSeedGeneratorHIC::forwardUpdateSeed::phipred::phihit::df "<<phi<<" "<<realhit.phi()<<" "<<df<<endl;
#endif
//
// calculate the new Pl
//
double delx = realhit.z() - theHICConst->zvert;
double delr = sqrt(realhit.y()*realhit.y()+realhit.x()*realhit.x());
double theta = atan2( delr, delx );
double ptmin = 0.;
double ptmax = 0.;
double ptnew = 0.;
double pznew = 0.;
// old ok double pznew = abs((aCharge*theHicConst->forwparam[1])/(df-theHicConst->forwparam[0]));
if( fabs(FTSOLD.parameters().momentum().eta()) > 1.9 )
{
#ifdef DEBUG
cout<<" First parametrization "<<df<<endl;
#endif
pznew = fabs(( df - 0.0191878 )/(-0.0015952))/3.;
if( df > 0.1 ) pznew = 5.;
if( fabs(pznew)<3.) pznew = 3.;
if( FTSOLD.parameters().position().z() < 0. ) pznew = (-1)*pznew;
ptnew = pznew * tan( theta );
}
if( fabs(FTSOLD.parameters().momentum().eta()) > 1.7 && fabs(FTSOLD.parameters().momentum().eta()) < 1.9 )
{
#ifdef DEBUG
cout<<" Second parametrization "<<df<<endl;
#endif
pznew = fabs(( df - 0.38 )/(-0.009))/3.;
if( fabs(pznew)<2.) pznew = 2.;
if( FTSOLD.parameters().position().z() < 0. ) pznew = (-1)*pznew;
ptnew = pznew * tan( theta );
}
if( fabs(FTSOLD.parameters().momentum().eta()) > 1.6 && fabs(FTSOLD.parameters().momentum().eta()) < 1.7 )
{
#ifdef DEBUG
cout<<" Third parametrization "<<df<<endl;
#endif
pznew = fabs(( df - 0.9 )/(-0.02))/3.;
if( fabs(pznew)<1.) pznew = 1.;
if( FTSOLD.parameters().position().z() < 0. ) pznew = (-1)*pznew;
ptnew = pznew * tan( theta );
}
if( fabs(FTSOLD.parameters().momentum().eta()) > 0.7 && fabs(FTSOLD.parameters().momentum().eta()) < 1.6 )
{
#ifdef DEBUG
cout<<" Forth parametrization "<<df<<endl;
#endif
double dfinv = 0.;
if( df < 0.0000001 ) {
dfinv = 1000000.;
}
else
{
dfinv = 1/df;
}
ptmin = (dfinv - 4.)/0.7 + 3.;
if( ptmin < 2. ) ptmin = 2.;
ptmax = (dfinv - 0.5)/0.3 + 3.;
ptnew = ( ptmin + ptmax )/2.;
pznew = ptnew/tan( theta );
// std::cout<<" Point 6 "<<std::endl;
#ifdef DEBUG
std::cout<<" Paramters of algorithm "<<df<<" "<<theHICConst->forwparam[1]<<" "<<theHICConst->forwparam[0]<<std::endl;
std::cout<<" dfinv "<<dfinv<<" ptmax "<<ptmax<<" ptmin "<<ptmin<<std::endl;
// std::cout<<" check "<<pt<<" "<<ptnew<<" "<<dpt<<" pz "<<pznew<<" "<<pz<<std::endl;
#endif
}
//
// Check if it is valid
//
if( (pt - ptnew)/pt < -2 || (pt - ptnew)/pt > 1 )
{
#ifdef DEBUG
cout<<" Return fake 0 pt::ptnew "<<pt<<" "<<ptnew<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good); // bad rhit
}
// cout<<" Start recalculation 0 "<<endl;
//
// Recalculate phi, and Z.
//
double alf = theHICConst->atra * ( realhit.z() - theHICConst->zvert )/fabs(pznew);
double alfnew = realhit.phi() + aCharge*alf;
GlobalPoint xnew0(realhit.x(), realhit.y(), realhit.z());
GlobalVector pnew0( ptnew*cos(alfnew), ptnew*sin(alfnew), pznew );
#ifdef DEBUG
cout<<" Start recalculation 1 FTSOLD eta, r hit, pt "<<FTSOLD.parameters().momentum().eta()<<" "<<realhit.perp()<<
" "<<FTSOLD.parameters().momentum().perp()<<endl;
#endif
if( fabs(FTSOLD.parameters().momentum().eta()) < 1.7 && fabs(FTSOLD.parameters().momentum().eta()) > 0.8 )
{
if( realhit.perp() < 80. ) {
// if( realhit.perp() < 72. ) {
#ifdef DEBUG
cout<<" Return fake 1 "<<realhit.perp()<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
}
// std::cout<<" Point 9 "<<std::endl;
if( FTSOLD.parameters().momentum().perp() > 2.0){
if( fabs(FTSOLD.parameters().momentum().eta()) < 2.0 && fabs(FTSOLD.parameters().momentum().eta()) >= 1.7 )
{
if( realhit.perp() > 100. || realhit.perp() < 60. ) {
#ifdef DEBUG
cout<<" Return fake 2 "<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
}
if( fabs(FTSOLD.parameters().momentum().eta()) < 2.4 && fabs(FTSOLD.parameters().momentum().eta()) >= 2.0 )
{
if( realhit.perp() > 75. || realhit.perp() < 40. ) {
// if( realhit.perp() > 82. || realhit.perp() < 40. ) {
#ifdef DEBUG
cout<<" Return fake 3 "<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
}
}
else // pt<2
{
if( fabs(FTSOLD.parameters().momentum().eta()) < 2.0 && fabs(FTSOLD.parameters().momentum().eta()) >= 1.7 )
{
if( realhit.perp() > 84. || realhit.perp() < 40. ) {
#ifdef DEBUG
cout<<" Return fake 4 "<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
}
if( fabs(FTSOLD.parameters().momentum().eta()) < 2.4 && fabs(FTSOLD.parameters().momentum().eta()) >= 2.0 )
{
if( realhit.perp() > 84. || realhit.perp() < 40. ) {
#ifdef DEBUG
cout<<" Return fake 5 "<<endl;
#endif
return pair<TrajectoryMeasurement,bool>(tm,good);
}
}
} // pt ><2
#ifdef DEBUG
cout<<" Create new TM "<<endl;
#endif
AlgebraicSymMatrix55 m;
m(0,0) = fabs(0.5*pznew);
m(1,1) = theHICConst->phiro[13];
m(2,2) = theHICConst->tetro[13];
m(3,3) = theHICConst->phiro[13];
m(4,4) = theHICConst->tetro[13];
TrajectoryStateOnSurface updatedTsosOnDet=TrajectoryStateOnSurface
(GlobalTrajectoryParameters( xnew0, pnew0, (int)aCharge, &(*magfield) ),
CurvilinearTrajectoryError(m), dl->surface() );
float estimate=1.;
TrajectoryMeasurement newtm(tm.forwardPredictedState(), updatedTsosOnDet, rh,estimate, dl);
good=true;
pair<TrajectoryMeasurement,bool> newtmr(newtm,good);
#ifdef DEBUG
std::cout<<" Endcap newtm estimate= "<<newtmr.first.estimate()<<" "<<newtmr.second<<" pt "<<pnew0.perp()<<" pz "<<pnew0.z()<<std::endl;
#endif
return newtmr;
}
| virtual std::map<DetLayer*,SeedContainer> cms::DiMuonSeedGeneratorHIC::produce | ( | const edm::Event & | e, |
| const edm::EventSetup & | c, | ||
| FreeTrajectoryState & | , | ||
| TrajectoryStateOnSurface & | , | ||
| FreeTrajectoryState & | , | ||
| const TransientTrackingRecHitBuilder * | RecHitBuilder, | ||
| const MeasurementTracker * | measurementTracker, | ||
| std::vector< DetLayer * > * | |||
| ) | [virtual] |
Referenced by cms::HITrackVertexMaker::produceTracks().
Definition at line 110 of file DiMuonSeedGeneratorHIC.h.
std::vector<LayerWithHits*> cms::DiMuonSeedGeneratorHIC::allLayersWithHits [private] |
Definition at line 108 of file DiMuonSeedGeneratorHIC.h.
std::vector<BarrelDetLayer*> cms::DiMuonSeedGeneratorHIC::bl [private] |
Definition at line 103 of file DiMuonSeedGeneratorHIC.h.
std::string cms::DiMuonSeedGeneratorHIC::builderName [private] |
Definition at line 100 of file DiMuonSeedGeneratorHIC.h.
std::vector<ForwardDetLayer*> cms::DiMuonSeedGeneratorHIC::fneg [private] |
Definition at line 105 of file DiMuonSeedGeneratorHIC.h.
std::vector<ForwardDetLayer*> cms::DiMuonSeedGeneratorHIC::fpos [private] |
Definition at line 104 of file DiMuonSeedGeneratorHIC.h.
bool cms::DiMuonSeedGeneratorHIC::isFirstCall [private] |
Definition at line 109 of file DiMuonSeedGeneratorHIC.h.
const MagneticField* cms::DiMuonSeedGeneratorHIC::magfield [private] |
Definition at line 123 of file DiMuonSeedGeneratorHIC.h.
Definition at line 121 of file DiMuonSeedGeneratorHIC.h.
std::vector<const DetLayer*> cms::DiMuonSeedGeneratorHIC::theDetLayer [private] |
Definition at line 106 of file DiMuonSeedGeneratorHIC.h.
Definition at line 97 of file DiMuonSeedGeneratorHIC.h.
Definition at line 96 of file DiMuonSeedGeneratorHIC.h.
const HICConst* cms::DiMuonSeedGeneratorHIC::theHICConst [private] |
Definition at line 98 of file DiMuonSeedGeneratorHIC.h.
const LayerMeasurements* cms::DiMuonSeedGeneratorHIC::theLayerMeasurements [private] |
Definition at line 126 of file DiMuonSeedGeneratorHIC.h.
int cms::DiMuonSeedGeneratorHIC::theLowMult [private] |
Definition at line 101 of file DiMuonSeedGeneratorHIC.h.
const MeasurementTracker* cms::DiMuonSeedGeneratorHIC::theMeasurementTracker [private] |
Definition at line 125 of file DiMuonSeedGeneratorHIC.h.
Definition at line 99 of file DiMuonSeedGeneratorHIC.h.
const GeometricSearchTracker* cms::DiMuonSeedGeneratorHIC::theTracker [private] |
Definition at line 124 of file DiMuonSeedGeneratorHIC.h.
const TransientTrackingRecHitBuilder* cms::DiMuonSeedGeneratorHIC::TTRHbuilder [private] |
Definition at line 122 of file DiMuonSeedGeneratorHIC.h.
1.7.3