New Track Event Model: TrackFirstCleverExtrapolator Class Reference

00039                                                                         :
00040   TrackExtrapolator(type, name, parent),
00041   m_tMax(10.),
00042   m_eMax(100.0*MeV),
00043   m_freeFieldExtrapolator(0),
00044   m_shortFieldExtrapolator(0),
00045   m_longFieldExtrapolator(0)
00046 {
00047   declareInterface<ITrackExtrapolator>( this );
00048 
00049   //job options
00050   declareProperty("fieldStart"    , m_zFieldStart = -500.*mm);
00051   declareProperty("fieldStop"     , m_zFieldStop= 21000.*mm );
00052   declareProperty("shortDist"     , m_shortDist = 100.*mm);
00053   declareProperty("freeFieldExtrapolatorName",
00054                   m_freeFieldExtrapolatorName = "TrLinearExtrapolator");
00055   declareProperty("shortFieldExtrapolatorName",
00056                   m_shortFieldExtrapolatorName ="TrParabolicExtrapolator");
00057   declareProperty("longFieldExtrapolatorName",
00058                   m_longFieldExtrapolatorName = "TrHerabExtrapolator");
00059   declareProperty("applyMultScattCorr"  , m_applyMultScattCorr  = true);
00060   declareProperty("fms2"                , m_fms2                = 1.0);
00061   declareProperty("thickWall"           , m_thickWall           = 0.*mm);
00062   declareProperty("applyEnergyLossCorr" , m_applyEnergyLossCorr = true);
00063   declareProperty("energyLoss"          , 
00064                   m_energyLoss = 354.1 * MeV*mm2/mole );
00065   declareProperty("maxStepSize"         , m_maxStepSize         = 1000. * mm );
00066   declareProperty("minRadThreshold"     , m_minRadThreshold     = 1.0e-4 );
00067 
00068   //for electrons
00069   declareProperty("applyElectronEnergyLossCorr",
00070                   m_applyElectronEnergyLossCorr = true);
00071   declareProperty("startElectronCorr",
00072                   m_startElectronCorr = 2500.*mm);
00073   declareProperty("stopElectronCorr",
00074                   m_stopElectronCorr = 9000.*mm);
00075 
00076   // TODO print here the parameters giving if debug mode!!
00077   
00078 
00079 }

00122 {
00123   //choose which extrapolator to use in magnetic field - simplifies the code
00124   return (fabs(zTarget-zStart)< m_shortDist ? m_shortFieldExtrapolator :  m_longFieldExtrapolator);
00125 }

