00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <cstdio>
00011 #include <fstream>
00012 #include <string.h>
00013 #include <map>
00014 #include <vector>
00015 #include <algorithm>
00016 #include "CLHEP/Geometry/Vector3D.h"
00017 #include "CLHEP/Geometry/Point3D.h"
00018 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00019 typedef HepGeom::Point3D<double> HepPoint3D;
00020 #endif
00021 #include "CLHEP/Matrix/Vector.h"
00022 #include "CLHEP/Matrix/SymMatrix.h"
00023 #include "CLHEP/Matrix/Matrix.h"
00024 #include "GaudiKernel/MsgStream.h"
00025 #include "GaudiKernel/AlgFactory.h"
00026 #include "GaudiKernel/ISvcLocator.h"
00027 #include "GaudiKernel/IDataProviderSvc.h"
00028 #include "GaudiKernel/PropertyMgr.h"
00029 #include "GaudiKernel/IMessageSvc.h"
00030 #include "GaudiKernel/IDataManagerSvc.h"
00031 #include "GaudiKernel/SmartDataPtr.h"
00032 #include "GaudiKernel/PropertyMgr.h"
00033 #include "GaudiKernel/IJobOptionsSvc.h"
00034 #include "GaudiKernel/Bootstrap.h"
00035 #include "GaudiKernel/StatusCode.h"
00036 #include "GaudiKernel/ContainedObject.h"
00037 #include "GaudiKernel/SmartRef.h"
00038 #include "GaudiKernel/SmartRefVector.h"
00039 #include "GaudiKernel/ObjectVector.h"
00040
00041
00042 #include "KalFitAlg/helix/Helix.h"
00043
00044 #include "KalFitAlg/lpav/Lpav.h"
00045
00046 #include "KalFitAlg/coil/Bfield.h"
00047 #include "KalFitAlg/KalFitTrack.h"
00048 #include "KalFitAlg/KalFitHitMdc.h"
00049 #include "KalFitAlg/KalFitHelixSeg.h"
00050 #include "KalFitAlg/KalFitElement.h"
00051 #include "KalFitAlg/KalFitAlg.h"
00052 #include "McTruth/McParticle.h"
00053 #include "EventModel/EventHeader.h"
00054 #include "EvTimeEvent/RecEsTime.h"
00055 #include "ReconEvent/ReconEvent.h"
00056 #include "MdcRawEvent/MdcDigi.h"
00057 #include "MdcRecEvent/RecMdcHit.h"
00058 #include "MdcRecEvent/RecMdcTrack.h"
00059 #include "MdcRecEvent/RecMdcKalHelixSeg.h"
00060 #include "MdcRecEvent/RecMdcKalTrack.h"
00061 #include "MdcGeomSvc/MdcGeomSvc.h"
00062 #include "MagneticField/IMagneticFieldSvc.h"
00063 #include "MagneticField/MagneticFieldSvc.h"
00064
00065 #include "VertexFit/IVertexDbSvc.h"
00066
00067 #include "Identifier/Identifier.h"
00068 #include "Identifier/MdcID.h"
00069 #include "GaudiKernel/IPartPropSvc.h"
00070 #include "GaudiKernel/INTupleSvc.h"
00071 using CLHEP::HepVector;
00072 using CLHEP::Hep3Vector;
00073 using CLHEP::HepMatrix;
00074 using CLHEP::HepSymMatrix;
00075
00076 using namespace Event;
00077 using namespace KalmanFit;
00078
00079
00080
00081
00082 const double KalFitAlg::RIW = 6.35;
00083
00085 KalFitAlg::KalFitAlg(const std::string& name, ISvcLocator* pSvcLocator):
00086 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc_(0),m_MFSvc_(0),
00087 _wire(0), _layer(0), _superLayer(0),
00088 pathl_(1), wsag_(4), back_(1), pT_(0.0), lead_(2), mhyp_(31),
00089 pe_cut_(2.0), pmu_cut_(2.0), ppi_cut_(2.0), pk_cut_(2.0), pp_cut_(2.0),
00090 muls_(1), loss_(1), enhance_(0),drifttime_choice_(0),choice_(0),
00091 fac_h1_(1),fac_h2_(1),fac_h3_(1),fac_h4_(1),fac_h5_(1),
00092 i_back_(-1), debug_kft_(0), debug_(0), ntuple_(0),eventno(-1),
00093 Tds_back_no(0),m_nt1(0),m_nt2(0),m_nt3(0),m_nt4(0),m_nt5(0),
00094 iqual_back_(1),tprop_(1),
00095 dchi2cut_inner_(0),dchi2cut_outer_(0),
00096 dchi2cut_mid1_(0),dchi2cut_mid2_(0),
00097 dchi2cut_layid2_(0),dchi2cut_layid3_(0),m_usevtxdb(0),m_dangcut(10),m_dphicut(10),
00098 nTotalTrks(0)
00099 {
00100 declareProperty("dchi2cut_layid2",dchi2cut_layid2_ = 10);
00101 declareProperty("dchi2cut_layid3",dchi2cut_layid3_ = 10);
00102 declareProperty("dchi2cut_inner",dchi2cut_inner_ = 10);
00103 declareProperty("dchi2cut_mid1",dchi2cut_mid1_ = 10);
00104 declareProperty("dchi2cut_mid2",dchi2cut_mid2_ = 10);
00105 declareProperty("dchi2cut_outer",dchi2cut_outer_ = 10);
00106 declareProperty("gain1",gain1_ = 1.);
00107 declareProperty("gain2",gain2_ = 1.);
00108 declareProperty("gain3",gain3_ = 1.);
00109 declareProperty("gain4",gain4_ = 1.);
00110 declareProperty("gain5",gain5_ = 1.);
00111 declareProperty("fitnocut",fitnocut_ = 5);
00112 declareProperty("inner_steps",inner_steps_ = 3);
00113 declareProperty("outer_steps",outer_steps_ = 3);
00114 declareProperty("choice",choice_ = 0);
00115 declareProperty("numfcor",numfcor_ = 0);
00116 declareProperty("numf",numf_ = 0);
00117 declareProperty("steplev",steplev_ = 2);
00118 declareProperty("usage",usage_=0);
00119 declareProperty("i_front",i_front_=1);
00120 declareProperty("lead",lead_=2);
00121 declareProperty("muls",muls_=1);
00122 declareProperty("loss",loss_=1);
00123 declareProperty("matrixg",matrixg_=100.0);
00124 declareProperty("lr",lr_=1);
00125 declareProperty("debug_kft",debug_kft_=0);
00126 declareProperty("debug",debug_=0);
00127 declareProperty("ntuple",ntuple_=0);
00128 declareProperty("activeonly",activeonly_=0);
00129 declareProperty("matfile",matfile_="geomdc_material.dat");
00130 declareProperty("cylfile",cylfile_="geomdc_cylinder.dat");
00131 declareProperty("dchi2cutf",dchi2cutf_=1000.0);
00132 declareProperty("dchi2cuts",dchi2cuts_=1000.0);
00133 declareProperty("pt",pT_=0.0);
00134 declareProperty("pe_cut",pe_cut_=2.0);
00135 declareProperty("pmu_cut",pmu_cut_=2.0);
00136 declareProperty("ppi_cut",ppi_cut_=2.0);
00137 declareProperty("pk_cut",pk_cut_=2.0);
00138 declareProperty("pp_cut",pp_cut_=2.0);
00139 declareProperty("wsag",wsag_=4);
00140 declareProperty("back",back_=1);
00141 declareProperty("i_back",i_back_=1);
00142 declareProperty("tofflag",tofflag_=1);
00143 declareProperty("tof_hyp",tof_hyp_=1);
00144 declareProperty("resolution",resolution_=1);
00145 declareProperty("fstrag",fstrag_=0.9);
00146 declareProperty("drifttime_choice",drifttime_choice_=0);
00147 declareProperty("tprop",tprop_=1);
00148 declareProperty("pt_cut",pt_cut_= 0.2);
00149 declareProperty("theta_cut",theta_cut_= 0.8);
00150 declareProperty("usevtxdb",m_usevtxdb= 0);
00151 declareProperty("cosmicflag",m_csmflag= 0);
00152 declareProperty("dangcut",m_dangcut=10.);
00153 declareProperty("dphicut",m_dphicut=10.);
00154
00155 for(int i=0; i<5; i++) nFailedTrks[i]=0;
00156 }
00157
00158
00159 KalFitAlg::~KalFitAlg(void)
00160 {
00161 MsgStream log(msgSvc(), name());
00162 log << MSG::INFO <<" Start cleaning, delete Mdc geometry objects" << endreq;
00163 clean();
00164 log << MSG::INFO << "End cleaning " << endreq;
00165 }
00166
00167 void KalFitAlg::clean(void)
00168 {
00169
00170 _BesKalmanFitWalls.clear();
00171 _BesKalmanFitMaterials.clear();
00172 if(_wire)free(_wire);
00173 if(_layer)free(_layer);
00174 if(_superLayer)free(_superLayer);
00175 }
00176
00177
00178 StatusCode KalFitAlg::initialize()
00179 {
00180 MsgStream log(msgSvc(), name());
00181 log << MSG::INFO << "in initize()"
00182 << "KalFit> Initialization for current run " << endreq;
00183 log << MSG::INFO << "Present Parameters: muls: " << muls_ <<" loss: "
00184 << loss_ <<" matrixg "<< matrixg_ <<" lr "<< lr_
00185 << " debug " << debug_ << " ntuple " << ntuple_
00186 << " activeonly "<< activeonly_ << endreq;
00187
00188 KalFitTrack::LR(lr_);
00189 KalFitTrack::resol(resolution_);
00190 KalFitTrack::Tof_correc_ = tofflag_;
00191 KalFitTrack::tofall_ = tof_hyp_;
00192 KalFitTrack::chi2_hitf_ = dchi2cutf_;
00193 KalFitTrack::chi2_hits_ = dchi2cuts_;
00194 KalFitTrack::factor_strag_ = fstrag_;
00195 KalFitTrack::debug_ = debug_kft_;
00196 KalFitTrack::drifttime_choice_ = drifttime_choice_;
00197 KalFitTrack::steplev_ = steplev_;
00198 KalFitTrack::numfcor_ = numfcor_;
00199 KalFitTrack::inner_steps_ = inner_steps_;
00200 KalFitTrack::outer_steps_ = outer_steps_;
00201 KalFitTrack::tprop_ = tprop_;
00202
00203 setDchisqCut();
00204
00205 clean();
00206 log << MSG::INFO << ".....building Mdc " << endreq;
00207
00208
00209
00210
00211
00212 setBesFromGdml();
00213
00214 setCalibSvc_init();
00215
00216
00217
00218
00219
00220
00221
00222 hist_def();
00223
00224
00225 IMagneticFieldSvc* IMFSvc;
00226 StatusCode sc = service ("MagneticFieldSvc",IMFSvc);
00227 if(sc!=StatusCode::SUCCESS) {
00228 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00229 }
00230 KalFitTrack::setMagneticFieldSvc(IMFSvc);
00231
00232
00233 if (KalFitTrack::numfcor_){
00234 KalFitTrack::Bznom_ = (IMFSvc->getReferField())*10000;
00235 if(0 == KalFitTrack::Bznom_) KalFitTrack::Bznom_ = -10;
00236
00237 if(4 == debug_){
00238 std::cout<<" initialize, referField from MagneticFieldSvc: "<< (IMFSvc->getReferField())*10000 <<std::endl;
00239 std::cout<<" magnetic field: "<<KalFitTrack::Bznom_<<std::endl;
00240 }
00241
00242 }
00243
00244
00245 if (ntuple_)
00246 log << MSG::INFO <<" ntuple out, the option is "<< ntuple_ <<endreq;
00247 if (debug_ >0 ) {
00248 cout << "KalFitAlg> DEBUG open,Here is the important Parameters :\n";
00249 cout << " Leading particule with mass hyp = " << lead_ << std::endl;
00250 cout << " mhyp = " << mhyp_ << std::endl;
00251 cout << "===== Effects taking into account : " << std::endl;
00252 cout << " - multiple scattering = " << muls_ << std::endl;
00253 cout << " - energy loss = " << loss_ << std::endl;
00254 if (KalFitTrack::strag_)
00255 cout << " - straggling for the energy loss " << std::endl;
00256 cout << " - nominal Bz value = " << KalFitTrack::Bznom_ << std::endl;
00257
00258 if (KalFitTrack::numf_ > 19)
00259 cout << " - non uniform magnetic field treatment "
00260 << KalFitTrack::numf_ << std::endl;
00261 cout << " - wire sag correction = " << wsag_ << std::endl;
00262 cout << " - Tof correction = " << KalFitTrack::Tof_correc_ << std::endl;
00263 cout << " - chi2 cut for each hit = " << KalFitTrack::chi2_hitf_
00264 << std::endl << " is applied after " << KalFitTrack::nmdc_hit2_
00265 << " hits included " << std::endl;
00266
00267 if (back_){
00268 cout << " Backward filter is on with a pT cut value = " << pT_ << endl;
00269 }
00270 if(debug_ == 4) cout << " pathl = " << pathl_ << std::endl;
00271
00272 if (KalFitTrack::LR_==1)
00273 cout << " Decision L/R from MdcRecHit " << std::endl;
00274 }
00275
00276 KalFitElement::muls(muls_);
00277 KalFitElement::loss(loss_);
00278 KalFitTrack::lead(lead_);
00279 KalFitTrack::back(back_);
00280 KalFitTrack::numf(numf_);
00281
00282 IPartPropSvc* p_PartPropSvc;
00283 static const bool CREATEIFNOTTHERE(true);
00284 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00285 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
00286 log << MSG::WARNING << " Could not initialize Particle Properties Service" << endreq;
00287 return StatusCode::SUCCESS;
00288 }
00289 m_particleTable = p_PartPropSvc->PDT();
00290
00291 return StatusCode::SUCCESS;
00292 }
00293
00294 StatusCode KalFitAlg::finalize()
00295 {
00296 MsgStream log(msgSvc(), name());
00297 log << MSG::DEBUG<<"KalFitAlg:: nTotalTrks = "<<nTotalTrks<<endreq;
00298 log << MSG::DEBUG<<" e: "<<nFailedTrks[0]<<" failed, "<<nTotalTrks-nFailedTrks[0]<<" successed"<<endreq;
00299 log << MSG::DEBUG<<" mu: "<<nFailedTrks[1]<<" failed, "<<nTotalTrks-nFailedTrks[1]<<" successed"<<endreq;
00300 log << MSG::DEBUG<<" pi: "<<nFailedTrks[2]<<" failed, "<<nTotalTrks-nFailedTrks[2]<<" successed"<<endreq;
00301 log << MSG::DEBUG<<" K: "<<nFailedTrks[3]<<" failed, "<<nTotalTrks-nFailedTrks[3]<<" successed"<<endreq;
00302 log << MSG::DEBUG<<" p: "<<nFailedTrks[4]<<" failed, "<<nTotalTrks-nFailedTrks[4]<<" successed"<<endreq;
00303 return StatusCode::SUCCESS;
00304 }
00305
00306
00307 StatusCode KalFitAlg::beginRun( )
00308 {
00309 MsgStream log(msgSvc(), name());
00310 log << MSG::INFO << "in beginRun()" << endreq;
00311 log << MSG::INFO << "Present Parameters: muls: " << muls_ <<" loss: "
00312 << " activeonly "<< activeonly_ << endreq;
00313
00314
00315 setGeomSvc_init();
00316
00317 set_Mdc();
00318
00319 IMagneticFieldSvc* IMFSvc;
00320 StatusCode sc = service ("MagneticFieldSvc",IMFSvc);
00321 if(sc!=StatusCode::SUCCESS) {
00322 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00323 }
00324
00325
00326 if (KalFitTrack::numfcor_){
00327 KalFitTrack::Bznom_ = (IMFSvc->getReferField())*10000;
00328 if(0 == KalFitTrack::Bznom_) KalFitTrack::Bznom_ = -10;
00329
00330 if(4 == debug_){
00331 std::cout<<" beginRun, referField from MagneticFieldSvc:"<< (IMFSvc->getReferField())*10000 <<std::endl;
00332 std::cout<<" magnetic field: "<<KalFitTrack::Bznom_<<std::endl;
00333 }
00334 }
00335
00336 return StatusCode::SUCCESS;
00337 }
00338
00339
00340
00341
00342
00343 void KalFitAlg::hist_def ( void )
00344 {
00345 if(ntuple_&1) {
00346 NTuplePtr nt1(ntupleSvc(),"FILE104/n101");
00347 StatusCode status;
00348 if ( nt1 ) m_nt1 = nt1;
00349 else {
00350 m_nt1= ntupleSvc()->book("FILE104/n101",CLID_ColumnWiseTuple,"KalFit");
00351 if ( m_nt1 ) {
00352
00353 status = m_nt1->addItem("trackid",m_trackid);
00354 status = m_nt1->addItem("stat",5,2,m_stat);
00355 status = m_nt1->addItem("ndf",5,2,m_ndf);
00356 status = m_nt1->addItem("chisq",5,2,m_chisq);
00357 status = m_nt1->addItem("length",5,m_length);
00358 status = m_nt1->addItem("tof",5,m_tof);
00359 status = m_nt1->addItem("nhits",5,m_nhits);
00360 status = m_nt1->addItem("zhelix",5,m_zhelix);
00361 status = m_nt1->addItem("zhelixe",5,m_zhelixe);
00362 status = m_nt1->addItem("zhelixmu",5,m_zhelixmu);
00363 status = m_nt1->addItem("zhelixk",5,m_zhelixk);
00364 status = m_nt1->addItem("zhelixp",5,m_zhelixp);
00365 status = m_nt1->addItem("zptot",m_zptot);
00366 status = m_nt1->addItem("zptote",m_zptote);
00367 status = m_nt1->addItem("zptotmu",m_zptotmu);
00368 status = m_nt1->addItem("zptotk",m_zptotk);
00369 status = m_nt1->addItem("zptotp",m_zptotp);
00370
00371 status = m_nt1->addItem("zpt",m_zpt);
00372 status = m_nt1->addItem("zpte",m_zpte);
00373 status = m_nt1->addItem("zptmu",m_zptmu);
00374 status = m_nt1->addItem("zptk",m_zptk);
00375 status = m_nt1->addItem("zptp",m_zptp);
00376
00377 status = m_nt1->addItem("fptot",m_fptot);
00378 status = m_nt1->addItem("fptote",m_fptote);
00379 status = m_nt1->addItem("fptotmu",m_fptotmu);
00380 status = m_nt1->addItem("fptotk",m_fptotk);
00381 status = m_nt1->addItem("fptotp",m_fptotp);
00382 status = m_nt1->addItem("fpt",m_fpt);
00383 status = m_nt1->addItem("fpte",m_fpte);
00384 status = m_nt1->addItem("fptmu",m_fptmu);
00385 status = m_nt1->addItem("fptk",m_fptk);
00386 status = m_nt1->addItem("fptp",m_fptp);
00387
00388 status = m_nt1->addItem("lptot",m_lptot);
00389 status = m_nt1->addItem("lptote",m_lptote);
00390 status = m_nt1->addItem("lptotmu",m_lptotmu);
00391 status = m_nt1->addItem("lptotk",m_lptotk);
00392 status = m_nt1->addItem("lptotp",m_lptotp);
00393 status = m_nt1->addItem("lpt",m_lpt);
00394 status = m_nt1->addItem("lpte",m_lpte);
00395 status = m_nt1->addItem("lptmu",m_lptmu);
00396 status = m_nt1->addItem("lptk",m_lptk);
00397 status = m_nt1->addItem("lptp",m_lptp);
00398
00399 status = m_nt1->addItem("zsigp",m_zsigp);
00400 status = m_nt1->addItem("zsigpe",m_zsigpe);
00401 status = m_nt1->addItem("zsigpmu",m_zsigpmu);
00402 status = m_nt1->addItem("zsigpk",m_zsigpk);
00403 status = m_nt1->addItem("zsigpp",m_zsigpp);
00404 status = m_nt1->addItem("fhelix",5,m_fhelix);
00405 status = m_nt1->addItem("fhelixe",5,m_fhelixe);
00406 status = m_nt1->addItem("fhelixmu",5,m_fhelixmu);
00407 status = m_nt1->addItem("fhelixk",5,m_fhelixk);
00408 status = m_nt1->addItem("fhelixp",5,m_fhelixp);
00409 status = m_nt1->addItem("lhelix",5,m_lhelix);
00410 status = m_nt1->addItem("lhelixe",5,m_lhelixe);
00411 status = m_nt1->addItem("lhelixmu",5,m_lhelixmu);
00412 status = m_nt1->addItem("lhelixk",5,m_lhelixk);
00413 status = m_nt1->addItem("lhelixp",5,m_lhelixp);
00414 if(ntuple_&32) {
00415 status = m_nt1->addItem("zerror",15,m_zerror);
00416 status = m_nt1->addItem("zerrore",15,m_zerrore);
00417 status = m_nt1->addItem("zerrormu",15,m_zerrormu);
00418 status = m_nt1->addItem("zerrork",15,m_zerrork);
00419 status = m_nt1->addItem("zerrorp",15,m_zerrorp);
00420 status = m_nt1->addItem("ferror",15,m_ferror);
00421 status = m_nt1->addItem("ferrore",15,m_ferrore);
00422 status = m_nt1->addItem("ferrormu",15,m_ferrormu);
00423 status = m_nt1->addItem("ferrork",15,m_ferrork);
00424 status = m_nt1->addItem("ferrorp",15,m_ferrorp);
00425 status = m_nt1->addItem("lerror",15,m_lerror);
00426 status = m_nt1->addItem("lerrore",15,m_lerrore);
00427 status = m_nt1->addItem("lerrormu",15,m_lerrormu);
00428 status = m_nt1->addItem("lerrork",15,m_lerrork);
00429 status = m_nt1->addItem("lerrorp",15,m_lerrorp);
00430 }
00431 if((ntuple_&16)&&(ntuple_&1)) {
00432 status = m_nt1->addItem("evtid",m_evtid);
00433 status = m_nt1->addItem("mchelix",5,m_mchelix);
00434 status = m_nt1->addItem("mcptot",m_mcptot);
00435 status = m_nt1->addItem("mcpid",m_mcpid);
00436 }
00437 if( status.isFailure() ) cout<<"Ntuple1 add item failed!"<<endl;
00438 }
00439 }
00440 }
00441
00442 if(ntuple_&4) {
00443 NTuplePtr nt2(ntupleSvc(),"FILE104/n102");
00444 StatusCode status2;
00445 if ( nt2 ) m_nt2 = nt2;
00446 else {
00447 m_nt2= ntupleSvc()->book("FILE104/n102",CLID_ColumnWiseTuple,"KalFitComp");
00448 if ( m_nt2 ) {
00449 status2 = m_nt2->addItem("delx",m_delx);
00450 status2 = m_nt2->addItem("dely",m_dely);
00451 status2 = m_nt2->addItem("delz",m_delz);
00452 status2 = m_nt2->addItem("delthe",m_delthe);
00453 status2 = m_nt2->addItem("delphi",m_delphi);
00454 status2 = m_nt2->addItem("delp",m_delp);
00455 status2 = m_nt2->addItem("delpx",m_delpx);
00456 status2 = m_nt2->addItem("delpy",m_delpy);
00457 status2 = m_nt2->addItem("delpz",m_delpz);
00458
00459 if( status2.isFailure() ) cout<<"Ntuple2 add item failed!"<<endl;
00460 }
00461 }
00462 }
00463
00464 if(ntuple_&2) {
00465 NTuplePtr nt3(ntupleSvc(),"FILE104/n103");
00466 StatusCode status3;
00467 if ( nt3 ) m_nt3 = nt3;
00468 else {
00469 m_nt3= ntupleSvc()->book("FILE104/n103",CLID_ColumnWiseTuple,"PatRec");
00470 if ( m_nt3 ) {
00471 status3 = m_nt3->addItem("trkhelix",5,m_trkhelix);
00472 status3 = m_nt3->addItem("trkptot",m_trkptot);
00473 if(ntuple_&32) {
00474 status3 = m_nt3->addItem("trkerror",15,m_trkerror);
00475 status3 = m_nt3->addItem("trksigp",m_trksigp);
00476 }
00477 status3 = m_nt3->addItem("trkndf",m_trkndf);
00478 status3 = m_nt3->addItem("trkchisq",m_trkchisq);
00479 if( status3.isFailure() ) cout<<"Ntuple3 add item failed!"<<endl;
00480 }
00481 }
00482 }
00483
00484 if(ntuple_&4) {
00485 NTuplePtr nt4(ntupleSvc(),"FILE104/n104");
00486 StatusCode status4;
00487 if ( nt4 ) m_nt4 = nt4;
00488 else {
00489 m_nt4= ntupleSvc()->book("FILE104/n104",CLID_ColumnWiseTuple,"PatRecComp");
00490 if ( m_nt4 ) {
00491 status4 = m_nt4->addItem("trkdelx",m_trkdelx);
00492 status4 = m_nt4->addItem("trkdely",m_trkdely);
00493 status4 = m_nt4->addItem("trkdelz",m_trkdelz);
00494 status4 = m_nt4->addItem("trkdelthe",m_trkdelthe);
00495 status4 = m_nt4->addItem("trkdelphi",m_trkdelphi);
00496 status4 = m_nt4->addItem("trkdelp",m_trkdelp);
00497 if( status4.isFailure() ) cout<<"Ntuple4 add item failed!"<<endl;
00498 }
00499 }
00500 }
00501 if(ntuple_&8) {
00502 NTuplePtr nt5(ntupleSvc(), "FILE104/n105");
00503 StatusCode status5;
00504 if ( nt5 ) m_nt5 = nt5;
00505 else {
00506 m_nt5= ntupleSvc()->book("FILE104/n105",CLID_ColumnWiseTuple,"KalFitdChisq");
00507 if ( m_nt5 ) {
00508 status5 = m_nt5->addItem("dchi2",m_dchi2);
00509 status5 = m_nt5->addItem("masshyp",m_masshyp);
00510 status5 = m_nt5->addItem("residual_estim",m_residest);
00511 status5 = m_nt5->addItem("residual",m_residnew);
00512 status5 = m_nt5->addItem("layer",m_layer);
00513 status5 = m_nt5->addItem("kaldr",m_anal_dr);
00514 status5 = m_nt5->addItem("kalphi0",m_anal_phi0);
00515 status5 = m_nt5->addItem("kalkappa",m_anal_kappa);
00516 status5 = m_nt5->addItem("kaldz",m_anal_dz);
00517 status5 = m_nt5->addItem("kaltanl",m_anal_tanl);
00518 status5 = m_nt5->addItem("dr_ea",m_anal_ea_dr);
00519 status5 = m_nt5->addItem("phi0_ea",m_anal_ea_phi0);
00520 status5 = m_nt5->addItem("kappa_ea",m_anal_ea_kappa);
00521 status5 = m_nt5->addItem("dz_ea",m_anal_ea_dz);
00522 status5 = m_nt5->addItem("tanl_ea",m_anal_ea_tanl);
00523 if( status5.isFailure() ) cout<<"Ntuple5 add item failed!"<<endl;
00524 }
00525 }
00526 NTuplePtr nt6(ntupleSvc(),"FILE104/n106");
00527 StatusCode status6;
00528 if ( nt6 ) m_nt6 = nt6;
00529 else {
00530 m_nt6= ntupleSvc()->book("FILE104/n106",CLID_ColumnWiseTuple,"kal seg");
00531 if ( m_nt6 ) {
00532 status6 = m_nt6->addItem("docaInc",m_docaInc);
00533 status6 = m_nt6->addItem("docaExc",m_docaExc);
00534 status6 = m_nt6->addItem("tdr",m_tdrift);
00535 status6 = m_nt6->addItem("layerid", m_layerid);
00536 status6 = m_nt6->addItem("event", m_eventNo);
00537 status6 = m_nt6->addItem("residualInc",m_residualInc);
00538 status6 = m_nt6->addItem("residualExc",m_residualExc);
00539 status6 = m_nt6->addItem("lr",m_lr);
00540 status6 = m_nt6->addItem("dd",m_dd);
00541 status6 = m_nt6->addItem("y",m_yposition);
00542
00543 if( status6.isFailure() ) cout<<"Ntuple6 add item failed!"<<endl;
00544 }
00545 }
00546 }
00547 }
00548
00549
00550
00551 void KalFitAlg::setDchisqCut()
00552 {
00553 int layid = 0;
00554
00556 for(layid = 0; layid<2; layid++) {
00557 KalFitTrack::dchi2cutf_anal[layid] = dchi2cut_inner_;
00558 }
00559 KalFitTrack::dchi2cutf_anal[2] = dchi2cut_layid2_;
00560 KalFitTrack::dchi2cutf_anal[3] = dchi2cut_layid3_;
00561 for(layid = 4; layid<12; layid++) {
00562 KalFitTrack::dchi2cutf_anal[layid] = dchi2cut_mid1_;
00563 }
00564 for(layid = 12; layid<20; layid++) {
00565 KalFitTrack::dchi2cutf_anal[layid] = dchi2cut_mid2_;
00566 }
00567 for(layid = 20; layid<43; layid++) {
00568 KalFitTrack::dchi2cutf_anal[layid] = dchi2cut_outer_;
00569 }
00570
00571
00573 for(layid = 0; layid<2; layid++) {
00574 KalFitTrack::dchi2cuts_anal[layid] = dchi2cut_inner_;
00575 }
00576
00577 KalFitTrack::dchi2cuts_anal[2] = dchi2cut_layid2_;
00578 KalFitTrack::dchi2cuts_anal[3] = dchi2cut_layid3_;
00579 for(layid = 4; layid<12; layid++) {
00580 KalFitTrack::dchi2cuts_anal[layid] = dchi2cut_mid1_;
00581 }
00582 for(layid = 12; layid<20; layid++) {
00583 KalFitTrack::dchi2cuts_anal[layid] = dchi2cut_mid2_;
00584 }
00585 for(layid = 20; layid<43; layid++) {
00586 KalFitTrack::dchi2cuts_anal[layid] = dchi2cut_outer_;
00587 }
00588
00590
00591
00592
00593
00594
00596 for(layid = 0; layid<2; layid++) {
00597 KalFitTrack::dchi2cutf_calib[layid] = dchi2cut_inner_;
00598 }
00599
00600 KalFitTrack::dchi2cutf_calib[2] = dchi2cut_layid2_;
00601 KalFitTrack::dchi2cutf_calib[3] = dchi2cut_layid3_;
00602
00603 for(layid = 4; layid<12; layid++) {
00604 KalFitTrack::dchi2cutf_calib[layid] = dchi2cut_mid1_;
00605 }
00606
00607 for(layid = 12; layid<20; layid++) {
00608 KalFitTrack::dchi2cutf_calib[layid] = dchi2cut_mid2_;
00609 }
00610
00611 for(layid = 20; layid<43; layid++) {
00612 KalFitTrack::dchi2cutf_calib[layid] = dchi2cut_outer_;
00613 }
00614
00616 if(usage_<2){
00617 for(layid = 0; layid<43; layid++) {
00618 KalFitTrack::dchi2cutf_calib[layid] = 10.0;
00619 }
00620 }
00621
00622
00624 for(layid = 0; layid<2; layid++) {
00625 KalFitTrack::dchi2cuts_calib[layid] = dchi2cut_inner_;
00626 }
00627
00628 KalFitTrack::dchi2cuts_calib[2] = dchi2cut_layid2_;
00629 KalFitTrack::dchi2cuts_calib[3] = dchi2cut_layid3_;
00630
00631 for(layid = 4; layid<12; layid++) {
00632 KalFitTrack::dchi2cuts_calib[layid] = dchi2cut_mid1_;
00633 }
00634
00635 for(layid = 12; layid<20; layid++) {
00636 KalFitTrack::dchi2cuts_calib[layid] = dchi2cut_mid2_;
00637 }
00638
00639 for(layid = 20; layid<43; layid++) {
00640 KalFitTrack::dchi2cuts_calib[layid] = dchi2cut_outer_;
00641 }
00642
00644 if(usage_<2){
00645 for(layid = 0; layid<43; layid++) {
00646 KalFitTrack::dchi2cuts_calib[layid] = 10.0;
00647 }
00648 }
00649 }
00650
00651
00652
00653
00654 StatusCode KalFitAlg::execute()
00655 {
00656 MsgStream log(msgSvc(), name());
00657 log << MSG::INFO << "in execute()" << endreq;
00658
00659 for(int i=0; i<5; i++) iqual_front_[i] = 1;
00660 iqual_back_ = 1;
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00703 IMagneticFieldSvc* IMFSvc;
00704 StatusCode sc = service ("MagneticFieldSvc",IMFSvc);
00705 if(sc!=StatusCode::SUCCESS) {
00706 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00707 }
00708
00709
00710 if (KalFitTrack::numfcor_){
00711 KalFitTrack::Bznom_ = (IMFSvc->getReferField())*10000;
00712 if(0 == KalFitTrack::Bznom_) KalFitTrack::Bznom_ = -10;
00713
00714 if(4 == debug_){
00715 std::cout<<" execute, referField from MagneticFieldSvc: "<< (IMFSvc->getReferField())*10000 <<std::endl;
00716 std::cout<<" magnetic field: "<<KalFitTrack::Bznom_<<std::endl;
00717 }
00718 }
00719
00720 RecMdcKalTrackCol* kalcol = new RecMdcKalTrackCol;
00721
00722 IDataProviderSvc* evtSvc = NULL;
00723 Gaudi::svcLocator()->service("EventDataSvc", evtSvc);
00724 if (evtSvc) {
00725 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
00726 } else {
00727 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
00728 return StatusCode::SUCCESS;
00729 }
00730
00731 StatusCode kalsc;
00732 IDataManagerSvc *dataManSvc;
00733 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
00734 DataObject *aKalTrackCol;
00735 evtSvc->findObject("/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
00736 if(aKalTrackCol != NULL) {
00737 dataManSvc->clearSubTree("/Event/Recon/RecMdcKalTrackCol");
00738 evtSvc->unregisterObject("/Event/Recon/RecMdcKalTrackCol");
00739 }
00740
00741 kalsc = evtSvc->registerObject("/Event/Recon/RecMdcKalTrackCol", kalcol);
00742 if( kalsc.isFailure() ) {
00743 log << MSG::FATAL << "Could not register RecMdcKalTrack" << endreq;
00744 return StatusCode::SUCCESS;
00745 }
00746 log << MSG::INFO << "RecMdcKalTrackCol registered successfully!" <<endreq;
00747
00748 DataObject *aKalHelixSegCol;
00749 evtSvc->findObject("/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
00750 if(aKalHelixSegCol != NULL){
00751 dataManSvc->clearSubTree("/Event/Recon/RecMdcKalHelixSegCol");
00752 evtSvc->unregisterObject("/Event/Recon/RecMdcKalHelixSegCol");
00753 }
00754 RecMdcKalHelixSegCol *helixsegcol = new RecMdcKalHelixSegCol;
00755 kalsc = evtSvc->registerObject("/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
00756 if( kalsc.isFailure()){
00757 log<< MSG::FATAL << "Could not register RecMdcKalHelixSeg" <<endreq;
00758 return StatusCode::SUCCESS;
00759 }
00760 log << MSG::INFO << "RecMdcKalHelixSegCol register successfully!" <<endreq;
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770 MdcGeomSvc* const geosvc = dynamic_cast<MdcGeomSvc*>(imdcGeomSvc_);
00771 if(!geosvc) {
00772 std::cout<<"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
00773 }
00774
00775 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00776 if (!eventHeader) {
00777 log << MSG::WARNING << "Could not find Event Header" << endreq;
00778 return StatusCode::FAILURE;
00779 }
00780 int eventNo = eventHeader->eventNumber();
00781 int runNo = eventHeader->runNumber();
00782 if(runNo>0) wsag_=4;
00783 else wsag_=0;
00784
00785 double t0=0.;
00786 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00787 if (estimeCol && estimeCol->size()) {
00788 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
00789 t0 = (*iter_evt)->getTest();
00790
00791 }else{
00792 log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
00793 return StatusCode::SUCCESS;
00794 }
00795
00796
00797 if(debug_==4) {
00798 std::cout<<"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
00799 }
00800 KalFitTrack::setT0(t0);
00801
00802 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,"/Event/Digi/MdcDigiCol");
00803 if (sc!=StatusCode::SUCCESS) {
00804 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
00805 return StatusCode::SUCCESS;
00806 }
00807 KalFitTrack::setMdcDigiCol(mdcDigiCol);
00808
00809
00810
00811 if((ntuple_&16)&&(ntuple_&1)) {
00812
00813
00814
00815
00816 m_evtid = eventHeader->eventNumber();
00817 bool mcstat = true;
00818
00819 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),"/Event/MC/McParticleCol");
00820 if (!mcPartCol) {
00821 log << MSG::WARNING << "Could not find McParticle" << endreq;
00822 mcstat = false;
00823 }
00824
00825 if(mcstat) {
00826 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
00827 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
00828 if(!(*i_mcTrk)->primaryParticle()) continue;
00829 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
00830 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
00831 log << MSG::DEBUG << "MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
00832 << " theta=" << mom.theta() <<" phi="<< mom.phi()
00833 <<" px="<< mom.px() <<" py="<< mom.py() <<" pz="<< mom.pz()
00834 << endreq;
00835 double charge = 0.0;
00836 int pid = (*i_mcTrk)->particleProperty();
00837 if( pid >0 ) {
00838 charge = m_particleTable->particle( pid )->charge();
00839 } else if ( pid <0 ) {
00840 charge = m_particleTable->particle( -pid )->charge();
00841 charge *= -1;
00842 } else {
00843 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00844 }
00845 HepPoint3D pos2(pos.x(),pos.y(),pos.z());
00846 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
00847
00848 Helix mchelix(pos2, mom2, charge);
00849 log << MSG::DEBUG << "charge of the track " << charge << endreq;
00850 if( debug_ == 4) cout<< "helix: "<<mchelix.a()<<endl;
00851 mchelix.pivot( HepPoint3D(0,0,0) );
00852 for( int j =0; j<5; j++) {
00853 m_mchelix[j] = mchelix.a()[j];
00854 }
00855 m_mcpid = pid;
00856 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
00857 }
00858 }
00859 }
00860
00861 Identifier mdcid;
00862
00863
00864 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00865 if (!newtrkCol) {
00866 log << MSG::FATAL << "Could not find RecMdcTrackCol" << endreq;
00867 return( StatusCode::SUCCESS);
00868 }
00869 log << MSG::INFO << "Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
00870
00871 vector<MdcRec_trk>* mtrk_mgr = MdcRecTrkCol::getMdcRecTrkCol();
00872 mtrk_mgr->clear();
00873 vector<MdcRec_trk_add>* mtrkadd_mgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
00874 mtrkadd_mgr->clear();
00875 vector<MdcRec_wirhit>* mhit_mgr = MdcRecWirhitCol::getMdcRecWirhitCol();
00876 mhit_mgr->clear();
00877
00878 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
00879 Hep3Vector csmp3[2];
00880 double csmphi[2];
00881 int status_temp=0;
00882 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
00883 for(int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
00884 if(kj<3){
00885 csmp3[kj-1]=(*iter_trk)->p3();
00886 csmphi[kj-1] = (*iter_trk)->phi();
00887 }
00888 if(ntuple_&2) {
00889
00890 for( int j = 0, ij = 0; j<5; j++) {
00891 m_trkhelix[j] = (*iter_trk)->helix()[j];
00892 if(ntuple_&32) {
00893 for(int k=0; k<=j; k++,ij++) {
00894 m_trkerror[ij] = (*iter_trk)->err()[j][k];
00895 }
00896 }
00897 }
00898 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
00899 if(ntuple_&32){
00900 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
00901 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
00902 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
00903 }
00904 m_trkndf = (*iter_trk)->ndof();
00905 m_trkchisq = (*iter_trk)->chi2();
00906
00907 if (debug_ == 4) cout<<"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
00908
00909 StatusCode sc3 = m_nt3->write();
00910 if( sc3.isFailure() ) cout<<"Ntuple3 filling failed!"<<endl;
00911 }
00912
00913 if(ntuple_&4) {
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937 }
00938
00939
00940 log << MSG::DEBUG << "retrieved MDC tracks:"
00941 << " Nhits " <<(*iter_trk)->getNhits()
00942 << " Nster " <<(*iter_trk)->nster() <<endreq;
00943
00944 HitRefVec gothits = (*iter_trk)->getVecHits();
00945
00946 MdcRec_trk* rectrk = new MdcRec_trk;
00947
00948 rectrk->id = (*iter_trk)->trackId();
00949 rectrk->chiSq = (*iter_trk)->chi2();
00950 rectrk->ndf = (*iter_trk)->ndof();
00951 rectrk->fiTerm = (*iter_trk)->getFiTerm();
00952 rectrk->nhits = (*iter_trk)->getNhits();
00953 rectrk->nster = (*iter_trk)->nster();
00954 rectrk->nclus = 0;
00955 rectrk->stat = (*iter_trk)->stat();
00956 status_temp = (*iter_trk)->stat();
00957 MdcRec_trk_add* trkadd = new MdcRec_trk_add;
00958 trkadd->id = (*iter_trk)->trackId();
00959 trkadd->quality = 0;
00960 trkadd->kind = 1;
00961 trkadd->decision = 0;
00962 trkadd->body = rectrk;
00963 rectrk->add = trkadd;
00964
00965 for ( int i=0; i<5; i++) {
00966 rectrk->helix[i] = (*iter_trk)->helix()[i];
00967 if( i<3 ) rectrk->pivot[i] = (*iter_trk)->getPivot()[i];
00968 for( int j = 1; j<i+2;j++) {
00969 rectrk->error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
00970 }
00971 }
00972 std::sort(gothits.begin(), gothits.end(), order_rechits);
00973 HitRefVec::iterator it_gothit = gothits.begin();
00974 for( ; it_gothit != gothits.end(); it_gothit++) {
00975
00976 if( (*it_gothit)->getStat() != 1 ) {
00977 if(activeonly_) {
00978 log<<MSG::WARNING<<"this hit is not used in helix fitting!"<<endreq;
00979 continue;
00980 }
00981 }
00982
00983 log << MSG::DEBUG << "retrieved hits in MDC tracks:"
00984 << " hits DDL " <<(*it_gothit)->getDriftDistLeft()
00985 << " hits DDR " <<(*it_gothit)->getDriftDistRight()
00986 << " error DDL " <<(*it_gothit)->getErrDriftDistLeft()
00987 << " error DDR " <<(*it_gothit)->getErrDriftDistRight()
00988 << " id of hit "<<(*it_gothit)->getId()
00989 << " track id of hit "<<(*it_gothit)->getTrkId()
00990 << " hits ADC " <<(*it_gothit)->getAdc() << endreq;
00991
00992 MdcRec_wirhit* whit = new MdcRec_wirhit;
00993 whit->id = (*it_gothit)->getId();
00994 whit->ddl = (*it_gothit)->getDriftDistLeft();
00995 whit->ddr = (*it_gothit)->getDriftDistRight();
00996 whit->erddl = (*it_gothit)->getErrDriftDistLeft();
00997 whit->erddr = (*it_gothit)->getErrDriftDistRight();
00998 whit->pChiSq = (*it_gothit)->getChisqAdd();
00999 whit->lr = (*it_gothit)->getFlagLR();
01000 whit->stat = (*it_gothit)->getStat();
01001 mdcid = (*it_gothit)->getMdcId();
01002 int layid = MdcID::layer(mdcid);
01003 int localwid = MdcID::wire(mdcid);
01004 int w0id = geosvc->Layer(layid)->Wirst();
01005 int wid = w0id + localwid;
01006 log << MSG::INFO
01007 << "lr from PR: "<<whit->lr
01008 << " layerId = " << layid
01009 << " wireId = " << localwid
01010 << endreq;
01011
01012 const MdcGeoWire * const wirgeo = geosvc->Wire(wid);
01013
01014
01015 whit->rechitptr = *it_gothit;
01016 whit->geo = wirgeo;
01017 whit->dat = 0;
01018 whit->trk = rectrk;
01019 whit->tdc = (*it_gothit)->getTdc();
01020 whit->adc= (*it_gothit)->getAdc();
01021 rectrk->hitcol.push_back(whit);
01022 mhit_mgr->push_back(*whit);
01023 }
01024 mtrk_mgr->push_back(*rectrk);
01025 mtrkadd_mgr->push_back(*trkadd);
01026
01027 delete rectrk;
01028 delete trkadd;
01029 }
01030
01031
01032 if(ntuple_&4) {
01033 m_trkdelx = trkx1 - trkx2;
01034 m_trkdely = trky1 - trky2;
01035 m_trkdelz = trkz1 - trkz2;
01036 m_trkdelthe = trkthe1 + trkthe2;
01037 m_trkdelphi = trkphi1- trkphi2;
01038 m_trkdelp = trkp1 - trkp2;
01039 StatusCode sc4 = m_nt4->write();
01040 if( sc4.isFailure() ) cout<<"Ntuple4 filling failed!"<<endl;
01041 }
01042
01043 if(debug_ == 4) { std::cout<<"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
01044 <<"********"<<mhit_mgr->size()<<"****"<<mtrkadd_mgr->size()<<endl;
01045 }
01046
01047
01048 if(usage_ == 0) kalman_fitting_anal();
01049 if(usage_ == 1) kalman_fitting_calib();
01050 double mdang = 180.0 - csmp3[0].angle(csmp3[1].unit())*180.0/M_PI;
01051 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/M_PI;
01052
01053 if(usage_ == 2 && (mtrk_mgr->size())==2 && fabs(mdang)<m_dangcut && fabs(mdphi)<m_dphicut) kalman_fitting_csmalign();
01054 if(usage_ == 3 && (mtrk_mgr->size())==1 && status_temp==-1) kalman_fitting_MdcxReco_Csmc_Sew();
01055
01056 log << MSG::DEBUG <<"after kalman_fitting(),but in execute...."<<endreq;
01057 clearTables();
01058
01059
01061
01062 MdcID mdcId;
01063 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
01064
01065
01066
01067 RecMdcKalTrack::setPidType(RecMdcKalTrack::electron);
01068 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
01069 int i_trk=0;
01070 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++){
01071
01072 for(int hypo=0; hypo<5; hypo++)
01073 {
01074 if((*KalTrk)->getStat(0,hypo)==1) nFailedTrks[hypo]++;
01075 }
01076 HelixSegRefVec gothelixsegs = (*KalTrk)->getVecHelixSegs();
01077 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
01078
01079 int nhitofthistrk=0;
01080 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
01081 nhitofthistrk++;
01082
01083
01084
01085
01086 }
01087 iter_hit=gothelixsegs.begin();
01088
01089
01090
01091 }
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135 return StatusCode::SUCCESS;
01136 }
01137
01138
01139 void KalFitAlg::fillTds(MdcRec_trk& TrasanTRK, KalFitTrack& track,
01140 RecMdcKalTrack* trk , int l_mass) {
01141
01142 HepPoint3D IP(0,0,0);
01143 track.pivot(IP);
01144
01145 int iqual(1);
01146 int trasster = TrasanTRK.nster, trakster = track.nster(),
01147 trasax(TrasanTRK.nhits-trasster), trakax(track.nchits()-trakster);
01148 if (TrasanTRK.nhits-track.nchits()>fitnocut_ ||
01149 TrasanTRK.helix[2]*track.a()[2]<0)
01150 iqual = 0;
01151
01152 if (debug_ == 4) {
01153 cout<< "trasster trakster trasax trakax TrasK trackK iqual"<<endl
01154 <<trasster<<" "<<trakster<<" "<<trasax<<" "<<trakax
01155 <<" "<<TrasanTRK.helix[2]<<" "<<track.a()[2]<<" "<<iqual<<endl;
01156 cout<<"FillTds> track.chiSq..."<<track.chiSq()<<" nchits "<<track.nchits()
01157 <<" nster "<<track.nster()<<" iqual "<<iqual<<" track.Ea "<< track.Ea()<<endl;
01158
01159 cout<<"fillTds>.....track.Ea[2][2] "<<track.Ea()[2][2]<<endl;
01160 cout << " TRASAN stereo = " << trasster
01161 << " and KalFitTrack = " << trakster << std::endl;
01162 cout << " TRASAN axial = " << trasax
01163 << " and KalFitTrack = " << trakax << std::endl;
01164
01165 if (!iqual) {
01166 cout << "...there is a problem during fit !! " << std::endl;
01167 if (trasster-trakster>5)
01168 cout << " because stereo " << trasster-trakster << std::endl;
01169 if (trasax-trakax >5)
01170 cout << " because axial " << std::endl;
01171 if (TrasanTRK.helix[2]*track.a()[2]<0)
01172 cout << " because kappa sign " << std::endl;
01173 }
01174 }
01175
01176 if (track.nchits() > 5 && track.nster() > 1 &&
01177 track.nchits()-track.nster() > 2 && track.chiSq() > 0 &&
01178 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01179 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01180 track.Ea()[4][4] > 0 && iqual) {
01181 if(debug_ == 4) cout<<"fillTds>.....going on "<<endl;
01182 trk->setStat(0,0,l_mass);
01183 trk->setMass(track.mass(),l_mass);
01184
01185
01186 trk->setChisq(track.chiSq(),0,l_mass);
01187 trk->setNdf(track.nchits()-5,0,l_mass);
01188 trk->setNhits(track.nchits(),l_mass);
01189
01190 trk->setFHelix(track.a(),l_mass);
01191 trk->setFError(track.Ea(),l_mass);
01192
01193 } else {
01194
01195 if(debug_) cout<<"ALARM: FillTds Not refit with KalFilter!!!"<<endl;
01196
01197 trk->setStat(1,0,l_mass);
01198 trk->setMass(KalFitTrack::mass(l_mass),l_mass);
01199
01200 trk->setChisq(TrasanTRK.chiSq,0,l_mass);
01201 trk->setNdf(TrasanTRK.nhits-5,0,l_mass);
01202
01203 trk->setNhits(TrasanTRK.nhits,l_mass);
01204 double a_trasan[5], ea_trasan[15];
01205 for( int i =0 ; i <5; i++){
01206 a_trasan[i] = TrasanTRK.helix[i];
01207 }
01208 for( int j =0 ; j <15; j++){
01209 ea_trasan[j] = TrasanTRK.error[j];
01210 }
01211 trk->setFHelix(a_trasan, l_mass);
01212 trk->setFError(ea_trasan,l_mass);
01213 }
01214 }
01215
01216
01217 void KalFitAlg::fillTds_lead(MdcRec_trk& TrasanTRK, KalFitTrack& track,
01218 RecMdcKalTrack* trk , int l_mass) {
01219
01220 HepPoint3D IP(0,0,0);
01221 track.pivot(IP);
01222
01223
01224 int trasster = TrasanTRK.nster, trakster = track.nster(),
01225 trasax(TrasanTRK.nhits-trasster), trakax(track.nchits()-trakster);
01226 if (TrasanTRK.nhits-track.nchits()>fitnocut_ ||
01227 TrasanTRK.helix[2]*track.a()[2]<0)
01228 iqual_front_[l_mass] = 0;
01229 if (debug_ == 4) {
01230
01231 cout<<"Nhit from PR "<<TrasanTRK.nhits<<" nhit "<<track.nchits()<<endl;
01232 cout<< "trasster trakster trasax trakax TrasK trackK iqual"<<endl
01233 <<trasster<<" "<<trakster<<" "<<trasax<<" "<<trakax
01234 <<" "<<TrasanTRK.helix[2]<<" "<<track.a()[2]<<" "<<iqual_front_[l_mass]<<endl;
01235 cout<<"FillTds_lead> track.chiSq..."<<track.chiSq()<<" nchits "<<track.nchits()
01236 <<" nster "<<track.nster()<<" iqual_front_[l_mass] "<<iqual_front_[l_mass]<<" track.Ea "<<track.Ea()<<endl;
01237
01238 cout << " TRASAN stereo = " << trasster
01239 << " and KalFitTrack = " << trakster << std::endl;
01240 cout << " TRASAN axial = " << trasax
01241 << " and KalFitTrack = " << trakax << std::endl;
01242
01243 if (!iqual_front_[l_mass]) {
01244 cout << "...there is a problem during fit !! " << std::endl;
01245 if (trasster-trakster>5)
01246 cout << " because stereo " << trasster-trakster << std::endl;
01247 if (trasax-trakax >5)
01248 cout << " because axial " << std::endl;
01249 if (TrasanTRK.helix[2]*track.a()[2]<0)
01250 cout << " because kappa sign " << std::endl;
01251 }
01252 }
01253
01254 if (track.nchits() > 5 && track.nster() > 1 &&
01255 track.nchits()-track.nster() > 2 && track.chiSq() > 0 &&
01256 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01257 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01258 track.Ea()[4][4] > 0 && iqual_front_[l_mass]) {
01259
01260 trk->setStat(0,0,l_mass);
01261 trk->setMass(track.mass(),l_mass);
01262 trk->setChisq(track.chiSq(),0,l_mass);
01263 trk->setNdf(track.nchits()-5,0,l_mass);
01264 trk->setNhits(track.nchits(),l_mass);
01265
01266 trk->setTrackId(TrasanTRK.id);
01267
01268 if (debug_ == 4) cout<<" trasan id...1 "<<TrasanTRK.id<<endl;
01269
01270 trk->setFHelix(track.a(),l_mass);
01271 trk->setFError(track.Ea(),l_mass);
01272
01273 } else {
01274
01275
01276
01277 if(debug_) cout<<"ALARM: FillTds_forMdc Not refit with KalFilter!!!"<<endl;
01278
01279 trk->setStat(1,0,l_mass);
01280 trk->setMass(KalFitTrack::mass(l_mass),l_mass);
01281
01282
01283 trk->setChisq(TrasanTRK.chiSq,0,l_mass);
01284 trk->setNdf(TrasanTRK.nhits-5,0,l_mass);
01285
01286 trk->setTrackId(TrasanTRK.id);
01287
01288 if (debug_ ==4) cout<<" trasan id...2 "<<TrasanTRK.id<<endl;
01289
01290
01291 trk->setNhits(TrasanTRK.nhits,l_mass);
01292 double a_trasan[5], ea_trasan[15];
01293 for( int i =0 ; i <5; i++){
01294 a_trasan[i] = TrasanTRK.helix[i];
01295 }
01296 for( int j =0 ; j <15; j++){
01297 ea_trasan[j] = TrasanTRK.error[j];
01298 }
01299 trk->setFHelix(a_trasan,l_mass);
01300 trk->setFError(ea_trasan,l_mass);
01301
01302
01303 }
01304 }
01305
01306
01307
01308
01309 void KalFitAlg::fillTds_ip(MdcRec_trk& TrasanTRK, KalFitTrack& track,
01310 RecMdcKalTrack* trk, int l_mass)
01311 {
01312 HepPoint3D IP(0,0,0);
01313 track.pivot(IP);
01314
01315 if (debug_ == 4&& l_mass==lead_) {
01316 cout << "fillTds_IP>......"<<endl;
01317 cout << " dr = " << track.a()[0]
01318 << ", Er_dr = " << sqrt(track.Ea()[0][0]) << std::endl;
01319 cout << " phi0 = " << track.a()[1]
01320 << ", Er_phi0 = " << sqrt(track.Ea()[1][1]) << std::endl;
01321 cout << " PT = " << 1/track.a()[2]
01322 << ", Er_kappa =" << sqrt(track.Ea()[2][2]) << std::endl;
01323 cout << " dz = " << track.a()[3]
01324 << ", Er_dz = " << sqrt(track.Ea()[3][3]) << std::endl;
01325 cout << " tanl = " << track.a()[4]
01326 << ", Er_tanl = " << sqrt(track.Ea()[4][4]) << std::endl;
01327 }
01328
01329 if (TrasanTRK.nhits-track.nchits()>fitnocut_ ||
01330 TrasanTRK.helix[2]*track.a()[2]<0)
01331 iqual_front_[l_mass] = 0;
01332
01333
01334 if (track.nchits() > 5 && track.nster() > 1 &&
01335 track.nchits()-track.nster() > 2 && track.chiSq() > 0 &&
01336 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01337 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01338 track.Ea()[4][4] > 0 && iqual_front_[l_mass]) {
01339
01340
01341 double dr = track.a()[0];
01342 double phi0 = track.a()[1];
01343 double kappa = track.a()[2];
01344 double dz = track.a()[3];
01345 double tanl = track.a()[4];
01346 int nLayer = track.nLayerUsed();
01347 trk->setNlayer(nLayer, l_mass);
01348
01349
01350 double vx = dr*cos(phi0);
01351 double vy = dr*sin(phi0);
01352 double vz = dz;
01353
01354
01355
01356 if(0==kappa) kappa = 10e-10;
01357 double px = -sin(phi0)/fabs(kappa);
01358 double py = cos(phi0)/fabs(kappa);
01359 double pz = tanl/fabs(kappa);
01360
01361 trk->setX(vx, l_mass);
01362 trk->setY(vy, l_mass);
01363 trk->setZ(vz, l_mass);
01364 trk->setPx(px, l_mass);
01365 trk->setPy(py, l_mass);
01366 trk->setPz(pz, l_mass);
01367
01368 const HepPoint3D poca(dr*cos(phi0),dr*sin(phi0),dz);
01369 trk->setPoca(poca,l_mass);
01370
01371 trk->setZHelix(track.a(),l_mass);
01372 trk->setZError(track.Ea(),l_mass);
01373
01374
01375
01376 int charge=0;
01377 if (kappa > 0.0000000001)
01378 charge = 1;
01379 else if (kappa < -0.0000000001)
01380 charge = -1;
01381 trk->setCharge(charge,l_mass);
01382
01383
01384 double ptot = sqrt(px*px+py*py+pz*pz);
01385 trk->setTheta(acos(pz/ptot),l_mass);
01386 }
01387
01388 else{
01389
01390
01391
01392 double dr = TrasanTRK.helix[0];
01393 double phi0 = TrasanTRK.helix[1];
01394 double kappa = TrasanTRK.helix[2];
01395 double dz = TrasanTRK.helix[3];
01396 double tanl = TrasanTRK.helix[4];
01397
01398 double vx = dr*cos(phi0);
01399 double vy = dr*sin(phi0);
01400 double vz = dz;
01401
01402 if(0==kappa) kappa = 10e-10;
01403 double px = -sin(phi0)/fabs(kappa);
01404 double py = cos(phi0)/fabs(kappa);
01405 double pz = tanl/fabs(kappa);
01406
01407 trk->setX(vx, l_mass);
01408 trk->setY(vy, l_mass);
01409 trk->setZ(vz, l_mass);
01410
01411 trk->setPx(px, l_mass);
01412 trk->setPy(py, l_mass);
01413 trk->setPz(pz, l_mass);
01414
01415 const HepPoint3D poca(dr*cos(phi0),dr*sin(phi0),dz);
01416
01417 trk->setPoca(poca,l_mass);
01418
01419
01420 double a_trasan[5], ea_trasan[15];
01421 for( int i =0 ; i <5; i++){
01422 a_trasan[i] = TrasanTRK.helix[i];
01423 }
01424 for( int j =0 ; j <15; j++){
01425 ea_trasan[j] = TrasanTRK.error[j];
01426 }
01427 trk->setZHelix(a_trasan,l_mass);
01428 trk->setZError(ea_trasan,l_mass);
01429
01430
01431 int charge=0;
01432 if (kappa > 0.0000000001)
01433 charge = 1;
01434 else if (kappa < -0.0000000001)
01435 charge = -1;
01436 trk->setCharge(charge,l_mass);
01437
01438
01439 double ptot = sqrt(px*px+py*py+pz*pz);
01440 trk->setTheta(acos(pz/ptot),l_mass);
01441
01442
01443
01444 SmartDataPtr<RecMdcTrackCol> mdcTrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
01445
01446
01447 RecMdcTrackCol::iterator iter_mdcTrk = mdcTrkCol->begin();
01448 bool findMdcTrk = false;
01449 for(; iter_mdcTrk != mdcTrkCol->end(); iter_mdcTrk++) {
01450 if(TrasanTRK.id==(*iter_mdcTrk)->trackId()) {
01451 findMdcTrk = true;
01452 break;
01453 }
01454 }
01455 int nLayer = (*iter_mdcTrk)->nlayer();
01456 trk->setNlayer(nLayer, l_mass);
01457
01458 }
01459
01460 if(4==debug_) {
01461 RecMdcKalTrack::setPidType(RecMdcKalTrack::electron);
01462 std::cout<<"px: "<<trk->px()<<" py: "<<trk->py()<<" pz: "<<trk->pz()<<std::endl;
01463 std::cout<<"vx: "<<trk->x()<<" vy: "<<trk->y()<<" vz: "<<trk->z()<<std::endl;
01464 }
01465 }
01466
01467
01468 void KalFitAlg::fillTds_back(KalFitTrack& track,
01469 RecMdcKalTrack* trk, MdcRec_trk& TrasanTRK, int l_mass)
01470 {
01471
01472 HepPoint3D IP(0,0,0);
01473
01474
01475
01476 int iqual(1);
01477
01478 if ((trk->getNdf(0,l_mass))-(track.ndf_back()-5)>5) iqual = 0;
01479
01480 if(debug_ == 4) cout<< "fillTds_back> mass "<<trk->getMass(2)<<" ndf[0] "<<trk->getNdf(0,2)<<endl;
01481 if(debug_ == 4) cout<<"ndf_back "<< track.ndf_back() << " chi2_back " << track.chiSq_back()<<endl;
01482
01483 if (track.ndf_back() > 5 && track.chiSq_back() > 0 &&
01484 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01485 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01486 track.Ea()[4][4] > 0 && fabs(track.a()[0]) < DBL_MAX &&
01487 fabs(track.a()[1]) < DBL_MAX && fabs(track.a()[2]) < DBL_MAX &&
01488 fabs(track.a()[3]) < DBL_MAX && fabs(track.a()[4]) < DBL_MAX &&
01489 iqual) {
01490
01491
01492
01493 trk->setStat(0,1,l_mass);
01494 trk->setChisq(track.chiSq_back(),1,l_mass);
01495 trk->setNdf(track.ndf_back()-5,1,l_mass);
01496 trk->setLength(track.pathip(),l_mass);
01497 if(debug_ == 4) cout<<"l_mass "<<l_mass<<" path set as "<<track.pathip()<<endl;
01498 trk->setTof(track.tof(),l_mass);
01499
01500 if (KalFitTrack::tofall_){
01501 if(l_mass == 3) trk->setTof(track.tof_kaon(),l_mass);
01502 if(l_mass == 4) trk->setTof(track.tof_proton(),l_mass);
01503 }
01504
01505
01506 if (pathl_)
01507 for (int i = 0; i<43; i++) {
01508 trk->setPathl(track.pathl()[i],i);
01509 }
01510
01511 trk->setLHelix(track.a(),l_mass);
01512 trk->setLError(track.Ea(),l_mass);
01513 trk->setLPivot(track.pivot(),l_mass);
01514
01515 trk->setLPoint(track.point_last(),l_mass);
01516 trk->setPathSM(track.getPathSM(),l_mass);
01517 trk->setTof(track.getTofSM(),l_mass);
01518 trk->setFiTerm(track.getFiTerm(),l_mass);
01519
01520 if(4 == debug_){
01521 std::cout<<" last pivot: "<< trk->getLPivot(0)<<std::endl;
01522 std::cout<<" pathl in SM: "<< trk->getPathSM(0)<<std::endl;
01523 std::cout<<" fiTerm: "<< trk->getFiTerm(0)<<std::endl;
01524 std::cout<<" last point: "<< trk->getLPoint(0)<<std::endl;
01525 }
01526
01527 } else {
01528 if(debug_) cout<<"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
01529
01530 trk->setStat(1,1,l_mass);
01531 HepPoint3D piv(TrasanTRK.pivot[0],
01532 TrasanTRK.pivot[1],
01533 TrasanTRK.pivot[2]);
01534
01535 HepVector a(5);
01536 for(int i = 0; i < 5; i++)
01537 a[i] = TrasanTRK.helix[i];
01538
01539 HepSymMatrix ea(5);
01540 for(int i = 0, k = 0; i < 5; i++) {
01541 for(int j = 0; j <= i; j++) {
01542 ea[i][j] = matrixg_*TrasanTRK.error[k++];
01543 ea[j][i] = ea[i][j];
01544 }
01545 }
01546
01547 KalFitTrack track_rep(piv, a, ea, lead_,
01548 TrasanTRK.chiSq, TrasanTRK.nhits);
01549 double fiTerm = TrasanTRK.fiTerm;
01550
01551 double fi0 = track_rep.phi0();
01552 HepPoint3D xc(track_rep.kappa()/fabs(track_rep.kappa())*
01553 track_rep.center() );
01554 double x = xc.x();
01555 double y = xc.y();
01556 double phi_x;
01557 if( fabs( x ) > 1.0e-10 ){
01558 phi_x = atan2( y, x );
01559 if( phi_x < 0 ) phi_x += 2*M_PI;
01560 } else {
01561 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
01562 }
01563 if(debug_ == 4) cout<<"fiterm "<<fiTerm<<" fi0 "<<fi0<<" phi_x "<<phi_x<<endl;
01564 double dphi = fabs( fiTerm + fi0 - phi_x );
01565 if( dphi >= 2*M_PI ) dphi -= 2*M_PI;
01566 double tanl = track_rep.tanl();
01567 double cosl_inv = sqrt( tanl*tanl + 1.0 );
01568 if(debug_ == 4) {
01569 cout<<"tanl= "<<tanl<<" radius "<<track_rep.radius()<<" dphi "<<dphi<<endl;
01570 cout<<" cosl_inv "<<cosl_inv<<" radius_numf "<<track_rep.radius_numf()<<endl;
01571 }
01572 double track_len(fabs( track_rep.radius() * dphi * cosl_inv ));
01573 double light_speed( 29.9792458 );
01574 double pt( 1.0 / track_rep.kappa() );
01575 double p( pt * sqrt( 1.0 + tanl*tanl ) );
01576
01577
01578 trk->setStat(1,1,l_mass);
01579 trk->setChisq(TrasanTRK.chiSq,1,l_mass);
01580 if(debug_ == 4) {
01581 std::cout<<".....fillTds_back...chiSq..."<< TrasanTRK.chiSq<<endl;
01582 std::cout<<"...track_len..."<<track_len<<" ndf[1] "<< trk->getNdf(0,l_mass)<<endl;
01583 }
01584 trk->setNdf(TrasanTRK.nhits-5,1,l_mass);
01585 trk->setLength(track_len,l_mass);
01586 double mass_over_p( KalFitTrack::mass(l_mass)/ p );
01587 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
01588 trk->setTof(track_len / ( light_speed * beta ), l_mass) ;
01589
01590 track_rep.pivot(IP);
01591
01592 trk->setLHelix(track_rep.a(),l_mass);
01593 trk->setLError(track_rep.Ea(),l_mass);
01594 trk->setLPivot(track.pivot(),l_mass);
01595
01597 trk->setLPoint(track.point_last(),l_mass);
01598 trk->setPathSM(track.getPathSM(),l_mass);
01599 trk->setTof(track.getTofSM(),l_mass);
01600 trk->setFiTerm(track.getFiTerm(),l_mass);
01601 }
01602
01603
01604 if(debug_ == 4) {
01605
01606 std::cout<<" last point: "<< trk->getLPoint(0)<<std::endl;
01607 std::cout<<" pathl in SM: "<< trk->getPathSM(0)<<std::endl;
01608 std::cout<<" fiTerm: "<< trk->getFiTerm(0)<<std::endl;
01609
01610 cout<<"Now let us see results after smoothering at IP:........."<<endl;
01611 cout << " dr = " << track.a()[0]
01612 << ", Er_dr = " << sqrt(track.Ea()[0][0]) << std::endl;
01613 cout<< " phi0 = " << track.a()[1]
01614 << ", Er_phi0 = " << sqrt(track.Ea()[1][1]) << std::endl;
01615 cout << " PT = " << 1/track.a()[2]
01616 << ", Er_kappa = " << sqrt(track.Ea()[2][2]) << std::endl;
01617 cout << " dz = " << track.a()[3]
01618 << ", Er_dz = " << sqrt(track.Ea()[3][3]) << std::endl;
01619 cout << " tanl = " << track.a()[4]
01620 << ", Er_tanl = " << sqrt(track.Ea()[4][4]) << std::endl;
01621 cout << " Ea = " << track.Ea() <<endl;
01622 }
01623
01624 }
01625
01626 void KalFitAlg::fillTds_back(KalFitTrack& track,
01627 RecMdcKalTrack* trk, MdcRec_trk& TrasanTRK, int l_mass,
01628 RecMdcKalHelixSegCol* segcol)
01629 {
01630
01631 HepPoint3D IP(0,0,0);
01632
01633
01634 track.pivot(IP);
01635
01636
01637
01638
01639
01640 if ((trk->getNdf(0,l_mass))-(track.ndf_back()-5)>5){
01641 iqual_back_ = 0;
01642 }
01643 if(usage_>1){
01644 for(int i=0; i<5; i++) iqual_front_[i] = 1;
01645 iqual_back_ = 1;
01646 }
01647 if(debug_ == 4){
01648 std::cout<< "fillTds_back> mass "<<trk->getMass(2)<<" ndf[0][l_mass] "<<trk->getNdf(0,l_mass)<<endl;
01649 std::cout<<"ndf_back "<< track.ndf_back() << " chi2_back " << track.chiSq_back()<<endl;
01650 std::cout<<"track.ndf_back(), track.chiSq_back(), track.Ea()[5][5], track.a()[5], iqual_front_, iqual_back_: "<<track.ndf_back()<<" , "<<track.chiSq_back()<<" , "<<track.Ea()<<" , "<<track.a()<<" , "<<iqual_front_[l_mass]<<" , "<<iqual_back_<<std::endl;
01651 }
01652
01653 if (track.ndf_back() > 5 && track.chiSq_back() > 0 &&
01654 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01655 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01656 track.Ea()[4][4] > 0 && fabs(track.a()[0]) < DBL_MAX &&
01657 fabs(track.a()[1]) < DBL_MAX && fabs(track.a()[2]) < DBL_MAX &&
01658 fabs(track.a()[3]) < DBL_MAX && fabs(track.a()[4]) < DBL_MAX && iqual_front_[l_mass] && iqual_back_){
01659
01660
01661
01662
01663
01664 HelixSegRefVec helixsegrefvec;
01665 for (vector<KalFitHelixSeg>::iterator it = track.HelixSegs().begin(); it!=track.HelixSegs().end();it ++)
01666 {
01667
01668
01669
01670
01671 it->pivot(IP);
01672
01673
01674
01675 RecMdcKalHelixSeg* helixseg = new RecMdcKalHelixSeg;
01676 helixseg->setResIncl(it->residual_include());
01677 helixseg->setResExcl(it->residual_exclude());
01678 if(debug_ == 4) {
01679 std::cout<<"helixseg->Res_inc ..."<<helixseg->getResIncl()<<std::endl;
01680 }
01681 helixseg->setDrIncl(it->a_include()[0]);
01682 helixseg->setFi0Incl(it->a_include()[1]);
01683 helixseg->setCpaIncl(it->a_include()[2]);
01684 helixseg->setDzIncl(it->a_include()[3]);
01685 helixseg->setTanlIncl(it->a_include()[4]);
01686
01687
01688 helixseg->setDrExcl(it->a_exclude()[0]);
01689 helixseg->setFi0Excl(it->a_exclude()[1]);
01690 helixseg->setCpaExcl(it->a_exclude()[2]);
01691 helixseg->setDzExcl(it->a_exclude()[3]);
01692 helixseg->setTanlExcl(it->a_exclude()[4]);
01693
01694 helixseg->setHelixIncl(it->a_include());
01695 helixseg->setErrorIncl(it->Ea_include());
01696
01697
01698
01699
01700
01701 helixseg->setHelixExcl(it->a_exclude());
01702 helixseg->setErrorExcl(it->Ea_exclude());
01703 helixseg->setLayerId(it->layer());
01704
01705 if(debug_ == 4) {
01706 std::cout<<"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
01707 std::cout<<"helixseg a: "<<it->a()<<std::endl;
01708 std::cout<<"helixseg a_excl: "<<helixseg->getHelixExcl()<<std::endl;
01709 std::cout<<"helixseg a_incl: "<<helixseg->getHelixIncl()<<std::endl;
01710
01711 std::cout<<"helixseg Ea: "<<it->Ea()<<std::endl;
01712 std::cout<<"helixseg Ea_excl: "<<helixseg->getErrorExcl()<<std::endl;
01713 std::cout<<"helixseg Ea_incl: "<<helixseg->getErrorIncl()<<std::endl;
01714
01715 std::cout<<"helixseg layer: "<<it->layer()<<std::endl;
01716 }
01717
01718
01719 helixseg->setTrackId(it->HitMdc()->rechitptr()->getTrkId());
01720 helixseg->setMdcId(it->HitMdc()->rechitptr()->getMdcId());
01721 helixseg->setFlagLR(it->HitMdc()->LR());
01722 helixseg->setTdc(it->HitMdc()->rechitptr()->getTdc());
01723 helixseg->setAdc(it->HitMdc()->rechitptr()->getAdc());
01724 helixseg->setZhit(it->HitMdc()->rechitptr()->getZhit());
01725 helixseg->setTof(it->tof());
01726 helixseg->setDocaIncl(it->doca_include());
01727 helixseg->setDocaExcl(it->doca_exclude());
01728 helixseg->setDD(it->dd());
01729 helixseg->setEntra(it->HitMdc()->rechitptr()->getEntra());
01730 helixseg->setDT(it->dt());
01731 segcol->push_back(helixseg);
01732 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
01733 helixsegrefvec.push_back(refhelixseg);
01734
01735 if(ntuple_&8){
01736 m_docaInc = helixseg -> getDocaIncl();
01737 m_docaExc = helixseg -> getDocaExcl();
01738 m_residualInc = helixseg -> getResIncl();
01739 m_residualExc = helixseg -> getResExcl();
01740 m_dd = helixseg -> getDD();
01741 m_lr = helixseg->getFlagLR();
01742 m_tdrift = helixseg -> getDT();
01743 m_layerid = helixseg -> getLayerId();
01744 m_yposition= it->HitMdc()->wire().fwd().y();
01745 m_eventNo = eventNo;
01746 StatusCode sc6 = m_nt6->write();
01747 if( sc6.isFailure() ) cout<<"Ntuple6 helixseg filling failed!"<<endl;
01748
01749 }
01750 }
01751
01752 trk->setVecHelixSegs(helixsegrefvec, l_mass);
01753 if(debug_ == 4) {
01754 std::cout<<"trk->getVecHelixSegs size..."<<(trk->getVecHelixSegs()).size()<<std::endl;
01755 }
01756 trk->setStat(0,1,l_mass);
01757 trk->setChisq(track.chiSq_back(),1,l_mass);
01758 trk->setNdf(track.ndf_back()-5,1,l_mass);
01759
01760 trk->setNhits(track.ndf_back(),l_mass);
01761 if(!(track.ndf_back()==track.HelixSegs().size())) {
01762 std::cout<<"THEY ARE NOT EQUALL!!!"<<std::endl;
01763 }
01764 trk->setLength(track.pathip(),l_mass);
01765 if(debug_ == 4) {
01766 std::cout<<"l_mass "<<l_mass<<" path set as "<<track.pathip()<<endl;
01767 }
01768 trk->setTof(track.tof(),l_mass);
01769 if (KalFitTrack::tofall_){
01770 if(l_mass == 3) trk->setTof(track.tof_kaon(),l_mass);
01771 if(l_mass == 4) trk->setTof(track.tof_proton(),l_mass);
01772 }
01773
01774 if (pathl_)
01775 for (int i = 0; i<43; i++) {
01776 trk->setPathl(track.pathl()[i],i);
01777 }
01778 trk->setLHelix(track.a(),l_mass);
01779 trk->setLError(track.Ea(),l_mass);
01780 trk->setLPivot(track.pivot(),l_mass);
01781
01782 trk->setLPoint(track.point_last(),l_mass);
01783 trk->setPathSM(track.getPathSM(),l_mass);
01784 trk->setTof(track.getTofSM(),l_mass);
01785 trk->setFiTerm(track.getFiTerm(),l_mass);
01786 double a_trasan[5], ea_trasan[15];
01787 for( int i =0 ; i <5; i++){
01788 a_trasan[i] = TrasanTRK.helix[i];
01789 }
01790 for( int j =0 ; j <15; j++){
01791 ea_trasan[j] = TrasanTRK.helix[j];
01792 }
01793 trk->setTHelix(a_trasan);
01794 trk->setTError(ea_trasan);
01795
01796 if(4 == debug_){
01797 std::cout<<" last pivot: "<< trk->getLPivot(0)<<std::endl;
01798 std::cout<<" pathl in SM: "<< trk->getPathSM(0)<<std::endl;
01799 std::cout<<" fiTerm: "<< trk->getFiTerm(0)<<std::endl;
01800 std::cout<<" last point: "<< trk->getLPoint(0)<<std::endl;
01801 }
01802
01803 } else {
01804
01805 if(debug_) cout<<"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
01806
01807 trk->setStat(1,1,l_mass);
01808
01809 HepPoint3D piv(TrasanTRK.pivot[0],
01810 TrasanTRK.pivot[1],
01811 TrasanTRK.pivot[2]);
01812
01813 HepVector a(5);
01814 for(int i = 0; i < 5; i++)
01815 a[i] = TrasanTRK.helix[i];
01816
01817 HepSymMatrix ea(5);
01818 for(int i = 0, k = 0; i < 5; i++) {
01819 for(int j = 0; j <= i; j++) {
01820 ea[i][j] = matrixg_*TrasanTRK.error[k++];
01821 ea[j][i] = ea[i][j];
01822 }
01823 }
01824
01825 KalFitTrack track_rep(piv, a, ea, lead_,
01826 TrasanTRK.chiSq, TrasanTRK.nhits);
01827 double fiTerm = TrasanTRK.fiTerm;
01828
01829 double fi0 = track_rep.phi0();
01830 HepPoint3D xc(track_rep.kappa()/fabs(track_rep.kappa())*
01831 track_rep.center() );
01832 double x = xc.x();
01833 double y = xc.y();
01834 double phi_x;
01835 if( fabs( x ) > 1.0e-10 ){
01836 phi_x = atan2( y, x );
01837 if( phi_x < 0 ) phi_x += 2*M_PI;
01838 } else {
01839 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
01840 }
01841 if(debug_ == 4) cout<<"fiterm "<<fiTerm<<" fi0 "<<fi0<<" phi_x "<<phi_x<<endl;
01842 double dphi = fabs( fiTerm + fi0 - phi_x );
01843 if( dphi >= 2*M_PI ) dphi -= 2*M_PI;
01844 double tanl = track_rep.tanl();
01845 double cosl_inv = sqrt( tanl*tanl + 1.0 );
01846 if(debug_ == 4) {
01847 cout<<"tanl= "<<tanl<<" radius "<<track_rep.radius()<<" dphi "<<dphi<<endl;
01848 cout<<" cosl_inv "<<cosl_inv<<" radius_numf "<<track_rep.radius_numf()<<endl;
01849 }
01850 double track_len(fabs( track_rep.radius() * dphi * cosl_inv ));
01851 double light_speed( 29.9792458 );
01852 double pt( 1.0 / track_rep.kappa() );
01853 double p( pt * sqrt( 1.0 + tanl*tanl ) );
01854
01855
01856
01857
01858 trk->setStat(1,1,l_mass);
01859 trk->setChisq(TrasanTRK.chiSq,1,l_mass);
01860 if(debug_ == 4) {
01861 std::cout<<".....fillTds_back...chiSq..."<< TrasanTRK.chiSq<<std::endl;
01862 std::cout<<"...track_len..."<<track_len<<" ndf[1] "<< trk->getNdf(0,l_mass)<<std::endl;
01863 }
01864 trk->setNdf(TrasanTRK.nhits-5,1,l_mass);
01865 trk->setLength(track_len,l_mass);
01866 double mass_over_p( KalFitTrack::mass(l_mass)/ p );
01867 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
01868 trk->setTof(track_len / ( light_speed * beta ), l_mass) ;
01869
01870 track_rep.pivot(IP);
01871
01872 trk->setLHelix(track_rep.a(),l_mass);
01873 trk->setLError(track_rep.Ea(),l_mass);
01874 trk->setLPivot(track.pivot(),l_mass);
01875
01877 trk->setLPoint(track.point_last(),l_mass);
01878 trk->setPathSM(track.getPathSM(),l_mass);
01879 trk->setTof(track.getTofSM(),l_mass);
01880 trk->setFiTerm(track.getFiTerm(),l_mass);
01881 trk->setTHelix(track_rep.a());
01882 trk->setTError(track_rep.Ea());
01883
01884 }
01885
01886
01887 if(debug_ == 4) {
01888 cout<<"Now let us see results after smoothering at IP:........."<<endl;
01889 cout << " dr = " << track.a()[0]
01890 << ", Er_dr = " << sqrt(track.Ea()[0][0]) << std::endl;
01891 cout<< " phi0 = " << track.a()[1]
01892 << ", Er_phi0 = " << sqrt(track.Ea()[1][1]) << std::endl;
01893 cout << " PT = " << 1/track.a()[2]
01894 << ", Er_kappa = " << sqrt(track.Ea()[2][2]) << std::endl;
01895 cout << " dz = " << track.a()[3]
01896 << ", Er_dz = " << sqrt(track.Ea()[3][3]) << std::endl;
01897 cout << " tanl = " << track.a()[4]
01898 << ", Er_tanl = " << sqrt(track.Ea()[4][4]) << std::endl;
01899 cout << " Ea = " << track.Ea() <<endl;
01900 }
01901
01902 }
01903
01904 void KalFitAlg::fillTds_back(KalFitTrack& track,
01905 RecMdcKalTrack* trk, MdcRec_trk& TrasanTRK, int l_mass,
01906 RecMdcKalHelixSegCol* segcol, int smoothflag)
01907 {
01908
01909 HepPoint3D IP(0,0,0);
01910
01911
01912
01913
01914
01915
01916
01917
01918 iqual_back_ = 1;
01919 if ((trk->getNdf(0,l_mass))-(track.ndf_back()-5)>5){
01920 iqual_back_ = 0;
01921 }
01922
01923 if(debug_ == 4){
01924 std::cout<< "fillTds_back> mass "<<trk->getMass(2)<<" ndf[0][l_mass] "<<trk->getNdf(0,l_mass)<<endl;
01925 std::cout<<"ndf_back "<< track.ndf_back() << " chi2_back " << track.chiSq_back()<<endl;
01926 }
01927
01928
01929 if (track.ndf_back() > 5 && track.chiSq_back() > 0 &&
01930 track.Ea()[0][0] > 0 && track.Ea()[1][1] > 0 &&
01931 track.Ea()[2][2] > 0 && track.Ea()[3][3] > 0 &&
01932 track.Ea()[4][4] > 0 && fabs(track.a()[0]) < DBL_MAX &&
01933 fabs(track.a()[1]) < DBL_MAX && fabs(track.a()[2]) < DBL_MAX &&
01934 fabs(track.a()[3]) < DBL_MAX && fabs(track.a()[4]) < DBL_MAX && iqual_front_[l_mass] && iqual_back_){
01935
01936
01937
01938
01939
01940 HelixSegRefVec helixsegrefvec;
01941 for (vector<KalFitHelixSeg>::iterator it = track.HelixSegs().begin(); it!=track.HelixSegs().end();it ++)
01942 {
01943
01944
01945
01946
01947 it->pivot(IP);
01948
01949
01950
01951 RecMdcKalHelixSeg* helixseg = new RecMdcKalHelixSeg;
01952 helixseg->setResIncl(it->residual_include());
01953 helixseg->setResExcl(it->residual_exclude());
01954 if(debug_ == 4) {
01955 std::cout<<"helixseg->Res_inc ..."<<helixseg->getResIncl()<<std::endl;
01956 }
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970 helixseg->setHelixIncl(it->a_include());
01971
01972
01973
01974
01975
01976
01977 helixseg->setHelixExcl(it->a_exclude());
01978
01979
01980
01981 if(debug_ == 4) {
01982 std::cout<<"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
01983 std::cout<<"helixseg a: "<<it->a()<<std::endl;
01984 std::cout<<"helixseg a_excl: "<<helixseg->getHelixExcl()<<std::endl;
01985 std::cout<<"helixseg a_incl: "<<helixseg->getHelixIncl()<<std::endl;
01986
01987 std::cout<<"helixseg Ea: "<<it->Ea()<<std::endl;
01988 std::cout<<"helixseg Ea_excl: "<<helixseg->getErrorExcl()<<std::endl;
01989 std::cout<<"helixseg Ea_incl: "<<helixseg->getErrorIncl()<<std::endl;
01990
01991 std::cout<<"helixseg layer: "<<it->layer()<<std::endl;
01992 }
01993
01994
01995 helixseg->setTrackId(it->HitMdc()->rechitptr()->getTrkId());
01996 helixseg->setMdcId(it->HitMdc()->rechitptr()->getMdcId());
01997 helixseg->setFlagLR(it->HitMdc()->LR());
01998 helixseg->setTdc(it->HitMdc()->rechitptr()->getTdc());
01999 helixseg->setAdc(it->HitMdc()->rechitptr()->getAdc());
02000 helixseg->setZhit(it->HitMdc()->rechitptr()->getZhit());
02001 helixseg->setTof(it->tof());
02002 helixseg->setDocaIncl(it->doca_include());
02003 helixseg->setDocaExcl(it->doca_exclude());
02004 helixseg->setDD(it->dd());
02005 helixseg->setEntra(it->HitMdc()->rechitptr()->getEntra());
02006 helixseg->setDT(it->dt());
02007
02008 segcol->push_back(helixseg);
02009 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
02010 helixsegrefvec.push_back(refhelixseg);
02011 if(ntuple_&8){
02012 m_docaInc = helixseg -> getDocaIncl();
02013 m_docaExc = helixseg -> getDocaExcl();
02014 m_residualInc = helixseg -> getResIncl();
02015 m_residualExc = helixseg -> getResExcl();
02016 m_dd = helixseg -> getDD();
02017 m_lr = helixseg->getFlagLR();
02018 m_tdrift = helixseg -> getDT();
02019 m_layerid = helixseg -> getLayerId();
02020 m_yposition= it->HitMdc()->wire().fwd().y();
02021 m_eventNo = eventNo;
02022 StatusCode sc6 = m_nt6->write();
02023 if( sc6.isFailure() ) cout<<"Ntuple6 helixseg filling failed!"<<endl;
02024
02025 }
02026 }
02027
02028 trk->setVecHelixSegs(helixsegrefvec, l_mass);
02029
02030 if(debug_ == 4) {
02031 std::cout<<"trk->getVecHelixSegs size..."<<(trk->getVecHelixSegs()).size()<<std::endl;
02032 }
02033 trk->setStat(0,1,l_mass);
02034 trk->setChisq(track.chiSq_back(),1,l_mass);
02035 trk->setNdf(track.ndf_back()-5,1,l_mass);
02036
02037 trk->setNhits(track.ndf_back(),l_mass);
02038 if(!(track.ndf_back()==track.HelixSegs().size())) {
02039 std::cout<<"THEY ARE NOT EQUALL!!!"<<std::endl;
02040 }
02041 trk->setLength(track.pathip(),l_mass);
02042 if(debug_ == 4) {
02043 std::cout<<"l_mass "<<l_mass<<" path set as "<<track.pathip()<<endl;
02044 }
02045 trk->setTof(track.tof(),l_mass);
02046 if (KalFitTrack::tofall_){
02047 if(l_mass == 3) trk->setTof(track.tof_kaon(),l_mass);
02048 if(l_mass == 4) trk->setTof(track.tof_proton(),l_mass);
02049 }
02050
02051 if (pathl_)
02052 for (int i = 0; i<43; i++) {
02053 trk->setPathl(track.pathl()[i],i);
02054 }
02055 trk->setLHelix(track.a(),l_mass);
02056 trk->setLError(track.Ea(),l_mass);
02057 trk->setLPivot(track.pivot(),l_mass);
02058
02059 trk->setLPoint(track.point_last(),l_mass);
02060 trk->setPathSM(track.getPathSM(),l_mass);
02061 trk->setTof(track.getTofSM(),l_mass);
02062 trk->setFiTerm(track.getFiTerm(),l_mass);
02063 double a_trasan[5], ea_trasan[15];
02064 for( int i =0 ; i <5; i++){
02065 a_trasan[i] = TrasanTRK.helix[i];
02066 }
02067 for( int j =0 ; j <15; j++){
02068 ea_trasan[j] = TrasanTRK.helix[j];
02069 }
02070 trk->setTHelix(a_trasan);
02071 trk->setTError(ea_trasan);
02072
02073 if(4 == debug_){
02074 std::cout<<" last pivot: "<< trk->getLPivot(0)<<std::endl;
02075 std::cout<<" pathl in SM: "<< trk->getPathSM(0)<<std::endl;
02076 std::cout<<" fiTerm: "<< trk->getFiTerm(0)<<std::endl;
02077 std::cout<<" last point: "<< trk->getLPoint(0)<<std::endl;
02078 }
02079
02080 } else {
02081
02082 if(debug_) cout<<"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
02083
02084 trk->setStat(1,1,l_mass);
02085
02086 HepPoint3D piv(TrasanTRK.pivot[0],
02087 TrasanTRK.pivot[1],
02088 TrasanTRK.pivot[2]);
02089
02090 HepVector a(5);
02091 for(int i = 0; i < 5; i++)
02092 a[i] = TrasanTRK.helix[i];
02093
02094 HepSymMatrix ea(5);
02095 for(int i = 0, k = 0; i < 5; i++) {
02096 for(int j = 0; j <= i; j++) {
02097 ea[i][j] = matrixg_*TrasanTRK.error[k++];
02098 ea[j][i] = ea[i][j];
02099 }
02100 }
02101
02102 KalFitTrack track_rep(piv, a, ea, lead_,
02103 TrasanTRK.chiSq, TrasanTRK.nhits);
02104 double fiTerm = TrasanTRK.fiTerm;
02105
02106 double fi0 = track_rep.phi0();
02107 HepPoint3D xc(track_rep.kappa()/fabs(track_rep.kappa())*
02108 track_rep.center() );
02109 double x = xc.x();
02110 double y = xc.y();
02111 double phi_x;
02112 if( fabs( x ) > 1.0e-10 ){
02113 phi_x = atan2( y, x );
02114 if( phi_x < 0 ) phi_x += 2*M_PI;
02115 } else {
02116 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
02117 }
02118 if(debug_ == 4) cout<<"fiterm "<<fiTerm<<" fi0 "<<fi0<<" phi_x "<<phi_x<<endl;
02119 double dphi = fabs( fiTerm + fi0 - phi_x );
02120 if( dphi >= 2*M_PI ) dphi -= 2*M_PI;
02121 double tanl = track_rep.tanl();
02122 double cosl_inv = sqrt( tanl*tanl + 1.0 );
02123 if(debug_ == 4) {
02124 cout<<"tanl= "<<tanl<<" radius "<<track_rep.radius()<<" dphi "<<dphi<<endl;
02125 cout<<" cosl_inv "<<cosl_inv<<" radius_numf "<<track_rep.radius_numf()<<endl;
02126 }
02127
02128 double track_len(fabs( track_rep.radius() * fiTerm * cosl_inv ));
02129
02130 double light_speed( 29.9792458 );
02131 double pt( 1.0 / track_rep.kappa() );
02132 double p( pt * sqrt( 1.0 + tanl*tanl ) );
02133
02134
02135
02136
02137 trk->setStat(1,1,l_mass);
02138 trk->setChisq(TrasanTRK.chiSq,1,l_mass);
02139 if(debug_ == 4) {
02140 std::cout<<".....fillTds_back...chiSq..."<< TrasanTRK.chiSq<<std::endl;
02141 std::cout<<"...track_len..."<<track_len<<" ndf[1] "<< trk->getNdf(0,l_mass)<<std::endl;
02142 }
02143 trk->setNdf(TrasanTRK.nhits-5,1,l_mass);
02144 trk->setLength(track_len,l_mass);
02145 double mass_over_p( KalFitTrack::mass(l_mass)/ p );
02146 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02147 trk->setTof(track_len / ( light_speed * beta ), l_mass) ;
02148
02149
02150 HepPoint3D LPiovt = track_rep.x(fiTerm);
02151 track_rep.pivot(LPiovt);
02152
02153 trk->setLHelix(track_rep.a(),l_mass);
02154 trk->setLError(track_rep.Ea(),l_mass);
02155
02156
02157 trk->setLPivot(LPiovt, l_mass);
02158
02160 trk->setLPoint(track.point_last(),l_mass);
02161
02162 trk->setPathSM(track_len,l_mass);
02163
02164
02165 trk->setFiTerm(fiTerm,l_mass);
02166 trk->setTHelix(track_rep.a());
02167 trk->setTError(track_rep.Ea());
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234 }
02235
02236
02237 if(debug_ == 4) {
02238 cout<<"Now let us see results after smoothering at IP:........."<<endl;
02239 cout << " dr = " << track.a()[0]
02240 << ", Er_dr = " << sqrt(track.Ea()[0][0]) << std::endl;
02241 cout<< " phi0 = " << track.a()[1]
02242 << ", Er_phi0 = " << sqrt(track.Ea()[1][1]) << std::endl;
02243 cout << " PT = " << 1/track.a()[2]
02244 << ", Er_kappa = " << sqrt(track.Ea()[2][2]) << std::endl;
02245 cout << " dz = " << track.a()[3]
02246 << ", Er_dz = " << sqrt(track.Ea()[3][3]) << std::endl;
02247 cout << " tanl = " << track.a()[4]
02248 << ", Er_tanl = " << sqrt(track.Ea()[4][4]) << std::endl;
02249 cout << " Ea = " << track.Ea() <<endl;
02250 }
02251
02252 }
02253
02254
02255
02256 void KalFitAlg::sameas(RecMdcKalTrack* trk, int l_mass, int imain)
02257 {
02258
02259
02260 trk->setMass(trk->getMass(imain), l_mass);
02261 trk->setLength(trk->getLength(imain), l_mass);
02262 trk->setTof(trk->getTof(imain), l_mass);
02263 trk->setNhits(trk->getNhits(imain), l_mass);
02264
02265 for(int jj = 0; jj<2; jj++) {
02266 trk->setStat(trk->getStat(jj,imain), jj, l_mass);
02267 trk->setChisq(trk->getChisq(jj, l_mass), jj, l_mass);
02268 trk->setNdf(trk->getChisq(jj, l_mass), jj, l_mass);
02269 }
02270 trk->setLHelix(trk->getFHelix(),l_mass);
02271 trk->setLError(trk->getFError(),l_mass);
02272 }
02273
02274
02275 void KalFitAlg::smoother_anal(KalFitTrack& track, int way)
02276 {
02277
02278 IMdcGeomSvc* igeomsvc;
02279 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", igeomsvc);
02280 if(sc==StatusCode::FAILURE) cout << "GeoSvc failing!!!!!!!SC=" << sc << endl;
02281 MdcGeomSvc* geomsvc = dynamic_cast<MdcGeomSvc*>(igeomsvc);
02282 if(!geomsvc){
02283 std::cout<<"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
02284 }
02285
02286 HepPoint3D ip(0,0,0);
02287 if(m_usevtxdb==1){
02288 Hep3Vector xorigin(0,0,0);
02289 IVertexDbSvc* vtxsvc;
02290 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
02291 if(vtxsvc->isVertexValid()){
02292 double* dbv = vtxsvc->PrimaryVertex();
02293 double* vv = vtxsvc->SigmaPrimaryVertex();
02294 xorigin.setX(dbv[0]);
02295 xorigin.setY(dbv[1]);
02296 xorigin.setZ(dbv[2]);
02297 }
02298 ip[0] = xorigin[0];
02299 ip[1] = xorigin[1];
02300 ip[2] = xorigin[2];
02301 }
02302
02303
02304
02305 Helix work = *(Helix*)&track;
02306 work.ignoreErrorMatrix();
02307 work.pivot(ip);
02308
02309
02310 double tanl = track.tanl();
02311 double phi_old = work.phi0();
02312 double phi = track.phi0();
02313
02314 if (fabs(phi - phi_old) > M_PI) {
02315 if (phi > phi_old) phi -= 2 * M_PI;
02316 else phi_old -= 2 * M_PI;
02317 }
02318
02319 double path_zero = fabs(track.radius() * (phi_old-phi)* sqrt(1 + tanl * tanl));
02320
02321
02322
02323 HepSymMatrix Eakal(5,0);
02324 track.pivot(ip);
02327 Eakal = track.Ea()*matrixg_;
02328 track.Ea(Eakal);
02329
02330
02331 unsigned int nhit = track.HitsMdc().size();
02332 int layer_prev = -1;
02333
02334 HepVector pos_old(3,0);
02335 double r0kal_prec(0);
02336 int nhits_read(0);
02337 for( unsigned i=0 ; i < nhit; i++ ) {
02338 int ihit = (nhit-1)-i;
02339 KalFitHitMdc& HitMdc = track.HitMdc(ihit);
02340 const KalFitWire& Wire = HitMdc.wire();
02341
02342 int wireid = Wire.geoID();
02343 nhits_read++;
02344
02345 int layer = Wire.layer().layerId();
02346 if (pathl_ && layer != layer_prev) {
02347
02348 if (debug_ == 4) cout<<"in smoother,layerid "<<layer<<" layer_prev "
02349 <<layer_prev <<" pathl_ "<<pathl_<<endl;
02350
02351
02352 layer_prev = layer;
02353 }
02354
02355 HepPoint3D fwd(Wire.fwd());
02356 HepPoint3D bck(Wire.bck());
02357 Hep3Vector wire = (CLHEP::Hep3Vector)fwd - (CLHEP::Hep3Vector)bck;
02358 Helix work = *(Helix*)&track;
02359 work.ignoreErrorMatrix();
02360 work.pivot((fwd + bck) * .5);
02361 HepPoint3D x0kal = (work.x(0).z() - bck.z()) / wire.z() * wire + bck;
02362
02363 if(4 == debug_) std::cout<<" x0kal before sag: "<<x0kal<<std::endl;
02364 if (wsag_ == 4){
02365 Hep3Vector result;
02366 const MdcGeoWire* geowire = geomsvc->Wire(wireid);
02367 double tension = geowire->Tension();
02368
02369
02370 double zinit(x0kal.z()), lzx(Wire.lzx());
02371
02372
02373 double A = 47.35E-6/tension;
02374 double Zp = (zinit - bck.z())*lzx/wire.z();
02375
02376 if(4 == debug_){
02377 std::cout<<" sag in smoother_anal: "<<std::endl;
02378 std::cout<<" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
02379 std::cout<<" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
02380 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
02381 std::cout<<"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
02382 std::cout<<" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
02383 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
02384 }
02385
02386 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
02387 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
02388 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
02389
02390 wire.setX(wire.x()/wire.z());
02391 wire.setY(result.z());
02392 wire.setZ(1);
02393
02394 x0kal.setX(result.x());
02395 x0kal.setY(result.y());
02396 }
02397
02398 if(4 == debug_) std::cout<<" x0kal after sag: "<<x0kal<<std::endl;
02399
02400
02401 double r0kal = x0kal.perp();
02402 if (debug_ == 4) {
02403 cout<<"wire direction "<<wire<<endl;
02404 cout<<"x0kal "<<x0kal<<endl;
02405 cout<<"smoother::r0kal "<<r0kal<<" r0kal_prec "<<r0kal_prec <<endl;
02406 }
02407
02408
02409
02410
02411
02412
02413 double pathl(0);
02414 track.pivot_numf(x0kal, pathl);
02415 track.addPathSM(pathl);
02416
02417
02418
02419
02420
02421
02422 double pmag( sqrt( 1.0 + track.a()[4]*track.a()[4]) / track.a()[2]);
02423 double mass_over_p( track.mass()/ pmag );
02424 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
02425 double tofest = pathl / ( 29.9792458 * beta );
02426 track.addTofSM(tofest);
02427
02428
02429
02430 if(KalFitElement::muls()) track.msgasmdc(pathl, way);
02431
02432
02433
02434
02435 if(!(way<0&&fabs(track.kappa())>1000.0)) {
02436 if(KalFitElement::loss()) track.eloss(pathl, _BesKalmanFitMaterials[0], way);
02437 }
02438
02439
02440
02441
02442
02443
02444
02445 if(fabs(track.kappa())>0&&fabs(track.kappa())<1000.0&&fabs(track.tanl())<7.02) {
02446 HepVector Va(5,0);
02447 HepSymMatrix Ma(5,0);
02448 KalFitHelixSeg HelixSeg(&HitMdc,x0kal,Va,Ma);
02449
02450 Hep3Vector meas = track.momentum(0).cross(wire).unit();
02451 double dchi2=-1;
02452 track.smoother_Mdc(HitMdc, meas, HelixSeg, dchi2, m_csmflag);
02453 if(dchi2>0.0) {
02454 track.HelixSegs().push_back(HelixSeg);
02455 }
02456 }
02457
02459
02460 if(i == nhit-1){
02461
02463 HepPoint3D point;
02464 point.setX(x0kal.x() + track.a()[0]*cos(track.a()[1]));
02465 point.setY(x0kal.y() + track.a()[0]*sin(track.a()[1]));
02466 point.setZ(x0kal.z() + track.a()[3]);
02467 track.point_last(point);
02468
02470 double phi_old = track.a()[1];
02471 KalFitTrack temp(x0kal, track.a(), track.Ea(), 0, 0, 0);
02472 temp.pivot(ip);
02473 double phi_new = temp.a()[1];
02474 double fi = phi_new - phi_old;
02476
02477
02478 if(fabs(fi) >= CLHEP::twopi) fi = fmod(fi+2*CLHEP::twopi,CLHEP::twopi);
02479
02480 track.fiTerm(fi);
02481 }
02482
02483 if (debug_) cout<<"track----7-----"<<track.a()<<endl;
02484 r0kal_prec = r0kal;
02485 }
02486 }
02487
02488
02489
02490 void KalFitAlg::smoother_calib(KalFitTrack& track, int way)
02491 {
02492
02493
02494 IMdcGeomSvc* igeomsvc;
02495 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", igeomsvc);
02496 if(sc==StatusCode::FAILURE) cout << "GeoSvc failing!!!!!!!SC=" << sc << endl;
02497 MdcGeomSvc* geomsvc = dynamic_cast<MdcGeomSvc*>(igeomsvc);
02498 if(!geomsvc){
02499 std::cout<<"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
02500 }
02501
02502 HepSymMatrix Eakal(5,0);
02503
02504 HepPoint3D ip(0,0,0);
02505 track.pivot(ip);
02506 Eakal = track.getInitMatrix();
02507 if( debug_ == 4) {
02508 std::cout<<"the initial error matrix in smoothing is "<<Eakal<<std::endl;
02509 }
02510 track.Ea(Eakal);
02511
02512
02513 unsigned int nseg = track.HelixSegs().size();
02514 int layer_prev = -1;
02515
02516 HepVector pos_old(3,0);
02517 double r0kal_prec(0);
02518 int nsegs_read(0);
02519 for( unsigned i=0 ; i < nseg; i++ ) {
02520
02521 int flag=0;
02522 int iseg = (nseg-1)-i;
02523 KalFitHelixSeg& HelixSeg = track.HelixSeg(iseg);
02524 const KalFitWire& Wire = HelixSeg.HitMdc()->wire();
02525
02526 int wireid = Wire.geoID();
02527 nsegs_read++;
02528
02529 int layer = Wire.layer().layerId();
02530 if (pathl_ && layer != layer_prev) {
02531
02532 if (debug_ == 4) cout<<"in smoother,layerid "<<layer<<" layer_prev "
02533 <<layer_prev <<" pathl_ "<<pathl_<<endl;
02534
02535
02536 layer_prev = layer;
02537 }
02538
02539 HepPoint3D fwd(Wire.fwd());
02540 HepPoint3D bck(Wire.bck());
02541 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
02542 Helix work = *(Helix*)&track;
02543 work.ignoreErrorMatrix();
02544 work.pivot((fwd + bck) * .5);
02545 HepPoint3D x0kal = (work.x(0).z() - bck.z()) / wire.z() * wire + bck;
02546
02547
02548 if(4 == debug_) std::cout<<" x0kal before sag: "<<x0kal<<std::endl;
02549
02550 if (wsag_ == 4){
02551
02552 Hep3Vector result;
02553 const MdcGeoWire* geowire = geomsvc->Wire(wireid);
02554 double tension = geowire->Tension();
02555
02556
02557 double zinit(x0kal.z()), lzx(Wire.lzx());
02558
02559
02560
02561 double A = 47.35E-6/tension;
02562 double Zp = (zinit - bck.z())*lzx/wire.z();
02563
02564 if(4 == debug_){
02565
02566 std::cout<<" sag in smoother_calib: "<<std::endl;
02567 std::cout<<" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
02568 std::cout<<" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
02569 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
02570 std::cout<<"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
02571 std::cout<<" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
02572 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
02573 }
02574
02575 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
02576 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
02577 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
02578
02579 wire.setX(wire.x()/wire.z());
02580 wire.setY(result.z());
02581 wire.setZ(1);
02582
02583 x0kal.setX(result.x());
02584 x0kal.setY(result.y());
02585 }
02586
02587 if(4 == debug_) std::cout<<" x0kal after sag: "<<x0kal<<std::endl;
02588
02589
02590
02591 double r0kal = x0kal.perp();
02592 if (debug_ == 4) {
02593 cout<<"wire direction "<<wire<<endl;
02594 cout<<"x0kal "<<x0kal<<endl;
02595 cout<<"smoother::r0kal "<<r0kal<<" r0kal_prec "<<r0kal_prec <<endl;
02596 }
02597
02598
02599 double pathl(0);
02600 track.pivot_numf(x0kal, pathl);
02601
02602 if (debug_ == 4) cout<<"track----6-----"<<track.a()<<" ...path..."<<pathl
02603 <<"momentum"<<track.momentum(0)<<endl;
02604 if(KalFitElement::muls()) track.msgasmdc(pathl, way);
02605 if(KalFitElement::loss()) track.eloss(pathl, _BesKalmanFitMaterials[0], way);
02606
02607
02608 if(fabs(track.kappa())>0&&fabs(track.kappa())<1000.0&&fabs(track.tanl())<7.02) {
02609
02610 Hep3Vector meas = track.momentum(0).cross(wire).unit();
02611
02612 if(usage_>1) track.smoother_Mdc_csmalign(HelixSeg, meas,flag, m_csmflag);
02613 else track.smoother_Mdc(HelixSeg, meas,flag, m_csmflag);
02614
02615 }
02616
02617 if (debug_) cout<<"track----7-----"<<track.a()<<endl;
02618 r0kal_prec = r0kal;
02619
02620 if(flag == 0) {
02621 track.HelixSegs().erase(track.HelixSegs().begin()+iseg);
02622 }
02623 }
02624 }
02625
02626
02627
02628 void KalFitAlg::filter_fwd_anal(KalFitTrack& track, int l_mass, int way, HepSymMatrix& Eakal)
02629 {
02630
02631
02632
02633
02634
02635 IMdcGeomSvc* igeomsvc;
02636 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", igeomsvc);
02637 if(sc==StatusCode::FAILURE) cout << "GeoSvc failing!!!!!!!SC=" << sc << endl;
02638 MdcGeomSvc* geomsvc = dynamic_cast<MdcGeomSvc*>(igeomsvc);
02639 if(!geomsvc){
02640 std::cout<<"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
02641 }
02642
02643 Hep3Vector x0inner(track.pivot());
02644 HepVector pos_old(3,0);
02645 double r0kal_prec(0);
02646 int nhits_read(0);
02647 int nhit = track.HitsMdc().size();
02648 if(debug_ == 4) cout<<"filter_fwd..........111 nhit="<<nhit<<endl;
02649 for( int i=0 ; i < nhit; i++ ) {
02650 KalFitHitMdc& HitMdc = track.HitMdc(i);
02651
02652 if (HitMdc.chi2()<0) continue;
02653 const KalFitWire& Wire = HitMdc.wire();
02654 int layerf = Wire.layer().layerId();
02655
02656
02657
02658 int wireid = Wire.geoID();
02659 nhits_read++;
02660 HepPoint3D fwd(Wire.fwd());
02661 HepPoint3D bck(Wire.bck());
02662 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
02663 Helix work = *(Helix*)&track;
02664 work.ignoreErrorMatrix();
02665 work.pivot((fwd + bck) * .5);
02666
02667
02668
02669
02670
02671
02672 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
02673
02674 if(4 == debug_) std::cout<<" x0kal before sag: "<<x0kal<<std::endl;
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715 if (wsag_ == 4){
02716 Hep3Vector result;
02717 const MdcGeoWire* geowire = geomsvc->Wire(wireid);
02718 double tension = geowire->Tension();
02719
02720 double zinit(x0kal.z()), lzx(Wire.lzx());
02721
02722 double A = 47.35E-6/tension;
02723 double Zp = (zinit - bck.z())*lzx/wire.z();
02724
02725 if(4 == debug_){
02726 std::cout<<" sag in filter_fwd_anal: "<<std::endl;
02727 std::cout<<" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
02728 std::cout<<"zinit: "<<zinit<<" bck.z(): "<<bck.z()<<std::endl;
02729 std::cout<<" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
02730 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
02731 std::cout<<"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
02732 std::cout<<" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
02733 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
02734 }
02735
02736 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
02737 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
02738 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
02739
02740 wire.setX(wire.x()/wire.z());
02741 wire.setY(result.z());
02742 wire.setZ(1);
02743 x0kal.setX(result.x());
02744 x0kal.setY(result.y());
02745 }
02746
02747 if(4 == debug_) std::cout<<" x0kal after sag: "<<x0kal<<std::endl;
02748
02749
02750 double r0kal = x0kal.perp();
02751
02752
02753 double pathl(0);
02754
02755 track.pivot_numf(x0kal, pathl);
02756
02757 if (nhits_read == 1) {
02758 track.Ea(Eakal);
02759 } else {
02760 if(KalFitElement::muls()) track.msgasmdc(pathl, way);
02761 if(KalFitElement::loss()) track.eloss(pathl, _BesKalmanFitMaterials[0], way);
02762 }
02763
02764 double dtracknew = 0.;
02765 double dtrack = 0.;
02766 double dtdc = 0.;
02767
02768 if(fabs(track.kappa())>0&&fabs(track.kappa())<1000.0&&fabs(track.tanl())<7.02) {
02769 Hep3Vector meas = track.momentum(0).cross(wire).unit();
02770 double diff_chi2 = track.chiSq();
02771 Hep3Vector IP(0,0,0);
02772 Helix work_bef = *(Helix*)&track;
02773 work_bef.ignoreErrorMatrix();
02774 work_bef.pivot(IP);
02775 int inext(-1);
02776 if (i+1<nhit)
02777 for( unsigned k=i+1 ; k < nhit; k++ )
02778 if (!(track.HitMdc(k).chi2()<0)) {
02779 inext = (signed) k;
02780 break;
02781 }
02782 double dchi2 = -1.0;
02783
02784 double chi2 = track.update_hits(HitMdc,inext,meas,way,dchi2,dtrack,dtracknew,dtdc,m_csmflag);
02785
02786
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807
02808
02809
02810
02811
02812
02813
02814
02815
02816
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832 if( dchi2 <0 ) {
02833 std::cout<<" ... ERROR OF dchi2... "<<std::endl;
02834 }
02835
02836 if (ntuple_&8) {
02837 m_dchi2 = dchi2;
02838 m_masshyp = l_mass;
02839 m_residest = dtrack-dtdc;
02840 m_residnew = dtracknew -dtdc;
02841 m_layer = Wire.layer().layerId();
02842 Helix worktemp = *(Helix*)&track;
02843 m_anal_dr = worktemp.a()[0];
02844 m_anal_phi0 = worktemp.a()[1];
02845 m_anal_kappa = worktemp.a()[2];
02846 m_anal_dz = worktemp.a()[3];
02847 m_anal_tanl = worktemp.a()[4];
02848 m_anal_ea_dr = worktemp.Ea()[0][0];
02849 m_anal_ea_phi0 = worktemp.Ea()[1][1];
02850 m_anal_ea_kappa = worktemp.Ea()[2][2];
02851 m_anal_ea_dz = worktemp.Ea()[3][3];
02852 m_anal_ea_tanl = worktemp.Ea()[4][4];
02853 StatusCode sc5 = m_nt5->write();
02854 if( sc5.isFailure() ) cout<<"Ntuple5 filling failed!"<<endl;
02855 }
02856 Helix work_aft = *(Helix*)&track;
02857 work_aft.ignoreErrorMatrix();
02858 work_aft.pivot(IP);
02859 diff_chi2 = chi2 - diff_chi2;
02860 HitMdc.chi2(diff_chi2);
02861 }
02862 r0kal_prec = r0kal;
02863 }
02864 }
02865
02866
02867
02868 void KalFitAlg::filter_fwd_calib(KalFitTrack& track, int l_mass, int way, HepSymMatrix& Eakal)
02869 {
02870
02871
02872 IMdcGeomSvc* igeomsvc;
02873 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", igeomsvc);
02874 if(sc==StatusCode::FAILURE) cout << "GeoSvc failing!!!!!!!SC=" << sc << endl;
02875 MdcGeomSvc* geomsvc = dynamic_cast<MdcGeomSvc*>(igeomsvc);
02876 if(!geomsvc){
02877 std::cout<<"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
02878 }
02879
02880 if(debug_ == 4) {
02881 std::cout<<"filter_fwd_calib starting ...the inital error Matirx is "<<track.Ea()<<std::endl;
02882 }
02883 Hep3Vector x0inner(track.pivot());
02884 HepVector pos_old(3,0);
02885 double r0kal_prec(0);
02886 int nhits_read(0);
02887
02888 unsigned int nhit = track.HitsMdc().size();
02889 if(debug_ == 4) cout<<"filter_fwd..........111 nhit="<<nhit<<endl;
02890 for( unsigned i=0 ; i < nhit; i++ ) {
02891
02892 KalFitHitMdc& HitMdc = track.HitMdc(i);
02893 if(debug_ == 4)
02894 cout<<"filter_fwd...........222 chi2="<<HitMdc.chi2()<<endl;
02895
02896 if (HitMdc.chi2()<0) continue;
02897 const KalFitWire& Wire = HitMdc.wire();
02898
02899 int wireid = Wire.geoID();
02900 nhits_read++;
02901
02902
02903 HepPoint3D fwd(Wire.fwd());
02904 HepPoint3D bck(Wire.bck());
02905 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
02906 Helix work = *(Helix*)&track;
02907 work.ignoreErrorMatrix();
02908 work.pivot((fwd + bck) * .5);
02909 HepPoint3D x0kal = (work.x(0).z() - bck.z())/ wire.z() * wire + bck;
02910
02911 if(4 == debug_)
02912 std::cout<<" x0kal before sag: "<<x0kal<<std::endl;
02913
02914 if (wsag_ == 4){
02915 Hep3Vector result;
02916 const MdcGeoWire* geowire = geomsvc->Wire(wireid);
02917 double tension = geowire->Tension();
02918
02919 double zinit(x0kal.z()), lzx(Wire.lzx());
02920
02921 double A = 47.35E-6/tension;
02922 double Zp = (zinit - bck.z())*lzx/wire.z();
02923
02924 if(4 == debug_){
02925
02926 std::cout<<" sag in filter_fwd_calib: "<<std::endl;
02927 std::cout<<" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
02928 std::cout<<" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
02929 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
02930 std::cout<<"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
02931 std::cout<<" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
02932 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
02933 }
02934
02935 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
02936 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
02937 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
02938 wire.setX(wire.x()/wire.z());
02939 wire.setY(result.z());
02940 wire.setZ(1);
02941 x0kal.setX(result.x());
02942 x0kal.setY(result.y());
02943 }
02944
02945 if(4 == debug_)
02946 std::cout<<" x0kal after sag: "<<x0kal<<std::endl;
02947
02948
02949 double r0kal = x0kal.perp();
02950
02951
02952 double pathl(0);
02953
02954
02955
02956
02957 if (debug_ == 4)
02958 cout << "*** move from " << track.pivot() << " ( Kappa = "
02959 << track.a()[2] << ")" << endl;
02960 track.pivot_numf(x0kal, pathl);
02961
02962
02963
02964
02965
02966 if (debug_ == 4)
02967 cout << "*** to " << track.pivot() << " ( Kappa = "
02968 << track.a()[2] << ")" << std::endl;
02969
02970 if (nhits_read == 1) {
02971 track.Ea(Eakal);
02972 if(debug_==4) cout << "filter_fwd::Ea = " << track.Ea()<<endl;
02973
02974
02975 } else {
02976 if(KalFitElement::muls()) track.msgasmdc(pathl, way);
02977 if(KalFitElement::loss()) track.eloss(pathl, _BesKalmanFitMaterials[0], way);
02978 }
02979
02980
02981
02982
02983
02984
02985 if(fabs(track.kappa())>0&&fabs(track.kappa())<1000.0&&fabs(track.tanl())<7.02) {
02986 Hep3Vector meas = track.momentum(0).cross(wire).unit();
02987
02988 double diff_chi2 = track.chiSq();
02989
02990 Hep3Vector IP(0,0,0);
02991 Helix work_bef = *(Helix*)&track;
02992 work_bef.ignoreErrorMatrix();
02993 work_bef.pivot(IP);
02994 int inext(-1);
02995 if (i+1<nhit)
02996 for( unsigned k=i+1 ; k < nhit; k++ )
02997 if (!(track.HitMdc(k).chi2()<0)) {
02998 inext = (signed) k;
02999 break;
03000 }
03001 double dchi2 = -1.0;
03002
03003 if (debug_ == 4) {
03004 std::cout<<"the value of x0kal is "<<x0kal<<std::endl;
03005 std::cout<<"the value of track.pivot() is "<<track.pivot()<<std::endl;
03006 }
03007
03008 HepVector Va(5,0);
03009 HepSymMatrix Ma(5,0);
03010 KalFitHelixSeg HelixSeg(&HitMdc,x0kal,Va,Ma);
03011
03012 if(debug_ == 4) {
03013 std::cout<<"HelixSeg.Ea_pre_fwd() ..."<<HelixSeg.Ea_pre_fwd()<<std::endl;
03014 std::cout<<"HelixSeg.a_pre_fwd() ..."<<HelixSeg.a_pre_fwd()<<std::endl;
03015 std::cout<<"HelixSeg.Ea_filt_fwd() ..."<<HelixSeg.Ea_filt_fwd()<<std::endl;
03016 }
03017
03018
03019
03020
03021 double chi2=0.;
03022 if(usage_>1) chi2=track.update_hits_csmalign(HelixSeg, inext, meas, way, dchi2, m_csmflag);
03023 else chi2 = track.update_hits(HelixSeg, inext, meas, way, dchi2, m_csmflag);
03024
03025
03026
03027
03028 if(debug_ ==4 ) cout<<"layer, cell, dchi2,chi2, a, p: "<<HitMdc.wire().layer().layerId()<<" , "<<HitMdc.wire().localId()<<" , "<<dchi2<<" , "<<chi2<<" , "<<track.a()<<" , "<<sqrt(1.0+track.a()[4]*track.a()[4])/track.a()[2]<<endl;
03029
03030 if(debug_ == 4) cout<< "****inext***"<<inext<<" !!!! dchi2= "<< dchi2
03031 << " chisq= "<< chi2<< endl;
03032
03033 if (ntuple_&8) {
03034 m_dchi2 = dchi2;
03035 m_masshyp = l_mass;
03036 StatusCode sc5 = m_nt5->write();
03037 if( sc5.isFailure() ) cout<<"Ntuple5 filling failed!"<<endl;
03038 }
03039
03040 Helix work_aft = *(Helix*)&track;
03041 work_aft.ignoreErrorMatrix();
03042 work_aft.pivot(IP);
03043 diff_chi2 = chi2 - diff_chi2;
03044 HitMdc.chi2(diff_chi2);
03045
03046 if(debug_ == 4) {
03047
03048 cout << " -> after include meas = " << meas
03049 << " at R = " << track.pivot().perp() << std::endl;
03050 cout << " chi2 = " << chi2 << ", diff_chi2 = "
03051 << diff_chi2 << ", LR = "
03052 << HitMdc.LR() << ", stereo = " << HitMdc.wire().stereo()
03053 << ", layer = " << HitMdc.wire().layer().layerId() << std::endl;
03054 cout<<"filter_fwd..........HitMdc.chi2... "<<HitMdc.chi2()<<endl;
03055 }
03056
03057 if(dchi2>0.0 && (way!=999)) {
03058 track.HelixSegs().push_back(HelixSeg);
03059 }
03060 }
03061 r0kal_prec = r0kal;
03062 }
03063 }
03064
03065
03066 void KalFitAlg::innerwall(KalFitTrack& track, int l_mass, int way){
03067
03068 if (debug_ ==4) cout<<"innerwall....."<<endl;
03069 for(int i = 0; i < _BesKalmanFitWalls.size(); i++) {
03070 _BesKalmanFitWalls[i].updateTrack(track, way);
03071 if (debug_ == 4) cout<<"Wall "<<i<<" update the track!(filter)"<<endl;
03072 }
03073 }
03074
03075
03076
03077
03078
03079
03080 void KalFitAlg::kalman_fitting_anal(void)
03081 {
03082
03083
03084
03085 MsgStream log(msgSvc(), name());
03086 double Pt_threshold(0.3);
03087 Hep3Vector IP(0,0,0);
03088 vector<MdcRec_trk>* mdcMgr = MdcRecTrkCol::getMdcRecTrkCol();
03089 vector<MdcRec_trk_add>* mdc_addMgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
03090 vector<MdcRec_wirhit>* whMgr = MdcRecWirhitCol::getMdcRecWirhitCol();
03091
03092
03093 if ( !&whMgr ) return;
03094
03095
03096 int ntrk = mdcMgr->size();
03097 double* rPt = new double[ntrk];
03098 int* rOM = new int[ntrk];
03099 unsigned int* order = new unsigned int[ntrk];
03100 unsigned int* rCont = new unsigned int[ntrk];
03101 unsigned int* rGen = new unsigned int[ntrk];
03102
03103 int index = 0;
03104 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
03105 end = mdcMgr->end(); it != end; it++) {
03106
03107 rPt[index] = 0;
03108 if (it->helix[2])
03109 rPt[index] = 1 / fabs(it->helix[2]);
03110 if(debug_ == 4) cout<<"rPt...[ "<<index<<" ]...."<< rPt[index] <<endl;
03111 if(rPt[index] < 0) rPt[index] = DBL_MAX;
03112
03113 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
03114 if(debug_ == 4) cout<<"ppt size: "<< pt.size()<<endl;
03115 int outermost(-1);
03116 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03117 ii !=pt.begin()-1; ii--) {
03118 int lyr((*ii)->geo->Lyr()->Id());
03119 if (outermost < lyr) outermost = lyr;
03120 if(debug_ == 4) cout<<"outmost: "<<outermost<<" lyr: "<<lyr<<endl;
03121 }
03122 rOM[index] = outermost;
03123 order[index] = index;
03124 ++index;
03125 }
03126
03127
03128 for (int j, k = ntrk - 1; k >= 0; k = j){
03129 j = -1;
03130 for(int i = 1; i <= k; i++)
03131 if(rPt[i - 1] < rPt[i]){
03132 j = i - 1;
03133 std::swap(order[i], order[j]);
03134 std::swap(rPt[i], rPt[j]);
03135 std::swap(rOM[i], rOM[j]);
03136 std::swap(rCont[i], rCont[j]);
03137 std::swap(rGen[i], rGen[j]);
03138 }
03139 }
03140 delete [] rPt;
03141
03142
03143 int newcount(0);
03144
03145 DataObject *aReconEvent;
03146 eventSvc()->findObject("/Event/Recon",aReconEvent);
03147 if(!aReconEvent) {
03148
03149 ReconEvent* recevt = new ReconEvent;
03150 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
03151 if(sc!=StatusCode::SUCCESS) {
03152 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
03153 return;
03154 }
03155 }
03156
03157 RecMdcKalTrackCol* kalcol = new RecMdcKalTrackCol;
03158 RecMdcKalHelixSegCol *segcol =new RecMdcKalHelixSegCol;
03159
03160 log << MSG::INFO << "beginning to make RecMdcKalTrackCol" <<endreq;
03161
03162 for(int l = 0; l < ntrk; l++) {
03163
03164
03165
03166
03167 nTotalTrks++;
03168
03169 for(int i=0; i<5; i++) iqual_front_[i] = 1;
03170
03171
03172 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
03173 MdcRec_trk_add& TrasanTRK_add = *(mdc_addMgr->begin()+order[l]);
03174
03175
03176 int trasqual = TrasanTRK_add.quality;
03177 if(debug_ == 4) cout<<"kalman_fitting>trasqual... "<<trasqual<<endl;
03178 if (trasqual) continue;
03179
03180 newcount++;
03181 if (debug_ == 4)
03182 cout << "******* KalFit NUMBER : " << newcount << std::endl;
03183
03184
03185 int type(0);
03186 if ((TrasanTRK_add.decision & 32) == 32 ||
03187 (TrasanTRK_add.decision & 64) == 64) type = 1;
03188
03189
03190 HepPoint3D x(TrasanTRK.pivot[0],
03191 TrasanTRK.pivot[1],
03192 TrasanTRK.pivot[2]);
03193
03194 HepVector a(5);
03195 for(int i = 0; i < 5; i++)
03196 a[i] = TrasanTRK.helix[i];
03197
03198 HepSymMatrix ea(5);
03199 for(int i = 0, k = 0; i < 5; i++) {
03200 for(int j = 0; j <= i; j++) {
03201 ea[i][j] = matrixg_*TrasanTRK.error[k++];
03202 ea[j][i] = ea[i][j];
03203 }
03204 }
03205 double fiTerm = TrasanTRK.fiTerm;
03206 int way(1);
03207
03208 KalFitTrack track_lead = KalFitTrack(x, a, ea, lead_, 0, 0);
03209
03210 track_lead.bFieldZ(KalFitTrack::Bznom_);
03211
03212
03213 int inlyr(999), outlyr(-1);
03214 int* rStat = new int[43];
03215 for(int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
03216 std::vector<MdcRec_wirhit*> pt=TrasanTRK.hitcol;
03217 int hit_in(0);
03218 if(debug_ == 4) cout<<"*********Pt size****"<< pt.size()<<endl;
03219
03220 int Num[43] = {0};
03221 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03222 ii != pt.begin()-1; ii--) {
03223 Num[(*ii)->geo->Lyr()->Id()]++;
03224 }
03225 int hit_asso(0);
03226 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03227 ii != pt.begin()-1; ii--) {
03228 hit_asso++;
03229 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
03230 if (debug_ >0)
03231 cout << "WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
03232 << " hits in the layer "
03233 << (*ii)->geo->Lyr()->Id() << std::endl;
03234 continue;
03235 }
03236
03237
03238
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249 hit_in++;
03250 MdcRec_wirhit & rechit = **ii;
03251 double dist[2] = {rechit.ddl, rechit.ddr};
03252 double erdist[2] = {rechit.erddl, rechit.erddr};
03253 const MdcGeoWire* geo = rechit.geo;
03254
03255 int lr_decision(0);
03256 if (KalFitTrack::LR_ == 1){
03257 if (rechit.lr==2 || rechit.lr==0) lr_decision=-1;
03258
03259 else if (rechit.lr==1) lr_decision=1;
03260 }
03261
03262 int ind(geo->Lyr()->Id());
03263
03264
03265 track_lead.appendHitsMdc( KalFitHitMdc(rechit.id,
03266 lr_decision, rechit.tdc,
03267 dist, erdist,
03268 _wire+(geo->Id()), rechit.rechitptr));
03269
03270 rStat[ind]++;
03271 if (inlyr>ind) inlyr = ind;
03272 if (outlyr<ind) outlyr = ind;
03273 }
03274 if (debug_ == 4)
03275 cout << "**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
03276
03277
03278 int empty_between(0), empty(0);
03279 for (int i= inlyr; i <= outlyr; i++)
03280 if (!rStat[i]) empty_between++;
03281 empty = empty_between+inlyr+(42-outlyr);
03282 delete [] rStat;
03283
03284
03285 track_lead.order_wirhit(1);
03286
03287
03288 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.HitsMdc().begin(); it_hit!=track_lead.HitsMdc().end(); it_hit++){
03289
03290
03291
03292 }
03293
03294 track_lead.type(type);
03295 unsigned int nhit = track_lead.HitsMdc().size();
03296 if (!nhit && debug_ == 4) {
03297 cout << " ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
03298 continue;
03299 }
03300
03301
03302 double KalFitst(0), KalFitax(0), KalFitschi2(0);
03303
03304
03305 Hep3Vector outer_pivot(track_lead.x(fiTerm));
03306
03307 track_lead.pivot(outer_pivot);
03308
03309 track_lead.bFieldZ(KalFitTrack::Bznom_);
03310
03311 if (nhit>=3 && !KalFitTrack::LR_)
03312 start_seed(track_lead, lead_, way, TrasanTRK);
03313 HepSymMatrix Eakal(5,0);
03314
03316 double costheta = track_lead.a()[4] / sqrt(1.0 + track_lead.a()[4]*track_lead.a()[4]);
03317 if( (1.0/fabs(track_lead.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
03318 choice_ = 6;
03319 }
03320
03322 init_matrix(choice_,TrasanTRK, Eakal);
03323
03324
03325 if (debug_ == 4){
03326 cout << "from Mdc Pattern Recognition: " << std::endl;
03327 HepPoint3D IP(0,0,0);
03328 Helix work(track_lead.pivot(),
03329 track_lead.a(),
03330 track_lead.Ea());
03331 work.pivot(IP);
03332 cout << " dr = " << work.a()[0]
03333 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
03334 cout << " phi0 = " << work.a()[1]
03335 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
03336 cout << " PT = " << 1/work.a()[2]
03337 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
03338 cout << " dz = " << work.a()[3]
03339 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
03340 cout << " tanl = " << work.a()[4]
03341 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
03342 }
03343
03344 filter_fwd_anal(track_lead, lead_, way, Eakal);
03345
03346
03347 track_lead.update_forMdc();
03348
03349 HepPoint3D IP(0,0,0);
03350 if (debug_ == 4) {
03351 cout << " Mdc FIRST KALMAN FIT " << std::endl;
03352 Helix work(track_lead.pivot(),
03353 track_lead.a(),
03354 track_lead.Ea());
03355 work.pivot(IP);
03356 cout << " dr = " << work.a()[0]
03357 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
03358 cout << " phi0 = " << work.a()[1]
03359 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
03360 cout << " PT = " << 1/work.a()[2]
03361 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
03362 cout << " dz = " << work.a()[3]
03363 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
03364 cout << " tanl = " << work.a()[4]
03365 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
03366 }
03367
03368
03369 RecMdcKalTrack* kaltrk = new RecMdcKalTrack;
03370
03371 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
03372 }
03373
03374
03375 IDataProviderSvc* eventSvc = NULL;
03376 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
03377 if (eventSvc) {
03378 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
03379 } else {
03380 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
03381 return ;
03382 }
03383
03384 StatusCode kalsc;
03385 IDataManagerSvc *dataManSvc;
03386 dataManSvc= dynamic_cast<IDataManagerSvc*>(eventSvc);
03387 DataObject *aKalTrackCol;
03388 eventSvc->findObject("/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
03389 if(aKalTrackCol != NULL) {
03390 dataManSvc->clearSubTree("/Event/Recon/RecMdcKalTrackCol");
03391 eventSvc->unregisterObject("/Event/Recon/RecMdcKalTrackCol");
03392 }
03393
03394 kalsc = eventSvc->registerObject("/Event/Recon/RecMdcKalTrackCol", kalcol);
03395 if( kalsc.isFailure() ) {
03396 log << MSG::FATAL << "Could not register RecMdcKalTrack" << endreq;
03397 return ;
03398 }
03399 log << MSG::INFO << "RecMdcKalTrackCol registered successfully!" <<endreq;
03400
03401 StatusCode segsc;
03402
03403 DataObject *aRecKalSegEvent;
03404 eventSvc->findObject("/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
03405 if(aRecKalSegEvent!=NULL) {
03406
03407 segsc = eventSvc->unregisterObject("/Event/Recon/RecMdcKalHelixSegCol");
03408 if(segsc != StatusCode::SUCCESS) {
03409 log << MSG::FATAL << "Could not unregister RecMdcKalHelixSegCol collection" << endreq;
03410 return;
03411 }
03412 }
03413
03414 segsc = eventSvc->registerObject("/Event/Recon/RecMdcKalHelixSegCol", segcol);
03415 if( segsc.isFailure() ) {
03416 log << MSG::FATAL << "Could not register RecMdcKalHelixSeg" << endreq;
03417 return;
03418 }
03419 log << MSG::INFO << "RecMdcKalHelixSegCol registered successfully!" <<endreq;
03420
03421
03422 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),phi1(0.),phi2(0.),p1(0.),p2(0.);
03423 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
03424 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
03425
03426
03427 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,"/Event/Recon/RecMdcKalTrackCol");
03428 if (!kaltrkCol) {
03429 log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endreq;
03430 return;
03431 }
03432 log << MSG::INFO << "Begin to check RecMdcKalTrackCol"<<endreq;
03433 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
03434 for( int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
03435 log << MSG::DEBUG << "retrieved MDC Kalmantrack:"
03436 << "Track Id: " << (*iter_trk)->getTrackId()
03437 << " Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
03438 << " Length of the track: "<< (*iter_trk)->getLength(2)
03439 << " Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
03440 << " Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
03441 <<" "<< (*iter_trk)->getChisq(1,2) << endreq
03442 << "Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
03443 <<" "<< (*iter_trk)->getNdf(1,2) << endreq
03444 << "Kappa " << (*iter_trk)->getZHelix()[2]
03445 << endreq;
03446 for( int i = 0; i<43; i++) {
03447 log << MSG::DEBUG << "retrieved pathl["<<i<<"]= "
03448 << (*iter_trk)->getPathl(i) <<endreq;
03449 }
03450
03451 if(ntuple_&1) {
03452 m_trackid = (*iter_trk)->getTrackId();
03453
03454 for( int jj =0, iii=0; jj<5; jj++) {
03455
03456 m_length[jj] = (*iter_trk)->getLength(jj);
03457 m_tof[jj] = (*iter_trk)->getTof(jj);
03458 m_nhits[jj] = (*iter_trk)->getNhits(jj);
03459 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
03460 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
03461 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
03462 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
03463 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
03464 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
03465 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
03466 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
03467 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
03468 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
03469 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
03470 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
03471 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
03472 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
03473 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
03474
03475 if(ntuple_&32) {
03476 for(int kk=0; kk<=jj; kk++,iii++) {
03477 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
03478 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
03479 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
03480 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
03481 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
03482 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
03483 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
03484 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
03485 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
03486 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
03487 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
03488 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
03489 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
03490 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
03491 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
03492 }
03493 }
03494 }
03495
03496
03497
03498
03499
03500
03501
03502
03503
03504
03505
03506
03507
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
03532 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
03533 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
03534 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
03535 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
03536 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
03537 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
03538 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
03539 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
03540 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
03541
03542 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
03543 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
03544 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
03545 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
03546 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
03547 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
03548 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
03549 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
03550 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
03551 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
03552
03553 m_stat[0][0] = (*iter_trk)->getStat(0,0);
03554 m_stat[1][0] = (*iter_trk)->getStat(0,1);
03555 m_stat[2][0] = (*iter_trk)->getStat(0,2);
03556 m_stat[3][0] = (*iter_trk)->getStat(0,3);
03557 m_stat[4][0] = (*iter_trk)->getStat(0,4);
03558 m_stat[0][1] = (*iter_trk)->getStat(1,0);
03559 m_stat[1][1] = (*iter_trk)->getStat(1,1);
03560 m_stat[2][1] = (*iter_trk)->getStat(1,2);
03561 m_stat[3][1] = (*iter_trk)->getStat(1,3);
03562 m_stat[4][1] = (*iter_trk)->getStat(1,4);
03563
03564 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
03565 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
03566 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
03567 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
03568 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
03569
03570 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
03571 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
03572 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
03573 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
03574 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
03575
03576 m_lpt = 1/m_lhelix[2];
03577 m_lpte = 1/m_lhelixe[2];
03578 m_lptmu = 1/m_lhelixmu[2];
03579 m_lptk = 1/m_lhelixk[2];
03580 m_lptp = 1/m_lhelixp[2];
03581
03582 m_fpt = 1/m_fhelix[2];
03583 m_fpte = 1/m_fhelixe[2];
03584 m_fptmu = 1/m_fhelixmu[2];
03585 m_fptk = 1/m_fhelixk[2];
03586 m_fptp = 1/m_fhelixp[2];
03587
03588 if(debug_ >= 3){
03589 std::cout<<" "<<std::endl;
03590 std::cout<<"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
03591 std::cout<<"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
03592 std::cout<<"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
03593 std::cout<<"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
03594 std::cout<<"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
03595 }
03596
03597 m_zpt = 1/m_zhelix[2];
03598 m_zpte = 1/m_zhelixe[2];
03599 m_zptmu = 1/m_zhelixmu[2];
03600 m_zptk = 1/m_zhelixk[2];
03601 m_zptp = 1/m_zhelixp[2];
03602
03603 if(debug_ >= 3) {
03604 std::cout<<"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
03605 std::cout<<"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
03606 std::cout<<"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
03607 std::cout<<"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
03608 std::cout<<"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
03609 }
03610 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
03611 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
03612 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
03613 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
03614 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
03615
03616 if(debug_ >= 3) {
03617 std::cout<<"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
03618 std::cout<<"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
03619 std::cout<<"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
03620 std::cout<<"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
03621 std::cout<<"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
03622 }
03623
03624 if(ntuple_&32) {
03625 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
03626 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
03627 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
03628 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
03629 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
03630 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
03631 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
03632 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
03633 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
03634 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
03635 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
03636 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
03637 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
03638 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
03639 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
03640 }
03641
03642 StatusCode sc1 = m_nt1->write();
03643 if( sc1.isFailure() ) cout<<"Ntuple1 filling failed!"<<endl;
03644 }
03645
03646 if(ntuple_&4) {
03647 if(jj == 1) {
03648 phi1 = (*iter_trk)->getFFi0();
03649 r1 = (*iter_trk)->getFDr();
03650 z1 = (*iter_trk)->getFDz();
03651 kap1 = (*iter_trk)->getFCpa();
03652 tanl1 = (*iter_trk)->getFTanl();
03653 charge1 = kap1/fabs(kap1);
03654 x1 = r1*cos(phi1);
03655 y1 = r1*sin(phi1);
03656 p1 = sqrt(1+tanl1*tanl1)/kap1;
03657 the1 = M_PI/2-atan(tanl1);
03658 px1 = -sin(phi1)/fabs(kap1);
03659 py1 = cos(phi1)/fabs(kap1);
03660 pz1= tanl1/fabs(kap1);
03661
03662 } else if(jj == 2) {
03663 phi2 = (*iter_trk)->getFFi0();
03664 r2 = (*iter_trk)->getFDr();
03665 z2 = (*iter_trk)->getFDz();
03666 kap2 = (*iter_trk)->getFCpa();
03667 tanl2 = (*iter_trk)->getFTanl();
03668 charge2 = kap2/fabs(kap2);
03669 x2 = r1*cos(phi2);
03670 y2 = r1*sin(phi2);
03671 p2 = sqrt(1+tanl2*tanl2)/kap1;
03672 the2 = M_PI/2-atan(tanl2);
03673 px2 = -sin(phi2)/fabs(kap2);
03674 py2 = cos(phi2)/fabs(kap2);
03675 pz2= tanl2/fabs(kap2);
03676 }
03677 }
03678 }
03679 if(ntuple_&4) {
03680 m_delx = x1 - x2;
03681 m_dely = y1 - y2;
03682 m_delz = z1 - z2;
03683 m_delthe = the1 + the2;
03684 m_delphi = phi1- phi2;
03685 m_delp = p1 - p2;
03686 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
03687 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
03688 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
03689
03690 StatusCode sc2 = m_nt2->write();
03691 if( sc2.isFailure() ) cout<<"Ntuple2 filling failed!"<<endl;
03692 }
03693
03694 delete [] order;
03695 delete [] rCont;
03696 delete [] rGen;
03697 delete [] rOM;
03698
03699 if (debug_ == 4)
03700 cout << "Kalfitting finished " << std::endl;
03701 }
03702
03703 void KalFitAlg::kalman_fitting_calib(void)
03704 {
03705
03706 MsgStream log(msgSvc(), name());
03707 double Pt_threshold(0.3);
03708 Hep3Vector IP(0,0,0);
03709
03710 vector<MdcRec_trk>* mdcMgr = MdcRecTrkCol::getMdcRecTrkCol();
03711 vector<MdcRec_trk_add>* mdc_addMgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
03712 vector<MdcRec_wirhit>* whMgr = MdcRecWirhitCol::getMdcRecWirhitCol();
03713
03714
03715 if ( !&whMgr ) return;
03716
03717
03718 int ntrk = mdcMgr->size();
03719 double* rPt = new double[ntrk];
03720 int* rOM = new int[ntrk];
03721 unsigned int* order = new unsigned int[ntrk];
03722 unsigned int* rCont = new unsigned int[ntrk];
03723 unsigned int* rGen = new unsigned int[ntrk];
03724
03725 int index = 0;
03726 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
03727 end = mdcMgr->end(); it != end; it++) {
03728
03729 rPt[index] = 0;
03730 if (it->helix[2])
03731 rPt[index] = 1 / fabs(it->helix[2]);
03732 if(debug_ == 4) cout<<"rPt...[ "<<index<<" ]...."<< rPt[index] <<endl;
03733 if(rPt[index] < 0) rPt[index] = DBL_MAX;
03734
03735 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
03736 if(debug_ == 4) cout<<"ppt size: "<< pt.size()<<endl;
03737 int outermost(-1);
03738 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03739 ii !=pt.begin()-1; ii--) {
03740 int lyr((*ii)->geo->Lyr()->Id());
03741 if (outermost < lyr) outermost = lyr;
03742 if(debug_ == 4) cout<<"outmost: "<<outermost<<" lyr: "<<lyr<<endl;
03743 }
03744 rOM[index] = outermost;
03745 order[index] = index;
03746 ++index;
03747 }
03748
03749
03750 for (int j, k = ntrk - 1; k >= 0; k = j){
03751 j = -1;
03752 for(int i = 1; i <= k; i++)
03753 if(rPt[i - 1] < rPt[i]){
03754 j = i - 1;
03755 std::swap(order[i], order[j]);
03756 std::swap(rPt[i], rPt[j]);
03757 std::swap(rOM[i], rOM[j]);
03758 std::swap(rCont[i], rCont[j]);
03759 std::swap(rGen[i], rGen[j]);
03760 }
03761 }
03762 delete [] rPt;
03763
03764 int newcount(0);
03765
03766 DataObject *aReconEvent;
03767 eventSvc()->findObject("/Event/Recon",aReconEvent);
03768 if(!aReconEvent) {
03769
03770 ReconEvent* recevt = new ReconEvent;
03771 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
03772 if(sc!=StatusCode::SUCCESS) {
03773 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
03774 return;
03775 }
03776 }
03777
03778 RecMdcKalTrackCol* kalcol = new RecMdcKalTrackCol;
03779 RecMdcKalHelixSegCol *segcol =new RecMdcKalHelixSegCol;
03780
03781 log << MSG::INFO << "beginning to make RecMdcKalTrackCol" <<endreq;
03782
03783
03784 for(int l = 0; l < ntrk; l++) {
03785
03786 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
03787 MdcRec_trk_add& TrasanTRK_add = *(mdc_addMgr->begin()+order[l]);
03788
03789
03790 int trasqual = TrasanTRK_add.quality;
03791 if(debug_ == 4) cout<<"kalman_fitting>trasqual... "<<trasqual<<endl;
03792 if (trasqual) continue;
03793
03794 newcount++;
03795 if (debug_ == 4)
03796 cout << "******* KalFit NUMBER : " << newcount << std::endl;
03797
03798
03799 int type(0);
03800 if ((TrasanTRK_add.decision & 32) == 32 ||
03801 (TrasanTRK_add.decision & 64) == 64) type = 1;
03802
03803
03804 HepPoint3D x(TrasanTRK.pivot[0],
03805 TrasanTRK.pivot[1],
03806 TrasanTRK.pivot[2]);
03807
03808 HepVector a(5);
03809 for(int i = 0; i < 5; i++)
03810 a[i] = TrasanTRK.helix[i];
03811
03812 HepSymMatrix ea(5);
03813 for(int i = 0, k = 0; i < 5; i++) {
03814 for(int j = 0; j <= i; j++) {
03815 ea[i][j] = matrixg_*TrasanTRK.error[k++];
03816 ea[j][i] = ea[i][j];
03817 }
03818 }
03819
03820 KalFitTrack::setInitMatrix(ea);
03821
03822 double fiTerm = TrasanTRK.fiTerm;
03823 int way(1);
03824
03825 KalFitTrack track_lead = KalFitTrack(x, a, ea, lead_, 0, 0);
03826 track_lead.bFieldZ(KalFitTrack::Bznom_);
03827
03828 int inlyr(999), outlyr(-1);
03829 int* rStat = new int[43];
03830 for(int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
03831 std::vector<MdcRec_wirhit*> pt=TrasanTRK.hitcol;
03832 int hit_in(0);
03833 if(debug_ == 4) cout<<"*********Pt size****"<< pt.size()<<endl;
03834
03835 int Num[43] = {0};
03836 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03837 ii != pt.begin()-1; ii--) {
03838 Num[(*ii)->geo->Lyr()->Id()]++;
03839 }
03840
03841 int hit_asso(0);
03842 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
03843 ii != pt.begin()-1; ii--) {
03844
03845 hit_asso++;
03846 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
03847 if (debug_ >0)
03848 cout << "WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
03849 << " hits in the layer "
03850 << (*ii)->geo->Lyr()->Id() << std::endl;
03851 continue;
03852 }
03853
03854
03855
03856
03857
03858
03859
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870
03871
03872
03873
03874
03875 hit_in++;
03876 MdcRec_wirhit & rechit = **ii;
03877 double dist[2] = {rechit.ddl, rechit.ddr};
03878 double erdist[2] = {rechit.erddl, rechit.erddr};
03879 const MdcGeoWire* geo = rechit.geo;
03880
03881 int lr_decision(0);
03882 if (KalFitTrack::LR_ == 1){
03883 if (rechit.lr==2 || rechit.lr==0) lr_decision=-1;
03884
03885 else if (rechit.lr==1) lr_decision=1;
03886 }
03887
03888 int ind(geo->Lyr()->Id());
03889 track_lead.appendHitsMdc( KalFitHitMdc(rechit.id,
03890 lr_decision, rechit.tdc,
03891 dist, erdist,
03892 _wire+(geo->Id()), rechit.rechitptr));
03893
03894 rStat[ind]++;
03895 if (inlyr>ind) inlyr = ind;
03896 if (outlyr<ind) outlyr = ind;
03897 }
03898 if (debug_ == 4)
03899 cout << "**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
03900
03901
03902 int empty_between(0), empty(0);
03903 for (int i= inlyr; i <= outlyr; i++)
03904 if (!rStat[i]) empty_between++;
03905 empty = empty_between+inlyr+(42-outlyr);
03906 delete [] rStat;
03907
03908
03909 track_lead.order_wirhit(1);
03910 track_lead.type(type);
03911 unsigned int nhit = track_lead.HitsMdc().size();
03912 if (!nhit && debug_ == 4) {
03913 cout << " ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
03914 continue;
03915 }
03916
03917
03918 double KalFitst(0), KalFitax(0), KalFitschi2(0);
03919
03920 Hep3Vector outer_pivot(track_lead.x(fiTerm));
03921
03922 if(debug_ == 4) {
03923 std::cout<<"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.Ea()<<std::endl;
03924 }
03925 track_lead.pivot(outer_pivot);
03926 track_lead.bFieldZ(KalFitTrack::Bznom_);
03927
03928 if (nhit>=3 && !KalFitTrack::LR_)
03929 start_seed(track_lead, lead_, way, TrasanTRK);
03930 HepSymMatrix Eakal(5,0);
03931
03932
03933
03934 double costheta = track_lead.a()[4] / sqrt(1.0 + track_lead.a()[4]*track_lead.a()[4]);
03935 if( (1.0/fabs(track_lead.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
03936 choice_ = 6;
03937 }
03938
03939 init_matrix(choice_,TrasanTRK, Eakal);
03940
03941
03942
03943 if (debug_ == 4){
03944 std::cout << "from Mdc Pattern Recognition: " << std::endl;
03945 HepPoint3D IP(0,0,0);
03946 Helix work(track_lead.pivot(),
03947 track_lead.a(),
03948 track_lead.Ea());
03949 work.pivot(IP);
03950 std::cout << " dr = " << work.a()[0]
03951 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
03952 std::cout << " phi0 = " << work.a()[1]
03953 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
03954 std::cout << " PT = " << 1/work.a()[2]
03955 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
03956 std::cout << " dz = " << work.a()[3]
03957 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
03958 std::cout << " tanl = " << work.a()[4]
03959 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
03960 }
03961
03962 filter_fwd_calib(track_lead, lead_, way, Eakal);
03963 track_lead.update_forMdc();
03964
03965 HepPoint3D IP(0,0,0);
03966 if (debug_ == 4) {
03967 cout << " Mdc FIRST KALMAN FIT " << std::endl;
03968 Helix work(track_lead.pivot(),
03969 track_lead.a(),
03970 track_lead.Ea());
03971 work.pivot(IP);
03972 cout << " dr = " << work.a()[0]
03973 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
03974 cout << " phi0 = " << work.a()[1]
03975 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
03976 cout << " PT = " << 1/work.a()[2]
03977 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
03978 cout << " dz = " << work.a()[3]
03979 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
03980 cout << " tanl = " << work.a()[4]
03981 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
03982 }
03983
03984
03985 RecMdcKalTrack* kaltrk = new RecMdcKalTrack;
03986
03987
03988 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
03989 }
03990
03991
03992 StatusCode kalsc;
03993
03994 DataObject *aRecKalEvent;
03995 eventSvc()->findObject("/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
03996 if(aRecKalEvent!=NULL) {
03997
03998 kalsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalTrackCol");
03999 if(kalsc != StatusCode::SUCCESS) {
04000 log << MSG::FATAL << "Could not unregister RecMdcKalTrack collection" << endreq;
04001 return;
04002 }
04003 }
04004
04005 kalsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalTrackCol", kalcol);
04006 if( kalsc.isFailure()) {
04007 log << MSG::FATAL << "Could not register RecMdcKalTrack" << endreq;
04008 return;
04009 }
04010 log << MSG::INFO << "RecMdcKalTrackCol registered successfully!" <<endreq;
04011
04012
04013
04014 StatusCode segsc;
04015
04016 DataObject *aRecKalSegEvent;
04017 eventSvc()->findObject("/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
04018 if(aRecKalSegEvent!=NULL) {
04019
04020 segsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalHelixSegCol");
04021 if(segsc != StatusCode::SUCCESS) {
04022 log << MSG::FATAL << "Could not unregister RecMdcKalHelixSegCol collection" << endreq;
04023 return;
04024 }
04025 }
04026
04027 segsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalHelixSegCol", segcol);
04028 if( segsc.isFailure() ) {
04029 log << MSG::FATAL << "Could not register RecMdcKalHelixSeg" << endreq;
04030 return;
04031 }
04032 log << MSG::INFO << "RecMdcKalHelixSegCol registered successfully!" <<endreq;
04033
04034
04035 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),phi1(0.),phi2(0.),p1(0.),p2(0.);
04036 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
04037
04038
04039 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
04040 if (!kaltrkCol) {
04041 log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endreq;
04042 return;
04043 }
04044 log << MSG::INFO << "Begin to check RecMdcKalTrackCol"<<endreq;
04045 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
04046 for( int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
04047 log << MSG::DEBUG << "retrieved MDC Kalmantrack:"
04048 << "Track Id: " << (*iter_trk)->getTrackId()
04049 << " Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
04050 << " Length of the track: "<< (*iter_trk)->getLength(2)
04051 << " Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
04052 << " Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
04053 <<" "<< (*iter_trk)->getChisq(1,2) << endreq
04054 << "Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
04055 <<" "<< (*iter_trk)->getNdf(1,2) << endreq
04056 << "Kappa " << (*iter_trk)->getZHelix()[2]
04057 << endreq;
04058
04059 HelixSegRefVec gothelixsegs = (*iter_trk)->getVecHelixSegs();
04060 if(debug_ == 4) {
04061 std::cout<<"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
04062 }
04063
04064 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
04065 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
04066 if(debug_ == 4) {
04067 std::cout<<"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
04068 std::cout<<"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
04069 std::cout<<"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
04070 std::cout<<"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
04071 std::cout<<"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
04072 std::cout<<"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
04073 }
04074 }
04075 for( int i = 0; i<43; i++) {
04076 log << MSG::DEBUG << "retrieved pathl["<<i<<"]= "
04077 << (*iter_trk)->getPathl(i) <<endreq;
04078 }
04079
04080 if(ntuple_&1) {
04081 m_trackid = (*iter_trk)->getTrackId();
04082 for( int jj =0, iii=0; jj<5; jj++) {
04083 m_length[jj] = (*iter_trk)->getLength(jj);
04084 m_tof[jj] = (*iter_trk)->getTof(jj);
04085 m_nhits[jj] = (*iter_trk)->getNhits(jj);
04086 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
04087 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
04088 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
04089 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
04090 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
04091 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
04092 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
04093 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
04094 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
04095 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
04096 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
04097 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
04098 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
04099 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
04100 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
04101 if(ntuple_&32) {
04102 for(int kk=0; kk<=jj; kk++,iii++) {
04103 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
04104 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
04105 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
04106 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
04107 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
04108 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
04109 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
04110 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
04111 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
04112 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
04113 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
04114 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
04115 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
04116 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
04117 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
04118 }
04119
04120 }
04121 }
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146
04147
04148
04149
04150
04151
04152
04153
04154
04155
04156
04157
04158 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
04159 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
04160 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
04161 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
04162 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
04163 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
04164 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
04165 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
04166 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
04167 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
04168
04169 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
04170 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
04171 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
04172 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
04173 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
04174 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
04175 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
04176 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
04177 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
04178 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
04179
04180 m_stat[0][0] = (*iter_trk)->getStat(0,0);
04181 m_stat[1][0] = (*iter_trk)->getStat(0,1);
04182 m_stat[2][0] = (*iter_trk)->getStat(0,2);
04183 m_stat[3][0] = (*iter_trk)->getStat(0,3);
04184 m_stat[4][0] = (*iter_trk)->getStat(0,4);
04185 m_stat[0][1] = (*iter_trk)->getStat(1,0);
04186 m_stat[1][1] = (*iter_trk)->getStat(1,1);
04187 m_stat[2][1] = (*iter_trk)->getStat(1,2);
04188 m_stat[3][1] = (*iter_trk)->getStat(1,3);
04189 m_stat[4][1] = (*iter_trk)->getStat(1,4);
04190
04191 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
04192 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
04193 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
04194 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
04195 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
04196
04197 m_zpt = 1/m_zhelix[2];
04198 m_zpte = 1/m_zhelixe[2];
04199 m_zptmu = 1/m_zhelixmu[2];
04200 m_zptk = 1/m_zhelixk[2];
04201 m_zptp = 1/m_zhelixp[2];
04202
04203 m_fpt = 1/m_fhelix[2];
04204 m_fpte = 1/m_fhelixe[2];
04205 m_fptmu = 1/m_fhelixmu[2];
04206 m_fptk = 1/m_fhelixk[2];
04207 m_fptp = 1/m_fhelixp[2];
04208
04209 m_lpt = 1/m_lhelix[2];
04210 m_lpte = 1/m_lhelixe[2];
04211 m_lptmu = 1/m_lhelixmu[2];
04212 m_lptk = 1/m_lhelixk[2];
04213 m_lptp = 1/m_lhelixp[2];
04214
04215 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
04216 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
04217 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
04218 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
04219 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
04220
04221 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
04222 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
04223 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
04224 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
04225 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
04226 if(ntuple_&32) {
04227 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
04228 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
04229 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
04230 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
04231 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
04232 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
04233 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
04234 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
04235 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
04236 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
04237 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
04238 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
04239 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
04240 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
04241 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
04242 }
04243
04244 StatusCode sc1 = m_nt1->write();
04245 if( sc1.isFailure() ) cout<<"Ntuple1 filling failed!"<<endl;
04246 }
04247
04248 if(ntuple_&4) {
04249 if(jj == 1) {
04250 phi1 = (*iter_trk)->getFFi0();
04251 r1 = (*iter_trk)->getFDr();
04252 z1 = (*iter_trk)->getFDz();
04253 kap1 = (*iter_trk)->getFCpa();
04254 tanl1 = (*iter_trk)->getFTanl();
04255 x1 = r1*cos(phi1);
04256 y1 = r1*sin(phi1);
04257 p1 = sqrt(1+tanl1*tanl1)/kap1;
04258 the1 = M_PI/2-atan(tanl1);
04259 } else if(jj == 2) {
04260 phi2 = (*iter_trk)->getFFi0();
04261 r2 = (*iter_trk)->getFDr();
04262 z2 = (*iter_trk)->getFDz();
04263 kap2 = (*iter_trk)->getFCpa();
04264 tanl2 = (*iter_trk)->getFTanl();
04265 x2 = r1*cos(phi2);
04266 y2 = r1*sin(phi2);
04267 p2 = sqrt(1+tanl2*tanl2)/kap1;
04268 the2 = M_PI/2-atan(tanl2);
04269 }
04270 }
04271 }
04272 if(ntuple_&4) {
04273 m_delx = x1 - x2;
04274 m_dely = y1 - y2;
04275 m_delz = z1 - z2;
04276 m_delthe = the1 + the2;
04277 m_delphi = phi1- phi2;
04278 m_delp = p1 - p2;
04279 StatusCode sc2 = m_nt2->write();
04280 if( sc2.isFailure() ) cout<<"Ntuple2 filling failed!"<<endl;
04281 }
04282 delete [] order;
04283 delete [] rCont;
04284 delete [] rGen;
04285 delete [] rOM;
04286
04287 if (debug_ == 4)
04288 cout << "Kalfitting finished " << std::endl;
04289 }
04290
04291
04292 void KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew(void)
04293 {
04294
04295 MsgStream log(msgSvc(), name());
04296 double Pt_threshold(0.3);
04297 Hep3Vector IP(0,0,0);
04298
04299 vector<MdcRec_trk>* mdcMgr = MdcRecTrkCol::getMdcRecTrkCol();
04300 vector<MdcRec_trk_add>* mdc_addMgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
04301 vector<MdcRec_wirhit>* whMgr = MdcRecWirhitCol::getMdcRecWirhitCol();
04302
04303
04304 if ( !&whMgr ) return;
04305
04306
04307 int ntrk = mdcMgr->size();
04308
04309
04310 int nhits = whMgr->size();
04311
04312
04313
04314
04315 DataObject *aReconEvent;
04316 eventSvc()->findObject("/Event/Recon",aReconEvent);
04317 if(!aReconEvent) {
04318
04319 ReconEvent* recevt = new ReconEvent;
04320 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
04321 if(sc!=StatusCode::SUCCESS) {
04322 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
04323 return;
04324 }
04325 }
04326
04327 RecMdcKalTrackCol* kalcol = new RecMdcKalTrackCol;
04328 RecMdcKalHelixSegCol *segcol =new RecMdcKalHelixSegCol;
04329
04330 log << MSG::INFO << "beginning to make RecMdcKalTrackCol" <<endreq;
04331
04332
04333 MdcRec_trk& TrasanTRK = *(mdcMgr->begin());
04334 MdcRec_trk_add& TrasanTRK_add = *(mdc_addMgr->begin());
04335
04336
04337
04338
04339
04340
04341 int type(0);
04342 if ((TrasanTRK_add.decision & 32) == 32 ||
04343 (TrasanTRK_add.decision & 64) == 64) type = 1;
04344
04345
04346 HepPoint3D x(TrasanTRK.pivot[0],
04347 TrasanTRK.pivot[1],
04348 TrasanTRK.pivot[2]);
04349
04350 HepVector a(5);
04351 for(int i = 0; i < 5; i++)
04352 a[i] = TrasanTRK.helix[i];
04353
04354 HepSymMatrix ea(5);
04355 for(int i = 0, k = 0; i < 5; i++) {
04356 for(int j = 0; j <= i; j++) {
04357 ea[i][j] = matrixg_*TrasanTRK.error[k++];
04358 ea[j][i] = ea[i][j];
04359 }
04360 }
04361
04362 KalFitTrack::setInitMatrix(ea);
04363
04364 double fiTerm = TrasanTRK.fiTerm;
04365 int way(1);
04366
04367 KalFitTrack track_lead = KalFitTrack(x, a, ea, lead_, 0, 0);
04368 track_lead.bFieldZ(KalFitTrack::Bznom_);
04369
04370 int hit_asso(0);
04371
04372 int trasqual = TrasanTRK_add.quality;
04373 if(debug_ == 4) cout<<"kalman_fitting>trasqual... "<<trasqual<<endl;
04374 if (trasqual) return;
04375
04376 int inlyr(999), outlyr(-1);
04377 int* rStat = new int[43];
04378 for(int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
04379 std::vector<MdcRec_wirhit*> pt=TrasanTRK.hitcol;
04380 int hit_in(0);
04381 if(debug_ == 4) cout<<"*********Pt size****"<< pt.size()<<endl;
04382
04383 int Num[43] = {0};
04384 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
04385 ii != pt.begin()-1; ii--) {
04386 Num[(*ii)->geo->Lyr()->Id()]++;
04387 }
04388
04389 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
04390 ii != pt.begin()-1; ii--) {
04391
04392 hit_asso++;
04393 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
04394 if (debug_ >0)
04395 cout << "WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
04396 << " hits in the layer "
04397 << (*ii)->geo->Lyr()->Id() << std::endl;
04398 continue;
04399 }
04400
04401 hit_in++;
04402 MdcRec_wirhit & rechit = **ii;
04403 double dist[2] = {rechit.ddl, rechit.ddr};
04404 double erdist[2] = {rechit.erddl, rechit.erddr};
04405 const MdcGeoWire* geo = rechit.geo;
04406
04407 int lr_decision(0);
04408 if (KalFitTrack::LR_ == 1){
04409 if (rechit.lr==2 || rechit.lr==0) lr_decision=-1;
04410
04411 else if (rechit.lr==1) lr_decision=1;
04412 }
04413
04414 int ind(geo->Lyr()->Id());
04415 track_lead.appendHitsMdc( KalFitHitMdc(rechit.id,
04416 lr_decision, rechit.tdc,
04417 dist, erdist,
04418 _wire+(geo->Id()), rechit.rechitptr));
04419
04420 rStat[ind]++;
04421 if (inlyr>ind) inlyr = ind;
04422 if (outlyr<ind) outlyr = ind;
04423 }
04424
04425 int empty_between(0), empty(0);
04426 for (int i= inlyr; i <= outlyr; i++)
04427 if (!rStat[i]) empty_between++;
04428 empty = empty_between+inlyr+(42-outlyr);
04429 delete [] rStat;
04430
04431 if (debug_ == 4)
04432 cout << "**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
04433
04434
04435
04436 track_lead.order_wirhit(0);
04437 track_lead.type(type);
04438 unsigned int nhit = track_lead.HitsMdc().size();
04439 if (nhit<70) {
04440 cout << " ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
04441 return;
04442 }
04443
04444
04445 double KalFitst(0), KalFitax(0), KalFitschi2(0);
04446
04447 Hep3Vector outer_pivot(track_lead.x(fiTerm));
04448
04449 if(debug_ == 4) {
04450 std::cout<<"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.Ea()<<std::endl;
04451 }
04452 track_lead.pivot(outer_pivot);
04453 track_lead.bFieldZ(KalFitTrack::Bznom_);
04454
04455 if (nhit>=3 && !KalFitTrack::LR_)
04456 start_seed(track_lead, lead_, way, TrasanTRK);
04457 HepSymMatrix Eakal(5,0);
04458
04459
04460
04461 double costheta = track_lead.a()[4] / sqrt(1.0 + track_lead.a()[4]*track_lead.a()[4]);
04462 if( (1.0/fabs(track_lead.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
04463 choice_ = 6;
04464 }
04465
04466 init_matrix(choice_,TrasanTRK, Eakal);
04467
04468
04469
04470 if (debug_ == 4){
04471 std::cout << "from Mdc Pattern Recognition: " << std::endl;
04472
04473 Helix work(track_lead.pivot(),
04474 track_lead.a(),
04475 track_lead.Ea());
04476 work.pivot(IP);
04477 std::cout << " dr = " << work.a()[0]
04478 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
04479 std::cout << " phi0 = " << work.a()[1]
04480 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
04481 std::cout << " PT = " << 1/work.a()[2]
04482 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
04483 std::cout << " dz = " << work.a()[3]
04484 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
04485 std::cout << " tanl = " << work.a()[4]
04486 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
04487 }
04488 filter_fwd_calib(track_lead, lead_, way, Eakal);
04489 track_lead.update_forMdc();
04490
04491
04492 if (debug_ == 4) {
04493 cout << " Mdc FIRST KALMAN FIT " << std::endl;
04494 Helix work1(track_lead.pivot(),
04495 track_lead.a(),
04496 track_lead.Ea());
04497 work1.pivot(IP);
04498 cout << " dr = " << work1.a()[0]
04499 << ", Er_dr = " << sqrt(work1.Ea()[0][0]) << std::endl;
04500 cout << " phi0 = " << work1.a()[1]
04501 << ", Er_phi0 = " << sqrt(work1.Ea()[1][1]) << std::endl;
04502 cout << " PT = " << 1/work1.a()[2]
04503 << ", Er_kappa = " << sqrt(work1.Ea()[2][2]) << std::endl;
04504 cout << " dz = " << work1.a()[3]
04505 << ", Er_dz = " << sqrt(work1.Ea()[3][3]) << std::endl;
04506 cout << " tanl = " << work1.a()[4]
04507 << ", Er_tanl = " << sqrt(work1.Ea()[4][4]) << std::endl;
04508 }
04509
04510
04511 RecMdcKalTrack* kaltrk = new RecMdcKalTrack;
04512
04513
04514 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
04515
04516
04517 StatusCode kalsc;
04518
04519 DataObject *aRecKalEvent;
04520 eventSvc()->findObject("/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
04521 if(aRecKalEvent!=NULL) {
04522
04523 kalsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalTrackCol");
04524 if(kalsc != StatusCode::SUCCESS) {
04525 log << MSG::FATAL << "Could not unregister RecMdcKalTrack collection" << endreq;
04526 return;
04527 }
04528 }
04529
04530 kalsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalTrackCol", kalcol);
04531 if( kalsc.isFailure()) {
04532 log << MSG::FATAL << "Could not register RecMdcKalTrack" << endreq;
04533 return;
04534 }
04535 log << MSG::INFO << "RecMdcKalTrackCol registered successfully!" <<endreq;
04536
04537
04538
04539 StatusCode segsc;
04540
04541 DataObject *aRecKalSegEvent;
04542 eventSvc()->findObject("/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
04543 if(aRecKalSegEvent!=NULL) {
04544
04545 segsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalHelixSegCol");
04546 if(segsc != StatusCode::SUCCESS) {
04547 log << MSG::FATAL << "Could not unregister RecMdcKalHelixSegCol collection" << endreq;
04548 return;
04549 }
04550 }
04551
04552 segsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalHelixSegCol", segcol);
04553 if( segsc.isFailure() ) {
04554 log << MSG::FATAL << "Could not register RecMdcKalHelixSeg" << endreq;
04555 return;
04556 }
04557 log << MSG::INFO << "RecMdcKalHelixSegCol registered successfully!" <<endreq;
04558
04559
04560 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),phi1(0.),phi2(0.),p1(0.),p2(0.);
04561 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
04562
04563
04564 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
04565 if (!kaltrkCol) {
04566 log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endreq;
04567 return;
04568 }
04569 log << MSG::INFO << "Begin to check RecMdcKalTrackCol"<<endreq;
04570 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
04571 for( int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
04572 log << MSG::DEBUG << "retrieved MDC Kalmantrack:"
04573 << "Track Id: " << (*iter_trk)->getTrackId()
04574 << " Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
04575 << " Length of the track: "<< (*iter_trk)->getLength(2)
04576 << " Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
04577 << " Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
04578 <<" "<< (*iter_trk)->getChisq(1,2) << endreq
04579 << "Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
04580 <<" "<< (*iter_trk)->getNdf(1,2) << endreq
04581 << "Kappa " << (*iter_trk)->getZHelix()[2]
04582 << "zhelixmu "<<(*iter_trk)->getZHelixMu()
04583 << endreq;
04584
04585 HelixSegRefVec gothelixsegs = (*iter_trk)->getVecHelixSegs();
04586 if(debug_ == 4) {
04587 std::cout<<"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
04588 }
04589
04590 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
04591 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
04592 if(debug_ == 4) {
04593 std::cout<<"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
04594 std::cout<<"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
04595 std::cout<<"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
04596 std::cout<<"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
04597 std::cout<<"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
04598 std::cout<<"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
04599 }
04600 }
04601 for( int i = 0; i<43; i++) {
04602 log << MSG::DEBUG << "retrieved pathl["<<i<<"]= "
04603 << (*iter_trk)->getPathl(i) <<endreq;
04604 }
04605
04606 if(ntuple_&1) {
04607 m_trackid = (*iter_trk)->getTrackId();
04608 for( int jj =0, iii=0; jj<5; jj++) {
04609 m_length[jj] = (*iter_trk)->getLength(jj);
04610 m_tof[jj] = (*iter_trk)->getTof(jj);
04611 m_nhits[jj] = (*iter_trk)->getNhits(jj);
04612 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
04613 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
04614 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
04615 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
04616 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
04617 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
04618 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
04619 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
04620 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
04621 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
04622 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
04623 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
04624 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
04625 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
04626 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
04627 if(ntuple_&32) {
04628 for(int kk=0; kk<=jj; kk++,iii++) {
04629 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
04630 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
04631 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
04632 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
04633 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
04634 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
04635 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
04636 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
04637 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
04638 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
04639 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
04640 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
04641 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
04642 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
04643 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
04644 }
04645
04646 }
04647 }
04648
04649
04650 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
04651 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
04652 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
04653 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
04654 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
04655 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
04656 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
04657 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
04658 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
04659 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
04660
04661 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
04662 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
04663 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
04664 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
04665 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
04666 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
04667 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
04668 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
04669 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
04670 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
04671
04672 m_stat[0][0] = (*iter_trk)->getStat(0,0);
04673 m_stat[1][0] = (*iter_trk)->getStat(0,1);
04674 m_stat[2][0] = (*iter_trk)->getStat(0,2);
04675 m_stat[3][0] = (*iter_trk)->getStat(0,3);
04676 m_stat[4][0] = (*iter_trk)->getStat(0,4);
04677 m_stat[0][1] = (*iter_trk)->getStat(1,0);
04678 m_stat[1][1] = (*iter_trk)->getStat(1,1);
04679 m_stat[2][1] = (*iter_trk)->getStat(1,2);
04680 m_stat[3][1] = (*iter_trk)->getStat(1,3);
04681 m_stat[4][1] = (*iter_trk)->getStat(1,4);
04682
04683 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
04684 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
04685 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
04686 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
04687 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
04688
04689 m_zpt = 1/m_zhelix[2];
04690 m_zpte = 1/m_zhelixe[2];
04691 m_zptmu = 1/m_zhelixmu[2];
04692 m_zptk = 1/m_zhelixk[2];
04693 m_zptp = 1/m_zhelixp[2];
04694
04695 m_fpt = 1/m_fhelix[2];
04696 m_fpte = 1/m_fhelixe[2];
04697 m_fptmu = 1/m_fhelixmu[2];
04698 m_fptk = 1/m_fhelixk[2];
04699 m_fptp = 1/m_fhelixp[2];
04700
04701 m_lpt = 1/m_lhelix[2];
04702 m_lpte = 1/m_lhelixe[2];
04703 m_lptmu = 1/m_lhelixmu[2];
04704 m_lptk = 1/m_lhelixk[2];
04705 m_lptp = 1/m_lhelixp[2];
04706
04707 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
04708 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
04709 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
04710 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
04711 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
04712
04713 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
04714 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
04715 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
04716 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
04717 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
04718 if(ntuple_&32) {
04719 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
04720 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
04721 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
04722 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
04723 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
04724 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
04725 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
04726 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
04727 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
04728 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
04729 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
04730 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
04731 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
04732 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
04733 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
04734 }
04735
04736 StatusCode sc1 = m_nt1->write();
04737 if( sc1.isFailure() ) cout<<"Ntuple1 filling failed!"<<endl;
04738 }
04739
04740 if(ntuple_&4) {
04741 if(jj == 1) {
04742 phi1 = (*iter_trk)->getFFi0();
04743 r1 = (*iter_trk)->getFDr();
04744 z1 = (*iter_trk)->getFDz();
04745 kap1 = (*iter_trk)->getFCpa();
04746 tanl1 = (*iter_trk)->getFTanl();
04747 x1 = r1*cos(phi1);
04748 y1 = r1*sin(phi1);
04749 p1 = sqrt(1+tanl1*tanl1)/kap1;
04750 the1 = M_PI/2-atan(tanl1);
04751 } else if(jj == 2) {
04752 phi2 = (*iter_trk)->getFFi0();
04753 r2 = (*iter_trk)->getFDr();
04754 z2 = (*iter_trk)->getFDz();
04755 kap2 = (*iter_trk)->getFCpa();
04756 tanl2 = (*iter_trk)->getFTanl();
04757 x2 = r1*cos(phi2);
04758 y2 = r1*sin(phi2);
04759 p2 = sqrt(1+tanl2*tanl2)/kap1;
04760 the2 = M_PI/2-atan(tanl2);
04761 }
04762 }
04763 }
04764 if(ntuple_&4) {
04765 m_delx = x1 - x2;
04766 m_dely = y1 - y2;
04767 m_delz = z1 - z2;
04768 m_delthe = the1 + the2;
04769 m_delphi = phi1- phi2;
04770 m_delp = p1 - p2;
04771 StatusCode sc2 = m_nt2->write();
04772 if( sc2.isFailure() ) cout<<"Ntuple2 filling failed!"<<endl;
04773 }
04774
04775 if (debug_ == 4)
04776 cout << "Kalfitting finished " << std::endl;
04777 }
04778
04779
04780
04781 void KalFitAlg::kalman_fitting_csmalign(void)
04782 {
04783
04784 MsgStream log(msgSvc(), name());
04785 double Pt_threshold(0.3);
04786 Hep3Vector IP(0,0,0);
04787
04788 vector<MdcRec_trk>* mdcMgr = MdcRecTrkCol::getMdcRecTrkCol();
04789 vector<MdcRec_trk_add>* mdc_addMgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
04790 vector<MdcRec_wirhit>* whMgr = MdcRecWirhitCol::getMdcRecWirhitCol();
04791
04792
04793 if ( !&whMgr ) return;
04794
04795
04796 int ntrk = mdcMgr->size();
04797
04798
04799 int nhits = whMgr->size();
04800
04801
04802
04803 double* rY = new double[ntrk];
04804 double* rfiTerm = new double[ntrk];
04805 double* rPt = new double[ntrk];
04806 int* rOM = new int[ntrk];
04807 unsigned int* order = new unsigned int[ntrk];
04808 unsigned int* rCont = new unsigned int[ntrk];
04809 unsigned int* rGen = new unsigned int[ntrk];
04810
04811 int index = 0;
04812 Hep3Vector csmp3[2];
04813 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
04814 end = mdcMgr->end(); it != end; it++) {
04815
04816 rfiTerm[index]=it->fiTerm;
04817
04818
04819 rPt[index] = 0;
04820 if (it->helix[2])
04821 rPt[index] = 1 / fabs(it->helix[2]);
04822 if(debug_ == 4) cout<<"rPt...[ "<<index<<" ]...."<< rPt[index] <<endl;
04823 if(rPt[index] < 0) rPt[index] = DBL_MAX;
04824
04825 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
04826 if(debug_ == 4) cout<<"ppt size: "<< pt.size()<<endl;
04827 int outermost(-1);
04828 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
04829 ii !=pt.begin()-1; ii--) {
04830 int lyr((*ii)->geo->Lyr()->Id());
04831 if (outermost < lyr) {
04832 outermost = lyr;
04833 rY[index] = (*ii)->geo->Forward().y();
04834 }
04835 if(debug_ == 4) cout<<"outmost: "<<outermost<<" lyr: "<<lyr<<endl;
04836 }
04837 rOM[index] = outermost;
04838 order[index] = index;
04839 ++index;
04840 }
04841
04842
04843 for (int j, k = ntrk - 1; k >= 0; k = j){
04844 j = -1;
04845 for(int i = 1; i <= k; i++)
04846 if(rY[i - 1] < rY[i]){
04847 j = i - 1;
04848 std::swap(order[i], order[j]);
04849 std::swap(rY[i], rY[j]);
04850 std::swap(rOM[i], rOM[j]);
04851 std::swap(rCont[i], rCont[j]);
04852 std::swap(rGen[i], rGen[j]);
04853 }
04854 }
04855 delete [] rPt;
04856 delete [] rY;
04857 delete [] rfiTerm;
04858
04859 int newcount(0);
04860
04861 DataObject *aReconEvent;
04862 eventSvc()->findObject("/Event/Recon",aReconEvent);
04863 if(!aReconEvent) {
04864
04865 ReconEvent* recevt = new ReconEvent;
04866 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
04867 if(sc!=StatusCode::SUCCESS) {
04868 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
04869 return;
04870 }
04871 }
04872
04873 RecMdcKalTrackCol* kalcol = new RecMdcKalTrackCol;
04874 RecMdcKalHelixSegCol *segcol =new RecMdcKalHelixSegCol;
04875
04876 log << MSG::INFO << "beginning to make RecMdcKalTrackCol" <<endreq;
04877
04878
04879
04880
04881
04882
04883 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
04884 MdcRec_trk_add& TrasanTRK_add = *(mdc_addMgr->begin()+order[1]);
04885
04886
04887
04888
04889
04890 newcount++;
04891 if (debug_ == 4)
04892 cout << "******* KalFit NUMBER : " << newcount << std::endl;
04893
04894
04895 int type(0);
04896 if ((TrasanTRK_add.decision & 32) == 32 ||
04897 (TrasanTRK_add.decision & 64) == 64) type = 1;
04898
04899
04900 HepPoint3D x(TrasanTRK.pivot[0],
04901 TrasanTRK.pivot[1],
04902 TrasanTRK.pivot[2]);
04903
04904 HepVector a(5);
04905 for(int i = 0; i < 5; i++)
04906 a[i] = TrasanTRK.helix[i];
04907
04908 HepSymMatrix ea(5);
04909 for(int i = 0, k = 0; i < 5; i++) {
04910 for(int j = 0; j <= i; j++) {
04911 ea[i][j] = matrixg_*TrasanTRK.error[k++];
04912 ea[j][i] = ea[i][j];
04913 }
04914 }
04915
04916 KalFitTrack::setInitMatrix(ea);
04917
04918 double fiTerm = TrasanTRK.fiTerm;
04919 int way(1);
04920
04921 KalFitTrack track_lead = KalFitTrack(x, a, ea, lead_, 0, 0);
04922 track_lead.bFieldZ(KalFitTrack::Bznom_);
04923
04924 int hit_asso(0);
04925 for(int l = 0; l < ntrk; l++) {
04926 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
04927 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
04928
04929 int trasqual = TrasanTRK_add1.quality;
04930 if(debug_ == 4) cout<<"kalman_fitting>trasqual... "<<trasqual<<endl;
04931 if (trasqual) continue;
04932
04933 int inlyr(999), outlyr(-1);
04934 int* rStat = new int[43];
04935 for(int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
04936 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.hitcol;
04937 int hit_in(0);
04938 if(debug_ == 4) cout<<"*********Pt size****"<< pt.size()<<endl;
04939
04940 int Num[43] = {0};
04941 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
04942 ii != pt.begin()-1; ii--) {
04943 Num[(*ii)->geo->Lyr()->Id()]++;
04944 }
04945
04946 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
04947 ii != pt.begin()-1; ii--) {
04948
04949 hit_asso++;
04950 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
04951 if (debug_ >0)
04952 cout << "WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
04953 << " hits in the layer "
04954 << (*ii)->geo->Lyr()->Id() << std::endl;
04955 continue;
04956 }
04957
04958 hit_in++;
04959 MdcRec_wirhit & rechit = **ii;
04960 double dist[2] = {rechit.ddl, rechit.ddr};
04961 double erdist[2] = {rechit.erddl, rechit.erddr};
04962 const MdcGeoWire* geo = rechit.geo;
04963
04964 int lr_decision(0);
04965 if (KalFitTrack::LR_ == 1){
04966 if (rechit.lr==2 || rechit.lr==0) lr_decision=-1;
04967
04968 else if (rechit.lr==1) lr_decision=1;
04969 }
04970
04971 int ind(geo->Lyr()->Id());
04972 track_lead.appendHitsMdc( KalFitHitMdc(rechit.id,
04973 lr_decision, rechit.tdc,
04974 dist, erdist,
04975 _wire+(geo->Id()), rechit.rechitptr));
04976
04977 rStat[ind]++;
04978 if (inlyr>ind) inlyr = ind;
04979 if (outlyr<ind) outlyr = ind;
04980 }
04981
04982 int empty_between(0), empty(0);
04983 for (int i= inlyr; i <= outlyr; i++)
04984 if (!rStat[i]) empty_between++;
04985 empty = empty_between+inlyr+(42-outlyr);
04986 delete [] rStat;
04987 }
04988 if (debug_ == 4)
04989 cout << "**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
04990
04991
04992
04993 track_lead.order_wirhit(0);
04994 track_lead.type(type);
04995 unsigned int nhit = track_lead.HitsMdc().size();
04996 if (nhit<70) {
04997 cout << " ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
04998 return;
04999 }
05000
05001
05002 double KalFitst(0), KalFitax(0), KalFitschi2(0);
05003
05004 Hep3Vector outer_pivot(track_lead.x(fiTerm));
05005
05006 if(debug_ == 4) {
05007 std::cout<<"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.Ea()<<std::endl;
05008 }
05009 track_lead.pivot(outer_pivot);
05010 track_lead.bFieldZ(KalFitTrack::Bznom_);
05011
05012 if (nhit>=3 && !KalFitTrack::LR_)
05013 start_seed(track_lead, lead_, way, TrasanTRK);
05014 HepSymMatrix Eakal(5,0);
05015
05016
05017
05018 double costheta = track_lead.a()[4] / sqrt(1.0 + track_lead.a()[4]*track_lead.a()[4]);
05019 if( (1.0/fabs(track_lead.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
05020 choice_ = 6;
05021 }
05022
05023 init_matrix(choice_,TrasanTRK, Eakal);
05024
05025
05026
05027 if (debug_ == 4){
05028 std::cout << "from Mdc Pattern Recognition: " << std::endl;
05029
05030 Helix work(track_lead.pivot(),
05031 track_lead.a(),
05032 track_lead.Ea());
05033 work.pivot(IP);
05034 std::cout << " dr = " << work.a()[0]
05035 << ", Er_dr = " << sqrt(work.Ea()[0][0]) << std::endl;
05036 std::cout << " phi0 = " << work.a()[1]
05037 << ", Er_phi0 = " << sqrt(work.Ea()[1][1]) << std::endl;
05038 std::cout << " PT = " << 1/work.a()[2]
05039 << ", Er_kappa = " << sqrt(work.Ea()[2][2]) << std::endl;
05040 std::cout << " dz = " << work.a()[3]
05041 << ", Er_dz = " << sqrt(work.Ea()[3][3]) << std::endl;
05042 std::cout << " tanl = " << work.a()[4]
05043 << ", Er_tanl = " << sqrt(work.Ea()[4][4]) << std::endl;
05044 }
05045 filter_fwd_calib(track_lead, lead_, way, Eakal);
05046 track_lead.update_forMdc();
05047
05048
05049 if (debug_ == 4) {
05050 cout << " Mdc FIRST KALMAN FIT " << std::endl;
05051 Helix work1(track_lead.pivot(),
05052 track_lead.a(),
05053 track_lead.Ea());
05054 work1.pivot(IP);
05055 cout << " dr = " << work1.a()[0]
05056 << ", Er_dr = " << sqrt(work1.Ea()[0][0]) << std::endl;
05057 cout << " phi0 = " << work1.a()[1]
05058 << ", Er_phi0 = " << sqrt(work1.Ea()[1][1]) << std::endl;
05059 cout << " PT = " << 1/work1.a()[2]
05060 << ", Er_kappa = " << sqrt(work1.Ea()[2][2]) << std::endl;
05061 cout << " dz = " << work1.a()[3]
05062 << ", Er_dz = " << sqrt(work1.Ea()[3][3]) << std::endl;
05063 cout << " tanl = " << work1.a()[4]
05064 << ", Er_tanl = " << sqrt(work1.Ea()[4][4]) << std::endl;
05065 }
05066
05067
05068 RecMdcKalTrack* kaltrk = new RecMdcKalTrack;
05069
05070
05071 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
05072
05073
05074
05075 StatusCode kalsc;
05076
05077 DataObject *aRecKalEvent;
05078 eventSvc()->findObject("/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
05079 if(aRecKalEvent!=NULL) {
05080
05081 kalsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalTrackCol");
05082 if(kalsc != StatusCode::SUCCESS) {
05083 log << MSG::FATAL << "Could not unregister RecMdcKalTrack collection" << endreq;
05084 return;
05085 }
05086 }
05087
05088 kalsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalTrackCol", kalcol);
05089 if( kalsc.isFailure()) {
05090 log << MSG::FATAL << "Could not register RecMdcKalTrack" << endreq;
05091 return;
05092 }
05093 log << MSG::INFO << "RecMdcKalTrackCol registered successfully!" <<endreq;
05094
05095
05096
05097 StatusCode segsc;
05098
05099 DataObject *aRecKalSegEvent;
05100 eventSvc()->findObject("/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
05101 if(aRecKalSegEvent!=NULL) {
05102
05103 segsc = eventSvc()->unregisterObject("/Event/Recon/RecMdcKalHelixSegCol");
05104 if(segsc != StatusCode::SUCCESS) {
05105 log << MSG::FATAL << "Could not unregister RecMdcKalHelixSegCol collection" << endreq;
05106 return;
05107 }
05108 }
05109
05110 segsc = eventSvc()->registerObject("/Event/Recon/RecMdcKalHelixSegCol", segcol);
05111 if( segsc.isFailure() ) {
05112 log << MSG::FATAL << "Could not register RecMdcKalHelixSeg" << endreq;
05113 return;
05114 }
05115 log << MSG::INFO << "RecMdcKalHelixSegCol registered successfully!" <<endreq;
05116
05117
05118 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),phi1(0.),phi2(0.),p1(0.),p2(0.);
05119 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
05120
05121
05122 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
05123 if (!kaltrkCol) {
05124 log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endreq;
05125 return;
05126 }
05127 log << MSG::INFO << "Begin to check RecMdcKalTrackCol"<<endreq;
05128 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
05129 for( int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
05130 log << MSG::DEBUG << "retrieved MDC Kalmantrack:"
05131 << "Track Id: " << (*iter_trk)->getTrackId()
05132 << " Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
05133 << " Length of the track: "<< (*iter_trk)->getLength(2)
05134 << " Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
05135 << " Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
05136 <<" "<< (*iter_trk)->getChisq(1,2) << endreq
05137 << "Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
05138 <<" "<< (*iter_trk)->getNdf(1,2) << endreq
05139 << "Kappa " << (*iter_trk)->getZHelix()[2]
05140 << "zhelixmu "<<(*iter_trk)->getZHelixMu()
05141 << endreq;
05142
05143 HelixSegRefVec gothelixsegs = (*iter_trk)->getVecHelixSegs();
05144 if(debug_ == 4) {
05145 std::cout<<"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
05146 }
05147
05148 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
05149 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
05150 if(debug_ == 4) {
05151 std::cout<<"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
05152 std::cout<<"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
05153 std::cout<<"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
05154 std::cout<<"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
05155 std::cout<<"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
05156 std::cout<<"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
05157 }
05158 }
05159 for( int i = 0; i<43; i++) {
05160 log << MSG::DEBUG << "retrieved pathl["<<i<<"]= "
05161 << (*iter_trk)->getPathl(i) <<endreq;
05162 }
05163
05164 if(ntuple_&1) {
05165 m_trackid = (*iter_trk)->getTrackId();
05166 for( int jj =0, iii=0; jj<5; jj++) {
05167 m_length[jj] = (*iter_trk)->getLength(jj);
05168 m_tof[jj] = (*iter_trk)->getTof(jj);
05169 m_nhits[jj] = (*iter_trk)->getNhits(jj);
05170 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
05171 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
05172 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
05173 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
05174 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
05175 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
05176 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
05177 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
05178 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
05179 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
05180 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
05181 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
05182 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
05183 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
05184 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
05185 if(ntuple_&32) {
05186 for(int kk=0; kk<=jj; kk++,iii++) {
05187 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
05188 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
05189 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
05190 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
05191 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
05192 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
05193 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
05194 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
05195 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
05196 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
05197 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
05198 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
05199 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
05200 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
05201 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
05202 }
05203
05204 }
05205 }
05206
05207
05208 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
05209 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
05210 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
05211 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
05212 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
05213 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
05214 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
05215 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
05216 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
05217 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
05218
05219 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
05220 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
05221 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
05222 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
05223 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
05224 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
05225 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
05226 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
05227 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
05228 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
05229
05230 m_stat[0][0] = (*iter_trk)->getStat(0,0);
05231 m_stat[1][0] = (*iter_trk)->getStat(0,1);
05232 m_stat[2][0] = (*iter_trk)->getStat(0,2);
05233 m_stat[3][0] = (*iter_trk)->getStat(0,3);
05234 m_stat[4][0] = (*iter_trk)->getStat(0,4);
05235 m_stat[0][1] = (*iter_trk)->getStat(1,0);
05236 m_stat[1][1] = (*iter_trk)->getStat(1,1);
05237 m_stat[2][1] = (*iter_trk)->getStat(1,2);
05238 m_stat[3][1] = (*iter_trk)->getStat(1,3);
05239 m_stat[4][1] = (*iter_trk)->getStat(1,4);
05240
05241 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
05242 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
05243 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
05244 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
05245 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
05246
05247 m_zpt = 1/m_zhelix[2];
05248 m_zpte = 1/m_zhelixe[2];
05249 m_zptmu = 1/m_zhelixmu[2];
05250 m_zptk = 1/m_zhelixk[2];
05251 m_zptp = 1/m_zhelixp[2];
05252
05253 m_fpt = 1/m_fhelix[2];
05254 m_fpte = 1/m_fhelixe[2];
05255 m_fptmu = 1/m_fhelixmu[2];
05256 m_fptk = 1/m_fhelixk[2];
05257 m_fptp = 1/m_fhelixp[2];
05258
05259 m_lpt = 1/m_lhelix[2];
05260 m_lpte = 1/m_lhelixe[2];
05261 m_lptmu = 1/m_lhelixmu[2];
05262 m_lptk = 1/m_lhelixk[2];
05263 m_lptp = 1/m_lhelixp[2];
05264
05265 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
05266 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
05267 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
05268 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
05269 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
05270
05271 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
05272 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
05273 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
05274 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
05275 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
05276 if(ntuple_&32) {
05277 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
05278 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
05279 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
05280 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
05281 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
05282 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
05283 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
05284 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
05285 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
05286 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
05287 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
05288 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
05289 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
05290 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
05291 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
05292 }
05293
05294 StatusCode sc1 = m_nt1->write();
05295 if( sc1.isFailure() ) cout<<"Ntuple1 filling failed!"<<endl;
05296 }
05297
05298 if(ntuple_&4) {
05299 if(jj == 1) {
05300 phi1 = (*iter_trk)->getFFi0();
05301 r1 = (*iter_trk)->getFDr();
05302 z1 = (*iter_trk)->getFDz();
05303 kap1 = (*iter_trk)->getFCpa();
05304 tanl1 = (*iter_trk)->getFTanl();
05305 x1 = r1*cos(phi1);
05306 y1 = r1*sin(phi1);
05307 p1 = sqrt(1+tanl1*tanl1)/kap1;
05308 the1 = M_PI/2-atan(tanl1);
05309 } else if(jj == 2) {
05310 phi2 = (*iter_trk)->getFFi0();
05311 r2 = (*iter_trk)->getFDr();
05312 z2 = (*iter_trk)->getFDz();
05313 kap2 = (*iter_trk)->getFCpa();
05314 tanl2 = (*iter_trk)->getFTanl();
05315 x2 = r1*cos(phi2);
05316 y2 = r1*sin(phi2);
05317 p2 = sqrt(1+tanl2*tanl2)/kap1;
05318 the2 = M_PI/2-atan(tanl2);
05319 }
05320 }
05321 }
05322 if(ntuple_&4) {
05323 m_delx = x1 - x2;
05324 m_dely = y1 - y2;
05325 m_delz = z1 - z2;
05326 m_delthe = the1 + the2;
05327 m_delphi = phi1- phi2;
05328 m_delp = p1 - p2;
05329 StatusCode sc2 = m_nt2->write();
05330 if( sc2.isFailure() ) cout<<"Ntuple2 filling failed!"<<endl;
05331 }
05332 delete [] order;
05333 delete [] rCont;
05334 delete [] rGen;
05335 delete [] rOM;
05336
05337 if (debug_ == 4)
05338 cout << "Kalfitting finished " << std::endl;
05339 }
05340
05341
05342
05343 void KalFitAlg::complete_track(MdcRec_trk& TrasanTRK,
05344 MdcRec_trk_add& TrasanTRK_add,
05345 KalFitTrack& track_lead,
05346 RecMdcKalTrack* kaltrk,
05347 RecMdcKalTrackCol* kalcol,RecMdcKalHelixSegCol *segcol,int flagsmooth)
05348 {
05349 static int nmass = KalFitTrack::nmass();
05350 int way(1);
05351 MsgStream log(msgSvc(), name());
05352 KalFitTrack track_first(track_lead);
05353 KalFitTrack track_ip(track_lead);
05354 innerwall(track_ip, lead_, way);
05355
05356 fillTds_ip(TrasanTRK, track_ip, kaltrk, lead_);
05357
05358
05359 fillTds_lead(TrasanTRK, track_first, kaltrk, lead_);
05360
05361 double pp = track_lead.momentum().mag();
05362
05363 if(!(i_front_<0)){
05364
05365 for(int l_mass = 0, flg = 1; l_mass < nmass;
05366 l_mass++, flg <<= 1) {
05367
05368 if (!(mhyp_ & flg)) continue;
05369 if (l_mass == lead_) continue;
05370
05371
05372 if ((lead_ != 0 && l_mass==0 && pp > pe_cut_) ||
05373 (lead_ != 1 && l_mass==1 && pp > pmu_cut_) ||
05374 (lead_ != 2 && l_mass==2 && pp > ppi_cut_) ||
05375 (lead_ != 3 && l_mass==3 && pp > pk_cut_) ||
05376 (lead_ != 4 && l_mass==4 && pp > pp_cut_))
05377 continue;
05378 if(debug_ == 4) cout<<"complete_track..REFIT ASSUMPION " << l_mass << endl;
05379
05380 double chiSq = 0;
05381 int nhits = 0;
05382
05383
05384 HepPoint3D x_trasan(TrasanTRK.pivot[0],
05385 TrasanTRK.pivot[1],
05386 TrasanTRK.pivot[2]);
05387 HepVector a_trasan(5);
05388 for(int i = 0; i < 5; i++)
05389 a_trasan[i] = TrasanTRK.helix[i];
05390
05391 HepSymMatrix ea_trasan(5);
05392 for(int i = 0, k = 0; i < 5; i++) {
05393 for(int j = 0; j <= i; j++) {
05394 ea_trasan[i][j] = matrixg_*TrasanTRK.error[k++];
05395 ea_trasan[j][i] = ea_trasan[i][j];
05396 }
05397 }
05398
05399 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
05400 track.HitsMdc(track_lead.HitsMdc());
05401 double fiTerm = TrasanTRK.fiTerm;
05402 track.pivot(track.x(fiTerm));
05403 HepSymMatrix Eakal(5,0);
05404
05405 double costheta = track.a()[4] / sqrt(1.0 + track.a()[4]*track.a()[4]);
05406 if( (1.0/fabs(track.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
05407 choice_ = 6;
05408 }
05409
05410 init_matrix(choice_,TrasanTRK, Eakal);
05411 filter_fwd_anal(track, l_mass, way, Eakal);
05412 KalFitTrack track_z(track);
05414 innerwall(track_z, l_mass, way);
05415 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
05416
05417 fillTds(TrasanTRK, track, kaltrk, l_mass);
05418 }
05419 }
05420
05421
05422
05423 if (enhance_) {
05424
05425 HepPoint3D x_first(0, 0, 0);
05426 HepVector a_first(kaltrk->getFHelix());
05427 HepSymMatrix ea_first(kaltrk->getFError());
05428 HepVector fac(5);
05429 fac[0]=fac_h1_; fac[1]=fac_h2_; fac[2]=fac_h3_; fac[3]=fac_h4_; fac[4]=fac_h5_;
05430 for(int i = 0; i < 5; i++)
05431 for(int j = 0; j <= i; j++)
05432 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
05433 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
05434 }
05435
05436
05437 KalFitTrack track_back(track_lead);
05438 if (debug_ == 4) {
05439 cout << " Backward fitting flag:" << back_<< endl;
05440 cout << "track_back pivot " << track_back.pivot()
05441 << " track_lead kappa " << track_lead.a()[2]
05442 <<endl;
05443 }
05444
05445 if (back_ && track_lead.a()[2] != 0 &&
05446 1/fabs(track_lead.a()[2]) > pT_) {
05447 track_back.HitsMdc(track_lead.HitsMdc());
05448
05449 if (KalFitTrack::tofall_) {
05450 double p_kaon(0), p_proton(0);
05451 if (!(kaltrk->getStat(0,3))) {
05452 p_kaon = 1 / fabs(kaltrk->getZHelixK()[2]) *
05453 sqrt(1 + kaltrk->getZHelixK()[4]*kaltrk->getZHelixK()[4]);
05454 track_back.p_kaon(p_kaon);
05455 } else {
05456 p_kaon = 1 / fabs(track_back.a()[2]) *
05457 sqrt(1 + track_back.a()[4]*track_back.a()[4]);
05458 track_back.p_kaon(p_kaon);
05459 }
05460 if (!(kaltrk->getStat(0,4))) {
05461 p_proton = 1 / fabs(kaltrk->getZHelixP()[2]) *
05462 sqrt(1 + kaltrk->getZHelixP()[4]*kaltrk->getZHelixP()[4]);
05463 track_back.p_proton(p_proton);
05464 } else {
05465 p_proton = 1 / fabs(track_back.a()[2]) *
05466 sqrt(1 + track_back.a()[4]*track_back.a()[4]);
05467 track_back.p_proton(p_proton);
05468 }
05469 }
05470
05471 if (!(i_back_<0)) {
05472
05473 for(int l_mass = 0; l_mass < nmass; l_mass++) {
05474
05475 KalFitTrack track_seed(track_back);
05476 track_seed.chgmass(l_mass);
05477
05478
05479
05480
05481
05482 smoother_anal(track_seed, -way);
05483
05484
05485
05486 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass, segcol, 1);
05487
05488 }
05489 } else {
05490 smoother_anal(track_back, -way);
05491
05492
05493
05494
05495 fillTds_back(track_back, kaltrk, TrasanTRK, lead_, segcol, 1);
05496 }
05497 }
05498
05499
05500
05501
05502
05503
05504
05505
05506
05507
05508
05509
05510 log << MSG::DEBUG << "registered MDC Kalmantrack:"
05511 << "Track Id: " << kaltrk->getTrackId()
05512 << " Mass of the fit: "<< kaltrk->getMass(2)<< endreq
05513 << "Length of the track: "<< kaltrk->getLength(2)
05514 << " Tof of the track: "<< kaltrk->getTof(2) << endreq
05515 << "Chisq of the fit: "<< kaltrk->getChisq(0,2)
05516 <<" "<< kaltrk->getChisq(1,2) << endreq
05517 << "Ndf of the fit: "<< kaltrk->getNdf(0,2)
05518 <<" "<< kaltrk->getNdf(1,2) << endreq
05519 << "Helix " << kaltrk->getZHelix()[2]
05520 <<endreq;
05521
05522 kalcol->push_back(kaltrk);
05523 track_lead.HitsMdc().clear();
05524 }
05525
05526
05527 void KalFitAlg::complete_track(MdcRec_trk& TrasanTRK,
05528 MdcRec_trk_add& TrasanTRK_add,
05529 KalFitTrack& track_lead,
05530 RecMdcKalTrack* kaltrk,
05531 RecMdcKalTrackCol* kalcol,RecMdcKalHelixSegCol *segcol)
05532 {
05533 static int nmass = KalFitTrack::nmass();
05534 int way(1);
05535 MsgStream log(msgSvc(), name());
05536 KalFitTrack track_first(track_lead);
05537 KalFitTrack track_ip(track_lead);
05538
05539 if (debug_ == 4){
05540 cout << "track_first pivot "<<track_first.pivot()<< " helix "<<track_first.a()<<endl;
05541 }
05542 if(usage_==1) innerwall(track_ip, lead_, way);
05543
05544 fillTds_ip(TrasanTRK, track_ip, kaltrk, lead_);
05545
05546 if (debug_ == 4) {
05547 cout << "after inner wall, track_ip pivot "<<track_first.pivot()<< " helix "<<track_first.a()<<endl;
05548 }
05549
05550 fillTds_lead(TrasanTRK, track_first, kaltrk, lead_);
05551
05552
05553 double pp = track_lead.momentum().mag();
05554
05555
05556 if(!(i_front_<0)){
05557
05558 for(int l_mass = 0, flg = 1; l_mass < nmass;
05559 l_mass++, flg <<= 1) {
05560
05561 if (!(mhyp_ & flg)) continue;
05562 if (l_mass == lead_) continue;
05563
05564
05565 if ((lead_ != 0 && l_mass==0 && pp > pe_cut_) ||
05566 (lead_ != 1 && l_mass==1 && pp > pmu_cut_) ||
05567 (lead_ != 2 && l_mass==2 && pp > ppi_cut_) ||
05568 (lead_ != 3 && l_mass==3 && pp > pk_cut_) ||
05569 (lead_ != 4 && l_mass==4 && pp > pp_cut_))
05570 continue;
05571
05572 if(debug_ == 4) {
05573 cout<<"complete_track..REFIT ASSUMPION " << l_mass << endl;
05574 }
05575
05576
05577 double chiSq = 0;
05578 int nhits = 0;
05579
05580
05581 HepPoint3D x_trasan(TrasanTRK.pivot[0],
05582 TrasanTRK.pivot[1],
05583 TrasanTRK.pivot[2]);
05584
05585 HepVector a_trasan(5);
05586 for(int i = 0; i < 5; i++){
05587 a_trasan[i] = TrasanTRK.helix[i];
05588 }
05589
05590 HepSymMatrix ea_trasan(5);
05591 for(int i = 0, k = 0; i < 5; i++) {
05592 for(int j = 0; j <= i; j++) {
05593 ea_trasan[i][j] = matrixg_*TrasanTRK.error[k++];
05594 ea_trasan[j][i] = ea_trasan[i][j];
05595 }
05596 }
05597
05598 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
05599 track.HitsMdc(track_lead.HitsMdc());
05600
05601 double fiTerm = TrasanTRK.fiTerm;
05602 track.pivot(track.x(fiTerm));
05603
05604 HepSymMatrix Eakal(5,0);
05605
05606 double costheta = track_lead.a()[4] / sqrt(1.0 + track_lead.a()[4]*track_lead.a()[4]);
05607 if( (1.0/fabs(track_lead.a()[2]) < pt_cut_ ) && (fabs(costheta)> theta_cut_) ) {
05608 choice_ = 6;
05609 }
05610
05611 init_matrix(choice_, TrasanTRK, Eakal);
05612
05613 filter_fwd_calib(track, l_mass, way, Eakal);
05614
05615 KalFitTrack track_z(track);
05617 innerwall(track_z, l_mass, way);
05618 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
05619
05620 fillTds(TrasanTRK, track, kaltrk, l_mass);
05621 }
05622 }
05623
05624
05625
05626 if (enhance_) {
05627 HepPoint3D x_first(0, 0, 0);
05628 HepVector a_first(kaltrk->getFHelix());
05629 HepSymMatrix ea_first(kaltrk->getFError());
05630 HepVector fac(5);
05631 fac[0]=fac_h1_; fac[1]=fac_h2_; fac[2]=fac_h3_; fac[3]=fac_h4_; fac[4]=fac_h5_;
05632
05633 for(int i = 0; i < 5; i++)
05634 for(int j = 0; j <= i; j++)
05635 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
05636 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
05637 }
05638
05639
05640
05641
05642
05643
05644
05645
05646 KalFitTrack track_back(track_lead);
05647
05648
05649
05650
05651 if (debug_ == 4) {
05652 cout << " Backward fitting flag:" << back_<< endl;
05653 cout << "track_back pivot " << track_back.pivot()
05654 << " track_lead kappa " << track_lead.a()[2]
05655 <<endl;
05656 }
05657
05658 if (back_ && track_lead.a()[2] != 0 &&
05659 1/fabs(track_lead.a()[2]) > pT_) {
05660 track_back.HitsMdc(track_lead.HitsMdc());
05661
05662 if (KalFitTrack::tofall_) {
05663
05664 double p_kaon(0), p_proton(0);
05665
05666 if (!(kaltrk->getStat(0,3))) {
05667 p_kaon = 1 / fabs(kaltrk->getZHelixK()[2]) *
05668 sqrt(1 + kaltrk->getZHelixK()[4]*kaltrk->getZHelixK()[4]);
05669 track_back.p_kaon(p_kaon);
05670 } else {
05671 p_kaon = 1 / fabs(track_back.a()[2]) *
05672 sqrt(1 + track_back.a()[4]*track_back.a()[4]);
05673 track_back.p_kaon(p_kaon);
05674 }
05675 if (!(kaltrk->getStat(0,4))) {
05676 p_proton = 1 / fabs(kaltrk->getZHelixP()[2]) *
05677 sqrt(1 + kaltrk->getZHelixP()[4]*kaltrk->getZHelixP()[4]);
05678 track_back.p_proton(p_proton);
05679 } else {
05680 p_proton = 1 / fabs(track_back.a()[2]) *
05681 sqrt(1 + track_back.a()[4]*track_back.a()[4]);
05682 track_back.p_proton(p_proton);
05683 }
05684
05685 }
05686
05687
05688 if (!(i_back_<0)) {
05689 for(int l_mass = 0; l_mass < nmass; l_mass++) {
05690 KalFitTrack track_seed(track_back);
05691 track_seed.chgmass(l_mass);
05692 smoother_calib(track_seed, -way);
05693
05694 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass,segcol);
05695 }
05696 } else {
05697
05698 smoother_calib(track_back, -way);
05699
05700 fillTds_back(track_back, kaltrk, TrasanTRK, lead_,segcol);
05701
05702 }
05703 }
05704
05705
05706
05707
05708
05709
05710
05711
05712
05713
05714
05715
05716
05717 log << MSG::DEBUG << "registered MDC Kalmantrack:"
05718 << "Track Id: " << kaltrk->getTrackId()
05719 << " Mass of the fit: "<< kaltrk->getMass(2)<< endreq
05720 << "Length of the track: "<< kaltrk->getLength(2)
05721 << " Tof of the track: "<< kaltrk->getTof(2) << endreq
05722 << "Chisq of the fit: "<< kaltrk->getChisq(0,2)
05723 <<" "<< kaltrk->getChisq(1,2) << endreq
05724 << "Ndf of the fit: "<< kaltrk->getNdf(0,2)
05725 <<" "<< kaltrk->getNdf(1,2) << endreq
05726 << "Helix " << kaltrk->getZHelix()[2]
05727 <<endreq;
05728
05729 kalcol->push_back(kaltrk);
05730 track_lead.HitsMdc().clear();
05731 track_back.HelixSegs().clear();
05732
05733 }
05734
05735
05736 void KalFitAlg::init_matrix(MdcRec_trk& trk, HepSymMatrix& Eakal )
05737 {
05738 for ( int i=0; i<5; i++) {
05739 for( int j = 1; j<i+2;j++) {
05740 Eakal(i+1,j) = matrixg_*(trk.error[i*(i+1)/2+j-1]);
05741 Eakal(j,i+1) = Eakal(i+1,j);
05742 }
05743 }
05744
05745 if (debug_ == 4) cout<<"initialised Ea.. "<<Eakal<<endl;
05746 }
05747
05748
05749
05750 void KalFitAlg::init_matrix(int k, MdcRec_trk& trk, HepSymMatrix& Eakal )
05751 {
05752 if(0==k){
05753 for ( int i=0; i<5; i++) {
05754 for( int j = 1; j<i+2;j++) {
05755 Eakal(i+1,j) = matrixg_*(trk.error[i*(i+1)/2+j-1]);
05756 Eakal(j,i+1) = Eakal(i+1,j);
05757 }
05758 Eakal(1,1) = Eakal(1,1)* gain1_;
05759 Eakal(2,2) = Eakal(2,2)* gain2_;
05760 Eakal(3,3) = Eakal(3,3)* gain3_;
05761 Eakal(4,4) = Eakal(4,4)* gain4_;
05762 Eakal(5,5) = Eakal(5,5)* gain5_;
05763 }
05764 }
05765
05766
05767
05768
05769
05770
05771
05772
05773 if(1==k){
05774 Eakal(1,1) = .2;
05775 Eakal(2,2) = .001;
05776 Eakal(3,3) = 1.0;
05777 Eakal(4,4) = 10.0;
05778 Eakal(5,5) = 0.002;
05779 }
05780
05781 if(2==k){
05782 Eakal(1,1) = .2;
05783 Eakal(2,2) = 0.1;
05784 Eakal(3,3) = 1.0;
05785 Eakal(4,4) = 25.0;
05786 Eakal(5,5) = 0.10;
05787 }
05788
05789
05790 if(3==k){
05791 Eakal(1,1) = .2;
05792 Eakal(2,2) = .001;
05793 Eakal(3,3) = 1.0;
05794 Eakal(4,4) = 25.0;
05795 Eakal(5,5) = 0.10;
05796 }
05797
05798 if(4==k){
05799 Eakal(1,1) = .2;
05800 Eakal(2,2) = .01;
05801 Eakal(3,3) = 0.01;
05802 Eakal(4,4) = 1.;
05803 Eakal(5,5) = .01;
05804 }
05805
05806 if(5==k) {
05807 Eakal(1,1) = 2.;
05808 Eakal(2,2) = 0.1;
05809 Eakal(3,3) = 1.;
05810 Eakal(4,4) = 20.;
05811 Eakal(5,5) = 0.1;
05812 }
05813
05814 if(6==k) {
05815 Eakal(1,1) = 0.01;
05816 Eakal(2,2) = 0.01;
05817 Eakal(3,3) = 0.01;
05818 Eakal(4,4) = 100.;
05819 Eakal(5,5) = 0.5;
05820 }
05821
05822 if(k!=0){
05823 Eakal(3,3) = 0.2;
05824 Eakal(1,1) = 1;
05825 Eakal(4,4) = 1;
05826 }
05827
05828 if (debug_ == 4) cout<<"initialised Eakal.. "<<Eakal<<endl;
05829 }
05830
05831
05832
05833
05834 void KalFitAlg::start_seed(KalFitTrack& track, int lead_, int way, MdcRec_trk& TrasanTRK)
05835 {
05836 if (debug_ == 4)
05837 cout << "start_seed begin... " << std::endl;
05838
05839 Hep3Vector x_init(track.pivot());
05840 HepSymMatrix Ea_init(5,0);
05841 Ea_init = track.Ea();
05842 HepVector a_init(5);
05843 a_init = track.a();
05844
05845
05846 unsigned int nhit_included(10);
05847 int LR[8][3] = {
05848 {1,1,1},
05849 {1,1,-1},
05850 {1,-1,1},
05851 {1,-1,-1},
05852 {-1,1,1},
05853 {-1,1,-1},
05854 {-1,-1,1},
05855 {-1,-1,-1}
05856 };
05857
05858 unsigned int nhit = track.HitsMdc().size();
05859 double chi2_min(DBL_MAX);
05860 int i_min(-1);
05861 for (int ktrial = 0; ktrial < 8; ktrial++) {
05862
05863
05864 track.pivot(x_init);
05865 track.a(a_init);
05866 track.Ea(Ea_init);
05867
05868 track.chiSq(0);
05869 track.chiSq_back(0);
05870 track.nchits(0);
05871 track.nster(0);
05872 track.ndf_back(0);
05873
05874 HepSymMatrix Eakal(5,0);
05875
05876 init_matrix(choice_,TrasanTRK, Eakal);
05877
05878 for( unsigned i=0 ; i < nhit; i++ ){
05879 KalFitHitMdc& HitMdc = track.HitMdc(i);
05880 int lr_decision(0);
05881 if (i<3) lr_decision = LR[ktrial][i];
05882 HitMdc.LR(lr_decision);
05883 if (i<nhit_included)
05884 HitMdc.chi2(0);
05885 else
05886 HitMdc.chi2(-1);
05887 }
05888
05889
05890 if(usage_==0) filter_fwd_anal(track, lead_, way, Eakal);
05891 way=999;
05892 if(usage_>0) filter_fwd_calib(track, lead_, way, Eakal);
05893
05894
05895 if (debug_ == 4)
05896 cout << "---- Result for " << ktrial << " case : chi2 = " << track.chiSq()
05897 << ", nhits included = " << track.nchits() << " for nhits available = "
05898 << nhit << std::endl;
05899
05900 if (track.chiSq() < chi2_min &&
05901 (track.nchits() == nhit_included || track.nchits() == nhit)){
05902 chi2_min = track.chiSq();
05903 i_min = ktrial;
05904 }
05905 }
05906
05907
05908 if (debug_ == 4)
05909 cout << "*** i_min = " << i_min << " with a chi2 = " << chi2_min << std::endl;
05910
05911 for( unsigned i=0 ; i < nhit; i++ ){
05912 KalFitHitMdc& HitMdc = track.HitMdc(i);
05913 int lr_decision(0);
05914 if (i_min >= 0 && i < 3)
05915 lr_decision = LR[i_min][i];
05916 HitMdc.LR(lr_decision);
05917 HitMdc.chi2(0);
05918 HitMdc.chi2_back(0);
05919 }
05920 track.pivot(x_init);
05921 track.a(a_init);
05922 track.Ea(Ea_init);
05923 track.chiSq(0);
05924 track.chiSq_back(0);
05925 track.nchits(0);
05926 track.nster(0);
05927 track.ndf_back(0);
05928
05929
05930 if (debug_ == 4) {
05931 for( unsigned i=0 ; i < 3; i++ ){
05932 KalFitHitMdc& HitMdc = track.HitMdc(i);
05933 cout << " LR(" << i << ") = " << HitMdc.LR()
05934 << ", stereo = " << HitMdc.wire().stereo()
05935 << ", layer = " << HitMdc.wire().layer().layerId() << std::endl;
05936 }
05937 }
05938 }
05939
05940
05941
05942
05943 void KalFitAlg::clearTables( ) {
05944
05945 if(debug_ == 4) cout<<"Begining to clear Tables ...."<<endl;
05946 vector<MdcRec_trk>* mdcMgr = MdcRecTrkCol::getMdcRecTrkCol();
05947 vector<MdcRec_trk_add>* mdc_addMgr = MdcRecTrkAddCol::getMdcRecTrkAddCol();
05948 vector<MdcRec_wirhit>* whMgr = MdcRecWirhitCol::getMdcRecWirhitCol();
05949 vector<MdcRec_trk>::iterator tit=mdcMgr->begin();
05950 for( ; tit != mdcMgr->end(); tit++) {
05951 vector<MdcRec_wirhit*>::iterator vit= tit->hitcol.begin() ;
05952 for(; vit != tit->hitcol.end(); vit++) {
05953 delete (*vit);
05954 }
05955 }
05956
05957 mdcMgr->clear();
05958 mdc_addMgr->clear();
05959 whMgr->clear();
05960
05961
05962
05963
05964
05965
05966
05967
05968 }
05969
05970 bool KalFitAlg::order_rechits(const SmartRef<RecMdcHit>& m1, const SmartRef<RecMdcHit>& m2) {
05971 return MdcID::layer(m1->getMdcId()) > MdcID::layer(m2->getMdcId());
05972 }
05973
