CMSSW Reference Manual

00001 // -------------------------------------------------------------------------
00002 // File and Version Information:
00003 //      $Id: DcxTrackCandidatesToTracks.cc,v 1.8 2007/03/07 22:04:03 gutsche Exp $
00004 //
00005 // Description:
00006 //      Class Implementation for |DcxTrackCandidatesToTracks|
00007 //
00008 // Environment:
00009 //      Software developed for the BaBar Detector at the SLAC B-Factory.
00010 //
00011 // Author List:
00012 //      S. Wagner
00013 //
00014 // Copyright Information:
00015 //      Copyright (C) 1995      SLAC
00016 //
00017 //------------------------------------------------------------------------
00018 #include <cmath>
00019 #include "RecoTracker/RoadSearchHelixMaker/interface/DcxTrackCandidatesToTracks.hh"
00020 #include "RecoTracker/RoadSearchHelixMaker/interface/Dcxmatinv.hh"
00021 #include "RecoTracker/RoadSearchHelixMaker/interface/DcxFittedHel.hh"
00022 #include "RecoTracker/RoadSearchHelixMaker/interface/DcxHit.hh"
00023 #include "RecoTracker/RoadSearchHelixMaker/interface/Dcxprobab.hh"
00024 
00025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00026 
00027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00028 
00029 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00030 
00031 // #include "TrackingTools/TrajectoryParametrization/interface/PerigeeTrajectoryParameters.h"
00032 // #include "TrackingTools/TrajectoryParametrization/interface/PerigeeTrajectoryError.h"
00033 
00034 #include "TrackingTools/PatternTools/interface/TrajectoryStateClosestToPointBuilder.h"
00035 
00036 #include "DataFormats/Math/interface/Vector3D.h"
00037 
00038 using std::cout;
00039 using std::endl;
00040 using std::ostream;
00041 
00042 double DcxTrackCandidatesToTracks::epsilon      =   0.000000001;
00043 double DcxTrackCandidatesToTracks::half_pi=1.570796327;
00044 
00045 //constructors
00046 DcxTrackCandidatesToTracks::DcxTrackCandidatesToTracks(){
00047   edm::LogInfo("RoadSearch") << "DcxTrackCandidatesToTracks null constructor - does nothing" ;}
00048 
00049 //points
00050 DcxTrackCandidatesToTracks::DcxTrackCandidatesToTracks(std::vector<DcxHit*> &listohits, reco::TrackCollection &output, const MagneticField *field)
00051 { 
00052   double rmin=1000.0; double phi_try=0.0; int ntrk=0;
00053   for (unsigned int i=0; i<listohits.size(); ++i) {
00054     if (!listohits[i]->stereo()){
00055       double rhit=sqrt(listohits[i]->x()*listohits[i]->x()+listohits[i]->y()*listohits[i]->y());
00056       if ( (rhit<rmin) ){phi_try=atan2(listohits[i]->y(),listohits[i]->x());rmin=rhit;}
00057     }
00058     for (unsigned int j=0; j<listohits.size(); ++j) {
00059       for (unsigned int k=0; k<listohits.size(); ++k) {
00060         if ( ( ( 1==listohits[i]->Layer())||( 3==listohits[i]->Layer()) ) 
00061              && ( ( 9==listohits[j]->Layer())||(11==listohits[j]->Layer()) )
00062              && ( (17==listohits[k]->Layer())||(19==listohits[k]->Layer()) ) ){
00063           makecircle(listohits[i]->x(),listohits[i]->y(),listohits[j]->x(),listohits[j]->y(),
00064                      listohits[k]->x(),listohits[k]->y());
00065           double xc=xc_cs; double yc=yc_cs; double rc=rc_cs;
00066           double d0=sqrt(xc*xc+yc*yc)-rc; 
00067 //        double dmax=sqrt(xc*xc+yc*yc)+rc; // not used
00068           double s3=s3_cs;
00069           if ( (fabs(s3)>0.1) ){
00070             double d0h=-s3*d0;
00071             double phi0h=atan2(yc,xc)+s3*half_pi;
00072             double omegah=-s3/rc;
00073             DcxHel make_a_circ(d0h,phi0h,omegah,0.0,0.0);// make_a_hel.print();
00074             std::vector<DcxHit*> axlist;
00075             check_axial( listohits, axlist, make_a_circ);
00076             int n_axial=axlist.size();
00077             for (unsigned int l=0; l<listohits.size(); ++l) {
00078               for (unsigned int m=0; m<listohits.size(); ++m) {
00079                 if ( ((2==listohits[l]->Layer())||(4==listohits[l]->Layer())) && 
00080                      ((12==listohits[m]->Layer())||(10==listohits[m]->Layer())) ){
00081                   double z1=find_z_at_cyl(make_a_circ,listohits[l]); 
00082                   double z2=find_z_at_cyl(make_a_circ,listohits[m]);
00083                   double l1=find_l_at_z(z1,make_a_circ,listohits[l]); 
00084                   double l2=find_l_at_z(z2,make_a_circ,listohits[m]);
00085                   double tanl=(z1-z2)/(l1-l2); double z0=z1-tanl*l1;
00086                   DcxHel make_a_hel(d0h,phi0h,omegah,z0,tanl);// make_a_hel.print();
00087                   std::vector<DcxHit*> outlist = axlist;
00088                   check_stereo( listohits, outlist, make_a_hel);
00089                   int n_stereo=outlist.size()-n_axial;
00090                   if ((n_stereo>2)&&(n_axial>6)){
00091                     DcxFittedHel real_fit(outlist,make_a_hel);// real_fit.FitPrint();
00092                     if (real_fit.Prob()>0.001){
00093                       ntrk++;
00094                       edm::LogInfo("RoadSearch") << "ntrk xprob pt nax nst d0 phi0 omega z0 tanl " << ntrk << " " << real_fit.Prob() 
00095                                                  << " " << real_fit.Pt() << " " << n_axial << " " << n_stereo
00096                                                  << " " << real_fit.D0() << " " << real_fit.Phi0() << " " << real_fit.Omega() 
00097                                                  << " " << real_fit.Z0() << " " << real_fit.Tanl() ;
00098 
00099 
00100                       AlgebraicVector paraVector(5);
00101                       paraVector[0] = real_fit.Omega();
00102                       paraVector[1] = half_pi - atan(real_fit.Tanl());
00103                       paraVector[2] = real_fit.Phi0();
00104                       paraVector[3] = - real_fit.D0();
00105                       paraVector[4] = real_fit.Z0();
00106 
00107                       PerigeeTrajectoryParameters params(paraVector);
00108 
00109                       PerigeeTrajectoryError error;
00110 
00111                       GlobalPoint reference(0,0,0);
00112 
00113                       TrajectoryStateClosestToPoint tscp(params, 
00114                                                          real_fit.Pt(),
00115                                                          error, 
00116                                                          reference,
00117                                                          field);
00118 
00119                       GlobalPoint v = tscp.theState().position();
00120                       math::XYZPoint  pos( v.x(), v.y(), v.z() );
00121                       GlobalVector p = tscp.theState().momentum();
00122                       math::XYZVector mom( p.x(), p.y(), p.z() );
00123 
00124                       output.push_back(reco::Track(real_fit.Chisq(),
00125                                                    outlist.size()-5,
00126                                                    pos, 
00127                                                    mom, 
00128                                                    tscp.charge(), 
00129                                                    tscp.theState().curvilinearError()));
00130 
00131                       for (unsigned int n=0; n<outlist.size(); ++n){outlist[n]->SetUsedOnHel(ntrk);}
00132                     }else{
00133                     }
00134                   }
00135                 }
00136               }
00137             }
00138           }
00139         }
00140       }
00141     }
00142   }
00143 }//endof DcxTrackCandidatesToTracks
00144 
00145 //destructor
00146 DcxTrackCandidatesToTracks::~DcxTrackCandidatesToTracks( ){ }//endof ~DcxTrackCandidatesToTracks
00147 
00148 void DcxTrackCandidatesToTracks::makecircle(double x1_cs, double y1_cs, double x2_cs,
00149                                    double y2_cs, double x3_cs, double y3_cs){
00150   x1t_cs=x1_cs-x3_cs; y1t_cs=y1_cs-y3_cs; r1s_cs=x1t_cs*x1t_cs+y1t_cs*y1t_cs;
00151   x2t_cs=x2_cs-x3_cs; y2t_cs=y2_cs-y3_cs; r2s_cs=x2t_cs*x2t_cs+y2t_cs*y2t_cs;
00152   double rho=x1t_cs*y2t_cs-x2t_cs*y1t_cs;
00153   if (fabs(rho)<DcxTrackCandidatesToTracks::epsilon){
00154     rc_cs=1.0/(DcxTrackCandidatesToTracks::epsilon);
00155     fac_cs=sqrt(x1t_cs*x1t_cs+y1t_cs*y1t_cs);
00156     xc_cs=x2_cs+y1t_cs*rc_cs/fac_cs;
00157     yc_cs=y2_cs-x1t_cs*rc_cs/fac_cs;
00158   }else{
00159     fac_cs=0.5/rho;
00160     xc_cs=fac_cs*(r1s_cs*y2t_cs-r2s_cs*y1t_cs);
00161     yc_cs=fac_cs*(r2s_cs*x1t_cs-r1s_cs*x2t_cs); 
00162     rc_cs=sqrt(xc_cs*xc_cs+yc_cs*yc_cs); xc_cs+=x3_cs; yc_cs+=y3_cs;
00163   }
00164   s3_cs=0.0;
00165   f1_cs=x1_cs*yc_cs-y1_cs*xc_cs; f2_cs=x2_cs*yc_cs-y2_cs*xc_cs; 
00166   f3_cs=x3_cs*yc_cs-y3_cs*xc_cs;
00167   if ((f1_cs<0.0)&&(f2_cs<0.0)&&(f3_cs<0.0))s3_cs=1.0;
00168   if ((f1_cs>0.0)&&(f2_cs>0.0)&&(f3_cs>0.0))s3_cs=-1.0;
00169 }
00170 
00171 void DcxTrackCandidatesToTracks::check_axial( std::vector<DcxHit*> &listohits, std::vector<DcxHit*> &outlist, DcxHel make_a_hel){
00172   for (unsigned int i=0; i<listohits.size(); ++i) {
00173     DcxHit* try_me = listohits[i];
00174     if ((!try_me->stereo())&&(!try_me->GetUsedOnHel())){
00175       double doca = try_me->residual(make_a_hel); 
00176 //      double len = make_a_hel.Doca_Len();
00177 //      edm::LogInfo("RoadSearch") << "In check_axial, doca(mmm) len xh yh " << 10000.0*doca << " " << len 
00178 //                            << " " << make_a_hel.Xh(len) << " " << make_a_hel.Yh() ;
00179       if (fabs(doca)<0.100)outlist.push_back(try_me);
00180     } 
00181   }
00182 }
00183 
00184 void DcxTrackCandidatesToTracks::check_stereo( std::vector<DcxHit*> &listohits, std::vector<DcxHit*> &outlist, DcxHel make_a_hel){
00185   for (unsigned int i=0; i<listohits.size(); ++i) {
00186     DcxHit* try_me = listohits[i];
00187     if ((try_me->stereo())&&(!try_me->GetUsedOnHel())){
00188       double doca = try_me->residual(make_a_hel);
00189 //      double len = make_a_hel.Doca_Len();
00190 //      edm::LogInfo("RoadSearch") << "In check_stereo, doca(mmm) len xh yh " << 10000.0*doca << " " << len 
00191 //                            << " " << make_a_hel.Xh(len) << " " << make_a_hel.Yh(len) ;
00192       if (fabs(doca)<0.100)outlist.push_back(try_me);
00193     } 
00194   }
00195 }
00196 
00197 double DcxTrackCandidatesToTracks::find_z_at_cyl(DcxHel he, DcxHit* hi){
00198   zint=-1000.0;
00199   x0_wr=hi->x(); y0_wr=hi->y(); sx_wr=hi->wx(); sy_wr=hi->wy();
00200   xc_cl=he.Xc(); yc_cl=he.Yc(); r_cl=fabs(1.0/he.Omega()); 
00201   double cx=x0_wr-xc_cl; double cy=y0_wr-yc_cl;
00202   double a=sx_wr*sx_wr+sy_wr*sy_wr; double b=2.0*(sx_wr*cx+sy_wr*cy); 
00203   double c=(cx*cx+cy*cy-r_cl*r_cl);
00204   double bsq=b*b; double fac=4.0*a*c; double ta=2.0*a;
00205   if (fac < bsq){
00206     double sfac=sqrt(bsq-fac); double z1=-(b-sfac)/ta; double z2=-(b+sfac)/ta;
00207     if (fabs(z1) < fabs(z2)){
00208       zint=z1;
00209     }else{
00210       zint=z2;
00211     }
00212   }
00213   return zint;
00214 }
00215 double DcxTrackCandidatesToTracks::find_l_at_z(double zi, DcxHel he, DcxHit* hi){
00216   double omega=he.Omega(), d0=he.D0();
00217   double xl = hi->x()+zi*hi->wx()/hi->wz();
00218   double yl = hi->y()+zi*hi->wy()/hi->wz();
00219   double rl = sqrt(xl*xl+yl*yl);
00220   double orcsq=(1.0+d0*omega)*(1.0+d0*omega);    
00221   double orlsq=(rl*omega)*(rl*omega);
00222   double cphil=(1.0+orcsq-orlsq)/(2.0*(1.0+d0*omega));
00223   double phil1=acos(cphil);
00224   double l1=fabs(phil1/omega); 
00225   return l1;
00226 }