00341                                                                           {
00342 
00343   //extrapolate the state to the target z
00344   // cases:  1.  Transport entirely outside of magnetic field
00345   //         2.  entirely inside magnetic field
00346   //         3.  from outside magnet field to in passing start of field
00347   //         4.  from inside to out passing end of field
00348   //         5.  traverse the whole magnetic field
00349 
00350  
00351   // case 1 transport entirely outside field
00352   if ((GSL_MAX(zTarget,zStart) <= m_zFieldStart) ||
00353       (GSL_MIN(zTarget,zStart) >= m_zFieldStop)){
00354     verbose() << "Field free Transport from "  << zStart  << " to "
00355         << zTarget << endreq;
00356     TrackTransportStep tStep = TrackTransportStep(m_freeFieldExtrapolator,zTarget);
00357     transStepList.push_back(tStep);
00358   }
00359 
00360   //case 2 transport entirely inside field
00361   else if ((GSL_MIN(zTarget,zStart) >= m_zFieldStart) &&
00362            (GSL_MAX(zTarget,zStart) <= m_zFieldStop)){
00363 
00364     ITrackExtrapolator* magExtrapolator= chooseMagFieldExtrapolator( zStart,
00365                                                                  zTarget);
00366 
00367     if (msgLevel(MSG::VERBOSE)) verbose() << "Magnet region Transport from "
00368               << zStart  << " to "  << zTarget
00369               << " using " <<magExtrapolator->name() << endreq;
00370 
00371     TrackTransportStep tStep = TrackTransportStep( magExtrapolator, zTarget);
00372     transStepList.push_back(tStep);
00373   }
00374 
00375   //case 3 - inside to out or vice versa passing m_ZFieldStart
00376   else if ((GSL_MIN(zTarget,zStart) < m_zFieldStart) &&
00377            (GSL_MAX(zTarget,zStart) <= m_zFieldStop)){
00378 
00379     verbose()
00380         << "Transporting from inside to outside - crossing m_zFieldStart" 
00381         << endreq;
00382 
00383     if ( zTarget < m_zFieldStart ){
00384 
00385       ITrackExtrapolator* magExtrapolator =
00386         chooseMagFieldExtrapolator(zStart,m_zFieldStart);
00387 
00388       verbose()
00389           << "Magnet region Transport from "  << zStart  << " to "  
00390           << m_zFieldStart
00391           << " using " << magExtrapolator->name() << endreq;
00392       TrackTransportStep tStep1 = TrackTransportStep(magExtrapolator,m_zFieldStart);
00393       transStepList.push_back(tStep1);
00394 
00395       verbose() << "Field free Transport from "  << m_zFieldStart  
00396           << " to " << zTarget << endreq;
00397       TrackTransportStep tStep2 = TrackTransportStep(m_freeFieldExtrapolator,zTarget);
00398       transStepList.push_back(tStep2);
00399     } else {
00400 
00401       verbose() << "Field free Transport from "  << zStart  << " to "
00402           << m_zFieldStart << endreq;
00403       TrackTransportStep tStep1 = TrackTransportStep(m_freeFieldExtrapolator, 
00404                                                m_zFieldStart);
00405       transStepList.push_back(tStep1);
00406 
00407       ITrackExtrapolator* magExtrapolator = chooseMagFieldExtrapolator(
00408                                                                  m_zFieldStart,
00409                                                                  zTarget);
00410 
00411       verbose() << "Magnet region Transport from "  << m_zFieldStart
00412           <<  " to "  << zTarget
00413                 << " using " <<magExtrapolator->name() << endreq;
00414       TrackTransportStep tStep2 = TrackTransportStep(magExtrapolator,zTarget);
00415       transStepList.push_back(tStep2);
00416     }
00417   }
00418 
00419 
00420   //case 4 - inside to out or vice versa passing m_zFieldStop
00421   else if ((GSL_MIN(zTarget,zStart) >= m_zFieldStart) &&
00422            (GSL_MAX(zTarget,zStart) > m_zFieldStop)){
00423 
00424     verbose()
00425          <<  "Transporting from inside to outside - crossing m_zFieldStop" 
00426          << endreq;
00427 
00428     if (zTarget>m_zFieldStop){
00429 
00430       ITrackExtrapolator* magExtrapolator = chooseMagFieldExtrapolator(zStart,
00431                                                                  m_zFieldStop);
00432 
00433       if (msgLevel(MSG::VERBOSE)) verbose() 
00434         << "Magnet region Transport from "  << zStart << " to "  << m_zFieldStop
00435         << " using " << magExtrapolator->name() << endreq;
00436       TrackTransportStep tStep1 = TrackTransportStep(magExtrapolator,m_zFieldStop);
00437       transStepList.push_back(tStep1);
00438 
00439       verbose() << "Field free Transport from "  << zStart  << " to "
00440           << zTarget << endreq;
00441       TrackTransportStep tStep2 = TrackTransportStep(m_freeFieldExtrapolator,zTarget);
00442       transStepList.push_back(tStep2);
00443 
00444     }
00445     else{
00446 
00447       verbose() << "Field free Transport from "  << zStart  << " to "
00448           << m_zFieldStop << endreq;
00449       TrackTransportStep tStep1 = TrackTransportStep(m_freeFieldExtrapolator,
00450                                                m_zFieldStop);
00451       transStepList.push_back(tStep1);
00452 
00453       ITrackExtrapolator* magExtrapolator = chooseMagFieldExtrapolator(
00454                                                                  m_zFieldStop,
00455                                                                  zTarget);
00456 
00457       if (msgLevel(MSG::VERBOSE)) verbose() 
00458         << "Magnet region Transport from " << zStart <<  " to "  << zTarget
00459         << " using "<< magExtrapolator->name() << endreq;
00460       TrackTransportStep tStep2 = TrackTransportStep(magExtrapolator,zTarget);
00461       transStepList.push_back(tStep2);
00462     }
00463   }
00464 
00465 
00466   // case 5 - whole magnet in-between
00467   else if ((GSL_MAX(zTarget,zStart) > m_zFieldStop) &&
00468            (GSL_MIN(zTarget,zStart) < m_zFieldStart)){
00469     verbose() << "Transport from outside magnet to outside magnet" 
00470         << endreq;
00471 
00472     //determine whether up or downstream case
00473     double z1 =0.;
00474     double z2 =0.;
00475     if ( zTarget < zStart ){
00476       // upstream
00477       z1 = m_zFieldStop;
00478       z2 = m_zFieldStart;
00479     }
00480     else{
00481       // down stream
00482       z1 = m_zFieldStart;
00483       z2 = m_zFieldStop;
00484     }
00485 
00486     verbose() << "Field free Transport from "  << zStart  << " to "
00487         << z1 << endreq;
00488     TrackTransportStep tStep1 = TrackTransportStep( m_freeFieldExtrapolator, z1);
00489     transStepList.push_back(tStep1);
00490 
00491     // seems reasonable to use long field extrapolator for whole field !!!
00492     verbose() << "Long Magnet region Transport from "  << zStart
00493         << " to "  << z2 << endreq;
00494     TrackTransportStep tStep2 = TrackTransportStep( m_longFieldExtrapolator, z2);
00495     transStepList.push_back(tStep2);
00496 
00497     // step out of magnetic field
00498     verbose() << "Field free Transport from "  << zStart  << " to "
00499         << zTarget << endreq;
00500     TrackTransportStep tStep3 = TrackTransportStep(m_freeFieldExtrapolator, zTarget);
00501     transStepList.push_back(tStep3);
00502   }
00503 
00504   // end
00505   return StatusCode::SUCCESS;
00506 }

