16#include "CLHEP/Geometry/Vector3D.h"
17#include "CLHEP/Geometry/Point3D.h"
18#ifndef ENABLE_BACKWARDS_COMPATIBILITY
21#include "CLHEP/Matrix/Vector.h"
22#include "CLHEP/Matrix/SymMatrix.h"
23#include "CLHEP/Matrix/Matrix.h"
24#include "GaudiKernel/MsgStream.h"
25#include "GaudiKernel/AlgFactory.h"
26#include "GaudiKernel/ISvcLocator.h"
27#include "GaudiKernel/IDataProviderSvc.h"
28#include "GaudiKernel/PropertyMgr.h"
29#include "GaudiKernel/IMessageSvc.h"
30#include "GaudiKernel/IDataManagerSvc.h"
31#include "GaudiKernel/SmartDataPtr.h"
32#include "GaudiKernel/PropertyMgr.h"
33#include "GaudiKernel/IJobOptionsSvc.h"
34#include "GaudiKernel/Bootstrap.h"
35#include "GaudiKernel/StatusCode.h"
36#include "GaudiKernel/ContainedObject.h"
37#include "GaudiKernel/SmartRef.h"
38#include "GaudiKernel/SmartRefVector.h"
39#include "GaudiKernel/ObjectVector.h"
42#include "KalFitAlg/helix/Helix.h"
44#include "KalFitAlg/lpav/Lpav.h"
46#include "KalFitAlg/coil/Bfield.h"
47#include "KalFitAlg/KalFitTrack.h"
48#include "KalFitAlg/KalFitHitMdc.h"
49#include "KalFitAlg/KalFitHelixSeg.h"
50#include "KalFitAlg/KalFitElement.h"
51#include "KalFitAlg/KalFitAlg.h"
52#include "McTruth/McParticle.h"
53#include "EventModel/EventHeader.h"
54#include "EvTimeEvent/RecEsTime.h"
55#include "ReconEvent/ReconEvent.h"
56#include "MdcRawEvent/MdcDigi.h"
57#include "MdcRecEvent/RecMdcHit.h"
58#include "MdcRecEvent/RecMdcTrack.h"
59#include "MdcRecEvent/RecMdcKalHelixSeg.h"
60#include "MdcRecEvent/RecMdcKalTrack.h"
61#include "MdcGeomSvc/MdcGeomSvc.h"
62#include "MagneticField/IMagneticFieldSvc.h"
63#include "MagneticField/MagneticFieldSvc.h"
65#include "VertexFit/IVertexDbSvc.h"
67#include "Identifier/Identifier.h"
68#include "Identifier/MdcID.h"
69#include "GaudiKernel/IPartPropSvc.h"
70#include "GaudiKernel/INTupleSvc.h"
71using CLHEP::HepVector;
72using CLHEP::Hep3Vector;
73using CLHEP::HepMatrix;
74using CLHEP::HepSymMatrix;
82const double KalFitAlg::RIW = 6.35;
86 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc_(0),m_MFSvc_(0),
87 _wire(0), _layer(0), _superLayer(0),
88 pathl_(1), wsag_(4), back_(1), pT_(0.0), lead_(2), mhyp_(31),
89 pe_cut_(2.0), pmu_cut_(2.0), ppi_cut_(2.0), pk_cut_(2.0), pp_cut_(2.0),
90 muls_(1), loss_(1), enhance_(0),drifttime_choice_(0),choice_(0),
91 fac_h1_(1),fac_h2_(1),fac_h3_(1),fac_h4_(1),fac_h5_(1),
92 i_back_(-1), debug_kft_(0), debug_(0), ntuple_(0),eventno(-1),
93 Tds_back_no(0),m_nt1(0),m_nt2(0),m_nt3(0),m_nt4(0),m_nt5(0),
94 iqual_back_(1),tprop_(1),
95 dchi2cut_inner_(0),dchi2cut_outer_(0),
96 dchi2cut_mid1_(0),dchi2cut_mid2_(0),
97 dchi2cut_layid2_(0),dchi2cut_layid3_(0),m_usevtxdb(0),m_dangcut(10),m_dphicut(10),
106 declareProperty(
"gain1",
gain1_ = 1.);
107 declareProperty(
"gain2",
gain2_ = 1.);
108 declareProperty(
"gain3",
gain3_ = 1.);
109 declareProperty(
"gain4",
gain4_ = 1.);
110 declareProperty(
"gain5",
gain5_ = 1.);
111 declareProperty(
"fitnocut",
fitnocut_ = 5);
114 declareProperty(
"choice",
choice_ = 0);
115 declareProperty(
"numfcor",
numfcor_ = 0);
116 declareProperty(
"numf",
numf_ = 0);
117 declareProperty(
"steplev",
steplev_ = 2);
118 declareProperty(
"usage",
usage_=0);
119 declareProperty(
"i_front",
i_front_=1);
120 declareProperty(
"lead",
lead_=2);
121 declareProperty(
"muls",
muls_=1);
122 declareProperty(
"loss",
loss_=1);
123 declareProperty(
"matrixg",
matrixg_=100.0);
124 declareProperty(
"lr",
lr_=1);
126 declareProperty(
"debug",
debug_=0);
127 declareProperty(
"ntuple",
ntuple_=0);
129 declareProperty(
"matfile",
matfile_=
"geomdc_material.dat");
130 declareProperty(
"cylfile",
cylfile_=
"geomdc_cylinder.dat");
131 declareProperty(
"dchi2cutf",
dchi2cutf_=1000.0);
132 declareProperty(
"dchi2cuts",
dchi2cuts_=1000.0);
133 declareProperty(
"pt",
pT_=0.0);
134 declareProperty(
"pe_cut",
pe_cut_=2.0);
135 declareProperty(
"pmu_cut",
pmu_cut_=2.0);
136 declareProperty(
"ppi_cut",
ppi_cut_=2.0);
137 declareProperty(
"pk_cut",
pk_cut_=2.0);
138 declareProperty(
"pp_cut",
pp_cut_=2.0);
139 declareProperty(
"wsag",
wsag_=4);
140 declareProperty(
"back",
back_=1);
141 declareProperty(
"i_back",
i_back_=1);
142 declareProperty(
"tofflag",
tofflag_=1);
143 declareProperty(
"tof_hyp",
tof_hyp_=1);
145 declareProperty(
"fstrag",
fstrag_=0.9);
147 declareProperty(
"tprop",
tprop_=1);
148 declareProperty(
"pt_cut",
pt_cut_= 0.2);
151 declareProperty(
"cosmicflag",
m_csmflag= 0);
152 declareProperty(
"dangcut",
m_dangcut=10.);
153 declareProperty(
"dphicut",
m_dphicut=10.);
155 for(
int i=0; i<5; i++) nFailedTrks[i]=0;
161 MsgStream log(
msgSvc(), name());
162 log << MSG::INFO <<
" Start cleaning, delete Mdc geometry objects" << endreq;
164 log << MSG::INFO <<
"End cleaning " << endreq;
170 _BesKalmanFitWalls.clear();
171 _BesKalmanFitMaterials.clear();
172 if(_wire)free(_wire);
173 if(_layer)free(_layer);
174 if(_superLayer)free(_superLayer);
180 MsgStream log(
msgSvc(), name());
181 log << MSG::INFO <<
"in initize()"
182 <<
"KalFit> Initialization for current run " << endreq;
183 log << MSG::INFO <<
"Present Parameters: muls: " <<
muls_ <<
" loss: "
206 log << MSG::INFO <<
".....building Mdc " << endreq;
226 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
227 if(sc!=StatusCode::SUCCESS) {
228 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
238 std::cout<<
" initialize, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
246 log << MSG::INFO <<
" ntuple out, the option is "<<
ntuple_ <<endreq;
248 cout <<
"KalFitAlg> DEBUG open,Here is the important Parameters :\n";
249 cout <<
" Leading particule with mass hyp = " <<
lead_ << std::endl;
250 cout <<
" mhyp = " <<
mhyp_ << std::endl;
251 cout <<
"===== Effects taking into account : " << std::endl;
252 cout <<
" - multiple scattering = " <<
muls_ << std::endl;
253 cout <<
" - energy loss = " <<
loss_ << std::endl;
255 cout <<
" - straggling for the energy loss " << std::endl;
259 cout <<
" - non uniform magnetic field treatment "
261 cout <<
" - wire sag correction = " <<
wsag_ << std::endl;
265 <<
" hits included " << std::endl;
268 cout <<
" Backward filter is on with a pT cut value = " <<
pT_ << endl;
270 if(
debug_ == 4) cout <<
" pathl = " <<
pathl_ << std::endl;
273 cout <<
" Decision L/R from MdcRecHit " << std::endl;
282 IPartPropSvc* p_PartPropSvc;
283 static const bool CREATEIFNOTTHERE(
true);
284 StatusCode PartPropStatus = service(
"PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
285 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
286 log << MSG::WARNING <<
" Could not initialize Particle Properties Service" << endreq;
287 return StatusCode::SUCCESS;
289 m_particleTable = p_PartPropSvc->PDT();
291 return StatusCode::SUCCESS;
296 MsgStream log(
msgSvc(), name());
297 log << MSG::DEBUG<<
"KalFitAlg:: nTotalTrks = "<<nTotalTrks<<endreq;
298 log << MSG::DEBUG<<
" e: "<<nFailedTrks[0]<<
" failed, "<<nTotalTrks-nFailedTrks[0]<<
" successed"<<endreq;
299 log << MSG::DEBUG<<
" mu: "<<nFailedTrks[1]<<
" failed, "<<nTotalTrks-nFailedTrks[1]<<
" successed"<<endreq;
300 log << MSG::DEBUG<<
" pi: "<<nFailedTrks[2]<<
" failed, "<<nTotalTrks-nFailedTrks[2]<<
" successed"<<endreq;
301 log << MSG::DEBUG<<
" K: "<<nFailedTrks[3]<<
" failed, "<<nTotalTrks-nFailedTrks[3]<<
" successed"<<endreq;
302 log << MSG::DEBUG<<
" p: "<<nFailedTrks[4]<<
" failed, "<<nTotalTrks-nFailedTrks[4]<<
" successed"<<endreq;
303 return StatusCode::SUCCESS;
309 MsgStream log(
msgSvc(), name());
310 log << MSG::INFO <<
"in beginRun()" << endreq;
311 log << MSG::INFO <<
"Present Parameters: muls: " <<
muls_ <<
" loss: "
320 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
321 if(sc!=StatusCode::SUCCESS) {
322 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
331 std::cout<<
" beginRun, referField from MagneticFieldSvc:"<< (IMFSvc->
getReferField())*10000 <<std::endl;
336 return StatusCode::SUCCESS;
346 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101");
348 if ( nt1 ) m_nt1 = nt1;
350 m_nt1=
ntupleSvc()->book(
"FILE104/n101",CLID_ColumnWiseTuple,
"KalFit");
353 status = m_nt1->addItem(
"trackid",m_trackid);
354 status = m_nt1->addItem(
"stat",5,2,m_stat);
355 status = m_nt1->addItem(
"ndf",5,2,m_ndf);
356 status = m_nt1->addItem(
"chisq",5,2,m_chisq);
357 status = m_nt1->addItem(
"length",5,m_length);
358 status = m_nt1->addItem(
"tof",5,m_tof);
359 status = m_nt1->addItem(
"nhits",5,m_nhits);
360 status = m_nt1->addItem(
"zhelix",5,m_zhelix);
361 status = m_nt1->addItem(
"zhelixe",5,m_zhelixe);
362 status = m_nt1->addItem(
"zhelixmu",5,m_zhelixmu);
363 status = m_nt1->addItem(
"zhelixk",5,m_zhelixk);
364 status = m_nt1->addItem(
"zhelixp",5,m_zhelixp);
365 status = m_nt1->addItem(
"zptot",m_zptot);
366 status = m_nt1->addItem(
"zptote",m_zptote);
367 status = m_nt1->addItem(
"zptotmu",m_zptotmu);
368 status = m_nt1->addItem(
"zptotk",m_zptotk);
369 status = m_nt1->addItem(
"zptotp",m_zptotp);
371 status = m_nt1->addItem(
"zpt",m_zpt);
372 status = m_nt1->addItem(
"zpte",m_zpte);
373 status = m_nt1->addItem(
"zptmu",m_zptmu);
374 status = m_nt1->addItem(
"zptk",m_zptk);
375 status = m_nt1->addItem(
"zptp",m_zptp);
377 status = m_nt1->addItem(
"fptot",m_fptot);
378 status = m_nt1->addItem(
"fptote",m_fptote);
379 status = m_nt1->addItem(
"fptotmu",m_fptotmu);
380 status = m_nt1->addItem(
"fptotk",m_fptotk);
381 status = m_nt1->addItem(
"fptotp",m_fptotp);
382 status = m_nt1->addItem(
"fpt",m_fpt);
383 status = m_nt1->addItem(
"fpte",m_fpte);
384 status = m_nt1->addItem(
"fptmu",m_fptmu);
385 status = m_nt1->addItem(
"fptk",m_fptk);
386 status = m_nt1->addItem(
"fptp",m_fptp);
388 status = m_nt1->addItem(
"lptot",m_lptot);
389 status = m_nt1->addItem(
"lptote",m_lptote);
390 status = m_nt1->addItem(
"lptotmu",m_lptotmu);
391 status = m_nt1->addItem(
"lptotk",m_lptotk);
392 status = m_nt1->addItem(
"lptotp",m_lptotp);
393 status = m_nt1->addItem(
"lpt",m_lpt);
394 status = m_nt1->addItem(
"lpte",m_lpte);
395 status = m_nt1->addItem(
"lptmu",m_lptmu);
396 status = m_nt1->addItem(
"lptk",m_lptk);
397 status = m_nt1->addItem(
"lptp",m_lptp);
399 status = m_nt1->addItem(
"zsigp",m_zsigp);
400 status = m_nt1->addItem(
"zsigpe",m_zsigpe);
401 status = m_nt1->addItem(
"zsigpmu",m_zsigpmu);
402 status = m_nt1->addItem(
"zsigpk",m_zsigpk);
403 status = m_nt1->addItem(
"zsigpp",m_zsigpp);
404 status = m_nt1->addItem(
"fhelix",5,m_fhelix);
405 status = m_nt1->addItem(
"fhelixe",5,m_fhelixe);
406 status = m_nt1->addItem(
"fhelixmu",5,m_fhelixmu);
407 status = m_nt1->addItem(
"fhelixk",5,m_fhelixk);
408 status = m_nt1->addItem(
"fhelixp",5,m_fhelixp);
409 status = m_nt1->addItem(
"lhelix",5,m_lhelix);
410 status = m_nt1->addItem(
"lhelixe",5,m_lhelixe);
411 status = m_nt1->addItem(
"lhelixmu",5,m_lhelixmu);
412 status = m_nt1->addItem(
"lhelixk",5,m_lhelixk);
413 status = m_nt1->addItem(
"lhelixp",5,m_lhelixp);
415 status = m_nt1->addItem(
"zerror",15,m_zerror);
416 status = m_nt1->addItem(
"zerrore",15,m_zerrore);
417 status = m_nt1->addItem(
"zerrormu",15,m_zerrormu);
418 status = m_nt1->addItem(
"zerrork",15,m_zerrork);
419 status = m_nt1->addItem(
"zerrorp",15,m_zerrorp);
420 status = m_nt1->addItem(
"ferror",15,m_ferror);
421 status = m_nt1->addItem(
"ferrore",15,m_ferrore);
422 status = m_nt1->addItem(
"ferrormu",15,m_ferrormu);
423 status = m_nt1->addItem(
"ferrork",15,m_ferrork);
424 status = m_nt1->addItem(
"ferrorp",15,m_ferrorp);
425 status = m_nt1->addItem(
"lerror",15,m_lerror);
426 status = m_nt1->addItem(
"lerrore",15,m_lerrore);
427 status = m_nt1->addItem(
"lerrormu",15,m_lerrormu);
428 status = m_nt1->addItem(
"lerrork",15,m_lerrork);
429 status = m_nt1->addItem(
"lerrorp",15,m_lerrorp);
432 status = m_nt1->addItem(
"evtid",m_evtid);
433 status = m_nt1->addItem(
"mchelix",5,m_mchelix);
434 status = m_nt1->addItem(
"mcptot",m_mcptot);
435 status = m_nt1->addItem(
"mcpid",m_mcpid);
437 if( status.isFailure() ) cout<<
"Ntuple1 add item failed!"<<endl;
443 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102");
445 if ( nt2 ) m_nt2 = nt2;
447 m_nt2=
ntupleSvc()->book(
"FILE104/n102",CLID_ColumnWiseTuple,
"KalFitComp");
449 status2 = m_nt2->addItem(
"delx",m_delx);
450 status2 = m_nt2->addItem(
"dely",m_dely);
451 status2 = m_nt2->addItem(
"delz",m_delz);
452 status2 = m_nt2->addItem(
"delthe",m_delthe);
453 status2 = m_nt2->addItem(
"delphi",m_delphi);
454 status2 = m_nt2->addItem(
"delp",m_delp);
455 status2 = m_nt2->addItem(
"delpx",m_delpx);
456 status2 = m_nt2->addItem(
"delpy",m_delpy);
457 status2 = m_nt2->addItem(
"delpz",m_delpz);
459 if( status2.isFailure() ) cout<<
"Ntuple2 add item failed!"<<endl;
465 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103");
467 if ( nt3 ) m_nt3 = nt3;
469 m_nt3=
ntupleSvc()->book(
"FILE104/n103",CLID_ColumnWiseTuple,
"PatRec");
471 status3 = m_nt3->addItem(
"trkhelix",5,m_trkhelix);
472 status3 = m_nt3->addItem(
"trkptot",m_trkptot);
474 status3 = m_nt3->addItem(
"trkerror",15,m_trkerror);
475 status3 = m_nt3->addItem(
"trksigp",m_trksigp);
477 status3 = m_nt3->addItem(
"trkndf",m_trkndf);
478 status3 = m_nt3->addItem(
"trkchisq",m_trkchisq);
479 if( status3.isFailure() ) cout<<
"Ntuple3 add item failed!"<<endl;
485 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104");
487 if ( nt4 ) m_nt4 = nt4;
489 m_nt4=
ntupleSvc()->book(
"FILE104/n104",CLID_ColumnWiseTuple,
"PatRecComp");
491 status4 = m_nt4->addItem(
"trkdelx",m_trkdelx);
492 status4 = m_nt4->addItem(
"trkdely",m_trkdely);
493 status4 = m_nt4->addItem(
"trkdelz",m_trkdelz);
494 status4 = m_nt4->addItem(
"trkdelthe",m_trkdelthe);
495 status4 = m_nt4->addItem(
"trkdelphi",m_trkdelphi);
496 status4 = m_nt4->addItem(
"trkdelp",m_trkdelp);
497 if( status4.isFailure() ) cout<<
"Ntuple4 add item failed!"<<endl;
502 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105");
504 if ( nt5 ) m_nt5 = nt5;
506 m_nt5=
ntupleSvc()->book(
"FILE104/n105",CLID_ColumnWiseTuple,
"KalFitdChisq");
508 status5 = m_nt5->addItem(
"dchi2",m_dchi2);
509 status5 = m_nt5->addItem(
"masshyp",m_masshyp);
510 status5 = m_nt5->addItem(
"residual_estim",m_residest);
511 status5 = m_nt5->addItem(
"residual",m_residnew);
512 status5 = m_nt5->addItem(
"layer",m_layer);
513 status5 = m_nt5->addItem(
"kaldr",m_anal_dr);
514 status5 = m_nt5->addItem(
"kalphi0",m_anal_phi0);
515 status5 = m_nt5->addItem(
"kalkappa",m_anal_kappa);
516 status5 = m_nt5->addItem(
"kaldz",m_anal_dz);
517 status5 = m_nt5->addItem(
"kaltanl",m_anal_tanl);
518 status5 = m_nt5->addItem(
"dr_ea",m_anal_ea_dr);
519 status5 = m_nt5->addItem(
"phi0_ea",m_anal_ea_phi0);
520 status5 = m_nt5->addItem(
"kappa_ea",m_anal_ea_kappa);
521 status5 = m_nt5->addItem(
"dz_ea",m_anal_ea_dz);
522 status5 = m_nt5->addItem(
"tanl_ea",m_anal_ea_tanl);
523 if( status5.isFailure() ) cout<<
"Ntuple5 add item failed!"<<endl;
526 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
528 if ( nt6 ) m_nt6 = nt6;
530 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
532 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
533 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
534 status6 = m_nt6->addItem(
"tdr",m_tdrift);
535 status6 = m_nt6->addItem(
"layerid", m_layerid);
536 status6 = m_nt6->addItem(
"event", m_eventNo);
537 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
538 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
539 status6 = m_nt6->addItem(
"lr",m_lr);
540 status6 = m_nt6->addItem(
"dd",m_dd);
541 status6 = m_nt6->addItem(
"y",m_yposition);
543 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
556 for(layid = 0; layid<2; layid++) {
561 for(layid = 4; layid<12; layid++) {
564 for(layid = 12; layid<20; layid++) {
567 for(layid = 20; layid<43; layid++) {
573 for(layid = 0; layid<2; layid++) {
579 for(layid = 4; layid<12; layid++) {
582 for(layid = 12; layid<20; layid++) {
585 for(layid = 20; layid<43; layid++) {
596 for(layid = 0; layid<2; layid++) {
603 for(layid = 4; layid<12; layid++) {
607 for(layid = 12; layid<20; layid++) {
611 for(layid = 20; layid<43; layid++) {
617 for(layid = 0; layid<43; layid++) {
624 for(layid = 0; layid<2; layid++) {
631 for(layid = 4; layid<12; layid++) {
635 for(layid = 12; layid<20; layid++) {
639 for(layid = 20; layid<43; layid++) {
645 for(layid = 0; layid<43; layid++) {
656 MsgStream log(
msgSvc(), name());
657 log << MSG::INFO <<
"in execute()" << endreq;
704 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
705 if(sc!=StatusCode::SUCCESS) {
706 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
715 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
722 IDataProviderSvc* evtSvc = NULL;
723 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
725 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
727 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
728 return StatusCode::SUCCESS;
732 IDataManagerSvc *dataManSvc;
733 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
734 DataObject *aKalTrackCol;
735 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
736 if(aKalTrackCol != NULL) {
737 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
738 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
741 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
742 if( kalsc.isFailure() ) {
743 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
744 return StatusCode::SUCCESS;
746 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
748 DataObject *aKalHelixSegCol;
749 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
750 if(aKalHelixSegCol != NULL){
751 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
752 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
755 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
756 if( kalsc.isFailure()){
757 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
758 return StatusCode::SUCCESS;
760 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
772 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
775 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
777 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
778 return StatusCode::FAILURE;
780 int eventNo = eventHeader->eventNumber();
781 int runNo = eventHeader->runNumber();
786 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
787 if (estimeCol && estimeCol->size()) {
788 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
789 t0 = (*iter_evt)->getTest();
792 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
793 return StatusCode::SUCCESS;
798 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
802 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
803 if (sc!=StatusCode::SUCCESS) {
804 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
805 return StatusCode::SUCCESS;
816 m_evtid = eventHeader->eventNumber();
819 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
821 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
826 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
827 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
828 if(!(*i_mcTrk)->primaryParticle())
continue;
829 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
830 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
831 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
832 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
833 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
836 int pid = (*i_mcTrk)->particleProperty();
838 charge = m_particleTable->particle( pid )->charge();
839 }
else if ( pid <0 ) {
840 charge = m_particleTable->particle( -pid )->charge();
843 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
846 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
849 log << MSG::DEBUG <<
"charge of the track " <<
charge << endreq;
850 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
852 for(
int j =0; j<5; j++) {
853 m_mchelix[j] = mchelix.
a()[j];
856 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
864 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
866 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
867 return( StatusCode::SUCCESS);
869 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
874 mtrkadd_mgr->clear();
878 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
882 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
883 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
885 csmp3[kj-1]=(*iter_trk)->p3();
886 csmphi[kj-1] = (*iter_trk)->phi();
890 for(
int j = 0, ij = 0; j<5; j++) {
891 m_trkhelix[j] = (*iter_trk)->helix()[j];
893 for(
int k=0; k<=j; k++,ij++) {
894 m_trkerror[ij] = (*iter_trk)->err()[j][k];
898 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
900 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
901 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
902 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
904 m_trkndf = (*iter_trk)->ndof();
905 m_trkchisq = (*iter_trk)->chi2();
907 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
909 StatusCode sc3 = m_nt3->write();
910 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
940 log << MSG::DEBUG <<
"retrieved MDC tracks:"
941 <<
" Nhits " <<(*iter_trk)->getNhits()
942 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
944 HitRefVec gothits = (*iter_trk)->getVecHits();
948 rectrk->
id = (*iter_trk)->trackId();
949 rectrk->
chiSq = (*iter_trk)->chi2();
950 rectrk->
ndf = (*iter_trk)->ndof();
951 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
952 rectrk->
nhits = (*iter_trk)->getNhits();
953 rectrk->
nster = (*iter_trk)->nster();
955 rectrk->
stat = (*iter_trk)->stat();
956 status_temp = (*iter_trk)->stat();
958 trkadd->
id = (*iter_trk)->trackId();
962 trkadd->
body = rectrk;
963 rectrk->
add = trkadd;
965 for (
int i=0; i<5; i++) {
966 rectrk->
helix[i] = (*iter_trk)->helix()[i];
967 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
968 for(
int j = 1; j<i+2;j++) {
969 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
972 std::sort(gothits.begin(), gothits.end(), order_rechits);
973 HitRefVec::iterator it_gothit = gothits.begin();
974 for( ; it_gothit != gothits.end(); it_gothit++) {
976 if( (*it_gothit)->getStat() != 1 ) {
978 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
983 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
984 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
985 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
986 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
987 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
988 <<
" id of hit "<<(*it_gothit)->getId()
989 <<
" track id of hit "<<(*it_gothit)->getTrkId()
990 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
993 whit->
id = (*it_gothit)->getId();
994 whit->
ddl = (*it_gothit)->getDriftDistLeft();
995 whit->
ddr = (*it_gothit)->getDriftDistRight();
996 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
997 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
998 whit->
pChiSq = (*it_gothit)->getChisqAdd();
999 whit->
lr = (*it_gothit)->getFlagLR();
1000 whit->
stat = (*it_gothit)->getStat();
1001 mdcid = (*it_gothit)->getMdcId();
1005 int wid = w0id + localwid;
1007 <<
"lr from PR: "<<whit->
lr
1008 <<
" layerId = " << layid
1009 <<
" wireId = " << localwid
1019 whit->
tdc = (*it_gothit)->getTdc();
1020 whit->
adc= (*it_gothit)->getAdc();
1021 rectrk->
hitcol.push_back(whit);
1022 mhit_mgr->push_back(*whit);
1024 mtrk_mgr->push_back(*rectrk);
1025 mtrkadd_mgr->push_back(*trkadd);
1033 m_trkdelx = trkx1 - trkx2;
1034 m_trkdely = trky1 - trky2;
1035 m_trkdelz = trkz1 - trkz2;
1036 m_trkdelthe = trkthe1 + trkthe2;
1037 m_trkdelphi = trkphi1- trkphi2;
1038 m_trkdelp = trkp1 - trkp2;
1039 StatusCode sc4 = m_nt4->write();
1040 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1043 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1044 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1050 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1051 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1056 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1063 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1068 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1070 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++){
1072 for(
int hypo=0; hypo<5; hypo++)
1074 if((*KalTrk)->getStat(0,hypo)==1) nFailedTrks[hypo]++;
1077 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1079 int nhitofthistrk=0;
1080 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1087 iter_hit=gothelixsegs.begin();
1135 return StatusCode::SUCCESS;
1146 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1147 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1149 TrasanTRK.
helix[2]*track.
a()[2]<0)
1153 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1154 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1155 <<
" "<<TrasanTRK.
helix[2]<<
" "<<track.
a()[2]<<
" "<<iqual<<endl;
1156 cout<<
"FillTds> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1157 <<
" nster "<<track.
nster()<<
" iqual "<<iqual<<
" track.Ea "<< track.
Ea()<<endl;
1159 cout<<
"fillTds>.....track.Ea[2][2] "<<track.
Ea()[2][2]<<endl;
1160 cout <<
" TRASAN stereo = " << trasster
1161 <<
" and KalFitTrack = " << trakster << std::endl;
1162 cout <<
" TRASAN axial = " << trasax
1163 <<
" and KalFitTrack = " << trakax << std::endl;
1166 cout <<
"...there is a problem during fit !! " << std::endl;
1167 if (trasster-trakster>5)
1168 cout <<
" because stereo " << trasster-trakster << std::endl;
1169 if (trasax-trakax >5)
1170 cout <<
" because axial " << std::endl;
1171 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1172 cout <<
" because kappa sign " << std::endl;
1178 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1179 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1180 track.
Ea()[4][4] > 0 && iqual) {
1181 if(
debug_ == 4) cout<<
"fillTds>.....going on "<<endl;
1195 if(
debug_) cout<<
"ALARM: FillTds Not refit with KalFilter!!!"<<endl;
1204 double a_trasan[5], ea_trasan[15];
1205 for(
int i =0 ; i <5; i++){
1206 a_trasan[i] = TrasanTRK.
helix[i];
1208 for(
int j =0 ; j <15; j++){
1209 ea_trasan[j] = TrasanTRK.
error[j];
1224 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1225 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1227 TrasanTRK.
helix[2]*track.
a()[2]<0)
1231 cout<<
"Nhit from PR "<<TrasanTRK.
nhits<<
" nhit "<<track.
nchits()<<endl;
1232 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1233 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1235 cout<<
"FillTds_lead> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1236 <<
" nster "<<track.
nster()<<
" iqual_front_[l_mass] "<<
iqual_front_[l_mass]<<
" track.Ea "<<track.
Ea()<<endl;
1238 cout <<
" TRASAN stereo = " << trasster
1239 <<
" and KalFitTrack = " << trakster << std::endl;
1240 cout <<
" TRASAN axial = " << trasax
1241 <<
" and KalFitTrack = " << trakax << std::endl;
1244 cout <<
"...there is a problem during fit !! " << std::endl;
1245 if (trasster-trakster>5)
1246 cout <<
" because stereo " << trasster-trakster << std::endl;
1247 if (trasax-trakax >5)
1248 cout <<
" because axial " << std::endl;
1249 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1250 cout <<
" because kappa sign " << std::endl;
1256 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1257 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1268 if (
debug_ == 4) cout<<
" trasan id...1 "<<TrasanTRK.
id<<endl;
1277 if(
debug_) cout<<
"ALARM: FillTds_forMdc Not refit with KalFilter!!!"<<endl;
1288 if (
debug_ ==4) cout<<
" trasan id...2 "<<TrasanTRK.
id<<endl;
1292 double a_trasan[5], ea_trasan[15];
1293 for(
int i =0 ; i <5; i++){
1294 a_trasan[i] = TrasanTRK.
helix[i];
1296 for(
int j =0 ; j <15; j++){
1297 ea_trasan[j] = TrasanTRK.
error[j];
1316 cout <<
"fillTds_IP>......"<<endl;
1317 cout <<
" dr = " << track.
a()[0]
1318 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1319 cout <<
" phi0 = " << track.
a()[1]
1320 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1321 cout <<
" PT = " << 1/track.
a()[2]
1322 <<
", Er_kappa =" << sqrt(track.
Ea()[2][2]) << std::endl;
1323 cout <<
" dz = " << track.
a()[3]
1324 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1325 cout <<
" tanl = " << track.
a()[4]
1326 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1330 TrasanTRK.
helix[2]*track.
a()[2]<0)
1336 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1337 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1341 double dr = track.
a()[0];
1342 double phi0 = track.
a()[1];
1343 double kappa = track.
a()[2];
1344 double dz = track.
a()[3];
1345 double tanl = track.
a()[4];
1350 double vx = dr*
cos(phi0);
1351 double vy = dr*
sin(phi0);
1356 if(0==kappa) kappa = 10e-10;
1357 double px = -
sin(phi0)/fabs(kappa);
1358 double py =
cos(phi0)/fabs(kappa);
1359 double pz = tanl/fabs(kappa);
1361 trk->
setX(vx, l_mass);
1362 trk->
setY(vy, l_mass);
1363 trk->
setZ(vz, l_mass);
1364 trk->
setPx(px, l_mass);
1365 trk->
setPy(py, l_mass);
1366 trk->
setPz(pz, l_mass);
1377 if (kappa > 0.0000000001)
1379 else if (kappa < -0.0000000001)
1384 double ptot = sqrt(px*px+py*py+pz*pz);
1385 trk->
setTheta(acos(pz/ptot),l_mass);
1392 double dr = TrasanTRK.
helix[0];
1393 double phi0 = TrasanTRK.
helix[1];
1394 double kappa = TrasanTRK.
helix[2];
1395 double dz = TrasanTRK.
helix[3];
1396 double tanl = TrasanTRK.
helix[4];
1398 double vx = dr*
cos(phi0);
1399 double vy = dr*
sin(phi0);
1402 if(0==kappa) kappa = 10e-10;
1403 double px = -
sin(phi0)/fabs(kappa);
1404 double py =
cos(phi0)/fabs(kappa);
1405 double pz = tanl/fabs(kappa);
1407 trk->
setX(vx, l_mass);
1408 trk->
setY(vy, l_mass);
1409 trk->
setZ(vz, l_mass);
1411 trk->
setPx(px, l_mass);
1412 trk->
setPy(py, l_mass);
1413 trk->
setPz(pz, l_mass);
1420 double a_trasan[5], ea_trasan[15];
1421 for(
int i =0 ; i <5; i++){
1422 a_trasan[i] = TrasanTRK.
helix[i];
1424 for(
int j =0 ; j <15; j++){
1425 ea_trasan[j] = TrasanTRK.
error[j];
1432 if (kappa > 0.0000000001)
1434 else if (kappa < -0.0000000001)
1439 double ptot = sqrt(px*px+py*py+pz*pz);
1440 trk->
setTheta(acos(pz/ptot),l_mass);
1444 SmartDataPtr<RecMdcTrackCol> mdcTrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
1447 RecMdcTrackCol::iterator iter_mdcTrk = mdcTrkCol->begin();
1448 bool findMdcTrk =
false;
1449 for(; iter_mdcTrk != mdcTrkCol->end(); iter_mdcTrk++) {
1450 if(TrasanTRK.
id==(*iter_mdcTrk)->trackId()) {
1455 int nLayer = (*iter_mdcTrk)->nlayer();
1462 std::cout<<
"px: "<<trk->
px()<<
" py: "<<trk->
py()<<
" pz: "<<trk->
pz()<<std::endl;
1463 std::cout<<
"vx: "<<trk->
x()<<
" vy: "<<trk->
y()<<
" vz: "<<trk->
z()<<std::endl;
1480 if(
debug_ == 4) cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0] "<<trk->
getNdf(0,2)<<endl;
1484 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1485 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1486 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
1497 if(
debug_ == 4) cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
1507 for (
int i = 0; i<43; i++) {
1521 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
1522 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1523 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1524 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1528 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
1533 TrasanTRK.
pivot[2]);
1536 for(
int i = 0; i < 5; i++)
1537 a[i] = TrasanTRK.
helix[i];
1540 for(
int i = 0, k = 0; i < 5; i++) {
1541 for(
int j = 0; j <= i; j++) {
1543 ea[j][i] = ea[i][j];
1549 double fiTerm = TrasanTRK.
fiTerm;
1551 double fi0 = track_rep.
phi0();
1557 if( fabs(
x ) > 1.0e-10 ){
1558 phi_x = atan2(
y,
x );
1559 if( phi_x < 0 ) phi_x += 2*
M_PI;
1561 phi_x = (
y > 0 ) ? M_PI_4: 3.0*M_PI_4;
1563 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
1564 double dphi = fabs( fiTerm + fi0 - phi_x );
1565 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
1566 double tanl = track_rep.
tanl();
1567 double cosl_inv = sqrt( tanl*tanl + 1.0 );
1569 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
1570 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
1572 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
1573 double light_speed( 29.9792458 );
1574 double pt( 1.0 / track_rep.
kappa() );
1575 double p( pt * sqrt( 1.0 + tanl*tanl ) );
1581 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<endl;
1582 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<endl;
1587 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
1588 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
1590 track_rep.
pivot(IP);
1606 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1607 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1608 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1610 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
1611 cout <<
" dr = " << track.
a()[0]
1612 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1613 cout<<
" phi0 = " << track.
a()[1]
1614 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1615 cout <<
" PT = " << 1/track.
a()[2]
1616 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
1617 cout <<
" dz = " << track.
a()[3]
1618 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1619 cout <<
" tanl = " << track.
a()[4]
1620 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1621 cout <<
" Ea = " << track.
Ea() <<endl;
1648 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
1649 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
1650 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;
1654 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1655 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1656 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
1665 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
1676 helixseg->
setResIncl(it->residual_include());
1677 helixseg->
setResExcl(it->residual_exclude());
1679 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
1681 helixseg->
setDrIncl(it->a_include()[0]);
1684 helixseg->
setDzIncl(it->a_include()[3]);
1688 helixseg->
setDrExcl(it->a_exclude()[0]);
1691 helixseg->
setDzExcl(it->a_exclude()[3]);
1706 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
1707 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
1708 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
1709 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
1711 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
1712 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
1713 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
1715 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
1719 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
1720 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
1721 helixseg->
setFlagLR(it->HitMdc()->LR());
1722 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
1723 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
1724 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
1725 helixseg->
setTof(it->tof());
1728 helixseg->
setDD(it->dd());
1729 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
1730 helixseg->
setDT(it->dt());
1731 segcol->push_back(helixseg);
1732 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
1733 helixsegrefvec.push_back(refhelixseg);
1736 m_docaInc = helixseg -> getDocaIncl();
1737 m_docaExc = helixseg -> getDocaExcl();
1738 m_residualInc = helixseg -> getResIncl();
1739 m_residualExc = helixseg -> getResExcl();
1740 m_dd = helixseg -> getDD();
1742 m_tdrift = helixseg -> getDT();
1743 m_layerid = helixseg -> getLayerId();
1744 m_yposition= it->HitMdc()->wire().fwd().y();
1746 StatusCode sc6 = m_nt6->write();
1747 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
1754 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
1762 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
1766 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
1775 for (
int i = 0; i<43; i++) {
1786 double a_trasan[5], ea_trasan[15];
1787 for(
int i =0 ; i <5; i++){
1788 a_trasan[i] = TrasanTRK.
helix[i];
1790 for(
int j =0 ; j <15; j++){
1791 ea_trasan[j] = TrasanTRK.
helix[j];
1797 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
1798 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1799 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1800 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1805 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
1811 TrasanTRK.
pivot[2]);
1814 for(
int i = 0; i < 5; i++)
1815 a[i] = TrasanTRK.
helix[i];
1818 for(
int i = 0, k = 0; i < 5; i++) {
1819 for(
int j = 0; j <= i; j++) {
1821 ea[j][i] = ea[i][j];
1827 double fiTerm = TrasanTRK.
fiTerm;
1829 double fi0 = track_rep.
phi0();
1835 if( fabs(
x ) > 1.0e-10 ){
1836 phi_x = atan2(
y,
x );
1837 if( phi_x < 0 ) phi_x += 2*
M_PI;
1839 phi_x = (
y > 0 ) ? M_PI_4: 3.0*M_PI_4;
1841 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
1842 double dphi = fabs( fiTerm + fi0 - phi_x );
1843 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
1844 double tanl = track_rep.
tanl();
1845 double cosl_inv = sqrt( tanl*tanl + 1.0 );
1847 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
1848 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
1850 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
1851 double light_speed( 29.9792458 );
1852 double pt( 1.0 / track_rep.
kappa() );
1853 double p( pt * sqrt( 1.0 + tanl*tanl ) );
1861 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
1862 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
1867 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
1868 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
1870 track_rep.
pivot(IP);
1888 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
1889 cout <<
" dr = " << track.
a()[0]
1890 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1891 cout<<
" phi0 = " << track.
a()[1]
1892 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1893 cout <<
" PT = " << 1/track.
a()[2]
1894 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
1895 cout <<
" dz = " << track.
a()[3]
1896 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1897 cout <<
" tanl = " << track.
a()[4]
1898 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1899 cout <<
" Ea = " << track.
Ea() <<endl;
1924 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
1925 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
1930 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1931 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1932 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
1941 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
1952 helixseg->
setResIncl(it->residual_include());
1953 helixseg->
setResExcl(it->residual_exclude());
1955 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
1982 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
1983 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
1984 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
1985 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
1987 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
1988 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
1989 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
1991 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
1995 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
1996 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
1997 helixseg->
setFlagLR(it->HitMdc()->LR());
1998 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
1999 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
2000 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
2001 helixseg->
setTof(it->tof());
2004 helixseg->
setDD(it->dd());
2005 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
2006 helixseg->
setDT(it->dt());
2008 segcol->push_back(helixseg);
2009 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
2010 helixsegrefvec.push_back(refhelixseg);
2012 m_docaInc = helixseg -> getDocaIncl();
2013 m_docaExc = helixseg -> getDocaExcl();
2014 m_residualInc = helixseg -> getResIncl();
2015 m_residualExc = helixseg -> getResExcl();
2016 m_dd = helixseg -> getDD();
2018 m_tdrift = helixseg -> getDT();
2019 m_layerid = helixseg -> getLayerId();
2020 m_yposition= it->HitMdc()->wire().fwd().y();
2022 StatusCode sc6 = m_nt6->write();
2023 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
2031 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
2039 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
2043 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
2052 for (
int i = 0; i<43; i++) {
2063 double a_trasan[5], ea_trasan[15];
2064 for(
int i =0 ; i <5; i++){
2065 a_trasan[i] = TrasanTRK.
helix[i];
2067 for(
int j =0 ; j <15; j++){
2068 ea_trasan[j] = TrasanTRK.
helix[j];
2074 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
2075 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2076 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2077 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2082 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
2088 TrasanTRK.
pivot[2]);
2091 for(
int i = 0; i < 5; i++)
2092 a[i] = TrasanTRK.
helix[i];
2095 for(
int i = 0, k = 0; i < 5; i++) {
2096 for(
int j = 0; j <= i; j++) {
2098 ea[j][i] = ea[i][j];
2104 double fiTerm = TrasanTRK.
fiTerm;
2106 double fi0 = track_rep.
phi0();
2112 if( fabs(
x ) > 1.0e-10 ){
2113 phi_x = atan2(
y,
x );
2114 if( phi_x < 0 ) phi_x += 2*
M_PI;
2116 phi_x = (
y > 0 ) ? M_PI_4: 3.0*M_PI_4;
2118 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
2119 double dphi = fabs( fiTerm + fi0 - phi_x );
2120 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
2121 double tanl = track_rep.
tanl();
2122 double cosl_inv = sqrt( tanl*tanl + 1.0 );
2124 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
2125 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
2128 double track_len(fabs( track_rep.
radius() * fiTerm * cosl_inv ));
2130 double light_speed( 29.9792458 );
2131 double pt( 1.0 / track_rep.
kappa() );
2132 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2140 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
2141 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
2146 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2147 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2151 track_rep.
pivot(LPiovt);
2238 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2239 cout <<
" dr = " << track.
a()[0]
2240 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2241 cout<<
" phi0 = " << track.
a()[1]
2242 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2243 cout <<
" PT = " << 1/track.
a()[2]
2244 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2245 cout <<
" dz = " << track.
a()[3]
2246 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2247 cout <<
" tanl = " << track.
a()[4]
2248 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2249 cout <<
" Ea = " << track.
Ea() <<endl;
2265 for(
int jj = 0; jj<2; jj++) {
2279 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2280 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2283 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2288 Hep3Vector xorigin(0,0,0);
2290 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
2294 xorigin.setX(dbv[0]);
2295 xorigin.setY(dbv[1]);
2296 xorigin.setZ(dbv[2]);
2310 double tanl = track.
tanl();
2311 double phi_old = work.
phi0();
2312 double phi = track.
phi0();
2314 if (fabs(phi - phi_old) >
M_PI) {
2315 if (phi > phi_old) phi -= 2 *
M_PI;
2316 else phi_old -= 2 *
M_PI;
2319 double path_zero = fabs(track.
radius() * (phi_old-phi)* sqrt(1 + tanl * tanl));
2323 HepSymMatrix Eakal(5,0);
2331 unsigned int nhit = track.
HitsMdc().size();
2332 int layer_prev = -1;
2334 HepVector pos_old(3,0);
2335 double r0kal_prec(0);
2337 for(
unsigned i=0 ; i < nhit; i++ ) {
2338 int ihit = (nhit-1)-i;
2342 int wireid = Wire.
geoID();
2346 if (
pathl_ && layer != layer_prev) {
2348 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
2349 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
2357 Hep3Vector wire = (CLHEP::Hep3Vector)fwd - (CLHEP::Hep3Vector)bck;
2360 work.
pivot((fwd + bck) * .5);
2361 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
2363 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2367 double tension = geowire->
Tension();
2370 double zinit(x0kal.z()), lzx(Wire.
lzx());
2373 double A = 47.35E-6/tension;
2374 double Zp = (zinit - bck.z())*lzx/wire.z();
2377 std::cout<<
" sag in smoother_anal: "<<std::endl;
2378 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2379 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2380 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2381 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2382 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2383 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2386 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2387 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2388 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2390 wire.setX(wire.x()/wire.z());
2391 wire.setY(result.z());
2394 x0kal.setX(result.x());
2395 x0kal.setY(result.y());
2398 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2401 double r0kal = x0kal.perp();
2403 cout<<
"wire direction "<<wire<<endl;
2404 cout<<
"x0kal "<<x0kal<<endl;
2405 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
2422 double pmag( sqrt( 1.0 + track.
a()[4]*track.
a()[4]) / track.
a()[2]);
2423 double mass_over_p( track.
mass()/ pmag );
2424 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2425 double tofest = pathl / ( 29.9792458 * beta );
2435 if(!(way<0&&fabs(track.
kappa())>1000.0)) {
2445 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2447 HepSymMatrix Ma(5,0);
2450 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2464 point.setX(x0kal.x() + track.
a()[0]*
cos(track.
a()[1]));
2465 point.setY(x0kal.y() + track.
a()[0]*
sin(track.
a()[1]));
2466 point.setZ(x0kal.z() + track.
a()[3]);
2470 double phi_old = track.
a()[1];
2473 double phi_new = temp.
a()[1];
2474 double fi = phi_new - phi_old;
2478 if(fabs(fi) >= CLHEP::twopi) fi = fmod(fi+2*CLHEP::twopi,CLHEP::twopi);
2483 if (
debug_) cout<<
"track----7-----"<<track.
a()<<endl;
2495 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2496 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2499 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2502 HepSymMatrix Eakal(5,0);
2508 std::cout<<
"the initial error matrix in smoothing is "<<Eakal<<std::endl;
2513 unsigned int nseg = track.
HelixSegs().size();
2514 int layer_prev = -1;
2516 HepVector pos_old(3,0);
2517 double r0kal_prec(0);
2519 for(
unsigned i=0 ; i < nseg; i++ ) {
2522 int iseg = (nseg-1)-i;
2526 int wireid = Wire.
geoID();
2530 if (
pathl_ && layer != layer_prev) {
2532 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
2533 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
2541 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
2544 work.
pivot((fwd + bck) * .5);
2545 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
2548 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2554 double tension = geowire->
Tension();
2557 double zinit(x0kal.z()), lzx(Wire.
lzx());
2561 double A = 47.35E-6/tension;
2562 double Zp = (zinit - bck.z())*lzx/wire.z();
2566 std::cout<<
" sag in smoother_calib: "<<std::endl;
2567 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2568 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2569 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2570 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2571 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2572 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2575 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2576 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2577 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2579 wire.setX(wire.x()/wire.z());
2580 wire.setY(result.z());
2583 x0kal.setX(result.x());
2584 x0kal.setY(result.y());
2587 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2591 double r0kal = x0kal.perp();
2593 cout<<
"wire direction "<<wire<<endl;
2594 cout<<
"x0kal "<<x0kal<<endl;
2595 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
2602 if (
debug_ == 4) cout<<
"track----6-----"<<track.
a()<<
" ...path..."<<pathl
2603 <<
"momentum"<<track.
momentum(0)<<endl;
2608 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2610 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2617 if (
debug_) cout<<
"track----7-----"<<track.
a()<<endl;
2636 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2637 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2640 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2643 Hep3Vector x0inner(track.
pivot());
2644 HepVector pos_old(3,0);
2645 double r0kal_prec(0);
2647 int nhit = track.
HitsMdc().size();
2648 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
2649 for(
int i=0 ; i < nhit; i++ ) {
2652 if (HitMdc.
chi2()<0)
continue;
2658 int wireid = Wire.
geoID();
2662 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
2665 work.
pivot((fwd + bck) * .5);
2672 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
2674 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2718 double tension = geowire->
Tension();
2720 double zinit(x0kal.z()), lzx(Wire.
lzx());
2722 double A = 47.35E-6/tension;
2723 double Zp = (zinit - bck.z())*lzx/wire.z();
2726 std::cout<<
" sag in filter_fwd_anal: "<<std::endl;
2727 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2728 std::cout<<
"zinit: "<<zinit<<
" bck.z(): "<<bck.z()<<std::endl;
2729 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2730 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2731 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2732 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2733 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2736 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2737 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2738 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2740 wire.setX(wire.x()/wire.z());
2741 wire.setY(result.z());
2743 x0kal.setX(result.x());
2744 x0kal.setY(result.y());
2747 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2750 double r0kal = x0kal.perp();
2757 if (nhits_read == 1) {
2764 double dtracknew = 0.;
2768 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2769 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2770 double diff_chi2 = track.
chiSq();
2771 Hep3Vector IP(0,0,0);
2777 for(
unsigned k=i+1 ; k < nhit; k++ )
2782 double dchi2 = -1.0;
2833 std::cout<<
" ... ERROR OF dchi2... "<<std::endl;
2839 m_residest = dtrack-dtdc;
2840 m_residnew = dtracknew -dtdc;
2843 m_anal_dr = worktemp.
a()[0];
2844 m_anal_phi0 = worktemp.
a()[1];
2845 m_anal_kappa = worktemp.
a()[2];
2846 m_anal_dz = worktemp.
a()[3];
2847 m_anal_tanl = worktemp.
a()[4];
2848 m_anal_ea_dr = worktemp.
Ea()[0][0];
2849 m_anal_ea_phi0 = worktemp.
Ea()[1][1];
2850 m_anal_ea_kappa = worktemp.
Ea()[2][2];
2851 m_anal_ea_dz = worktemp.
Ea()[3][3];
2852 m_anal_ea_tanl = worktemp.
Ea()[4][4];
2853 StatusCode sc5 = m_nt5->write();
2854 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
2859 diff_chi2 = chi2 - diff_chi2;
2860 HitMdc.
chi2(diff_chi2);
2873 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2874 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2877 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2881 std::cout<<
"filter_fwd_calib starting ...the inital error Matirx is "<<track.
Ea()<<std::endl;
2883 Hep3Vector x0inner(track.
pivot());
2884 HepVector pos_old(3,0);
2885 double r0kal_prec(0);
2888 unsigned int nhit = track.
HitsMdc().size();
2889 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
2890 for(
unsigned i=0 ; i < nhit; i++ ) {
2894 cout<<
"filter_fwd...........222 chi2="<<HitMdc.
chi2()<<endl;
2896 if (HitMdc.
chi2()<0)
continue;
2899 int wireid = Wire.
geoID();
2905 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
2908 work.
pivot((fwd + bck) * .5);
2909 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
2912 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2917 double tension = geowire->
Tension();
2919 double zinit(x0kal.z()), lzx(Wire.
lzx());
2921 double A = 47.35E-6/tension;
2922 double Zp = (zinit - bck.z())*lzx/wire.z();
2926 std::cout<<
" sag in filter_fwd_calib: "<<std::endl;
2927 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2928 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2929 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2930 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2931 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2932 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2935 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2936 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2937 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2938 wire.setX(wire.x()/wire.z());
2939 wire.setY(result.z());
2941 x0kal.setX(result.x());
2942 x0kal.setY(result.y());
2946 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2949 double r0kal = x0kal.perp();
2958 cout <<
"*** move from " << track.
pivot() <<
" ( Kappa = "
2959 << track.
a()[2] <<
")" << endl;
2967 cout <<
"*** to " << track.
pivot() <<
" ( Kappa = "
2968 << track.
a()[2] <<
")" << std::endl;
2970 if (nhits_read == 1) {
2972 if(
debug_==4) cout <<
"filter_fwd::Ea = " << track.
Ea()<<endl;
2985 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2986 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2988 double diff_chi2 = track.
chiSq();
2990 Hep3Vector IP(0,0,0);
2996 for(
unsigned k=i+1 ; k < nhit; k++ )
3001 double dchi2 = -1.0;
3004 std::cout<<
"the value of x0kal is "<<x0kal<<std::endl;
3005 std::cout<<
"the value of track.pivot() is "<<track.
pivot()<<std::endl;
3009 HepSymMatrix Ma(5,0);
3013 std::cout<<
"HelixSeg.Ea_pre_fwd() ..."<<HelixSeg.
Ea_pre_fwd()<<std::endl;
3014 std::cout<<
"HelixSeg.a_pre_fwd() ..."<<HelixSeg.
a_pre_fwd()<<std::endl;
3015 std::cout<<
"HelixSeg.Ea_filt_fwd() ..."<<HelixSeg.
Ea_filt_fwd()<<std::endl;
3028 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;
3030 if(
debug_ == 4) cout<<
"****inext***"<<inext<<
" !!!! dchi2= "<< dchi2
3031 <<
" chisq= "<< chi2<< endl;
3036 StatusCode sc5 = m_nt5->write();
3037 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3043 diff_chi2 = chi2 - diff_chi2;
3044 HitMdc.
chi2(diff_chi2);
3048 cout <<
" -> after include meas = " << meas
3049 <<
" at R = " << track.
pivot().perp() << std::endl;
3050 cout <<
" chi2 = " << chi2 <<
", diff_chi2 = "
3051 << diff_chi2 <<
", LR = "
3052 << HitMdc.
LR() <<
", stereo = " << HitMdc.
wire().
stereo()
3054 cout<<
"filter_fwd..........HitMdc.chi2... "<<HitMdc.
chi2()<<endl;
3057 if(dchi2>0.0 && (way!=999)) {
3068 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3069 for(
int i = 0; i < _BesKalmanFitWalls.size(); i++) {
3070 _BesKalmanFitWalls[i].updateTrack(track, way);
3071 if (
debug_ == 4) cout<<
"Wall "<<i<<
" update the track!(filter)"<<endl;
3085 MsgStream log(
msgSvc(), name());
3086 double Pt_threshold(0.3);
3087 Hep3Vector IP(0,0,0);
3093 if ( !&whMgr )
return;
3096 int ntrk = mdcMgr->size();
3097 double* rPt =
new double[ntrk];
3098 int* rOM =
new int[ntrk];
3099 unsigned int* order =
new unsigned int[ntrk];
3100 unsigned int* rCont =
new unsigned int[ntrk];
3101 unsigned int* rGen =
new unsigned int[ntrk];
3104 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3105 end = mdcMgr->end(); it != end; it++) {
3109 rPt[index] = 1 / fabs(it->helix[2]);
3110 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3111 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3113 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3114 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3116 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3117 ii !=pt.begin()-1; ii--) {
3118 int lyr((*ii)->geo->Lyr()->Id());
3119 if (outermost < lyr) outermost = lyr;
3120 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3122 rOM[index] = outermost;
3123 order[index] = index;
3128 for (
int j, k = ntrk - 1; k >= 0; k = j){
3130 for(
int i = 1; i <= k; i++)
3131 if(rPt[i - 1] < rPt[i]){
3133 std::swap(order[i], order[j]);
3134 std::swap(rPt[i], rPt[j]);
3135 std::swap(rOM[i], rOM[j]);
3136 std::swap(rCont[i], rCont[j]);
3137 std::swap(rGen[i], rGen[j]);
3145 DataObject *aReconEvent;
3146 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3150 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3151 if(sc!=StatusCode::SUCCESS) {
3152 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3160 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3162 for(
int l = 0; l < ntrk; l++) {
3172 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3176 int trasqual = TrasanTRK_add.
quality;
3177 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3178 if (trasqual)
continue;
3182 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3186 if ((TrasanTRK_add.
decision & 32) == 32 ||
3187 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3192 TrasanTRK.
pivot[2]);
3195 for(
int i = 0; i < 5; i++)
3196 a[i] = TrasanTRK.
helix[i];
3199 for(
int i = 0, k = 0; i < 5; i++) {
3200 for(
int j = 0; j <= i; j++) {
3202 ea[j][i] = ea[i][j];
3205 double fiTerm = TrasanTRK.
fiTerm;
3213 int inlyr(999), outlyr(-1);
3214 int* rStat =
new int[43];
3215 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3216 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3218 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3221 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3222 ii != pt.begin()-1; ii--) {
3223 Num[(*ii)->geo->Lyr()->Id()]++;
3226 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3227 ii != pt.begin()-1; ii--) {
3229 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
3231 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3232 <<
" hits in the layer "
3233 << (*ii)->geo->Lyr()->Id() << std::endl;
3251 double dist[2] = {rechit.
ddl, rechit.
ddr};
3252 double erdist[2] = {rechit.
erddl, rechit.
erddr};
3257 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
3259 else if (rechit.
lr==1) lr_decision=1;
3262 int ind(geo->
Lyr()->
Id());
3266 lr_decision, rechit.
tdc,
3271 if (inlyr>ind) inlyr = ind;
3272 if (outlyr<ind) outlyr = ind;
3275 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3278 int empty_between(0), empty(0);
3279 for (
int i= inlyr; i <= outlyr; i++)
3280 if (!rStat[i]) empty_between++;
3281 empty = empty_between+inlyr+(42-outlyr);
3288 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
3294 track_lead.
type(type);
3295 unsigned int nhit = track_lead.
HitsMdc().size();
3296 if (!nhit &&
debug_ == 4) {
3297 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3302 double KalFitst(0), KalFitax(0), KalFitschi2(0);
3305 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
3307 track_lead.
pivot(outer_pivot);
3313 HepSymMatrix Eakal(5,0);
3316 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
3326 cout <<
"from Mdc Pattern Recognition: " << std::endl;
3332 cout <<
" dr = " << work.
a()[0]
3333 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3334 cout <<
" phi0 = " << work.
a()[1]
3335 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3336 cout <<
" PT = " << 1/work.
a()[2]
3337 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3338 cout <<
" dz = " << work.
a()[3]
3339 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3340 cout <<
" tanl = " << work.
a()[4]
3341 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3351 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3356 cout <<
" dr = " << work.
a()[0]
3357 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3358 cout <<
" phi0 = " << work.
a()[1]
3359 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3360 cout <<
" PT = " << 1/work.
a()[2]
3361 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3362 cout <<
" dz = " << work.
a()[3]
3363 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3364 cout <<
" tanl = " << work.
a()[4]
3365 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3371 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
3375 IDataProviderSvc* eventSvc = NULL;
3376 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
3378 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
3380 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
3385 IDataManagerSvc *dataManSvc;
3386 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
3387 DataObject *aKalTrackCol;
3388 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
3389 if(aKalTrackCol != NULL) {
3390 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
3391 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
3394 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
3395 if( kalsc.isFailure() ) {
3396 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
3399 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
3403 DataObject *aRecKalSegEvent;
3404 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
3405 if(aRecKalSegEvent!=NULL) {
3407 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
3408 if(segsc != StatusCode::SUCCESS) {
3409 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
3414 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
3415 if( segsc.isFailure() ) {
3416 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
3419 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
3422 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.);
3423 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
3424 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
3427 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
3429 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
3432 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
3433 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
3434 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
3435 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
3436 <<
"Track Id: " << (*iter_trk)->getTrackId()
3437 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
3438 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
3439 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
3440 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
3441 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
3442 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
3443 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
3444 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
3446 for(
int i = 0; i<43; i++) {
3447 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
3448 << (*iter_trk)->getPathl(i) <<endreq;
3452 m_trackid = (*iter_trk)->getTrackId();
3454 for(
int jj =0, iii=0; jj<5; jj++) {
3456 m_length[jj] = (*iter_trk)->getLength(jj);
3457 m_tof[jj] = (*iter_trk)->getTof(jj);
3458 m_nhits[jj] = (*iter_trk)->getNhits(jj);
3459 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
3460 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
3461 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
3462 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
3463 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
3464 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
3465 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
3466 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
3467 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
3468 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
3469 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
3470 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
3471 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
3472 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
3473 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
3476 for(
int kk=0; kk<=jj; kk++,iii++) {
3477 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
3478 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
3479 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
3480 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
3481 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
3482 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
3483 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
3484 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
3485 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
3486 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
3487 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
3488 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
3489 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
3490 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
3491 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
3531 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
3532 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
3533 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
3534 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
3535 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
3536 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
3537 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
3538 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
3539 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
3540 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
3542 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
3543 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
3544 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
3545 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
3546 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
3547 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
3548 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
3549 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
3550 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
3551 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
3553 m_stat[0][0] = (*iter_trk)->getStat(0,0);
3554 m_stat[1][0] = (*iter_trk)->getStat(0,1);
3555 m_stat[2][0] = (*iter_trk)->getStat(0,2);
3556 m_stat[3][0] = (*iter_trk)->getStat(0,3);
3557 m_stat[4][0] = (*iter_trk)->getStat(0,4);
3558 m_stat[0][1] = (*iter_trk)->getStat(1,0);
3559 m_stat[1][1] = (*iter_trk)->getStat(1,1);
3560 m_stat[2][1] = (*iter_trk)->getStat(1,2);
3561 m_stat[3][1] = (*iter_trk)->getStat(1,3);
3562 m_stat[4][1] = (*iter_trk)->getStat(1,4);
3564 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
3565 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
3566 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
3567 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
3568 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
3570 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
3571 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
3572 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
3573 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
3574 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
3576 m_lpt = 1/m_lhelix[2];
3577 m_lpte = 1/m_lhelixe[2];
3578 m_lptmu = 1/m_lhelixmu[2];
3579 m_lptk = 1/m_lhelixk[2];
3580 m_lptp = 1/m_lhelixp[2];
3582 m_fpt = 1/m_fhelix[2];
3583 m_fpte = 1/m_fhelixe[2];
3584 m_fptmu = 1/m_fhelixmu[2];
3585 m_fptk = 1/m_fhelixk[2];
3586 m_fptp = 1/m_fhelixp[2];
3589 std::cout<<
" "<<std::endl;
3590 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
3591 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
3592 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
3593 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
3594 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
3597 m_zpt = 1/m_zhelix[2];
3598 m_zpte = 1/m_zhelixe[2];
3599 m_zptmu = 1/m_zhelixmu[2];
3600 m_zptk = 1/m_zhelixk[2];
3601 m_zptp = 1/m_zhelixp[2];
3604 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
3605 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
3606 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
3607 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
3608 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
3610 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
3611 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
3612 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
3613 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
3614 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
3617 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
3618 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
3619 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
3620 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
3621 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
3625 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
3626 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
3627 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
3628 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
3629 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
3630 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
3631 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
3632 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
3633 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
3634 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
3635 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
3636 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
3637 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
3638 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
3639 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
3642 StatusCode sc1 = m_nt1->write();
3643 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
3648 phi1 = (*iter_trk)->getFFi0();
3649 r1 = (*iter_trk)->getFDr();
3650 z1 = (*iter_trk)->getFDz();
3651 kap1 = (*iter_trk)->getFCpa();
3652 tanl1 = (*iter_trk)->getFTanl();
3653 charge1 = kap1/fabs(kap1);
3656 p1 = sqrt(1+tanl1*tanl1)/kap1;
3657 the1 =
M_PI/2-atan(tanl1);
3660 pz1= tanl1/fabs(kap1);
3662 }
else if(jj == 2) {
3663 phi2 = (*iter_trk)->getFFi0();
3664 r2 = (*iter_trk)->getFDr();
3665 z2 = (*iter_trk)->getFDz();
3666 kap2 = (*iter_trk)->getFCpa();
3667 tanl2 = (*iter_trk)->getFTanl();
3668 charge2 = kap2/fabs(kap2);
3671 p2 = sqrt(1+tanl2*tanl2)/kap1;
3672 the2 =
M_PI/2-atan(tanl2);
3675 pz2= tanl2/fabs(kap2);
3683 m_delthe = the1 + the2;
3686 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
3687 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
3688 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
3690 StatusCode sc2 = m_nt2->write();
3691 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
3700 cout <<
"Kalfitting finished " << std::endl;
3706 MsgStream log(
msgSvc(), name());
3707 double Pt_threshold(0.3);
3708 Hep3Vector IP(0,0,0);
3715 if ( !&whMgr )
return;
3718 int ntrk = mdcMgr->size();
3719 double* rPt =
new double[ntrk];
3720 int* rOM =
new int[ntrk];
3721 unsigned int* order =
new unsigned int[ntrk];
3722 unsigned int* rCont =
new unsigned int[ntrk];
3723 unsigned int* rGen =
new unsigned int[ntrk];
3726 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3727 end = mdcMgr->end(); it != end; it++) {
3731 rPt[index] = 1 / fabs(it->helix[2]);
3732 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3733 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3735 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3736 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3738 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3739 ii !=pt.begin()-1; ii--) {
3740 int lyr((*ii)->geo->Lyr()->Id());
3741 if (outermost < lyr) outermost = lyr;
3742 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3744 rOM[index] = outermost;
3745 order[index] = index;
3750 for (
int j, k = ntrk - 1; k >= 0; k = j){
3752 for(
int i = 1; i <= k; i++)
3753 if(rPt[i - 1] < rPt[i]){
3755 std::swap(order[i], order[j]);
3756 std::swap(rPt[i], rPt[j]);
3757 std::swap(rOM[i], rOM[j]);
3758 std::swap(rCont[i], rCont[j]);
3759 std::swap(rGen[i], rGen[j]);
3766 DataObject *aReconEvent;
3767 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3771 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3772 if(sc!=StatusCode::SUCCESS) {
3773 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3781 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3784 for(
int l = 0; l < ntrk; l++) {
3786 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3790 int trasqual = TrasanTRK_add.
quality;
3791 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3792 if (trasqual)
continue;
3796 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3800 if ((TrasanTRK_add.
decision & 32) == 32 ||
3801 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3806 TrasanTRK.
pivot[2]);
3809 for(
int i = 0; i < 5; i++)
3810 a[i] = TrasanTRK.
helix[i];
3813 for(
int i = 0, k = 0; i < 5; i++) {
3814 for(
int j = 0; j <= i; j++) {
3816 ea[j][i] = ea[i][j];
3822 double fiTerm = TrasanTRK.
fiTerm;
3828 int inlyr(999), outlyr(-1);
3829 int* rStat =
new int[43];
3830 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3831 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3833 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3836 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3837 ii != pt.begin()-1; ii--) {
3838 Num[(*ii)->geo->Lyr()->Id()]++;
3842 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3843 ii != pt.begin()-1; ii--) {
3846 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
3848 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3849 <<
" hits in the layer "
3850 << (*ii)->geo->Lyr()->Id() << std::endl;
3877 double dist[2] = {rechit.
ddl, rechit.
ddr};
3878 double erdist[2] = {rechit.
erddl, rechit.
erddr};
3883 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
3885 else if (rechit.
lr==1) lr_decision=1;
3888 int ind(geo->
Lyr()->
Id());
3890 lr_decision, rechit.
tdc,
3895 if (inlyr>ind) inlyr = ind;
3896 if (outlyr<ind) outlyr = ind;
3899 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3902 int empty_between(0), empty(0);
3903 for (
int i= inlyr; i <= outlyr; i++)
3904 if (!rStat[i]) empty_between++;
3905 empty = empty_between+inlyr+(42-outlyr);
3910 track_lead.
type(type);
3911 unsigned int nhit = track_lead.
HitsMdc().size();
3912 if (!nhit &&
debug_ == 4) {
3913 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3918 double KalFitst(0), KalFitax(0), KalFitschi2(0);
3920 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
3923 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
3925 track_lead.
pivot(outer_pivot);
3930 HepSymMatrix Eakal(5,0);
3934 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
3944 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
3950 std::cout <<
" dr = " << work.
a()[0]
3951 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3952 std::cout <<
" phi0 = " << work.
a()[1]
3953 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3954 std::cout <<
" PT = " << 1/work.
a()[2]
3955 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3956 std::cout <<
" dz = " << work.
a()[3]
3957 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3958 std::cout <<
" tanl = " << work.
a()[4]
3959 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3967 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3972 cout <<
" dr = " << work.
a()[0]
3973 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3974 cout <<
" phi0 = " << work.
a()[1]
3975 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3976 cout <<
" PT = " << 1/work.
a()[2]
3977 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3978 cout <<
" dz = " << work.
a()[3]
3979 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3980 cout <<
" tanl = " << work.
a()[4]
3981 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3988 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
3994 DataObject *aRecKalEvent;
3995 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
3996 if(aRecKalEvent!=NULL) {
3998 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
3999 if(kalsc != StatusCode::SUCCESS) {
4000 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4005 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4006 if( kalsc.isFailure()) {
4007 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4010 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4016 DataObject *aRecKalSegEvent;
4017 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4018 if(aRecKalSegEvent!=NULL) {
4020 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4021 if(segsc != StatusCode::SUCCESS) {
4022 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4027 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4028 if( segsc.isFailure() ) {
4029 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4032 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4035 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.);
4036 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4039 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4041 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4044 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4045 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4046 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4047 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4048 <<
"Track Id: " << (*iter_trk)->getTrackId()
4049 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4050 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4051 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4052 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4053 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4054 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4055 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4056 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4061 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4064 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4065 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4067 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4068 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4069 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4070 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4071 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4072 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4075 for(
int i = 0; i<43; i++) {
4076 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4077 << (*iter_trk)->getPathl(i) <<endreq;
4081 m_trackid = (*iter_trk)->getTrackId();
4082 for(
int jj =0, iii=0; jj<5; jj++) {
4083 m_length[jj] = (*iter_trk)->getLength(jj);
4084 m_tof[jj] = (*iter_trk)->getTof(jj);
4085 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4086 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4087 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4088 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4089 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4090 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4091 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4092 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4093 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4094 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4095 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4096 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4097 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4098 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4099 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4100 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4102 for(
int kk=0; kk<=jj; kk++,iii++) {
4103 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4104 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4105 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4106 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4107 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4108 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4109 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4110 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4111 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4112 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4113 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4114 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4115 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4116 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4117 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4158 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4159 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4160 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4161 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4162 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4163 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4164 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4165 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4166 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4167 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4169 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4170 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4171 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4172 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4173 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4174 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4175 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4176 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4177 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4178 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4180 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4181 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4182 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4183 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4184 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4185 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4186 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4187 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4188 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4189 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4191 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4192 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4193 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4194 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4195 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4197 m_zpt = 1/m_zhelix[2];
4198 m_zpte = 1/m_zhelixe[2];
4199 m_zptmu = 1/m_zhelixmu[2];
4200 m_zptk = 1/m_zhelixk[2];
4201 m_zptp = 1/m_zhelixp[2];
4203 m_fpt = 1/m_fhelix[2];
4204 m_fpte = 1/m_fhelixe[2];
4205 m_fptmu = 1/m_fhelixmu[2];
4206 m_fptk = 1/m_fhelixk[2];
4207 m_fptp = 1/m_fhelixp[2];
4209 m_lpt = 1/m_lhelix[2];
4210 m_lpte = 1/m_lhelixe[2];
4211 m_lptmu = 1/m_lhelixmu[2];
4212 m_lptk = 1/m_lhelixk[2];
4213 m_lptp = 1/m_lhelixp[2];
4215 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4216 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4217 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4218 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4219 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4221 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4222 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4223 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4224 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4225 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4227 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4228 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4229 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4230 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4231 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4232 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4233 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4234 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4235 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4236 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4237 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4238 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4239 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4240 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4241 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4244 StatusCode sc1 = m_nt1->write();
4245 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4250 phi1 = (*iter_trk)->getFFi0();
4251 r1 = (*iter_trk)->getFDr();
4252 z1 = (*iter_trk)->getFDz();
4253 kap1 = (*iter_trk)->getFCpa();
4254 tanl1 = (*iter_trk)->getFTanl();
4257 p1 = sqrt(1+tanl1*tanl1)/kap1;
4258 the1 =
M_PI/2-atan(tanl1);
4259 }
else if(jj == 2) {
4260 phi2 = (*iter_trk)->getFFi0();
4261 r2 = (*iter_trk)->getFDr();
4262 z2 = (*iter_trk)->getFDz();
4263 kap2 = (*iter_trk)->getFCpa();
4264 tanl2 = (*iter_trk)->getFTanl();
4267 p2 = sqrt(1+tanl2*tanl2)/kap1;
4268 the2 =
M_PI/2-atan(tanl2);
4276 m_delthe = the1 + the2;
4279 StatusCode sc2 = m_nt2->write();
4280 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4288 cout <<
"Kalfitting finished " << std::endl;
4295 MsgStream log(
msgSvc(), name());
4296 double Pt_threshold(0.3);
4297 Hep3Vector IP(0,0,0);
4304 if ( !&whMgr )
return;
4307 int ntrk = mdcMgr->size();
4310 int nhits = whMgr->size();
4315 DataObject *aReconEvent;
4316 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4320 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4321 if(sc!=StatusCode::SUCCESS) {
4322 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4330 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4342 if ((TrasanTRK_add.
decision & 32) == 32 ||
4343 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4348 TrasanTRK.
pivot[2]);
4351 for(
int i = 0; i < 5; i++)
4352 a[i] = TrasanTRK.
helix[i];
4355 for(
int i = 0, k = 0; i < 5; i++) {
4356 for(
int j = 0; j <= i; j++) {
4358 ea[j][i] = ea[i][j];
4364 double fiTerm = TrasanTRK.
fiTerm;
4372 int trasqual = TrasanTRK_add.
quality;
4373 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4374 if (trasqual)
return;
4376 int inlyr(999), outlyr(-1);
4377 int* rStat =
new int[43];
4378 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4379 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4381 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4384 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4385 ii != pt.begin()-1; ii--) {
4386 Num[(*ii)->geo->Lyr()->Id()]++;
4389 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4390 ii != pt.begin()-1; ii--) {
4393 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4395 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4396 <<
" hits in the layer "
4397 << (*ii)->geo->Lyr()->Id() << std::endl;
4403 double dist[2] = {rechit.
ddl, rechit.
ddr};
4404 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4409 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4411 else if (rechit.
lr==1) lr_decision=1;
4414 int ind(geo->
Lyr()->
Id());
4416 lr_decision, rechit.
tdc,
4421 if (inlyr>ind) inlyr = ind;
4422 if (outlyr<ind) outlyr = ind;
4425 int empty_between(0), empty(0);
4426 for (
int i= inlyr; i <= outlyr; i++)
4427 if (!rStat[i]) empty_between++;
4428 empty = empty_between+inlyr+(42-outlyr);
4432 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4437 track_lead.
type(type);
4438 unsigned int nhit = track_lead.
HitsMdc().size();
4440 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4445 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4447 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4450 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4452 track_lead.
pivot(outer_pivot);
4457 HepSymMatrix Eakal(5,0);
4461 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4471 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4477 std::cout <<
" dr = " << work.
a()[0]
4478 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4479 std::cout <<
" phi0 = " << work.
a()[1]
4480 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4481 std::cout <<
" PT = " << 1/work.
a()[2]
4482 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4483 std::cout <<
" dz = " << work.
a()[3]
4484 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4485 std::cout <<
" tanl = " << work.
a()[4]
4486 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4493 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4498 cout <<
" dr = " << work1.
a()[0]
4499 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
4500 cout <<
" phi0 = " << work1.
a()[1]
4501 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
4502 cout <<
" PT = " << 1/work1.
a()[2]
4503 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
4504 cout <<
" dz = " << work1.
a()[3]
4505 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
4506 cout <<
" tanl = " << work1.
a()[4]
4507 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
4514 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4519 DataObject *aRecKalEvent;
4520 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4521 if(aRecKalEvent!=NULL) {
4523 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4524 if(kalsc != StatusCode::SUCCESS) {
4525 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4530 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4531 if( kalsc.isFailure()) {
4532 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4535 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4541 DataObject *aRecKalSegEvent;
4542 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4543 if(aRecKalSegEvent!=NULL) {
4545 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4546 if(segsc != StatusCode::SUCCESS) {
4547 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4552 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4553 if( segsc.isFailure() ) {
4554 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4557 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4560 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.);
4561 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4564 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4566 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4569 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4570 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4571 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4572 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4573 <<
"Track Id: " << (*iter_trk)->getTrackId()
4574 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4575 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4576 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4577 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4578 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4579 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4580 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4581 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4582 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
4587 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4590 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4591 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4593 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4594 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4595 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4596 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4597 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4598 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4601 for(
int i = 0; i<43; i++) {
4602 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4603 << (*iter_trk)->getPathl(i) <<endreq;
4607 m_trackid = (*iter_trk)->getTrackId();
4608 for(
int jj =0, iii=0; jj<5; jj++) {
4609 m_length[jj] = (*iter_trk)->getLength(jj);
4610 m_tof[jj] = (*iter_trk)->getTof(jj);
4611 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4612 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4613 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4614 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4615 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4616 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4617 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4618 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4619 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4620 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4621 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4622 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4623 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4624 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4625 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4626 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4628 for(
int kk=0; kk<=jj; kk++,iii++) {
4629 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4630 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4631 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4632 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4633 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4634 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4635 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4636 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4637 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4638 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4639 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4640 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4641 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4642 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4643 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4650 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4651 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4652 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4653 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4654 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4655 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4656 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4657 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4658 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4659 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4661 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4662 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4663 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4664 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4665 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4666 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4667 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4668 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4669 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4670 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4672 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4673 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4674 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4675 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4676 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4677 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4678 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4679 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4680 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4681 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4683 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4684 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4685 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4686 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4687 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4689 m_zpt = 1/m_zhelix[2];
4690 m_zpte = 1/m_zhelixe[2];
4691 m_zptmu = 1/m_zhelixmu[2];
4692 m_zptk = 1/m_zhelixk[2];
4693 m_zptp = 1/m_zhelixp[2];
4695 m_fpt = 1/m_fhelix[2];
4696 m_fpte = 1/m_fhelixe[2];
4697 m_fptmu = 1/m_fhelixmu[2];
4698 m_fptk = 1/m_fhelixk[2];
4699 m_fptp = 1/m_fhelixp[2];
4701 m_lpt = 1/m_lhelix[2];
4702 m_lpte = 1/m_lhelixe[2];
4703 m_lptmu = 1/m_lhelixmu[2];
4704 m_lptk = 1/m_lhelixk[2];
4705 m_lptp = 1/m_lhelixp[2];
4707 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4708 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4709 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4710 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4711 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4713 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4714 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4715 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4716 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4717 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4719 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4720 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4721 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4722 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4723 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4724 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4725 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4726 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4727 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4728 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4729 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4730 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4731 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4732 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4733 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4736 StatusCode sc1 = m_nt1->write();
4737 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4742 phi1 = (*iter_trk)->getFFi0();
4743 r1 = (*iter_trk)->getFDr();
4744 z1 = (*iter_trk)->getFDz();
4745 kap1 = (*iter_trk)->getFCpa();
4746 tanl1 = (*iter_trk)->getFTanl();
4749 p1 = sqrt(1+tanl1*tanl1)/kap1;
4750 the1 =
M_PI/2-atan(tanl1);
4751 }
else if(jj == 2) {
4752 phi2 = (*iter_trk)->getFFi0();
4753 r2 = (*iter_trk)->getFDr();
4754 z2 = (*iter_trk)->getFDz();
4755 kap2 = (*iter_trk)->getFCpa();
4756 tanl2 = (*iter_trk)->getFTanl();
4759 p2 = sqrt(1+tanl2*tanl2)/kap1;
4760 the2 =
M_PI/2-atan(tanl2);
4768 m_delthe = the1 + the2;
4771 StatusCode sc2 = m_nt2->write();
4772 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4776 cout <<
"Kalfitting finished " << std::endl;
4784 MsgStream log(
msgSvc(), name());
4785 double Pt_threshold(0.3);
4786 Hep3Vector IP(0,0,0);
4793 if ( !&whMgr )
return;
4796 int ntrk = mdcMgr->size();
4799 int nhits = whMgr->size();
4803 double* rY =
new double[ntrk];
4804 double* rfiTerm =
new double[ntrk];
4805 double* rPt =
new double[ntrk];
4806 int* rOM =
new int[ntrk];
4807 unsigned int* order =
new unsigned int[ntrk];
4808 unsigned int* rCont =
new unsigned int[ntrk];
4809 unsigned int* rGen =
new unsigned int[ntrk];
4812 Hep3Vector csmp3[2];
4813 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
4814 end = mdcMgr->end(); it != end; it++) {
4816 rfiTerm[index]=it->fiTerm;
4821 rPt[index] = 1 / fabs(it->helix[2]);
4822 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4823 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4825 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4826 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4828 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4829 ii !=pt.begin()-1; ii--) {
4830 int lyr((*ii)->geo->Lyr()->Id());
4831 if (outermost < lyr) {
4833 rY[index] = (*ii)->geo->Forward().y();
4835 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4837 rOM[index] = outermost;
4838 order[index] = index;
4843 for (
int j, k = ntrk - 1; k >= 0; k = j){
4845 for(
int i = 1; i <= k; i++)
4846 if(rY[i - 1] < rY[i]){
4848 std::swap(order[i], order[j]);
4849 std::swap(rY[i], rY[j]);
4850 std::swap(rOM[i], rOM[j]);
4851 std::swap(rCont[i], rCont[j]);
4852 std::swap(rGen[i], rGen[j]);
4861 DataObject *aReconEvent;
4862 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4866 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4867 if(sc!=StatusCode::SUCCESS) {
4868 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4876 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4883 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
4892 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4896 if ((TrasanTRK_add.
decision & 32) == 32 ||
4897 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4902 TrasanTRK.
pivot[2]);
4905 for(
int i = 0; i < 5; i++)
4906 a[i] = TrasanTRK.
helix[i];
4909 for(
int i = 0, k = 0; i < 5; i++) {
4910 for(
int j = 0; j <= i; j++) {
4912 ea[j][i] = ea[i][j];
4918 double fiTerm = TrasanTRK.
fiTerm;
4925 for(
int l = 0; l < ntrk; l++) {
4926 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
4927 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
4929 int trasqual = TrasanTRK_add1.
quality;
4930 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4931 if (trasqual)
continue;
4933 int inlyr(999), outlyr(-1);
4934 int* rStat =
new int[43];
4935 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4936 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
4938 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4941 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4942 ii != pt.begin()-1; ii--) {
4943 Num[(*ii)->geo->Lyr()->Id()]++;
4946 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4947 ii != pt.begin()-1; ii--) {
4950 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4952 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4953 <<
" hits in the layer "
4954 << (*ii)->geo->Lyr()->Id() << std::endl;
4960 double dist[2] = {rechit.
ddl, rechit.
ddr};
4961 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4966 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4968 else if (rechit.
lr==1) lr_decision=1;
4971 int ind(geo->
Lyr()->
Id());
4973 lr_decision, rechit.
tdc,
4978 if (inlyr>ind) inlyr = ind;
4979 if (outlyr<ind) outlyr = ind;
4982 int empty_between(0), empty(0);
4983 for (
int i= inlyr; i <= outlyr; i++)
4984 if (!rStat[i]) empty_between++;
4985 empty = empty_between+inlyr+(42-outlyr);
4989 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4994 track_lead.
type(type);
4995 unsigned int nhit = track_lead.
HitsMdc().size();
4997 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5002 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5004 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5007 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5009 track_lead.
pivot(outer_pivot);
5014 HepSymMatrix Eakal(5,0);
5018 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5028 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5034 std::cout <<
" dr = " << work.
a()[0]
5035 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5036 std::cout <<
" phi0 = " << work.
a()[1]
5037 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5038 std::cout <<
" PT = " << 1/work.
a()[2]
5039 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5040 std::cout <<
" dz = " << work.
a()[3]
5041 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5042 std::cout <<
" tanl = " << work.
a()[4]
5043 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5050 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5055 cout <<
" dr = " << work1.
a()[0]
5056 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5057 cout <<
" phi0 = " << work1.
a()[1]
5058 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5059 cout <<
" PT = " << 1/work1.
a()[2]
5060 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5061 cout <<
" dz = " << work1.
a()[3]
5062 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5063 cout <<
" tanl = " << work1.
a()[4]
5064 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5071 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5077 DataObject *aRecKalEvent;
5078 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5079 if(aRecKalEvent!=NULL) {
5081 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5082 if(kalsc != StatusCode::SUCCESS) {
5083 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5088 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5089 if( kalsc.isFailure()) {
5090 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5093 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5099 DataObject *aRecKalSegEvent;
5100 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5101 if(aRecKalSegEvent!=NULL) {
5103 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5104 if(segsc != StatusCode::SUCCESS) {
5105 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5110 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5111 if( segsc.isFailure() ) {
5112 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5115 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5118 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.);
5119 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5122 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5124 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5127 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5128 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5129 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5130 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5131 <<
"Track Id: " << (*iter_trk)->getTrackId()
5132 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5133 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5134 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5135 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5136 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5137 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5138 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5139 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5140 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5145 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5148 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5149 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5151 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5152 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5153 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5154 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5155 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5156 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5159 for(
int i = 0; i<43; i++) {
5160 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5161 << (*iter_trk)->getPathl(i) <<endreq;
5165 m_trackid = (*iter_trk)->getTrackId();
5166 for(
int jj =0, iii=0; jj<5; jj++) {
5167 m_length[jj] = (*iter_trk)->getLength(jj);
5168 m_tof[jj] = (*iter_trk)->getTof(jj);
5169 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5170 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5171 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5172 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5173 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5174 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5175 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5176 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5177 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5178 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5179 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5180 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5181 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5182 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5183 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5184 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5186 for(
int kk=0; kk<=jj; kk++,iii++) {
5187 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5188 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5189 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5190 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5191 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5192 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5193 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5194 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5195 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5196 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5197 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5198 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5199 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5200 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5201 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5208 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5209 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5210 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5211 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5212 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5213 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5214 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5215 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5216 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5217 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5219 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5220 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5221 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5222 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5223 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5224 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5225 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5226 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5227 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5228 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5230 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5231 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5232 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5233 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5234 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5235 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5236 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5237 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5238 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5239 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5241 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5242 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5243 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5244 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5245 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5247 m_zpt = 1/m_zhelix[2];
5248 m_zpte = 1/m_zhelixe[2];
5249 m_zptmu = 1/m_zhelixmu[2];
5250 m_zptk = 1/m_zhelixk[2];
5251 m_zptp = 1/m_zhelixp[2];
5253 m_fpt = 1/m_fhelix[2];
5254 m_fpte = 1/m_fhelixe[2];
5255 m_fptmu = 1/m_fhelixmu[2];
5256 m_fptk = 1/m_fhelixk[2];
5257 m_fptp = 1/m_fhelixp[2];
5259 m_lpt = 1/m_lhelix[2];
5260 m_lpte = 1/m_lhelixe[2];
5261 m_lptmu = 1/m_lhelixmu[2];
5262 m_lptk = 1/m_lhelixk[2];
5263 m_lptp = 1/m_lhelixp[2];
5265 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5266 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5267 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5268 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5269 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5271 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5272 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5273 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5274 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5275 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5277 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5278 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5279 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5280 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5281 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5282 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5283 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5284 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5285 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5286 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5287 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5288 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5289 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5290 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5291 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5294 StatusCode sc1 = m_nt1->write();
5295 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5300 phi1 = (*iter_trk)->getFFi0();
5301 r1 = (*iter_trk)->getFDr();
5302 z1 = (*iter_trk)->getFDz();
5303 kap1 = (*iter_trk)->getFCpa();
5304 tanl1 = (*iter_trk)->getFTanl();
5307 p1 = sqrt(1+tanl1*tanl1)/kap1;
5308 the1 =
M_PI/2-atan(tanl1);
5309 }
else if(jj == 2) {
5310 phi2 = (*iter_trk)->getFFi0();
5311 r2 = (*iter_trk)->getFDr();
5312 z2 = (*iter_trk)->getFDz();
5313 kap2 = (*iter_trk)->getFCpa();
5314 tanl2 = (*iter_trk)->getFTanl();
5317 p2 = sqrt(1+tanl2*tanl2)/kap1;
5318 the2 =
M_PI/2-atan(tanl2);
5326 m_delthe = the1 + the2;
5329 StatusCode sc2 = m_nt2->write();
5330 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5338 cout <<
"Kalfitting finished " << std::endl;
5351 MsgStream log(
msgSvc(), name());
5361 double pp = track_lead.
momentum().mag();
5365 for(
int l_mass = 0, flg = 1; l_mass < nmass;
5366 l_mass++, flg <<= 1) {
5368 if (!(
mhyp_ & flg))
continue;
5369 if (l_mass ==
lead_)
continue;
5378 if(
debug_ == 4) cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
5386 TrasanTRK.
pivot[2]);
5387 HepVector a_trasan(5);
5388 for(
int i = 0; i < 5; i++)
5389 a_trasan[i] = TrasanTRK.
helix[i];
5391 HepSymMatrix ea_trasan(5);
5392 for(
int i = 0, k = 0; i < 5; i++) {
5393 for(
int j = 0; j <= i; j++) {
5395 ea_trasan[j][i] = ea_trasan[i][j];
5401 double fiTerm = TrasanTRK.
fiTerm;
5402 track.
pivot(track.
x(fiTerm));
5403 HepSymMatrix Eakal(5,0);
5405 double costheta = track.
a()[4] / sqrt(1.0 + track.
a()[4]*track.
a()[4]);
5415 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
5417 fillTds(TrasanTRK, track, kaltrk, l_mass);
5427 HepSymMatrix ea_first(kaltrk->
getFError());
5430 for(
int i = 0; i < 5; i++)
5431 for(
int j = 0; j <= i; j++)
5432 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
5433 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
5439 cout <<
" Backward fitting flag:" <<
back_<< endl;
5440 cout <<
"track_back pivot " << track_back.
pivot()
5441 <<
" track_lead kappa " << track_lead.
a()[2]
5445 if (
back_ && track_lead.
a()[2] != 0 &&
5446 1/fabs(track_lead.
a()[2]) >
pT_) {
5450 double p_kaon(0), p_proton(0);
5451 if (!(kaltrk->
getStat(0,3))) {
5454 track_back.
p_kaon(p_kaon);
5456 p_kaon = 1 / fabs(track_back.
a()[2]) *
5457 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5458 track_back.
p_kaon(p_kaon);
5460 if (!(kaltrk->
getStat(0,4))) {
5461 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
5465 p_proton = 1 / fabs(track_back.
a()[2]) *
5466 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5473 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
5486 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass, segcol, 1);
5510 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
5512 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
5513 <<
"Length of the track: "<< kaltrk->
getLength(2)
5514 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
5515 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
5516 <<
" "<< kaltrk->
getChisq(1,2) << endreq
5517 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
5518 <<
" "<< kaltrk->
getNdf(1,2) << endreq
5522 kalcol->push_back(kaltrk);
5535 MsgStream log(
msgSvc(), name());
5540 cout <<
"track_first pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
5547 cout <<
"after inner wall, track_ip pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
5553 double pp = track_lead.
momentum().mag();
5558 for(
int l_mass = 0, flg = 1; l_mass < nmass;
5559 l_mass++, flg <<= 1) {
5561 if (!(
mhyp_ & flg))
continue;
5562 if (l_mass ==
lead_)
continue;
5573 cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
5583 TrasanTRK.
pivot[2]);
5585 HepVector a_trasan(5);
5586 for(
int i = 0; i < 5; i++){
5587 a_trasan[i] = TrasanTRK.
helix[i];
5590 HepSymMatrix ea_trasan(5);
5591 for(
int i = 0, k = 0; i < 5; i++) {
5592 for(
int j = 0; j <= i; j++) {
5594 ea_trasan[j][i] = ea_trasan[i][j];
5601 double fiTerm = TrasanTRK.
fiTerm;
5602 track.
pivot(track.
x(fiTerm));
5604 HepSymMatrix Eakal(5,0);
5606 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5618 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
5620 fillTds(TrasanTRK, track, kaltrk, l_mass);
5629 HepSymMatrix ea_first(kaltrk->
getFError());
5633 for(
int i = 0; i < 5; i++)
5634 for(
int j = 0; j <= i; j++)
5635 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
5636 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
5652 cout <<
" Backward fitting flag:" <<
back_<< endl;
5653 cout <<
"track_back pivot " << track_back.
pivot()
5654 <<
" track_lead kappa " << track_lead.
a()[2]
5658 if (
back_ && track_lead.
a()[2] != 0 &&
5659 1/fabs(track_lead.
a()[2]) >
pT_) {
5664 double p_kaon(0), p_proton(0);
5666 if (!(kaltrk->
getStat(0,3))) {
5669 track_back.
p_kaon(p_kaon);
5671 p_kaon = 1 / fabs(track_back.
a()[2]) *
5672 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5673 track_back.
p_kaon(p_kaon);
5675 if (!(kaltrk->
getStat(0,4))) {
5676 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
5680 p_proton = 1 / fabs(track_back.
a()[2]) *
5681 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5689 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
5694 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass,segcol);
5717 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
5719 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
5720 <<
"Length of the track: "<< kaltrk->
getLength(2)
5721 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
5722 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
5723 <<
" "<< kaltrk->
getChisq(1,2) << endreq
5724 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
5725 <<
" "<< kaltrk->
getNdf(1,2) << endreq
5729 kalcol->push_back(kaltrk);
5738 for (
int i=0; i<5; i++) {
5739 for(
int j = 1; j<i+2;j++) {
5741 Eakal(j,i+1) = Eakal(i+1,j);
5745 if (
debug_ == 4) cout<<
"initialised Ea.. "<<Eakal<<endl;
5753 for (
int i=0; i<5; i++) {
5754 for(
int j = 1; j<i+2;j++) {
5756 Eakal(j,i+1) = Eakal(i+1,j);
5758 Eakal(1,1) = Eakal(1,1)*
gain1_;
5759 Eakal(2,2) = Eakal(2,2)*
gain2_;
5760 Eakal(3,3) = Eakal(3,3)*
gain3_;
5761 Eakal(4,4) = Eakal(4,4)*
gain4_;
5762 Eakal(5,5) = Eakal(5,5)*
gain5_;
5828 if (
debug_ == 4) cout<<
"initialised Eakal.. "<<Eakal<<endl;
5837 cout <<
"start_seed begin... " << std::endl;
5839 Hep3Vector x_init(track.
pivot());
5840 HepSymMatrix Ea_init(5,0);
5841 Ea_init = track.
Ea();
5842 HepVector a_init(5);
5846 unsigned int nhit_included(10);
5858 unsigned int nhit = track.
HitsMdc().size();
5861 for (
int ktrial = 0; ktrial < 8; ktrial++) {
5864 track.
pivot(x_init);
5874 HepSymMatrix Eakal(5,0);
5878 for(
unsigned i=0 ; i < nhit; i++ ){
5881 if (i<3) lr_decision = LR[ktrial][i];
5882 HitMdc.
LR(lr_decision);
5883 if (i<nhit_included)
5896 cout <<
"---- Result for " << ktrial <<
" case : chi2 = " << track.
chiSq()
5897 <<
", nhits included = " << track.
nchits() <<
" for nhits available = "
5898 << nhit << std::endl;
5900 if (track.
chiSq() < chi2_min &&
5901 (track.
nchits() == nhit_included || track.
nchits() == nhit)){
5902 chi2_min = track.
chiSq();
5909 cout <<
"*** i_min = " << i_min <<
" with a chi2 = " << chi2_min << std::endl;
5911 for(
unsigned i=0 ; i < nhit; i++ ){
5914 if (i_min >= 0 && i < 3)
5915 lr_decision = LR[i_min][i];
5916 HitMdc.
LR(lr_decision);
5920 track.
pivot(x_init);
5931 for(
unsigned i=0 ; i < 3; i++ ){
5933 cout <<
" LR(" << i <<
") = " << HitMdc.
LR()
5945 if(
debug_ == 4) cout<<
"Begining to clear Tables ...."<<endl;
5949 vector<MdcRec_trk>::iterator tit=mdcMgr->begin();
5950 for( ; tit != mdcMgr->end(); tit++) {
5951 vector<MdcRec_wirhit*>::iterator vit= tit->hitcol.begin() ;
5952 for(; vit != tit->hitcol.end(); vit++) {
5958 mdc_addMgr->clear();
5970bool KalFitAlg::order_rechits(
const SmartRef<RecMdcHit>& m1,
const SmartRef<RecMdcHit>& m2) {
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
double sin(const BesAngle a)
double cos(const BesAngle a)
ObjectVector< RecMdcKalHelixSeg > RecMdcKalHelixSegCol
SmartRefVector< RecMdcKalHelixSeg > HelixSegRefVec
ObjectVector< RecMdcKalTrack > RecMdcKalTrackCol
SmartRefVector< RecMdcHit > HitRefVec
HepGeom::Point3D< double > HepPoint3D
HepGeom::Point3D< double > HepPoint3D
void setFHelix(const HepVector &fhelix, const int pid)
void setStat(int stat, int i, int pid)
void setY(const double y, const int pid)
void setPz(const double pz, const int pid)
void setChisq(double chisq, int i, int pid)
void setZError(const HepSymMatrix &error, const int pid)
void setX(const double x, const int pid)
void setPoca(const HepPoint3D &poca, const int pid)
void setPx(const double px, const int pid)
void setZ(const double z, const int pid)
void setTheta(const double theta, const int pid)
void setPy(const double py, const int pid)
void setFError(const HepSymMatrix &ferror, const int pid)
static void setPidType(PidType pidType)
void setTrackId(int trackId)
void setZHelix(const HepVector &helix, const int pid)
void setCharge(const int charge, const int pid)
void setNdf(int ndf, int i, int pid)
void setNlayer(int nlayer, int pid)
virtual double getReferField()=0
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
double pe_cut_
value of the momentum cut to decide refit
void filter_fwd_anal(KalFitTrack &trk, int l_mass, int way, HepSymMatrix &Eakal)
Kalman filter (forward) in Mdc.
void filter_fwd_calib(KalFitTrack &trk, int l_mass, int way, HepSymMatrix &Eakal)
double fstrag_
factor of energy loss straggling for electron
int back_
flag to perform smoothing
void fillTds_lead(MdcRec_trk &TrasanTRK, KalFitTrack &track, RecMdcKalTrack *trk, int l_mass)
~KalFitAlg(void)
destructor
void kalman_fitting_MdcxReco_Csmc_Sew(void)
void setCalibSvc_init(void)
initialize for the services
void smoother_anal(KalFitTrack &trk, int way)
Kalman filter (smoothing or backward part)
int i_back_
mass assumption for backward filter (if <0 means use leading mass)
void complete_track(MdcRec_trk &TrasanTRK, MdcRec_trk_add &TrasanTRK_add, KalFitTrack &track_lead, RecMdcKalTrack *kaltrk, RecMdcKalTrackCol *kalcol, RecMdcKalHelixSegCol *segcol, int flagsmooth)
int wsag_
flag to take account the wire sag into account
void kalman_fitting_calib(void)
void sameas(RecMdcKalTrack *trk, int l_mass, int imain)
complete the RecMdcKalTrackCol
void fillTds_ip(MdcRec_trk &TrasanTRK, KalFitTrack &track, RecMdcKalTrack *trk, int l_mass)
with results got at (0,0,0)
void fillTds(MdcRec_trk &TrasanTRK, KalFitTrack &track, RecMdcKalTrack *trk, int l_mass)
with results got at the inner Mdc hit
void innerwall(KalFitTrack &trk, int l_mass, int way)
Take the inner walls (eloss and mult scat) into account.
int debug_
Debug flag for the track finder part.
void start_seed(KalFitTrack &track, int lead_, int way, MdcRec_trk &trk)
void setBesFromGdml(void)
void setGeomSvc_init(void)
KalFitAlg(const std::string &name, ISvcLocator *pSvcLocator)
constructor
void kalman_fitting_anal(void)
int tprop_
propagation correction
void set_Mdc(void)
Set up the wires, layers and superlayers...
int enhance_
flag to enhance the error matrix at the inner hit of Mdc (cosmic)
StatusCode execute()
event function
int lead_
leading mass assumption
void smoother_calib(KalFitTrack &trk, int way)
void init_matrix(MdcRec_trk &trk, HepSymMatrix &Ea)
double pT_
value of the pT cut for backward filter
StatusCode initialize()
initialize
void kalman_fitting_csmalign(void)
int ntuple_
Fill ntuples of KalFit.
void hist_def(void)
hist definition
void fillTds_back(KalFitTrack &track, RecMdcKalTrack *trk, MdcRec_trk &TrasanTRK, int l_mass)
with results got at the outer Mdc hit
CLHEP::HepVector a_pre_fwd(void)
CLHEP::HepSymMatrix Ea_filt_fwd(void)
CLHEP::HepSymMatrix & Ea_pre_fwd(void)
KalFitHitMdc * HitMdc(void)
Description of a Hit in Mdc.
double chi2_back(void) const
const KalFitWire & wire(void) const
const int layerId(void) const
returns layer ID
Description of a track class (<- Helix.cc)
static void setMdcDigiCol(MdcDigiCol *digicol)
static int lead(void)
Magnetic field map.
static double chi2_hitf_
Cut chi2 for each hit.
static double dchi2cuts_anal[43]
unsigned int nchits(void) const
static double dchi2cuts_calib[43]
double tof_proton(void) const
static int debug_
for debug
double tof_kaon(void) const
void order_wirhit(int index)
HepSymMatrix getInitMatrix(void) const
static void setInitMatrix(HepSymMatrix m)
static int nmdc_hit2_
Cut chi2 for each hit.
KalFitHelixSeg & HelixSeg(int i)
void addTofSM(double time)
static double factor_strag_
factor of energy loss straggling for electron
const HepPoint3D & pivot_numf(const HepPoint3D &newPivot)
Sets pivot position in a given mag field.
static int Tof_correc_
Flag for TOF correction.
double chiSq_back(void) const
void HitsMdc(vector< KalFitHitMdc > &lh)
double pathip(void) const
void point_last(const HepPoint3D &point)
set and give out the last point of the track
void HelixSegs(vector< KalFitHelixSeg > &vs)
static int numf_
Flag for treatment of non-uniform mag field.
unsigned int nster(void) const
double p_proton(void) const
double update_hits_csmalign(KalFitHelixSeg &HelixSeg, int inext, CLHEP::Hep3Vector &meas, int way, double &dchi2, int csmflag)
void tof(double path)
Update the tof estimation.
static int drifttime_choice_
the drifttime choice
double smoother_Mdc(KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag)
Kalman smoother for Mdc.
void msgasmdc(double path, int index)
Calculate multiple scattering angle.
static int tprop_
for signal propagation correction
void appendHitsMdc(KalFitHitMdc h)
Functions for Mdc hits list.
void addPathSM(double path)
static double dchi2cutf_calib[43]
void eloss(double path, const KalFitMaterial &m, int index)
Calculate total energy lost in material.
static void setT0(double t0)
static void setMagneticFieldSvc(IMagneticFieldSvc *)
static int numfcor_
NUMF treatment improved.
double radius_numf(void) const
Estimation of the radius in a given mag field.
double smoother_Mdc_csmalign(KalFitHelixSeg &seg, CLHEP::Hep3Vector &meas, int &flg, int csmflag)
double update_hits(KalFitHitMdc &HitMdc, int inext, CLHEP::Hep3Vector &meas, int way, double &dchi2, double &dtrack, double &dtracknew, double &dtdc, int csmflag)
Include the Mdc wire hits.
static double dchi2cutf_anal[43]
static int LR_
Use L/R decision from MdcRecHit information :
double p_kaon(void) const
KalFitHitMdc & HitMdc(int i)
static int strag_
Flag to take account of energy loss straggling :
Description of a Wire class.
unsigned int geoID(void) const
unsigned int localId(void) const
Extractor :
HepPoint3D bck(void) const
HepPoint3D fwd(void) const
Geometry :
const KalFitLayer_Mdc & layer(void) const
unsigned int stereo(void) const
double bFieldZ(double)
sets/returns z componet of the magnetic field.
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
void ignoreErrorMatrix(void)
unsets error matrix. Error calculations will be ignored after this function call until an error matri...
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
const HepSymMatrix & Ea(void) const
returns error matrix.
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
double radius(void) const
returns radious of helix.
const HepVector & a(void) const
returns helix parameters.
const HepPoint3D & pivot(void) const
returns pivot position.
MdcGeoLayer * Lyr(void) const
double Tension(void) const
const MdcGeoWire *const Wire(unsigned id)
const MdcGeoLayer *const Layer(unsigned id)
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
static vector< MdcRec_trk_add > * getMdcRecTrkAddCol(void)
static vector< MdcRec_trk > * getMdcRecTrkCol(void)
static vector< MdcRec_wirhit > * getMdcRecWirhitCol(void)
vector< MdcRec_wirhit * > hitcol
void setTrackId(int trackid)
void setResIncl(double res)
void setDocaExcl(double doca)
HepVector & getHelixExcl(void)
void setHelixIncl(const HepVector &helix)
void setLayerId(int layerId)
int getFlagLR(void) const
void setDrIncl(double dr)
HepVector & getHelixIncl(void)
void setTanlIncl(double tanl)
void setDzIncl(double dz)
void setEntra(double entra)
void setDzExcl(double dz)
void setFi0Excl(double fi0)
void setFi0Incl(double fi0)
void setFlagLR(int flagLR)
HepSymMatrix & getErrorExcl(void)
void setCpaIncl(double cpa)
void setResExcl(double res)
void setCpaExcl(double cpa)
void setDrExcl(double dr)
void setHelixExcl(const HepVector &helix)
void setZhit(double zhit)
HepSymMatrix & getErrorIncl(void)
void setErrorExcl(const HepSymMatrix &error)
void setTanlExcl(double tanl)
void setMdcId(Identifier mdcid)
void setDocaIncl(double doca)
void setErrorIncl(const HepSymMatrix &error)
const HepVector & getZHelix() const
void setLHelix(const HepVector &helix, const int pid)
void setFiTerm(double fi, const int pid)
double getPathSM(int pid) const
void setPathSM(double length, int pid)
const HepPoint3D & getLPoint() const
const double getFiTerm(const int pid) const
double getTof(int pid) const
const HepSymMatrix & getFError() const
double getChisq(int i, int pid) const
void setLength(double length, int pid)
void setTError(const HepSymMatrix &error)
void setMass(double mass, int pid)
HelixSegRefVec getVecHelixSegs(int pid=-1) const
const HepVector & getFHelix() const
void setTHelix(const HepVector &helix)
void setPathl(double pathl, int i)
const HepPoint3D & getLPivot() const
int getTrackId(void) const
void setVecHelixSegs(const HelixSegRefVec &vechelixsegs, int pid=-1)
void setLError(const HepSymMatrix &error, const int pid)
int getNhits(int pid) const
double getMass(int pid) const
void setLPoint(const HepPoint3D &point, const int pid)
int getStat(int i, int pid) const
void setLPivot(const HepPoint3D &pivot, const int pid)
double getLength(int pid) const
void setTof(double tof, int pid)
int getNdf(int i, int pid) const
void setNhits(int nhits, int pid)