00625 {
00626   //hard energy loss for electrons
00627 
00628   double t;
00629   double norm = sqrt(1.0+gsl_pow_2(state.tx())+gsl_pow_2(state.ty()));
00630 
00631   if (m_upStream){
00632     t = radLength*norm;
00633   }
00634   else{
00635     t = -radLength*norm;
00636   }
00637 
00638   // protect against t too big
00639   if (fabs(t)>m_tMax) t = GSL_SIGN(t)*m_tMax;
00640 
00641   // apply correction
00642   HepVector& tX = state.state();
00643   HepSymMatrix& tC = state.covariance();
00644 
00645   tC.fast(5,5) += gsl_pow_2(tX[4]) * (exp(-t*log(3.0)/log(2.0))-exp(-2.0*t));
00646   tX(4) *= exp(-t);
00647 
00648   return StatusCode::SUCCESS;
00649 }

00599 {
00600   // Apply correction for dE/dx energy loss (Bethe-Block)
00601   double norm = sqrt(1.0+gsl_pow_2(state.tx())+gsl_pow_2(state.ty()));
00602   double bbLoss = tWall * norm * aMaterial->density() * m_energyLoss * 
00603     aMaterial->Z() / aMaterial->A();
00604   bbLoss = GSL_MIN(m_eMax, bbLoss);
00605   if (m_upStream == false){
00606     bbLoss *= -1.0;
00607   }
00608 
00609   // apply correction - note for now only correct the state vector
00610   HepVector& tX = state.state();
00611 
00612   if (tX[4]>0.){
00613     tX[4] = 1./(1./tX[4]+(bbLoss));
00614   }
00615   else{
00616     tX[4] = 1./(1./tX[4]-(bbLoss));
00617   }
00618 
00619   return StatusCode::SUCCESS;
00620 
00621 }

00088 {
00089 
00090   // initialize
00091   StatusCode sc = GaudiTool::initialize();
00092   if (sc.isFailure()){
00093     return Error("Failed to initialize", sc);
00094   }
00095 
00096   // request field free region extrapolator
00097   m_freeFieldExtrapolator=tool<ITrackExtrapolator>(m_freeFieldExtrapolatorName);
00098  
00099   // request a short distance magnetic field extrapolator
00100   m_shortFieldExtrapolator=tool<ITrackExtrapolator>(m_shortFieldExtrapolatorName);
00101 
00102   // request extrapolator for going short distances in magnetic field
00103   m_longFieldExtrapolator=tool<ITrackExtrapolator>(m_longFieldExtrapolatorName);
00104   
00105   // initialize transport service
00106   m_transportSvc = svc<ITransportSvc>( "TransportSvc",true);
00107 
00108   return StatusCode::SUCCESS;
00109 }

00116                                                                        {
00117   // Extrapolate the TrStateP 'state' to z=zNew
00118   StatusCode sc;
00119 
00120   size_t ndim = state.nParameters();
00121 
00122   // create transport matrix
00123   m_F = HepMatrix(ndim, ndim, 1);
00124 
00125   //check not already at current
00126   double zStart = state.z();
00127   if (fabs(zNew-zStart) < TrackParameters::hiTolerance) {
00128     // already at required z position
00129     debug() << "already at required z position" << endreq;
00130     return StatusCode::SUCCESS;
00131   }
00132 
00133   //check whether upstream or downstream
00134   if (zStart>zNew) {
00135     m_upStream = true;
00136   } else{
00137     m_upStream = false;
00138   }
00139 
00140   int    nbStep = (int)( fabs( zNew-zStart ) / m_maxStepSize ) + 1;
00141   double zStep  = ( zNew - zStart ) / nbStep;
00142   int    nWall;
00143   double maxSlope      = 5.;
00144   double maxTransverse = 10 * m;
00145   
00146   for ( int step=0 ; nbStep > step ; ++step ) {
00147 
00148     ILVolume::Intersections intersept;
00149     HepVector& tX = state.state();
00150     HepPoint3D start( tX[0], tX[1], state.z() );  // Initial position
00151     HepVector3D vect( tX[2]*zStep, tX[3]*zStep, zStep );
00152 
00153     //protect against vertical or looping tracks
00154 
00155     if ( fabs(tX[0]) > maxTransverse ) {
00156       double pos = GSL_SIGN(tX[0])* maxTransverse ;
00157       warning() << "Protect against absurd tracks: X " 
00158           << state.x() << " set to " << pos << ", abort. " << endreq;
00159       state.setX(pos);
00160       return StatusCode::FAILURE;
00161     }
00162     
00163     if ( fabs(tX[1]) > maxTransverse ) {
00164       double pos = GSL_SIGN(tX[1])*maxTransverse;
00165       warning() << "Protect against absurd tracks: Y " 
00166           << state.y() << " set to " << pos << ", abort. " << endreq;
00167       state.setY(pos);
00168       return StatusCode::FAILURE;
00169     }
00170 
00171     if (fabs(state.tx()) > maxSlope){
00172       double slopeX = GSL_SIGN(state.tx())* maxSlope;
00173       warning() << "Protect against looping tracks: Tx " 
00174           << state.tx() << " set to " << slopeX << ", abort. " << endreq;
00175       state.setTx(slopeX);
00176       return StatusCode::FAILURE;
00177     }
00178     
00179     if (fabs(state.ty()) > maxSlope){
00180       double slopeY = GSL_SIGN(state.ty())*maxSlope;
00181       warning() << "Protect against looping tracks: Ty " 
00182           << state.ty() << " set to " << slopeY << ", abort. " << endreq;
00183       state.setTy(slopeY);
00184       return StatusCode::FAILURE;
00185     }
00186 
00187     // short cut if inside the same plane of the OT stations
00188     if ( (7.0*mm > fabs(zStep)) && (5000.0*mm < start.z() ) ) {
00189       nWall = 0;
00190       debug() << "No transport between z= " << start.z() << " and " 
00191           << start.z() + vect.z() << endreq;
00192     } 
00193     // check if transport is within LHCb
00194     else if (fabs(start.x()) > 25.*m || fabs(start.y()) > 25.*m ||
00195              fabs(start.z()) > 25.*m ||
00196              fabs(start.x()+vect.x()) > 25.*m ||
00197              fabs(start.y()+vect.y()) > 25.*m ||
00198              fabs(start.z()+vect.z()) > 25.*m ) {
00199       nWall = 0;
00200       debug() << "No transport between z= " << start.z() << " and " 
00201           << start.z() + vect.z() 
00202           << ", since it reaches outside LHCb" << endreq;
00203     } 
00204     else {
00205       // chronoSvc()->chronoStart("TransportSvcT");
00206       nWall = m_transportSvc->intersections( start, vect, 0., 1., 
00207                                              intersept, m_minRadThreshold );
00208 
00209       // if (nWall == 0) warning() << " zero walls !" << endmsg; 
00210       
00211       // chronoSvc()->chronoStop("TransportSvcT");
00212       if (msgLevel(MSG::DEBUG)) {
00213         debug() << "Transport : " << nWall
00214                 << " intersepts between z= " << start.z() << " and " 
00215                 << start.z() + vect.z()                 << endreq;
00216           //       << " time " << chronoSvc()->chronoDelta("TransportSvcT" , 
00217           // IChronoStatSvc::ELAPSED )
00218       }
00219     }
00220   
00221     std::vector<double> zWall;          // z positions of walls
00222     std::vector<double> tWall;          //thickness of the wall
00223     std::vector<double> radLengthWall;  // rad length of walls
00224     std::vector<const Material*> matWall;     // Materials of the walls
00225 
00226     double t1, t2, dz;
00227     double pos, thick, radl;
00228     for ( int ii = 0 ; nWall > ii ; ++ii ) {
00229       t1 = intersept[ii].first.first  * zStep;
00230       t2 = intersept[ii].first.second * zStep;
00231       dz = zScatter(t1,t2);
00232       pos   = start.z() + dz;
00233       thick = fabs( t2 - t1 );
00234       radl  = thick / intersept[ii].second->radiationLength();
00235       if ( m_minRadThreshold < radl  ) {
00236         zWall.push_back( pos );
00237         tWall.push_back( thick);
00238         radLengthWall.push_back( radl );
00239         matWall.push_back( intersept[ii].second );
00240       }
00241 
00242       if ((msgLevel(MSG::DEBUG))&& (m_minRadThreshold < radl)){ 
00243         debug() << format(" %2i x%8.2f y%8.2f z%10.3f thick%8.3f radl%7.3f %%",
00244                   ii , 
00245                   tX[0] + dz * tX[2], 
00246                   tX[1] + dz * tX[3], 
00247                   pos, thick, 100.*radl) <<endreq ;
00248       }
00249     } // ii
00250 
00251     zWall.push_back( start.z() + zStep );
00252     radLengthWall.push_back( 0. );
00253     tWall.push_back( 0. );
00254     matWall.push_back( 0);
00255     nWall = tWall.size();
00256 
00257     // loop over the walls - last wall is `virtual one at target z'
00258     for (int iStep=0; iStep<nWall; ++iStep){
00259 
00260       //break transport into steps using different extrapolators
00261       std::list<TrackTransportStep> transStepList;
00262       sc = createTransportSteps( state.z(), zWall[iStep], transStepList );
00263       
00264       if (transStepList.empty()){
00265         return Error( "no transport steps !!!!", StatusCode::FAILURE);
00266       }
00267       
00268       // loop over the tracking steps
00269       std::list<TrackTransportStep>::iterator iTransStep = transStepList.begin();
00270       bool stopTransport = false;
00271 
00272       while ((iTransStep!=transStepList.end())&&(stopTransport == false)){
00273         
00274         //extrapolate
00275         ITrackExtrapolator* thisExtrapolator = (*iTransStep).extrapolator();      
00276         sc = thisExtrapolator->propagate( state, (*iTransStep).z() );
00277         
00278         //update f
00279         updateTransportMatrix( thisExtrapolator->transportMatrix() );
00280         
00281         // check for success
00282         if (sc.isFailure()) {
00283           stopTransport = true;
00284           warning() << "Transport to " << (*iTransStep).z()
00285               << "using "+thisExtrapolator->name() << " FAILED" << endreq;
00286         }
00287         ++iTransStep;
00288       }
00289       
00290       if (stopTransport == true) return sc;
00291       
00292       //protect against vertical or looping tracks
00293       if (fabs(state.tx()) > maxSlope){
00294         double slopeX = GSL_SIGN(  state.tx() )*maxSlope;
00295         warning() << "Protect against looping tracks: Tx " 
00296             << state.tx() << " set to " << slopeX << ", abort. " << endreq;
00297         state.setTx(slopeX);
00298         return StatusCode::FAILURE;
00299       }
00300       
00301       if (fabs(state.ty()) > maxSlope){
00302         double slopeY = GSL_SIGN(state.ty())*maxSlope;
00303         warning() << "Protect against looping tracks: Ty " 
00304             << state.ty() << " set to " << slopeY << ", abort. " << endreq;
00305         state.setTy(slopeY);
00306         return StatusCode::FAILURE;
00307       }
00308       
00309       //multiple scattering
00310       if (m_applyMultScattCorr == true){
00311         if (tWall[iStep] < m_thickWall){
00312           sc = thinScatter( state, radLengthWall[iStep] );
00313         }
00314         else{
00315           sc = thickScatter( state, tWall[iStep], radLengthWall[iStep] );
00316         }
00317       }
00318       
00319       // dE_dx energy loss
00320       if (( m_applyEnergyLossCorr == true)&&(matWall[iStep] != 0)){
00321         sc = energyLoss( state, tWall[iStep], matWall[iStep] );
00322       }
00323       
00324       //electron energy loss
00325       if ((m_applyElectronEnergyLossCorr == true) &&(11 == partId.abspid())){
00326         if ((state.z() > m_startElectronCorr) &&
00327                (state.z() < m_stopElectronCorr)){
00328           sc = electronEnergyLoss(state,radLengthWall[iStep]);
00329         }
00330       }
00331       
00332     } // loop over walls
00333   } // loop over steps
00334   
00335   verbose() << "TrStateP extrapolated succesfully" << endreq;
00336 
00337   return StatusCode::SUCCESS;
00338 }

00549 {
00550   // apply - thick scatter multiple scattering
00551   double scatLength = 0.;
00552   double norm2 = 1.0 + gsl_pow_2(state.tx()) + gsl_pow_2(state.ty());
00553   if (radLength > TrackParameters::lowTolerance){
00554     double radThick = sqrt(norm2)*radLength;
00555     //assumed GEV
00556     scatLength = radThick * gsl_pow_2(TrackParameters::moliereFactor *
00557                                 (1. + 0.038*log(radThick)));
00558   }
00559 
00560   // protect zero momentum
00561   double p = GSL_MAX( state.p(), 1.0 * MeV );
00562   double cnoise = m_fms2 * scatLength / gsl_pow_2(p);
00563 
00564   // slope covariances
00565   double covTxTx = norm2 * (1. + gsl_pow_2(state.tx())) * cnoise;
00566   double covTyTy = norm2 * (1. + gsl_pow_2(state.ty())) * cnoise;
00567   double covTxTy = norm2 * state.tx() * state.ty() * cnoise;
00568 
00569   //D - depends on whether up or downstream
00570   double D = (m_upStream) ? -1. : 1. ;
00571 
00572   HepSymMatrix Q(5,0);
00573 
00574   Q.fast(1,1) = covTxTx * gsl_pow_2(tWall) / 3.;
00575   Q.fast(2,1) = covTxTy * gsl_pow_2(tWall) / 3.;
00576   Q.fast(3,1) = 0.5*covTxTx * D * tWall;
00577   Q.fast(4,1) = 0.5*covTxTy * D * tWall;
00578 
00579   Q.fast(2,2) = covTyTy * gsl_pow_2(tWall) / 3.;
00580   Q.fast(3,2) = 0.5*covTxTy * D * tWall;
00581   Q.fast(4,2) = 0.5*covTyTy * D * tWall;
00582 
00583   Q.fast(3,3) = covTxTx;
00584   Q.fast(4,3) = covTxTy;
00585 
00586   Q.fast(4,4) = covTyTy;
00587 
00588   // update covariance matrix C = C + Q
00589   HepSymMatrix& tC = state.covariance();
00590   tC += Q;
00591 
00592   return StatusCode::SUCCESS;
00593 
00594 }

00511                                                                       {
00512 
00513   // apply multiple scattering - thin scatter
00514 
00515   double scatLength = 0.;
00516   double norm2 = 1.0+gsl_pow_2(state.tx())+gsl_pow_2(state.ty());
00517 
00518   if (radLength >TrackParameters::lowTolerance){
00519     double radThick = sqrt(norm2)*radLength;
00520     //assumed GEV
00521     scatLength = radThick*gsl_pow_2(TrackParameters::moliereFactor*
00522                               (1.+0.038*log(radThick)));
00523   }
00524 
00525   // protect 0 momentum
00526   double p = GSL_MAX(state.p(),1.0*MeV);
00527 
00528   double cnoise = m_fms2 * scatLength/gsl_pow_2(p);
00529 
00530   // multiple scattering covariance matrix - initialized to 0
00531   HepSymMatrix Q = HepSymMatrix(5,0);
00532 
00533   Q.fast(3,3) = norm2 * (1. + gsl_pow_2(state.tx())) * cnoise;
00534   Q.fast(4,4) = norm2 * (1. + gsl_pow_2(state.ty())) * cnoise;
00535   Q.fast(4,3) = norm2 * state.tx() * state.ty() * cnoise;
00536 
00537   // update covariance matrix C = C + Q
00538   HepSymMatrix& tC = state.covariance();
00539   tC += Q;
00540 
00541   return StatusCode::SUCCESS;
00542 
00543 }

00128                             {
00129   //update F - after transport step
00130   m_F = newStepF* m_F;
00131 }

00653 {
00654   double zS;
00655   if (fabs(z1-z2) <= m_thickWall){
00656     // thin scatter - at centre in z
00657     zS = 0.5*(z1+z2); 
00658   }
00659   else {
00660     // thick scatter
00661     if (m_upStream == true){
00662       zS = GSL_MIN(z1,z2);
00663     }
00664     else {
00665       zS = GSL_MAX(z1,z2);
00666     }
00667   }
00668   return zS;
00669 }


Public Member Functions
	TrackFirstCleverExtrapolator (const std::string &type, const std::string &name, const IInterface *parent)
	Constructor.
virtual	~TrackFirstCleverExtrapolator ()
	destructor
virtual StatusCode	initialize ()
	intialize and finalize
virtual StatusCode	propagate (State &state, double zNew=0, ParticleID partId=ParticleID(211))
	propagate with Linear state
Private Member Functions
StatusCode	createTransportSteps (double zStart, double zTarget, std::list< TrackTransportStep > &transList)
	make transport steps
StatusCode	thinScatter (State &state, double radLength)
	apply thick scatter Q/p state
StatusCode	thickScatter (State &state, double tWall, double radLength)
	apply thick scatter Q/p state
StatusCode	energyLoss (State &state, double tWall, const Material *aMaterial)
	apply energy loss P state
StatusCode	electronEnergyLoss (State &state, double radLength)
	electron energy loss Q/p state
ITrackExtrapolator *	chooseMagFieldExtrapolator (const double zStart, const double zTarget)
	choose Extrapolator to use in field free region
void	updateTransportMatrix (const HepMatrix &newStepF)
	update transport matrix
double	zScatter (const double z1, const double z2) const
	z scatter
Private Attributes
bool	m_upStream
int	m_particleType
double	m_tMax
	max radiation length - avoid underflow on NT
double	m_eMax
double	m_zFieldStart
	start of field
double	m_zFieldStop
	end of field
double	m_shortDist
	min distance to use RungaKutta
std::string	m_freeFieldExtrapolatorName
	extrapolator to use in field free region
std::string	m_shortFieldExtrapolatorName
	extrapolator to use for short transport in mag field
std::string	m_longFieldExtrapolatorName
	extrapolator to use for long transport in mag field
bool	m_applyMultScattCorr
	turn on/off multiple scattering correction
double	m_fms2
	factor for inflating scattering errors
double	m_thickWall
	thick wall
bool	m_applyEnergyLossCorr
	turn on/off multiple scattering correction
double	m_energyLoss
	tuneable energy loss correction
double	m_maxStepSize
	maximum length of a step
double	m_minRadThreshold
	minimal thickness of a wall
bool	m_applyElectronEnergyLossCorr
	turn on/off electron energy loss correction
double	m_startElectronCorr
	z start for electron energy loss
double	m_stopElectronCorr
	z start for electron energy loss
ITrackExtrapolator *	m_freeFieldExtrapolator
	extrapolators
ITrackExtrapolator *	m_shortFieldExtrapolator
ITrackExtrapolator *	m_longFieldExtrapolator
ITransportSvc *	m_transportSvc
	Pointer to the transport service.

TrackFirstCleverExtrapolator Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation