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"
40#include "GaudiKernel/RndmGenerators.h"
71#include "GaudiKernel/IPartPropSvc.h"
72#include "GaudiKernel/INTupleSvc.h"
79using CLHEP::HepVector;
80using CLHEP::Hep3Vector;
81using CLHEP::HepMatrix;
82using CLHEP::HepSymMatrix;
86const double newr[3]={7.8338,12.3604,16.6104};
87const double dmr_layer[3] = {79.838,125.104,167.604};
88const double dor_layer[3] = {81.338,126.604,169.104};
89const double dr_layer[3] = {78.338,123.604,166.104};
94const double KalFitAlg::RIW = 6.35;
100 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc_(0),m_MFSvc_(0),
101 _wire(0), _layer(0), _superLayer(0),
102 pathl_(1), wsag_(4), back_(1), pT_(0.0), lead_(2), mhyp_(31),
103 pe_cut_(2.0), pmu_cut_(2.0), ppi_cut_(2.0), pk_cut_(2.0), pp_cut_(2.0),
104 muls_(1), loss_(1), enhance_(0),drifttime_choice_(0),choice_(0),
105 fac_h1_(1),fac_h2_(1),fac_h3_(1),fac_h4_(1),fac_h5_(1),
106 i_back_(-1), debug_kft_(0), debug_(0), ntuple_(0),eventno(-1),
107 Tds_back_no(0),m_nt1(0),m_nt2(0),m_nt3(0),m_nt4(0),m_nt5(0),
108 iqual_back_(1),tprop_(1),
109 dchi2cut_inner_(0),dchi2cut_outer_(0),
110 dchi2cut_mid1_(0),dchi2cut_mid2_(0),
111 dchi2cut_layid2_(0),dchi2cut_layid3_(0),m_usevtxdb(0),m_dangcut(10),m_dphicut(10),
120 declareProperty(
"gain1",
gain1_ = 1.);
121 declareProperty(
"gain2",
gain2_ = 1.);
122 declareProperty(
"gain3",
gain3_ = 1.);
123 declareProperty(
"gain4",
gain4_ = 1.);
124 declareProperty(
"gain5",
gain5_ = 1.);
125 declareProperty(
"fitnocut",
fitnocut_ = 5);
128 declareProperty(
"choice",
choice_ = 0);
129 declareProperty(
"numfcor",
numfcor_ = 0);
130 declareProperty(
"numf",
numf_ = 0);
131 declareProperty(
"steplev",
steplev_ = 2);
132 declareProperty(
"usage",
usage_=0);
133 declareProperty(
"i_front",
i_front_=1);
134 declareProperty(
"lead",
lead_=2);
135 declareProperty(
"muls",
muls_=1);
136 declareProperty(
"loss",
loss_=1);
137 declareProperty(
"matrixg",
matrixg_=100.0);
138 declareProperty(
"lr",
lr_=1);
140 declareProperty(
"debug",
debug_=0);
141 declareProperty(
"ntuple",
ntuple_=0);
143 declareProperty(
"matfile",
matfile_=
"geomdc_material.dat");
144 declareProperty(
"cylfile",
cylfile_=
"geomdc_cylinder.dat");
145 declareProperty(
"dchi2cutf",
dchi2cutf_=1000.0);
146 declareProperty(
"dchi2cuts",
dchi2cuts_=1000.0);
147 declareProperty(
"pt",
pT_=0.0);
148 declareProperty(
"pe_cut",
pe_cut_=2.0);
149 declareProperty(
"pmu_cut",
pmu_cut_=2.0);
150 declareProperty(
"ppi_cut",
ppi_cut_=2.0);
151 declareProperty(
"pk_cut",
pk_cut_=2.0);
152 declareProperty(
"pp_cut",
pp_cut_=2.0);
153 declareProperty(
"wsag",
wsag_=4);
154 declareProperty(
"back",
back_=1);
155 declareProperty(
"i_back",
i_back_=1);
156 declareProperty(
"tofflag",
tofflag_=1);
157 declareProperty(
"tof_hyp",
tof_hyp_=1);
159 declareProperty(
"fstrag",
fstrag_=0.9);
161 declareProperty(
"tprop",
tprop_=1);
162 declareProperty(
"pt_cut",
pt_cut_= 0.2);
165 declareProperty(
"cosmicflag",
m_csmflag= 0);
166 declareProperty(
"dangcut",
m_dangcut=10.);
167 declareProperty(
"dphicut",
m_dphicut=10.);
168 declareProperty(
"useNCGem",useNCGem_=0);
169 declareProperty(
"nInnerCFiber",nInnerCFiber_=0);
170 declareProperty(
"SigRphiGem",myRphiResGem=0.033);
171 declareProperty(
"SigZGem", myZResGem=0.04);
172 declareProperty(
"UseGemPos3D", myUseGemPos3D=
false);
173 declareProperty(
"testCDC", testCDC=0);
174 declareProperty(
"ifProdNt10", ifProdNt10=
false);
175 declareProperty(
"ifProdNt11", ifProdNt11=
false);
176 declareProperty(
"ifProdNt12", ifProdNt12=
false);
178 for(
int i=0; i<5; i++) nFailedTrks[i]=0;
184 MsgStream log(
msgSvc(), name());
185 log << MSG::INFO <<
" Start cleaning, delete Mdc geometry objects" << endreq;
187 log << MSG::INFO <<
"End cleaning " << endreq;
193 _BesKalmanFitWalls.clear();
194 _BesKalmanFitMaterials.clear();
197 m_CGEMMaterials.clear();
198 if(_wire)free(_wire);
199 if(_layer)free(_layer);
200 if(_superLayer)free(_superLayer);
206 MsgStream log(
msgSvc(), name());
207 log << MSG::INFO <<
"in initize()"
208 <<
"KalFit> Initialization for current run " << endreq;
209 log << MSG::INFO <<
"Present Parameters: muls: " <<
muls_ <<
" loss: "
232 log << MSG::INFO <<
".....building Mdc " << endreq;
252 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
253 if(sc!=StatusCode::SUCCESS) {
254 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
264 std::cout<<
" initialize, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
272 log << MSG::INFO <<
" ntuple out, the option is "<<
ntuple_ <<endreq;
274 cout <<
"KalFitAlg> DEBUG open,Here is the important Parameters :\n";
275 cout <<
" Leading particule with mass hyp = " <<
lead_ << std::endl;
276 cout <<
" mhyp = " <<
mhyp_ << std::endl;
277 cout <<
"===== Effects taking into account : " << std::endl;
278 cout <<
" - multiple scattering = " <<
muls_ << std::endl;
279 cout <<
" - energy loss = " <<
loss_ << std::endl;
281 cout <<
" - straggling for the energy loss " << std::endl;
285 cout <<
" - non uniform magnetic field treatment "
287 cout <<
" - wire sag correction = " <<
wsag_ << std::endl;
291 <<
" hits included " << std::endl;
294 cout <<
" Backward filter is on with a pT cut value = " <<
pT_ << endl;
296 if(
debug_ == 4) cout <<
" pathl = " <<
pathl_ << std::endl;
299 cout <<
" Decision L/R from MdcRecHit " << std::endl;
308 IPartPropSvc* p_PartPropSvc;
309 static const bool CREATEIFNOTTHERE(
true);
310 StatusCode PartPropStatus = service(
"PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
311 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
312 log << MSG::WARNING <<
" Could not initialize Particle Properties Service" << endreq;
313 return StatusCode::SUCCESS;
315 m_particleTable = p_PartPropSvc->PDT();
317 return StatusCode::SUCCESS;
322 MsgStream log(
msgSvc(), name());
323 log << MSG::DEBUG<<
"KalFitAlg:: nTotalTrks = "<<nTotalTrks<<endreq;
324 log << MSG::DEBUG<<
" e: "<<nFailedTrks[0]<<
" failed, "<<nTotalTrks-nFailedTrks[0]<<
" successed"<<endreq;
325 log << MSG::DEBUG<<
" mu: "<<nFailedTrks[1]<<
" failed, "<<nTotalTrks-nFailedTrks[1]<<
" successed"<<endreq;
326 log << MSG::DEBUG<<
" pi: "<<nFailedTrks[2]<<
" failed, "<<nTotalTrks-nFailedTrks[2]<<
" successed"<<endreq;
327 log << MSG::DEBUG<<
" K: "<<nFailedTrks[3]<<
" failed, "<<nTotalTrks-nFailedTrks[3]<<
" successed"<<endreq;
328 log << MSG::DEBUG<<
" p: "<<nFailedTrks[4]<<
" failed, "<<nTotalTrks-nFailedTrks[4]<<
" successed"<<endreq;
331 cout<<
"KalFitAlg:: nTotalTrks = "<<nTotalTrks<<endl;
332 cout<<
" e: "<<nFailedTrks[0]<<
" failed, "<<nTotalTrks-nFailedTrks[0]<<
" successed"<<endl;
333 cout<<
" mu: "<<nFailedTrks[1]<<
" failed, "<<nTotalTrks-nFailedTrks[1]<<
" successed"<<endl;
334 cout<<
" pi: "<<nFailedTrks[2]<<
" failed, "<<nTotalTrks-nFailedTrks[2]<<
" successed"<<endl;
335 cout<<
" K: "<<nFailedTrks[3]<<
" failed, "<<nTotalTrks-nFailedTrks[3]<<
" successed"<<endl;
336 cout<<
" p: "<<nFailedTrks[4]<<
" failed, "<<nTotalTrks-nFailedTrks[4]<<
" successed"<<endl;
338 return StatusCode::SUCCESS;
344 MsgStream log(
msgSvc(), name());
345 log << MSG::INFO <<
"in beginRun()" << endreq;
346 log << MSG::INFO <<
"Present Parameters: muls: " <<
muls_ <<
" loss: "
355 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
356 if(sc!=StatusCode::SUCCESS) {
357 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
366 std::cout<<
" beginRun, referField from MagneticFieldSvc:"<< (IMFSvc->
getReferField())*10000 <<std::endl;
371 return StatusCode::SUCCESS;
381 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101");
383 if ( nt1 ) m_nt1 = nt1;
385 m_nt1=
ntupleSvc()->book(
"FILE104/n101",CLID_ColumnWiseTuple,
"KalFit");
388 status = m_nt1->addItem(
"trackid",m_trackid);
389 status = m_nt1->addItem(
"stat",5,2,m_stat);
390 status = m_nt1->addItem(
"ndf",5,2,m_ndf);
391 status = m_nt1->addItem(
"chisq",5,2,m_chisq);
392 status = m_nt1->addItem(
"length",5,m_length);
393 status = m_nt1->addItem(
"tof",5,m_tof);
394 status = m_nt1->addItem(
"nhits",5,m_nhits);
395 status = m_nt1->addItem(
"zhelix",5,m_zhelix);
396 status = m_nt1->addItem(
"zhelixe",5,m_zhelixe);
397 status = m_nt1->addItem(
"zhelixmu",5,m_zhelixmu);
398 status = m_nt1->addItem(
"zhelixk",5,m_zhelixk);
399 status = m_nt1->addItem(
"zhelixp",5,m_zhelixp);
400 status = m_nt1->addItem(
"zptot",m_zptot);
401 status = m_nt1->addItem(
"zptote",m_zptote);
402 status = m_nt1->addItem(
"zptotmu",m_zptotmu);
403 status = m_nt1->addItem(
"zptotk",m_zptotk);
404 status = m_nt1->addItem(
"zptotp",m_zptotp);
406 status = m_nt1->addItem(
"zpt",m_zpt);
407 status = m_nt1->addItem(
"zpte",m_zpte);
408 status = m_nt1->addItem(
"zptmu",m_zptmu);
409 status = m_nt1->addItem(
"zptk",m_zptk);
410 status = m_nt1->addItem(
"zptp",m_zptp);
412 status = m_nt1->addItem(
"fptot",m_fptot);
413 status = m_nt1->addItem(
"fptote",m_fptote);
414 status = m_nt1->addItem(
"fptotmu",m_fptotmu);
415 status = m_nt1->addItem(
"fptotk",m_fptotk);
416 status = m_nt1->addItem(
"fptotp",m_fptotp);
417 status = m_nt1->addItem(
"fpt",m_fpt);
418 status = m_nt1->addItem(
"fpte",m_fpte);
419 status = m_nt1->addItem(
"fptmu",m_fptmu);
420 status = m_nt1->addItem(
"fptk",m_fptk);
421 status = m_nt1->addItem(
"fptp",m_fptp);
423 status = m_nt1->addItem(
"lptot",m_lptot);
424 status = m_nt1->addItem(
"lptote",m_lptote);
425 status = m_nt1->addItem(
"lptotmu",m_lptotmu);
426 status = m_nt1->addItem(
"lptotk",m_lptotk);
427 status = m_nt1->addItem(
"lptotp",m_lptotp);
428 status = m_nt1->addItem(
"lpt",m_lpt);
429 status = m_nt1->addItem(
"lpte",m_lpte);
430 status = m_nt1->addItem(
"lptmu",m_lptmu);
431 status = m_nt1->addItem(
"lptk",m_lptk);
432 status = m_nt1->addItem(
"lptp",m_lptp);
434 status = m_nt1->addItem(
"zsigp",m_zsigp);
435 status = m_nt1->addItem(
"zsigpe",m_zsigpe);
436 status = m_nt1->addItem(
"zsigpmu",m_zsigpmu);
437 status = m_nt1->addItem(
"zsigpk",m_zsigpk);
438 status = m_nt1->addItem(
"zsigpp",m_zsigpp);
439 status = m_nt1->addItem(
"fhelix",5,m_fhelix);
440 status = m_nt1->addItem(
"fhelixe",5,m_fhelixe);
441 status = m_nt1->addItem(
"fhelixmu",5,m_fhelixmu);
442 status = m_nt1->addItem(
"fhelixk",5,m_fhelixk);
443 status = m_nt1->addItem(
"fhelixp",5,m_fhelixp);
444 status = m_nt1->addItem(
"lhelix",5,m_lhelix);
445 status = m_nt1->addItem(
"lhelixe",5,m_lhelixe);
446 status = m_nt1->addItem(
"lhelixmu",5,m_lhelixmu);
447 status = m_nt1->addItem(
"lhelixk",5,m_lhelixk);
448 status = m_nt1->addItem(
"lhelixp",5,m_lhelixp);
449 status = m_nt1->addItem(
"rGem",4,m_rGem);
450 status = m_nt1->addItem(
"chi2Gem",4,m_chi2Gem);
451 status = m_nt1->addItem(
"phiGemExp",4,m_phiGemExp);
452 status = m_nt1->addItem(
"phiGemHit",4,m_phiGemHit);
453 status = m_nt1->addItem(
"zGemExp",4,m_zGemExp);
454 status = m_nt1->addItem(
"zGemHit",4,m_zGemHit);
456 status = m_nt1->addItem(
"zerror",15,m_zerror);
457 status = m_nt1->addItem(
"zerrore",15,m_zerrore);
458 status = m_nt1->addItem(
"zerrormu",15,m_zerrormu);
459 status = m_nt1->addItem(
"zerrork",15,m_zerrork);
460 status = m_nt1->addItem(
"zerrorp",15,m_zerrorp);
461 status = m_nt1->addItem(
"ferror",15,m_ferror);
462 status = m_nt1->addItem(
"ferrore",15,m_ferrore);
463 status = m_nt1->addItem(
"ferrormu",15,m_ferrormu);
464 status = m_nt1->addItem(
"ferrork",15,m_ferrork);
465 status = m_nt1->addItem(
"ferrorp",15,m_ferrorp);
466 status = m_nt1->addItem(
"lerror",15,m_lerror);
467 status = m_nt1->addItem(
"lerrore",15,m_lerrore);
468 status = m_nt1->addItem(
"lerrormu",15,m_lerrormu);
469 status = m_nt1->addItem(
"lerrork",15,m_lerrork);
470 status = m_nt1->addItem(
"lerrorp",15,m_lerrorp);
473 status = m_nt1->addItem(
"evtid",m_evtid);
474 status = m_nt1->addItem(
"mchelix",5,m_mchelix);
475 status = m_nt1->addItem(
"mcptot",m_mcptot);
476 status = m_nt1->addItem(
"mcpid",m_mcpid);
478 if( status.isFailure() ) cout<<
"Ntuple1 add item failed!"<<endl;
484 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102");
486 if ( nt2 ) m_nt2 = nt2;
488 m_nt2=
ntupleSvc()->book(
"FILE104/n102",CLID_ColumnWiseTuple,
"KalFitComp");
490 status2 = m_nt2->addItem(
"delx",m_delx);
491 status2 = m_nt2->addItem(
"dely",m_dely);
492 status2 = m_nt2->addItem(
"delz",m_delz);
493 status2 = m_nt2->addItem(
"delthe",m_delthe);
494 status2 = m_nt2->addItem(
"delphi",m_delphi);
495 status2 = m_nt2->addItem(
"delp",m_delp);
496 status2 = m_nt2->addItem(
"delpx",m_delpx);
497 status2 = m_nt2->addItem(
"delpy",m_delpy);
498 status2 = m_nt2->addItem(
"delpz",m_delpz);
500 if( status2.isFailure() ) cout<<
"Ntuple2 add item failed!"<<endl;
506 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103");
508 if ( nt3 ) m_nt3 = nt3;
510 m_nt3=
ntupleSvc()->book(
"FILE104/n103",CLID_ColumnWiseTuple,
"PatRec");
512 status3 = m_nt3->addItem(
"trkhelix",5,m_trkhelix);
513 status3 = m_nt3->addItem(
"trkptot",m_trkptot);
515 status3 = m_nt3->addItem(
"trkerror",15,m_trkerror);
516 status3 = m_nt3->addItem(
"trksigp",m_trksigp);
518 status3 = m_nt3->addItem(
"trkndf",m_trkndf);
519 status3 = m_nt3->addItem(
"trkchisq",m_trkchisq);
520 if( status3.isFailure() ) cout<<
"Ntuple3 add item failed!"<<endl;
526 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104");
528 if ( nt4 ) m_nt4 = nt4;
530 m_nt4=
ntupleSvc()->book(
"FILE104/n104",CLID_ColumnWiseTuple,
"PatRecComp");
532 status4 = m_nt4->addItem(
"trkdelx",m_trkdelx);
533 status4 = m_nt4->addItem(
"trkdely",m_trkdely);
534 status4 = m_nt4->addItem(
"trkdelz",m_trkdelz);
535 status4 = m_nt4->addItem(
"trkdelthe",m_trkdelthe);
536 status4 = m_nt4->addItem(
"trkdelphi",m_trkdelphi);
537 status4 = m_nt4->addItem(
"trkdelp",m_trkdelp);
538 if( status4.isFailure() ) cout<<
"Ntuple4 add item failed!"<<endl;
543 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105");
545 if ( nt5 ) m_nt5 = nt5;
547 m_nt5=
ntupleSvc()->book(
"FILE104/n105",CLID_ColumnWiseTuple,
"KalFitdChisq");
549 status5 = m_nt5->addItem(
"dchi2",m_dchi2);
550 status5 = m_nt5->addItem(
"masshyp",m_masshyp);
551 status5 = m_nt5->addItem(
"residual_estim",m_residest);
552 status5 = m_nt5->addItem(
"residual",m_residnew);
553 status5 = m_nt5->addItem(
"layer",m_layer);
554 status5 = m_nt5->addItem(
"kaldr",m_anal_dr);
555 status5 = m_nt5->addItem(
"kalphi0",m_anal_phi0);
556 status5 = m_nt5->addItem(
"kalkappa",m_anal_kappa);
557 status5 = m_nt5->addItem(
"kaldz",m_anal_dz);
558 status5 = m_nt5->addItem(
"kaltanl",m_anal_tanl);
559 status5 = m_nt5->addItem(
"dr_ea",m_anal_ea_dr);
560 status5 = m_nt5->addItem(
"phi0_ea",m_anal_ea_phi0);
561 status5 = m_nt5->addItem(
"kappa_ea",m_anal_ea_kappa);
562 status5 = m_nt5->addItem(
"dz_ea",m_anal_ea_dz);
563 status5 = m_nt5->addItem(
"tanl_ea",m_anal_ea_tanl);
564 if( status5.isFailure() ) cout<<
"Ntuple5 add item failed!"<<endl;
571 NTuplePtr nt7(
ntupleSvc(),
"FILE104/n120");
575 m_nt7 =
ntupleSvc()-> book(
"FILE104/n120",CLID_ColumnWiseTuple,
"kalfit_failure");
577 status7 = m_nt7 ->addItem(
"run_kal",m_run_kal);
578 status7 = m_nt7 ->addItem(
"event_kal",m_event_kal);
579 status7 = m_nt7 ->addItem(
"trkid_kal",m_trkid_kal);
580 status7 = m_nt7 ->addItem(
"dropedHits_kal_e",m_dropedHits_kal_e);
581 status7 = m_nt7 ->addItem(
"kappa2_kal_e",m_kappa2_kal_e);
582 status7 = m_nt7 ->addItem(
"trackNhits_kal_e",m_trackNhits_kal_e);
583 status7 = m_nt7 ->addItem(
"trackNster_kal_e",m_trackNster_kal_e);
584 status7 = m_nt7 ->addItem(
"trackNaxis_kal_e",m_trackNaxis_kal_e);
585 status7 = m_nt7 ->addItem(
"chi2_kal_e",m_chi2_kal_e);
586 status7 = m_nt7 ->addItem(
"Ea00_kal_e",m_Ea00_kal_e);
587 status7 = m_nt7 ->addItem(
"Ea11_kal_e",m_Ea11_kal_e);
588 status7 = m_nt7 ->addItem(
"Ea22_kal_e",m_Ea22_kal_e);
589 status7 = m_nt7 ->addItem(
"Ea33_kal_e",m_Ea33_kal_e);
590 status7 = m_nt7 ->addItem(
"Ea44_kal_e",m_Ea44_kal_e);
591 status7 = m_nt7 ->addItem(
"dropedHits_kal_mu",m_dropedHits_kal_mu);
592 status7 = m_nt7 ->addItem(
"kappa2_kal_mu",m_kappa2_kal_mu);
593 status7 = m_nt7 ->addItem(
"trackNhits_kal_mu",m_trackNhits_kal_mu);
594 status7 = m_nt7 ->addItem(
"trackNster_kal_mu",m_trackNster_kal_mu);
595 status7 = m_nt7 ->addItem(
"trackNaxis_kal_mu",m_trackNaxis_kal_mu);
596 status7 = m_nt7 ->addItem(
"chi2_kal_mu",m_chi2_kal_mu);
597 status7 = m_nt7 ->addItem(
"Ea00_kal_mu",m_Ea00_kal_mu);
598 status7 = m_nt7 ->addItem(
"Ea11_kal_mu",m_Ea11_kal_mu);
599 status7 = m_nt7 ->addItem(
"Ea22_kal_mu",m_Ea22_kal_mu);
600 status7 = m_nt7 ->addItem(
"Ea33_kal_mu",m_Ea33_kal_mu);
601 status7 = m_nt7 ->addItem(
"Ea44_kal_mu",m_Ea44_kal_mu);
602 status7 = m_nt7 ->addItem(
"iqual_front_kal_mu",m_iqual_front_kal_mu);
603 status7 = m_nt7 ->addItem(
"dropedHits_kal_pi",m_dropedHits_kal_pi);
604 status7 = m_nt7 ->addItem(
"kappa2_kal_pi",m_kappa2_kal_pi);
605 status7 = m_nt7 ->addItem(
"trackNhits_kal_pi",m_trackNhits_kal_pi);
606 status7 = m_nt7 ->addItem(
"trackNster_kal_pi",m_trackNster_kal_pi);
607 status7 = m_nt7 ->addItem(
"trackNaxis_kal_pi",m_trackNaxis_kal_pi);
608 status7 = m_nt7 ->addItem(
"chi2_kal_pi",m_chi2_kal_pi);
609 status7 = m_nt7 ->addItem(
"Ea00_kal_pi",m_Ea00_kal_pi);
610 status7 = m_nt7 ->addItem(
"Ea11_kal_pi",m_Ea11_kal_pi);
611 status7 = m_nt7 ->addItem(
"Ea22_kal_pi",m_Ea22_kal_pi);
612 status7 = m_nt7 ->addItem(
"Ea33_kal_pi",m_Ea33_kal_pi);
613 status7 = m_nt7 ->addItem(
"Ea44_kal_pi",m_Ea44_kal_pi);
614 status7 = m_nt7 ->addItem(
"dropedHits_kal_k",m_dropedHits_kal_k);
615 status7 = m_nt7 ->addItem(
"kappa2_kal_k",m_kappa2_kal_k);
616 status7 = m_nt7 ->addItem(
"trackNhits_kal_k",m_trackNhits_kal_k);
617 status7 = m_nt7 ->addItem(
"trackNster_kal_k",m_trackNster_kal_k);
618 status7 = m_nt7 ->addItem(
"trackNaxis_kal_k",m_trackNaxis_kal_k);
619 status7 = m_nt7 ->addItem(
"chi2_kal_k",m_chi2_kal_k);
620 status7 = m_nt7 ->addItem(
"Ea00_kal_k",m_Ea00_kal_k);
621 status7 = m_nt7 ->addItem(
"Ea11_kal_k",m_Ea11_kal_k);
622 status7 = m_nt7 ->addItem(
"Ea22_kal_k",m_Ea22_kal_k);
623 status7 = m_nt7 ->addItem(
"Ea33_kal_k",m_Ea33_kal_k);
624 status7 = m_nt7 ->addItem(
"Ea44_kal_k",m_Ea44_kal_k);
625 status7 = m_nt7 ->addItem(
"dropedHits_kal_p",m_dropedHits_kal_p);
626 status7 = m_nt7 ->addItem(
"kappa2_kal_p",m_kappa2_kal_p);
627 status7 = m_nt7 ->addItem(
"trackNhits_kal_p",m_trackNhits_kal_p);
628 status7 = m_nt7 ->addItem(
"trackNster_kal_p",m_trackNster_kal_p);
629 status7 = m_nt7 ->addItem(
"trackNaxis_kal_p",m_trackNaxis_kal_p);
630 status7 = m_nt7 ->addItem(
"chi2_kal_p",m_chi2_kal_p);
631 status7 = m_nt7 ->addItem(
"Ea00_kal_p",m_Ea00_kal_p);
632 status7 = m_nt7 ->addItem(
"Ea11_kal_p",m_Ea11_kal_p);
633 status7 = m_nt7 ->addItem(
"Ea22_kal_p",m_Ea22_kal_p);
634 status7 = m_nt7 ->addItem(
"Ea33_kal_p",m_Ea33_kal_p);
635 status7 = m_nt7 ->addItem(
"Ea44_kal_p",m_Ea44_kal_p);
637 status7 = m_nt7 ->addItem(
"hit_number",m_hit_no,0, 1000);
638 status7 = m_nt7 ->addItem(
"nCluster",m_nCluster,0, 1000);
639 status7 = m_nt7->addIndexedItem(
"dchi2_hit_e",m_hit_no,m_dchi2_hit_e);
640 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_e",m_hit_no,m_residest_hit_e);
641 status7 = m_nt7->addIndexedItem(
"residual_hit_e",m_hit_no,m_residnew_hit_e);
642 status7 = m_nt7->addIndexedItem(
"layer_hit_e",m_hit_no,m_layer_hit_e);
643 status7 = m_nt7->addIndexedItem(
"kaldr_hit_e",m_hit_no,m_anal_dr_hit_e);
644 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_e",m_hit_no,m_anal_phi0_hit_e);
645 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_e",m_hit_no,m_anal_kappa_hit_e);
646 status7 = m_nt7->addIndexedItem(
"kaldz_hit_e",m_hit_no,m_anal_dz_hit_e);
647 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_e",m_hit_no,m_anal_tanl_hit_e);
648 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_e",m_hit_no,m_anal_ea_dr_hit_e);
649 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_e",m_hit_no,m_anal_ea_phi0_hit_e);
650 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_e",m_hit_no,m_anal_ea_kappa_hit_e);
651 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_e",m_hit_no,m_anal_ea_dz_hit_e);
652 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_e",m_hit_no,m_anal_ea_tanl_hit_e);
653 status7 = m_nt7->addIndexedItem(
"dchi2_hit_mu",m_hit_no,m_dchi2_hit_mu);
654 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_mu",m_hit_no,m_residest_hit_mu);
655 status7 = m_nt7->addIndexedItem(
"residual_hit_mu",m_hit_no,m_residnew_hit_mu);
656 status7 = m_nt7->addIndexedItem(
"layer_hit_mu",m_hit_no,m_layer_hit_mu);
657 status7 = m_nt7->addIndexedItem(
"kaldr_hit_mu",m_hit_no,m_anal_dr_hit_mu);
658 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_mu",m_hit_no,m_anal_phi0_hit_mu);
659 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_mu",m_hit_no,m_anal_kappa_hit_mu);
660 status7 = m_nt7->addIndexedItem(
"kaldz_hit_mu",m_hit_no,m_anal_dz_hit_mu);
661 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_mu",m_hit_no,m_anal_tanl_hit_mu);
662 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_mu",m_hit_no,m_anal_ea_dr_hit_mu);
663 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_mu",m_hit_no,m_anal_ea_phi0_hit_mu);
664 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_mu",m_hit_no,m_anal_ea_kappa_hit_mu);
665 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_mu",m_hit_no,m_anal_ea_dz_hit_mu);
666 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_mu",m_hit_no,m_anal_ea_tanl_hit_mu);
667 status7 = m_nt7->addIndexedItem(
"dchi2_hit_pi",m_hit_no,m_dchi2_hit_pi);
668 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_pi",m_hit_no,m_residest_hit_pi);
669 status7 = m_nt7->addIndexedItem(
"residual_hit_pi",m_hit_no,m_residnew_hit_pi);
670 status7 = m_nt7->addIndexedItem(
"layer_hit_pi",m_hit_no,m_layer_hit_pi);
671 status7 = m_nt7->addIndexedItem(
"kaldr_hit_pi",m_hit_no,m_anal_dr_hit_pi);
672 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_pi",m_hit_no,m_anal_phi0_hit_pi);
673 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_pi",m_hit_no,m_anal_kappa_hit_pi);
674 status7 = m_nt7->addIndexedItem(
"kaldz_hit_pi",m_hit_no,m_anal_dz_hit_pi);
675 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_pi",m_hit_no,m_anal_tanl_hit_pi);
676 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_pi",m_hit_no,m_anal_ea_dr_hit_pi);
677 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_pi",m_hit_no,m_anal_ea_phi0_hit_pi);
678 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_pi",m_hit_no,m_anal_ea_kappa_hit_pi);
679 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_pi",m_hit_no,m_anal_ea_dz_hit_pi);
680 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_pi",m_hit_no,m_anal_ea_tanl_hit_pi);
681 status7 = m_nt7->addIndexedItem(
"dchi2_hit_k",m_hit_no,m_dchi2_hit_k);
682 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_k",m_hit_no,m_residest_hit_k);
683 status7 = m_nt7->addIndexedItem(
"residual_hit_k",m_hit_no,m_residnew_hit_k);
684 status7 = m_nt7->addIndexedItem(
"layer_hit_k",m_hit_no,m_layer_hit_k);
685 status7 = m_nt7->addIndexedItem(
"kaldr_hit_k",m_hit_no,m_anal_dr_hit_k);
686 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_k",m_hit_no,m_anal_phi0_hit_k);
687 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_k",m_hit_no,m_anal_kappa_hit_k);
688 status7 = m_nt7->addIndexedItem(
"kaldz_hit_k",m_hit_no,m_anal_dz_hit_k);
689 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_k",m_hit_no,m_anal_tanl_hit_k);
690 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_k",m_hit_no,m_anal_ea_dr_hit_k);
691 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_k",m_hit_no,m_anal_ea_phi0_hit_k);
692 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_k",m_hit_no,m_anal_ea_kappa_hit_k);
693 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_k",m_hit_no,m_anal_ea_dz_hit_k);
694 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_k",m_hit_no,m_anal_ea_tanl_hit_k);
695 status7 = m_nt7->addIndexedItem(
"dchi2_hit_p",m_hit_no,m_dchi2_hit_p);
696 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_p",m_hit_no,m_residest_hit_p);
697 status7 = m_nt7->addIndexedItem(
"residual_hit_p",m_hit_no,m_residnew_hit_p);
698 status7 = m_nt7->addIndexedItem(
"layer_hit_p",m_hit_no,m_layer_hit_p);
699 status7 = m_nt7->addIndexedItem(
"kaldr_hit_p",m_hit_no,m_anal_dr_hit_p);
700 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_p",m_hit_no,m_anal_phi0_hit_p);
701 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_p",m_hit_no,m_anal_kappa_hit_p);
702 status7 = m_nt7->addIndexedItem(
"kaldz_hit_p",m_hit_no,m_anal_dz_hit_p);
703 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_p",m_hit_no,m_anal_tanl_hit_p);
704 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_p",m_hit_no,m_anal_ea_dr_hit_p);
705 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_p",m_hit_no,m_anal_ea_phi0_hit_p);
706 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_p",m_hit_no,m_anal_ea_kappa_hit_p);
707 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_p",m_hit_no,m_anal_ea_dz_hit_p);
708 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_p",m_hit_no,m_anal_ea_tanl_hit_p);
710 if( status7.isFailure() ) cout<<
"Ntuple7 add item failed!"<<endl;
717 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
719 if ( nt6 ) m_nt6 = nt6;
721 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
723 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
724 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
725 status6 = m_nt6->addItem(
"tdr",m_tdrift);
726 status6 = m_nt6->addItem(
"layerid", m_layerid);
727 status6 = m_nt6->addItem(
"event", m_eventNo);
728 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
729 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
730 status6 = m_nt6->addItem(
"lr",m_lr);
731 status6 = m_nt6->addItem(
"dd",m_dd);
732 status6 = m_nt6->addItem(
"y",m_yposition);
734 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
739 NTuplePtr nt10(
ntupleSvc(),
"FILE104/n110");
741 if ( nt10 ) m_nt10 = nt10;
743 m_nt10=
ntupleSvc()->book(
"FILE104/n110",CLID_ColumnWiseTuple,
"test");
745 status10 = m_nt10->addItem(
"evt",m_evt3);
746 status10 = m_nt10->addItem(
"qua",m_qua);
747 status10 = m_nt10->addItem(
"nhit",m_nGemHits,0,4000000);
748 status10 = m_nt10->addIndexedItem(
"meas_r",m_nGemHits,m_meas_r);
749 status10 = m_nt10->addIndexedItem(
"meas_phi",m_nGemHits,m_meas_phi);
750 status10 = m_nt10->addIndexedItem(
"meas_z",m_nGemHits,m_meas_z);
751 status10 = m_nt10->addIndexedItem(
"esti1_r",m_nGemHits,m_esti1_r);
752 status10 = m_nt10->addIndexedItem(
"esti1_phi",m_nGemHits,m_esti1_phi);
753 status10 = m_nt10->addIndexedItem(
"esti1_z",m_nGemHits,m_esti1_z);
754 status10 = m_nt10->addIndexedItem(
"esti2_r",m_nGemHits,m_esti2_r);
755 status10 = m_nt10->addIndexedItem(
"esti2_phi",m_nGemHits,m_esti2_phi);
756 status10 = m_nt10->addIndexedItem(
"esti2_z",m_nGemHits,m_esti2_z);
757 status10 = m_nt10->addIndexedItem(
"diff1_phi",m_nGemHits,m_diff1_phi);
758 status10 = m_nt10->addIndexedItem(
"diff1_z",m_nGemHits,m_diff1_z);
759 status10 = m_nt10->addIndexedItem(
"diff2_phi",m_nGemHits,m_diff2_phi);
760 status10 = m_nt10->addIndexedItem(
"diff2_z",m_nGemHits,m_diff2_z);
761 status10 = m_nt10->addIndexedItem(
"layer",m_nGemHits,m_GemLayer);
762 status10 = m_nt10->addIndexedItem(
"mass",m_nGemHits,m_mass);
763 status10 = m_nt10->addIndexedItem(
"dchi2",m_nGemHits,m_Gchi2);
764 status10 = m_nt10->addIndexedItem(
"meas_phierr",m_nGemHits,m_meas_phierr);
765 status10 = m_nt10->addIndexedItem(
"meas_zerr",m_nGemHits,m_meas_zerr);
766 status10 = m_nt10->addIndexedItem(
"esti_phierr",m_nGemHits,m_esti_phierr);
767 status10 = m_nt10->addIndexedItem(
"esti_zerr",m_nGemHits,m_esti_zerr);
768 if( status10.isFailure() ) cout<<
"Ntuple10 add item failed!"<<endl;
774 NTuplePtr nt11(
ntupleSvc(),
"FILE104/n111");
776 if ( nt11 ) m_nt11 = nt11;
778 m_nt11=
ntupleSvc()->book(
"FILE104/n111",CLID_ColumnWiseTuple,
"truth");
780 status11 = m_nt11->addItem(
"evt",m_evt4);
781 status11 = m_nt11->addItem(
"ntruth",m_ntruth,0,400000);
782 status11 = m_nt11->addIndexedItem(
"phi",m_ntruth,m_dtphi);
783 status11 = m_nt11->addIndexedItem(
"z",m_ntruth,m_dtv);
784 status11 = m_nt11->addIndexedItem(
"postphi",m_ntruth,m_dtpostphi);
785 status11 = m_nt11->addIndexedItem(
"postz",m_ntruth,m_dtpostz);
786 status11 = m_nt11->addIndexedItem(
"layer",m_ntruth,m_tlayer);
787 if( status11.isFailure() ) cout<<
"Ntuple11 add item failed!"<<endl;
793 NTuplePtr nt12(
ntupleSvc(),
"FILE104/n112");
795 if ( nt12 ) m_nt12 = nt12;
797 m_nt12=
ntupleSvc()->book(
"FILE104/n112",CLID_ColumnWiseTuple,
"PatRecComp");
799 status12 = m_nt12->addItem(
"evt",m_evt);
800 status12 = m_nt12->addItem(
"track",m_track);
801 status12 = m_nt12->addItem(
"diff_dr",m_diff_dr);
802 status12 = m_nt12->addItem(
"diff_phi0",m_diff_phi0);
803 status12 = m_nt12->addItem(
"diff_kappa",m_diff_kappa);
804 status12 = m_nt12->addItem(
"diff_dz",m_diff_dz);
805 status12 = m_nt12->addItem(
"diff_tanl",m_diff_tanl);
806 status12 = m_nt12->addItem(
"diff_p",m_diff_p);
807 if( status12.isFailure() ) cout<<
"Ntuple12 add item failed!"<<endl;
820 for(layid = 0; layid<2; layid++) {
825 for(layid = 4; layid<12; layid++) {
828 for(layid = 12; layid<20; layid++) {
831 for(layid = 20; layid<43; layid++) {
837 for(layid = 0; layid<2; layid++) {
843 for(layid = 4; layid<12; layid++) {
846 for(layid = 12; layid<20; layid++) {
849 for(layid = 20; layid<43; layid++) {
860 for(layid = 0; layid<2; layid++) {
867 for(layid = 4; layid<12; layid++) {
871 for(layid = 12; layid<20; layid++) {
875 for(layid = 20; layid<43; layid++) {
881 for(layid = 0; layid<43; layid++) {
888 for(layid = 0; layid<2; layid++) {
895 for(layid = 4; layid<12; layid++) {
899 for(layid = 12; layid<20; layid++) {
903 for(layid = 20; layid<43; layid++) {
909 for(layid = 0; layid<43; layid++) {
920 MsgStream log(
msgSvc(), name());
921 log << MSG::INFO <<
"in execute()" << endreq;
922 if(
testOutput) std::cout<<
"begin to deal with EVENT ..."<<(++
eventno)<<std::endl;
968 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
969 if(sc!=StatusCode::SUCCESS) {
970 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
979 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
986 IDataProviderSvc* evtSvc = NULL;
987 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
989 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
991 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
992 return StatusCode::SUCCESS;
996 IDataManagerSvc *dataManSvc;
997 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
998 DataObject *aKalTrackCol;
999 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
1000 if(aKalTrackCol != NULL) {
1001 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
1002 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
1005 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
1006 if( kalsc.isFailure() ) {
1007 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
1008 return StatusCode::SUCCESS;
1010 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
1012 DataObject *aKalHelixSegCol;
1013 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
1014 if(aKalHelixSegCol != NULL){
1015 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
1016 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
1019 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
1020 if( kalsc.isFailure()){
1021 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
1022 return StatusCode::SUCCESS;
1024 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
1034 ISvcLocator* svcLocator = Gaudi::svcLocator();
1036 StatusCode Cgem_sc=svcLocator->service(
"CgemGeomSvc", ISvc);
1038 if (!Cgem_sc.isSuccess()) log<< MSG::INFO <<
"KalFitAlg::execute(): Could not open CGEM geometry file" << endreq;
1044 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
1047 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
1049 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
1050 return StatusCode::FAILURE;
1052 int eventNo = eventHeader->eventNumber();
1054 int runNo = eventHeader->runNumber();
1058 if(
testOutput) cout<<endl<<
"$$$$$$$$$$$ run="<<
runNo<<
", evt="<<
eventNo<<
" $$$$$$$$$$$$$$$$$"<<endl<<endl;
1061 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
1062 if (estimeCol && estimeCol->size()) {
1063 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
1064 t0 = (*iter_evt)->getTest();
1067 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
1068 return StatusCode::SUCCESS;
1073 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
1077 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
1078 if (sc!=StatusCode::SUCCESS) {
1079 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
1080 return StatusCode::SUCCESS;
1091 m_evtid = eventHeader->eventNumber();
1094 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
1096 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
1101 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
1102 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
1103 if(!(*i_mcTrk)->primaryParticle())
continue;
1104 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
1105 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
1106 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
1107 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
1108 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
1110 double charge = 0.0;
1111 int pid = (*i_mcTrk)->particleProperty();
1112 int mother_id = ((*i_mcTrk)->mother()).particleProperty();
1114 charge = m_particleTable->particle( pid )->charge();
1115 }
else if ( pid <0 ) {
1116 charge = m_particleTable->particle( -pid )->charge();
1119 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
1122 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
1125 log << MSG::DEBUG <<
"charge of the track " << charge << endreq;
1126 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
1128 for(
int j =0; j<5; j++) {
1129 m_mchelix[j] = mchelix.
a()[j];
1132 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
1135 cout<<
"MC pid, mother_id = "<<pid<<
", "<<mother_id<<endl;
1136 cout<<
" p4 = "<<(*i_mcTrk)->initialFourMomentum()<<endl;
1137 cout<<
" start from "<<(*i_mcTrk)->initialPosition()<<endl;
1138 cout<<
" end at "<<(*i_mcTrk)->finalPosition()<<endl;
1139 cout<<
" Helix: "<<mchelix.
a()<<endl;
1140 cout<<
"mc ptot, theta, phi, R = "<<m_mcptot<<
", "<<mom2.theta()/acos(-1.)*180<<
", "<<mom2.phi()/acos(-1.)*180<<
", "<<mchelix.
radius()<<endl;
1149 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
1151 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
1152 return( StatusCode::SUCCESS);
1154 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
1160 mtrkadd_mgr->clear();
1164 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
1165 Hep3Vector csmp3[2];
1168 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
1169 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
1171 csmp3[kj-1]=(*iter_trk)->p3();
1172 csmphi[kj-1] = (*iter_trk)->phi();
1176 for(
int j = 0, ij = 0; j<5; j++) {
1177 m_trkhelix[j] = (*iter_trk)->helix()[j];
1179 for(
int k=0; k<=j; k++,ij++) {
1180 m_trkerror[ij] = (*iter_trk)->err()[j][k];
1184 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
1186 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
1187 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
1188 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
1190 m_trkndf = (*iter_trk)->ndof();
1191 m_trkchisq = (*iter_trk)->chi2();
1193 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
1195 StatusCode sc3 = m_nt3->write();
1196 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
1226 log << MSG::DEBUG <<
"retrieved MDC tracks:"
1227 <<
" Nhits " <<(*iter_trk)->getNhits()
1228 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
1230 HitRefVec gothits = (*iter_trk)->getVecHits();
1234 rectrk->
id = (*iter_trk)->trackId();
1235 rectrk->
chiSq = (*iter_trk)->chi2();
1236 rectrk->
ndf = (*iter_trk)->ndof();
1237 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
1238 rectrk->
nhits = (*iter_trk)->getNhits();
1239 rectrk->
nster = (*iter_trk)->nster();
1240 rectrk->
nclus = (*iter_trk)->getNcluster();
1241 rectrk->
stat = (*iter_trk)->stat();
1242 status_temp = (*iter_trk)->stat();
1244 trkadd->
id = (*iter_trk)->trackId();
1248 trkadd->
body = rectrk;
1249 rectrk->
add = trkadd;
1251 for (
int i=0; i<5; i++) {
1252 rectrk->
helix[i] = (*iter_trk)->helix()[i];
1253 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
1254 for(
int j = 1; j<i+2;j++) {
1255 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
1258 std::sort(gothits.begin(), gothits.end(), order_rechits);
1259 HitRefVec::iterator it_gothit = gothits.begin();
1260 for( ; it_gothit != gothits.end(); it_gothit++) {
1262 if( (*it_gothit)->getStat() != 1 ) {
1264 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
1269 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
1270 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
1271 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
1272 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
1273 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
1274 <<
" id of hit "<<(*it_gothit)->getId()
1275 <<
" track id of hit "<<(*it_gothit)->getTrkId()
1276 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
1279 whit->
id = (*it_gothit)->getId();
1280 whit->
ddl = (*it_gothit)->getDriftDistLeft();
1281 whit->
ddr = (*it_gothit)->getDriftDistRight();
1282 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
1283 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
1284 whit->
pChiSq = (*it_gothit)->getChisqAdd();
1285 whit->
lr = (*it_gothit)->getFlagLR();
1286 whit->
stat = (*it_gothit)->getStat();
1287 mdcid = (*it_gothit)->getMdcId();
1291 int wid = w0id + localwid;
1293 <<
"lr from PR: "<<whit->
lr
1294 <<
" layerId = " << layid
1295 <<
" wireId = " << localwid
1305 whit->
tdc = (*it_gothit)->getTdc();
1306 whit->
adc= (*it_gothit)->getAdc();
1307 rectrk->
hitcol.push_back(whit);
1308 mhit_mgr->push_back(*whit);
1310 mtrk_mgr->push_back(*rectrk);
1311 mtrkadd_mgr->push_back(*trkadd);
1319 m_trkdelx = trkx1 - trkx2;
1320 m_trkdely = trky1 - trky2;
1321 m_trkdelz = trkz1 - trkz2;
1322 m_trkdelthe = trkthe1 + trkthe2;
1323 m_trkdelphi = trkphi1- trkphi2;
1324 m_trkdelp = trkp1 - trkp2;
1325 StatusCode sc4 = m_nt4->write();
1326 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1329 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1330 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1341 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1342 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1347 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1354 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1359 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1361 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++){
1363 for(
int hypo=0; hypo<5; hypo++)
1365 if((*KalTrk)->getStat(0,hypo)==1) nFailedTrks[hypo]++;
1369 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1371 int nhitofthistrk=0;
1372 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1379 iter_hit=gothelixsegs.begin();
1427 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
1429 log << MSG::WARNING <<
"Could not retrieve Cgem MC truth" << endreq;
1430 return StatusCode::FAILURE;
1432 m_evt4=eventHeader->eventNumber();
1433 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
1435 double dprex,dprey,dprez,dpostx,dposty,dpostz;
1437 for(; iter_truth != cgemMcHitCol->end(); ++iter_truth){
1438 layer = (*iter_truth)->GetLayerID();
1439 dprex = (*iter_truth)->GetPositionXOfPrePoint();
1440 dprey = (*iter_truth)->GetPositionYOfPrePoint();
1441 dprez = (*iter_truth)->GetPositionZOfPrePoint();
1442 dpostx = (*iter_truth)->GetPositionXOfPostPoint();
1443 dposty = (*iter_truth)->GetPositionYOfPostPoint();
1444 dpostz = (*iter_truth)->GetPositionZOfPostPoint();
1452 if((dprex >=0&&dprey>=0)||(dprex <0&&dprey>=0)){
1454 m_dtphi[jj] = acos(dprex/diR);
1457 m_dtpostphi[jj] = acos(dpostx/doR);
1458 m_dtpostz[jj]= dpostz;
1460 if((dprex <0&&dprey<0)||(dprex >=0&&dprey<0)){
1462 m_dtphi[jj] = -acos(dprex/diR);
1465 m_dtpostphi[jj] = -acos(dpostx/doR);
1466 m_dtpostz[jj]= dpostz;
1468 midx=(dprex+dpostx)/2;
1469 midy=(dprey+dposty)/2;
1470 if((midx>=0&&midy>=0)||(midx<0&&midy>=0)){
1472 m_dtphi[jj]=acos(midx/dmR);
1474 if((midx<0&&midy<0)||(midx>=0&&midy<0)){
1476 m_dtphi[jj]=-acos(midx/dmR);
1478 m_dtv[jj] = (m_dtv[jj]+m_dtpostz[jj])/20;
1479 m_tlayer[jj] = layer;
1489 return StatusCode::SUCCESS;
1500 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1501 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1505 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
1507 cout<<
"Open event header file failed!"<<endl;
1509 m_run_kal = eventHeader->runNumber();
1510 m_event_kal = eventHeader->eventNumber();
1515 m_dropedHits_kal_e = TrasanTRK.
ndf+5-track.
nchits();
1516 m_kappa2_kal_e = TrasanTRK.
helix[2]*track.
a()[2];
1517 m_trackNhits_kal_e = track.
nchits();
1518 m_trackNster_kal_e = track.
nster();
1519 m_trackNaxis_kal_e = track.
nchits()-track.
nster();
1520 m_chi2_kal_e = track.
chiSq();
1521 m_Ea00_kal_e = track.
Ea()[0][0];
1522 m_Ea11_kal_e = track.
Ea()[1][1];
1523 m_Ea22_kal_e = track.
Ea()[2][2];
1524 m_Ea33_kal_e = track.
Ea()[3][3];
1525 m_Ea44_kal_e = track.
Ea()[4][4];
1529 m_dropedHits_kal_mu = TrasanTRK.
ndf+5-track.
nchits();
1530 m_kappa2_kal_mu = TrasanTRK.
helix[2]*track.
a()[2];
1531 m_trackNhits_kal_mu = track.
nchits();
1532 m_trackNster_kal_mu = track.
nster();
1533 m_trackNaxis_kal_mu = track.
nchits()-track.
nster();
1534 m_chi2_kal_mu = track.
chiSq();
1535 m_Ea00_kal_mu = track.
Ea()[0][0];
1536 m_Ea11_kal_mu = track.
Ea()[1][1];
1537 m_Ea22_kal_mu = track.
Ea()[2][2];
1538 m_Ea33_kal_mu = track.
Ea()[3][3];
1539 m_Ea44_kal_mu = track.
Ea()[4][4];
1544 m_dropedHits_kal_pi = TrasanTRK.
ndf+5-track.
nchits();
1545 m_kappa2_kal_pi = TrasanTRK.
helix[2]*track.
a()[2];
1546 m_trackNhits_kal_pi = track.
nchits();
1547 m_trackNster_kal_pi = track.
nster();
1548 m_trackNaxis_kal_pi = track.
nchits()-track.
nster();
1549 m_chi2_kal_pi = track.
chiSq();
1550 m_Ea00_kal_pi = track.
Ea()[0][0];
1551 m_Ea11_kal_pi = track.
Ea()[1][1];
1552 m_Ea22_kal_pi = track.
Ea()[2][2];
1553 m_Ea33_kal_pi = track.
Ea()[3][3];
1554 m_Ea44_kal_pi = track.
Ea()[4][4];
1559 m_dropedHits_kal_k = TrasanTRK.
ndf+5-track.
nchits();
1560 m_kappa2_kal_k = TrasanTRK.
helix[2]*track.
a()[2];
1561 m_trackNhits_kal_k = track.
nchits();
1562 m_trackNster_kal_k = track.
nster();
1563 m_trackNaxis_kal_k = track.
nchits()-track.
nster();
1564 m_chi2_kal_k = track.
chiSq();
1565 m_Ea00_kal_k = track.
Ea()[0][0];
1566 m_Ea11_kal_k = track.
Ea()[1][1];
1567 m_Ea22_kal_k = track.
Ea()[2][2];
1568 m_Ea33_kal_k = track.
Ea()[3][3];
1569 m_Ea44_kal_k = track.
Ea()[4][4];
1574 m_dropedHits_kal_p = TrasanTRK.
ndf+5-track.
nchits();
1575 m_kappa2_kal_p = TrasanTRK.
helix[2]*track.
a()[2];
1576 m_trackNhits_kal_p = track.
nchits();
1577 m_trackNster_kal_p = track.
nster();
1578 m_trackNaxis_kal_p = track.
nchits()-track.
nster();
1579 m_chi2_kal_p = track.
chiSq();
1580 m_Ea00_kal_p = track.
Ea()[0][0];
1581 m_Ea11_kal_p = track.
Ea()[1][1];
1582 m_Ea22_kal_p = track.
Ea()[2][2];
1583 m_Ea33_kal_p = track.
Ea()[3][3];
1584 m_Ea44_kal_p = track.
Ea()[4][4];
1592 || TrasanTRK.
helix[2]*track.
a()[2]<0)
1596 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1597 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1598 <<
" "<<TrasanTRK.
helix[2]<<
" "<<track.
a()[2]<<
" "<<iqual<<endl;
1599 cout<<
"FillTds> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1600 <<
" nster "<<track.
nster()<<
" iqual "<<iqual<<
" track.Ea "<< track.
Ea()<<endl;
1602 cout<<
"fillTds>.....track.Ea[2][2] "<<track.
Ea()[2][2]<<endl;
1603 cout <<
" TRASAN stereo = " << trasster
1604 <<
" and KalFitTrack = " << trakster << std::endl;
1605 cout <<
" TRASAN axial = " << trasax
1606 <<
" and KalFitTrack = " << trakax << std::endl;
1609 cout <<
"...there is a problem during fit !! " << std::endl;
1610 if (trasster-trakster>5)
1611 cout <<
" because stereo " << trasster-trakster << std::endl;
1612 if (trasax-trakax >5)
1613 cout <<
" because axial " << std::endl;
1614 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1615 cout <<
" because kappa sign " << std::endl;
1621 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1622 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1623 track.
Ea()[4][4] > 0 && iqual) {
1624 if(
debug_ == 4) cout<<
"fillTds>.....going on "<<endl;
1638 if(
debug_==4) cout<<
"ALARM: FillTds Not refit with KalFilter!!! lmass="<<l_mass<<endl;
1649 double a_trasan[5], ea_trasan[15];
1650 for(
int i =0 ; i <5; i++){
1651 a_trasan[i] = TrasanTRK.
helix[i];
1653 for(
int j =0 ; j <15; j++){
1654 ea_trasan[j] = TrasanTRK.
error[j];
1670 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1671 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1674 TrasanTRK.
helix[2]*track.
a()[2]<0)
1678 cout<<
"Nhit from PR "<<TrasanTRK.
nhits<<
" nhit "<<track.
nchits()<<endl;
1679 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1680 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1682 cout<<
"FillTds_lead> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1683 <<
" nster "<<track.
nster()<<
" iqual_front_[l_mass] "<<
iqual_front_[l_mass]<<
" track.Ea "<<track.
Ea()<<endl;
1685 cout <<
" TRASAN stereo = " << trasster
1686 <<
" and KalFitTrack = " << trakster << std::endl;
1687 cout <<
" TRASAN axial = " << trasax
1688 <<
" and KalFitTrack = " << trakax << std::endl;
1691 cout <<
"...there is a problem during fit !! " << std::endl;
1692 if (trasster-trakster>5)
1693 cout <<
" because stereo " << trasster-trakster << std::endl;
1694 if (trasax-trakax >5)
1695 cout <<
" because axial " << std::endl;
1696 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1697 cout <<
" because kappa sign " << std::endl;
1703 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1704 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1715 if (
debug_ == 4) cout<<
" trasan id...1 "<<TrasanTRK.
id<<endl;
1724 if(
debug_==4) cout<<
"ALARM: FillTds_forMdc Not refit with KalFilter!!! lmass="<<l_mass<<endl;
1736 if (
debug_ ==4) cout<<
" trasan id...2 "<<TrasanTRK.
id<<endl;
1741 double a_trasan[5], ea_trasan[15];
1742 for(
int i =0 ; i <5; i++){
1743 a_trasan[i] = TrasanTRK.
helix[i];
1745 for(
int j =0 ; j <15; j++){
1746 ea_trasan[j] = TrasanTRK.
error[j];
1766 cout <<
"fillTds_IP>......"<<endl;
1767 cout <<
" dr = " << track.
a()[0]
1768 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1769 cout <<
" phi0 = " << track.
a()[1]
1770 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1771 cout <<
" PT = " << 1/track.
a()[2]
1772 <<
", Er_kappa =" << sqrt(track.
Ea()[2][2]) << std::endl;
1773 cout <<
" dz = " << track.
a()[3]
1774 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1775 cout <<
" tanl = " << track.
a()[4]
1776 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1902 TrasanTRK.
helix[2]*track.
a()[2]<0)
1908 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1909 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1913 double dr = track.
a()[0];
1914 double phi0 = track.
a()[1];
1915 double kappa = track.
a()[2];
1916 double dz = track.
a()[3];
1917 double tanl = track.
a()[4];
1920 double vx = dr*
cos(phi0);
1921 double vy = dr*
sin(phi0);
1926 if(0==kappa) kappa = 10e-10;
1927 double px = -
sin(phi0)/fabs(kappa);
1928 double py =
cos(phi0)/fabs(kappa);
1929 double pz = tanl/fabs(kappa);
1931 trk->
setX(vx, l_mass);
1932 trk->
setY(vy, l_mass);
1933 trk->
setZ(vz, l_mass);
1934 trk->
setPx(px, l_mass);
1935 trk->
setPy(py, l_mass);
1936 trk->
setPz(pz, l_mass);
1947 if (kappa > 0.0000000001)
1949 else if (kappa < -0.0000000001)
1954 double ptot = sqrt(px*px+py*py+pz*pz);
1955 trk->
setTheta(acos(pz/ptot),l_mass);
1961 if(4==
debug_) cout<<
"track "<< track.
trasan_id() <<
", pid "<< l_mass <<
" fails"<<endl;
1964 double dr = TrasanTRK.
helix[0];
1965 double phi0 = TrasanTRK.
helix[1];
1966 double kappa = TrasanTRK.
helix[2];
1967 double dz = TrasanTRK.
helix[3];
1968 double tanl = TrasanTRK.
helix[4];
1970 double vx = dr*
cos(phi0);
1971 double vy = dr*
sin(phi0);
1974 if(0==kappa) kappa = 10e-10;
1975 double px = -
sin(phi0)/fabs(kappa);
1976 double py =
cos(phi0)/fabs(kappa);
1977 double pz = tanl/fabs(kappa);
1979 trk->
setX(vx, l_mass);
1980 trk->
setY(vy, l_mass);
1981 trk->
setZ(vz, l_mass);
1983 trk->
setPx(px, l_mass);
1984 trk->
setPy(py, l_mass);
1985 trk->
setPz(pz, l_mass);
1992 double a_trasan[5], ea_trasan[15];
1993 for(
int i =0 ; i <5; i++){
1994 a_trasan[i] = TrasanTRK.
helix[i];
1996 for(
int j =0 ; j <15; j++){
1997 ea_trasan[j] = TrasanTRK.
error[j];
2004 if (kappa > 0.0000000001)
2006 else if (kappa < -0.0000000001)
2011 double ptot = sqrt(px*px+py*py+pz*pz);
2012 trk->
setTheta(acos(pz/ptot),l_mass);
2017 std::cout<<
"px: "<<trk->
px()<<
" py: "<<trk->
py()<<
" pz: "<<trk->
pz()<<std::endl;
2018 std::cout<<
"vx: "<<trk->
x()<<
" vy: "<<trk->
y()<<
" vz: "<<trk->
z()<<std::endl;
2035 if(
debug_ == 4) cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0] "<<trk->
getNdf(0,2)<<endl;
2039 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
2040 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
2041 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
2052 if(
debug_ == 4) cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
2062 for (
int i = 0; i<43; i++) {
2076 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
2077 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2078 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2079 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2083 if(
debug_==4) 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;
2127 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
2128 double light_speed( 29.9792458 );
2129 double pt( 1.0 / track_rep.
kappa() );
2130 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2136 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<endl;
2137 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<endl;
2142 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2143 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2145 track_rep.
pivot(IP);
2161 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2162 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2163 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2165 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2166 cout <<
" dr = " << track.
a()[0]
2167 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2168 cout<<
" phi0 = " << track.
a()[1]
2169 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2170 cout <<
" PT = " << 1/track.
a()[2]
2171 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2172 cout <<
" dz = " << track.
a()[3]
2173 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2174 cout <<
" tanl = " << track.
a()[4]
2175 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2176 cout <<
" Ea = " << track.
Ea() <<endl;
2203 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
2204 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
2205 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;
2209 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
2210 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
2211 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
2220 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
2231 helixseg->
setResIncl(it->residual_include());
2232 helixseg->
setResExcl(it->residual_exclude());
2234 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
2236 helixseg->
setDrIncl(it->a_include()[0]);
2239 helixseg->
setDzIncl(it->a_include()[3]);
2243 helixseg->
setDrExcl(it->a_exclude()[0]);
2246 helixseg->
setDzExcl(it->a_exclude()[3]);
2261 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
2262 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
2263 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
2264 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
2266 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
2267 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
2268 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
2270 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
2274 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
2275 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
2276 helixseg->
setFlagLR(it->HitMdc()->LR());
2277 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
2278 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
2279 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
2280 helixseg->
setTof(it->tof());
2283 helixseg->
setDD(it->dd());
2284 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
2285 helixseg->
setDT(it->dt());
2286 segcol->push_back(helixseg);
2287 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
2288 helixsegrefvec.push_back(refhelixseg);
2291 m_docaInc = helixseg -> getDocaIncl();
2292 m_docaExc = helixseg -> getDocaExcl();
2293 m_residualInc = helixseg -> getResIncl();
2294 m_residualExc = helixseg -> getResExcl();
2295 m_dd = helixseg -> getDD();
2297 m_tdrift = helixseg -> getDT();
2298 m_layerid = helixseg -> getLayerId();
2299 m_yposition= it->HitMdc()->wire().fwd().y();
2301 StatusCode sc6 = m_nt6->write();
2302 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
2309 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
2317 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
2321 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
2330 for (
int i = 0; i<43; i++) {
2341 double a_trasan[5], ea_trasan[15];
2342 for(
int i =0 ; i <5; i++){
2343 a_trasan[i] = TrasanTRK.
helix[i];
2345 for(
int j =0 ; j <15; j++){
2346 ea_trasan[j] = TrasanTRK.
helix[j];
2352 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
2353 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2354 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2355 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2360 if(
debug_==4) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
2366 TrasanTRK.
pivot[2]);
2369 for(
int i = 0; i < 5; i++)
2370 a[i] = TrasanTRK.
helix[i];
2373 for(
int i = 0, k = 0; i < 5; i++) {
2374 for(
int j = 0; j <= i; j++) {
2376 ea[j][i] = ea[i][j];
2382 double fiTerm = TrasanTRK.
fiTerm;
2384 double fi0 = track_rep.
phi0();
2390 if( fabs(
x ) > 1.0e-10 ){
2391 phi_x = atan2( y,
x );
2392 if( phi_x < 0 ) phi_x += 2*
M_PI;
2394 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
2396 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
2397 double dphi = fabs( fiTerm + fi0 - phi_x );
2398 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
2399 double tanl = track_rep.
tanl();
2400 double cosl_inv = sqrt( tanl*tanl + 1.0 );
2402 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
2403 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
2405 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
2406 double light_speed( 29.9792458 );
2407 double pt( 1.0 / track_rep.
kappa() );
2408 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2416 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
2417 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
2422 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2423 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2425 track_rep.
pivot(IP);
2443 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2444 cout <<
" dr = " << track.
a()[0]
2445 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2446 cout<<
" phi0 = " << track.
a()[1]
2447 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2448 cout <<
" PT = " << 1/track.
a()[2]
2449 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2450 cout <<
" dz = " << track.
a()[3]
2451 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2452 cout <<
" tanl = " << track.
a()[4]
2453 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2454 cout <<
" Ea = " << track.
Ea() <<endl;
2479 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
2480 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
2485 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
2486 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
2487 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
2496 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
2507 helixseg->
setResIncl(it->residual_include());
2508 helixseg->
setResExcl(it->residual_exclude());
2510 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
2537 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
2538 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
2539 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
2540 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
2542 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
2543 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
2544 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
2546 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
2550 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
2551 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
2552 helixseg->
setFlagLR(it->HitMdc()->LR());
2553 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
2554 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
2555 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
2556 helixseg->
setTof(it->tof());
2559 helixseg->
setDD(it->dd());
2560 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
2561 helixseg->
setDT(it->dt());
2563 segcol->push_back(helixseg);
2564 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
2565 helixsegrefvec.push_back(refhelixseg);
2567 m_docaInc = helixseg -> getDocaIncl();
2568 m_docaExc = helixseg -> getDocaExcl();
2569 m_residualInc = helixseg -> getResIncl();
2570 m_residualExc = helixseg -> getResExcl();
2571 m_dd = helixseg -> getDD();
2573 m_tdrift = helixseg -> getDT();
2574 m_layerid = helixseg -> getLayerId();
2575 m_yposition= it->HitMdc()->wire().fwd().y();
2577 StatusCode sc6 = m_nt6->write();
2578 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
2586 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
2594 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
2598 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
2607 for (
int i = 0; i<43; i++) {
2618 double a_trasan[5], ea_trasan[15];
2619 for(
int i =0 ; i <5; i++){
2620 a_trasan[i] = TrasanTRK.
helix[i];
2622 for(
int j =0 ; j <15; j++){
2623 ea_trasan[j] = TrasanTRK.
helix[j];
2629 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
2630 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2631 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2632 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2637 if(
debug_==4) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
2643 TrasanTRK.
pivot[2]);
2646 for(
int i = 0; i < 5; i++)
2647 a[i] = TrasanTRK.
helix[i];
2650 for(
int i = 0, k = 0; i < 5; i++) {
2651 for(
int j = 0; j <= i; j++) {
2653 ea[j][i] = ea[i][j];
2659 double fiTerm = TrasanTRK.
fiTerm;
2661 double fi0 = track_rep.
phi0();
2667 if( fabs(
x ) > 1.0e-10 ){
2668 phi_x = atan2( y,
x );
2669 if( phi_x < 0 ) phi_x += 2*
M_PI;
2671 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
2673 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
2674 double dphi = fabs( fiTerm + fi0 - phi_x );
2675 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
2676 double tanl = track_rep.
tanl();
2677 double cosl_inv = sqrt( tanl*tanl + 1.0 );
2679 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
2680 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
2683 double track_len(fabs( track_rep.
radius() * fiTerm * cosl_inv ));
2685 double light_speed( 29.9792458 );
2686 double pt( 1.0 / track_rep.
kappa() );
2687 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2695 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
2696 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
2701 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2702 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2706 track_rep.
pivot(LPiovt);
2793 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2794 cout <<
" dr = " << track.
a()[0]
2795 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2796 cout<<
" phi0 = " << track.
a()[1]
2797 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2798 cout <<
" PT = " << 1/track.
a()[2]
2799 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2800 cout <<
" dz = " << track.
a()[3]
2801 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2802 cout <<
" tanl = " << track.
a()[4]
2803 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2804 cout <<
" Ea = " << track.
Ea() <<endl;
2820 for(
int jj = 0; jj<2; jj++) {
2834 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2835 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2838 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2843 Hep3Vector xorigin(0,0,0);
2845 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
2849 xorigin.setX(dbv[0]);
2850 xorigin.setY(dbv[1]);
2851 xorigin.setZ(dbv[2]);
2865 double tanl = track.
tanl();
2866 double phi_old = work.
phi0();
2867 double phi = track.
phi0();
2869 if (fabs(phi - phi_old) >
M_PI) {
2870 if (phi > phi_old) phi -= 2 *
M_PI;
2871 else phi_old -= 2 *
M_PI;
2874 double path_zero = fabs(track.
radius() * (phi_old-phi)* sqrt(1 + tanl * tanl));
2878 HepSymMatrix Eakal(5,0);
2886 unsigned int nhit = track.
HitsMdc().size();
2887 int layer_prev = -1;
2889 HepVector pos_old(3,0);
2890 double r0kal_prec(0);
2892 for(
unsigned i=0 ; i < nhit; i++ ) {
2893 int ihit = (nhit-1)-i;
2897 int wireid = Wire.
geoID();
2901 if (
pathl_ && layer != layer_prev) {
2903 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
2904 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
2912 Hep3Vector wire = (CLHEP::Hep3Vector)fwd - (CLHEP::Hep3Vector)bck;
2915 work.
pivot((fwd + bck) * .5);
2916 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
2918 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2922 double tension = geowire->
Tension();
2925 double zinit(x0kal.z()), lzx(Wire.
lzx());
2928 double A = 47.35E-6/tension;
2929 double Zp = (zinit - bck.z())*lzx/wire.z();
2932 std::cout<<
" sag in smoother_anal: "<<std::endl;
2933 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2934 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2935 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2936 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2937 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2938 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2941 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2942 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2943 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2945 wire.setX(wire.x()/wire.z());
2946 wire.setY(result.z());
2949 x0kal.setX(result.x());
2950 x0kal.setY(result.y());
2953 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2956 double r0kal = x0kal.perp();
2958 cout<<
"wire direction "<<wire<<endl;
2959 cout<<
"x0kal "<<x0kal<<endl;
2960 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
2977 double pmag( sqrt( 1.0 + track.
a()[4]*track.
a()[4]) / track.
a()[2]);
2978 double mass_over_p( track.
mass()/ pmag );
2979 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2980 double tofest = pathl / ( 29.9792458 * beta );
2990 if(!(way<0&&fabs(track.
kappa())>1000.0)) {
3000 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3002 HepSymMatrix Ma(5,0);
3005 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3019 point.setX(x0kal.x() + track.
a()[0]*
cos(track.
a()[1]));
3020 point.setY(x0kal.y() + track.
a()[0]*
sin(track.
a()[1]));
3021 point.setZ(x0kal.z() + track.
a()[3]);
3025 double phi_old = track.
a()[1];
3028 double phi_new = temp.
a()[1];
3029 double fi = phi_new - phi_old;
3033 if(fabs(fi) >= CLHEP::twopi) fi = fmod(fi+2*CLHEP::twopi,CLHEP::twopi);
3038 if (
debug_==4) cout<<
"track----7-----"<<track.
a()<<endl;
3050 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3051 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3054 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3057 HepSymMatrix Eakal(5,0);
3063 std::cout<<
"the initial error matrix in smoothing is "<<Eakal<<std::endl;
3068 unsigned int nseg = track.
HelixSegs().size();
3069 int layer_prev = -1;
3071 HepVector pos_old(3,0);
3072 double r0kal_prec(0);
3074 for(
unsigned i=0 ; i < nseg; i++ ) {
3077 int iseg = (nseg-1)-i;
3081 int wireid = Wire.
geoID();
3085 if (
pathl_ && layer != layer_prev) {
3087 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
3088 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
3096 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3099 work.
pivot((fwd + bck) * .5);
3100 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
3103 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3109 double tension = geowire->
Tension();
3112 double zinit(x0kal.z()), lzx(Wire.
lzx());
3116 double A = 47.35E-6/tension;
3117 double Zp = (zinit - bck.z())*lzx/wire.z();
3121 std::cout<<
" sag in smoother_calib: "<<std::endl;
3122 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3123 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3124 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3125 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3126 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3127 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3130 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3131 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3132 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3134 wire.setX(wire.x()/wire.z());
3135 wire.setY(result.z());
3138 x0kal.setX(result.x());
3139 x0kal.setY(result.y());
3142 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3146 double r0kal = x0kal.perp();
3148 cout<<
"wire direction "<<wire<<endl;
3149 cout<<
"x0kal "<<x0kal<<endl;
3150 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
3157 if (
debug_ == 4) cout<<
"track----6-----"<<track.
a()<<
" ...path..."<<pathl
3158 <<
"momentum"<<track.
momentum(0)<<endl;
3163 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3165 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3172 if (
debug_==4) cout<<
"track----7-----"<<track.
a()<<endl;
3187 cout<<
"**********************"<<endl;
3188 cout<<
"filter pid type "<<l_mass<<endl;
3193 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3194 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3197 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3201 HepVector pos_old(3,0);
3202 double r0kal_prec(0);
3204 int nhit = track.
HitsMdc().size();
3205 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
3206 for(
int i=0 ; i < nhit; i++ ) {
3209 if (HitMdc.
chi2()<0)
continue;
3215 int wireid = Wire.
geoID();
3219 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3222 work.
pivot((fwd + bck) * .5);
3229 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
3232 std::cout<<
" a="<<track.
a()<<std::endl;
3233 std::cout<<
" pivot="<<track.
pivot()<<std::endl;
3234 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3279 double tension = geowire->
Tension();
3281 double zinit(x0kal.z()), lzx(Wire.
lzx());
3283 double A = 47.35E-6/tension;
3284 double Zp = (zinit - bck.z())*lzx/wire.z();
3287 std::cout<<
" sag in filter_fwd_anal: "<<std::endl;
3288 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3289 std::cout<<
"zinit: "<<zinit<<
" bck.z(): "<<bck.z()<<std::endl;
3290 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3291 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3292 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3293 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3294 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3297 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3298 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3299 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3301 wire.setX(wire.x()/wire.z());
3302 wire.setY(result.z());
3304 x0kal.setX(result.x());
3305 x0kal.setY(result.y());
3310 double r0kal = x0kal.perp();
3315 if(i==0) track.
pivot(x0kal);
3319 if (nhits_read == 1) {
3330 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3331 std::cout<<
" a="<<track.
a()<<std::endl;
3332 std::cout<<
" Ea="<<track.
Ea()<<std::endl;
3335 double dtracknew = 0.;
3339 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3340 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3341 double diff_chi2 = track.
chiSq();
3342 Hep3Vector IP(0,0,0);
3348 for(
unsigned k=i+1 ; k < nhit; k++ )
3353 double dchi2 = -1.0;
3408 m_residest = dtrack-dtdc;
3409 m_residnew = dtracknew -dtdc;
3412 m_anal_dr = worktemp.
a()[0];
3413 m_anal_phi0 = worktemp.
a()[1];
3414 m_anal_kappa = worktemp.
a()[2];
3415 m_anal_dz = worktemp.
a()[3];
3416 m_anal_tanl = worktemp.
a()[4];
3417 m_anal_ea_dr = worktemp.
Ea()[0][0];
3418 m_anal_ea_phi0 = worktemp.
Ea()[1][1];
3419 m_anal_ea_kappa = worktemp.
Ea()[2][2];
3420 m_anal_ea_dz = worktemp.
Ea()[3][3];
3421 m_anal_ea_tanl = worktemp.
Ea()[4][4];
3422 StatusCode sc5 = m_nt5->write();
3423 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3431 m_dchi2_hit_e[m_hit_no] = dchi2;
3432 m_residest_hit_e[m_hit_no] = dtrack-dtdc;
3433 m_residnew_hit_e[m_hit_no] = dtracknew -dtdc;
3436 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
3437 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
3438 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
3439 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
3440 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
3441 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
3442 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
3443 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
3444 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
3445 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
3449 else if(l_mass == 1){
3450 m_dchi2_hit_mu[m_hit_no] = dchi2;
3451 m_residest_hit_mu[m_hit_no] = dtrack-dtdc;
3452 m_residnew_hit_mu[m_hit_no] = dtracknew -dtdc;
3455 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
3456 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
3457 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
3458 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
3459 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
3460 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
3461 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
3462 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
3463 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
3464 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
3469 else if(l_mass == 2){
3470 m_dchi2_hit_pi[m_hit_no] = dchi2;
3471 m_residest_hit_pi[m_hit_no] = dtrack-dtdc;
3472 m_residnew_hit_pi[m_hit_no] = dtracknew -dtdc;
3475 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
3476 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
3477 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
3478 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
3479 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
3480 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
3481 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
3482 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
3483 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
3484 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
3489 else if(l_mass == 3){
3490 m_dchi2_hit_k[m_hit_no] = dchi2;
3491 m_residest_hit_k[m_hit_no] = dtrack-dtdc;
3492 m_residnew_hit_k[m_hit_no] = dtracknew -dtdc;
3495 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
3496 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
3497 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
3498 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
3499 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
3500 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
3501 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
3502 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
3503 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
3504 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
3509 else if(l_mass == 4){
3510 m_dchi2_hit_p[m_hit_no] = dchi2;
3511 m_residest_hit_p[m_hit_no] = dtrack-dtdc;
3512 m_residnew_hit_p[m_hit_no] = dtracknew -dtdc;
3515 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
3516 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
3517 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
3518 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
3519 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
3520 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
3521 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
3522 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
3523 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
3524 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
3535 diff_chi2 = chi2 - diff_chi2;
3536 HitMdc.
chi2(diff_chi2);
3550 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3551 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3554 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3558 std::cout<<
"filter_fwd_calib starting ...the inital error Matirx is "<<track.
Ea()<<std::endl;
3560 Hep3Vector x0inner(track.
pivot());
3561 HepVector pos_old(3,0);
3562 double r0kal_prec(0);
3565 unsigned int nhit = track.
HitsMdc().size();
3566 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
3567 for(
unsigned i=0 ; i < nhit; i++ ) {
3571 cout<<
"filter_fwd...........222 chi2="<<HitMdc.
chi2()<<endl;
3573 if (HitMdc.
chi2()<0)
continue;
3576 int wireid = Wire.
geoID();
3582 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3585 work.
pivot((fwd + bck) * .5);
3586 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
3589 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3594 double tension = geowire->
Tension();
3596 double zinit(x0kal.z()), lzx(Wire.
lzx());
3598 double A = 47.35E-6/tension;
3599 double Zp = (zinit - bck.z())*lzx/wire.z();
3603 std::cout<<
" sag in filter_fwd_calib: "<<std::endl;
3604 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3605 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3606 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3607 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3608 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3609 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3612 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3613 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3614 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3615 wire.setX(wire.x()/wire.z());
3616 wire.setY(result.z());
3618 x0kal.setX(result.x());
3619 x0kal.setY(result.y());
3623 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3626 double r0kal = x0kal.perp();
3635 cout <<
"*** move from " << track.
pivot() <<
" ( Kappa = "
3636 << track.
a()[2] <<
")" << endl;
3644 cout <<
"*** to " << track.
pivot() <<
" ( Kappa = "
3645 << track.
a()[2] <<
")" << std::endl;
3647 if (nhits_read == 1) {
3649 if(
debug_==4) cout <<
"filter_fwd::Ea = " << track.
Ea()<<endl;
3662 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3663 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3665 double diff_chi2 = track.
chiSq();
3667 Hep3Vector IP(0,0,0);
3673 for(
unsigned k=i+1 ; k < nhit; k++ )
3678 double dchi2 = -1.0;
3681 std::cout<<
"the value of x0kal is "<<x0kal<<std::endl;
3682 std::cout<<
"the value of track.pivot() is "<<track.
pivot()<<std::endl;
3686 HepSymMatrix Ma(5,0);
3690 std::cout<<
"HelixSeg.Ea_pre_fwd() ..."<<HelixSeg.
Ea_pre_fwd()<<std::endl;
3691 std::cout<<
"HelixSeg.a_pre_fwd() ..."<<HelixSeg.
a_pre_fwd()<<std::endl;
3692 std::cout<<
"HelixSeg.Ea_filt_fwd() ..."<<HelixSeg.
Ea_filt_fwd()<<std::endl;
3705 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;
3707 if(
debug_ == 4) cout<<
"****inext***"<<inext<<
" !!!! dchi2= "<< dchi2
3708 <<
" chisq= "<< chi2<< endl;
3713 StatusCode sc5 = m_nt5->write();
3714 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3720 diff_chi2 = chi2 - diff_chi2;
3721 HitMdc.
chi2(diff_chi2);
3725 cout <<
" -> after include meas = " << meas
3726 <<
" at R = " << track.
pivot().perp() << std::endl;
3727 cout <<
" chi2 = " << chi2 <<
", diff_chi2 = "
3728 << diff_chi2 <<
", LR = "
3729 << HitMdc.
LR() <<
", stereo = " << HitMdc.
wire().
stereo()
3731 cout<<
"filter_fwd..........HitMdc.chi2... "<<HitMdc.
chi2()<<endl;
3734 if(dchi2>0.0 && (way!=999)) {
3745 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3746 for(
int i = 0; i < _BesKalmanFitWalls.size(); i++) {
3747 _BesKalmanFitWalls[i].updateTrack(track, way);
3748 if (
debug_ == 4) cout<<
"Wall "<<i<<
" update the track!(filter)"<<endl;
3753 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3754 for(
int i = begin; i <= end; i++) {
3755 _BesKalmanFitWalls[i].updateTrack(track, way);
3756 if(
debug_ == 4){cout <<
"wall" << i << endl;
3757 cout <<
"track" <<
"pivot=" << track.
pivot() <<
"a=" << track.
a() <<
"ea=" <<track.
Ea() << endl;}
3758 if (
debug_ == 4) cout<<
"Wall "<<i<<
" update the track!(filter)"<<endl;
3763 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3764 double ri=(r1<r2)?r1:r2;
3765 double ro=(r1<r2)?r2:r1;
3768 vector<KalFitCylinder>::iterator it;
3769 vector<KalFitCylinder>::iterator itStart;
3770 vector<KalFitCylinder>::iterator itStop;
3773 itStart=_BesKalmanFitWalls.end()-index-1;
3774 itStop =_BesKalmanFitWalls.begin()-1;
3778 itStart=_BesKalmanFitWalls.begin()+index;
3779 itStop =_BesKalmanFitWalls.end();
3781 for(it=itStart; it != itStop; it=it+step) {
3782 if(ri < it->rmin() && it->radius() < ro) {
3783 it->updateTrack(track,way);
3788 else if(it->rmin() <= ri && ri <= it->radius() && it->radius() < ro){
3789 it->updateTrack(track,way,ri,it->
radius());
3798 else if(ri < it->rmin() && it->rmin() <= ro && ro <= it->radius()){
3799 it->updateTrack(track,way,it->rmin(),ro);
3808 else if(it->radius()<ri || ro < it->rmin()){
3812 it->updateTrack(track,way,ri,ro);
3824 index = _BesKalmanFitWalls.end() - it-1;
3831 index = it - _BesKalmanFitWalls.begin();
3846 MsgStream log(
msgSvc(), name());
3847 double Pt_threshold(0.3);
3848 Hep3Vector IP(0,0,0);
3854 if ( !&whMgr )
return;
3857 int ntrk = mdcMgr->size();
3858 double* rPt =
new double[ntrk];
3859 int* rOM =
new int[ntrk];
3860 unsigned int* order =
new unsigned int[ntrk];
3861 unsigned int* rCont =
new unsigned int[ntrk];
3862 unsigned int* rGen =
new unsigned int[ntrk];
3864 if(
testOutput) cout<<
"nMdcTrk = "<<ntrk<<endl;
3867 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3868 end = mdcMgr->end(); it != end; it++) {
3872 rPt[index] = 1 / fabs(it->helix[2]);
3873 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3874 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3876 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3877 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3879 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3880 ii !=pt.begin()-1; ii--) {
3881 int lyr((*ii)->geo->Lyr()->Id());
3882 if (outermost < lyr) outermost = lyr;
3883 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3885 rOM[index] = outermost;
3886 order[index] = index;
3891 for (
int j, k = ntrk - 1; k >= 0; k = j){
3893 for(
int i = 1; i <= k; i++)
3894 if(rPt[i - 1] < rPt[i]){
3896 std::swap(order[i], order[j]);
3897 std::swap(rPt[i], rPt[j]);
3898 std::swap(rOM[i], rOM[j]);
3899 std::swap(rCont[i], rCont[j]);
3900 std::swap(rGen[i], rGen[j]);
3908 DataObject *aReconEvent;
3909 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3913 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3914 if(sc!=StatusCode::SUCCESS) {
3915 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3923 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3925 for(
int l = 0; l < ntrk; l++) {
3934 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3938 int trasqual = TrasanTRK_add.
quality;
3939 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3940 if (trasqual)
continue;
3944 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3948 if ((TrasanTRK_add.
decision & 32) == 32 ||
3949 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3954 TrasanTRK.
pivot[2]);
3957 for(
int i = 0; i < 5; i++)
3958 a[i] = TrasanTRK.
helix[i];
3961 for(
int i = 0, k = 0; i < 5; i++) {
3962 for(
int j = 0; j <= i; j++) {
3964 ea[j][i] = ea[i][j];
3967 double fiTerm = TrasanTRK.
fiTerm;
3973 double pT_trk = fabs(track_lead.
pt());
3976 int inlyr(999), outlyr(-1);
3977 int* rStat =
new int[43];
3978 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3979 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3981 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3984 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3985 ii != pt.begin()-1; ii--) {
3986 Num[(*ii)->geo->Lyr()->Id()]++;
3989 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3990 ii != pt.begin()-1; ii--) {
3993 if (pT_trk>0.12 && Num[(*ii)->geo->Lyr()->Id()]>3) {
3995 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3996 <<
" hits in the layer "
3997 << (*ii)->geo->Lyr()->Id() << std::endl;
4017 double dist[2] = {rechit.
ddl, rechit.
ddr};
4018 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4023 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4025 else if (rechit.
lr==1) lr_decision=1;
4028 int ind(geo->
Lyr()->
Id());
4032 lr_decision, rechit.
tdc,
4037 if (inlyr>ind) inlyr = ind;
4038 if (outlyr<ind) outlyr = ind;
4041 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4044 int empty_between(0), empty(0);
4045 for (
int i= inlyr; i <= outlyr; i++)
4046 if (!rStat[i]) empty_between++;
4047 empty = empty_between+inlyr+(42-outlyr);
4054 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
4060 track_lead.
type(type);
4062 unsigned int nhit = track_lead.
HitsMdc().size();
4063 if (!nhit &&
debug_ == 4) {
4064 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4069 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4072 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4073 track_lead.
pivot(outer_pivot);
4079 HepSymMatrix Eakal(5,0);
4082 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4092 cout <<
"from Mdc Pattern Recognition: " << std::endl;
4098 cout <<
" dr = " << work.
a()[0]
4099 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4100 cout <<
" phi0 = " << work.
a()[1]
4101 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4102 cout <<
" PT = " << 1/work.
a()[2]
4103 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4104 cout <<
" dz = " << work.
a()[3]
4105 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4106 cout <<
" tanl = " << work.
a()[4]
4107 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4111 cout<<
"----------------------"<<endl;
4112 cout<<
"track "<<l<<
", pid = "<<
lead_<<
": "<<endl;
4129 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4134 cout <<
" dr = " << work.
a()[0]
4135 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4136 cout <<
" phi0 = " << work.
a()[1]
4137 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4138 cout <<
" PT = " << 1/work.
a()[2]
4139 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4140 cout <<
" dz = " << work.
a()[3]
4141 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4142 cout <<
" tanl = " << work.
a()[4]
4143 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4150 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
4154 IDataProviderSvc* eventSvc = NULL;
4155 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
4157 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
4159 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
4163 IDataManagerSvc *dataManSvc;
4164 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
4165 DataObject *aKalTrackCol;
4166 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
4167 if(aKalTrackCol != NULL) {
4168 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
4169 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4171 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4172 if( kalsc.isFailure() ) {
4173 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4176 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4179 DataObject *aRecKalSegEvent;
4180 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4181 if(aRecKalSegEvent!=NULL) {
4183 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4184 if(segsc != StatusCode::SUCCESS) {
4185 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4189 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4190 if( segsc.isFailure() ) {
4191 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4194 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4196 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.);
4197 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
4198 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
4201 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
4203 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4206 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4207 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4208 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4209 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4210 <<
"Track Id: " << (*iter_trk)->getTrackId()
4211 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4212 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4213 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4214 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4215 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4216 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
4217 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4218 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4220 for(
int i = 0; i<43; i++) {
4221 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4222 << (*iter_trk)->getPathl(i) <<endreq;
4225 m_trackid = (*iter_trk)->getTrackId();
4227 for(
int jj =0, iii=0; jj<5; jj++) {
4229 m_length[jj] = (*iter_trk)->getLength(jj);
4230 m_tof[jj] = (*iter_trk)->getTof(jj);
4231 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4232 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4233 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4234 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4235 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4236 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4237 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4238 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4239 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4240 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4241 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4242 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4243 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4244 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4245 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4246 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4249 for(
int kk=0; kk<=jj; kk++,iii++) {
4250 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4251 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4252 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4253 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4254 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4255 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4256 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4257 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4258 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4259 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4260 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4261 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4262 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4263 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4264 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4304 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4305 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4306 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4307 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4308 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4309 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4310 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4311 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4312 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4313 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4315 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4316 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4317 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4318 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4319 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4320 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4321 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4322 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4323 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4324 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4326 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4327 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4328 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4329 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4330 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4331 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4332 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4333 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4334 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4335 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4337 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4338 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4339 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4340 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4341 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4343 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4344 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4345 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4346 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4347 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4349 m_lpt = 1/m_lhelix[2];
4350 m_lpte = 1/m_lhelixe[2];
4351 m_lptmu = 1/m_lhelixmu[2];
4352 m_lptk = 1/m_lhelixk[2];
4353 m_lptp = 1/m_lhelixp[2];
4355 m_fpt = 1/m_fhelix[2];
4356 m_fpte = 1/m_fhelixe[2];
4357 m_fptmu = 1/m_fhelixmu[2];
4358 m_fptk = 1/m_fhelixk[2];
4359 m_fptp = 1/m_fhelixp[2];
4362 std::cout<<
" "<<std::endl;
4363 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
4364 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
4365 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
4366 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
4367 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
4370 m_zpt = 1/m_zhelix[2];
4371 m_zpte = 1/m_zhelixe[2];
4372 m_zptmu = 1/m_zhelixmu[2];
4373 m_zptk = 1/m_zhelixk[2];
4374 m_zptp = 1/m_zhelixp[2];
4377 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
4378 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
4379 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
4380 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
4381 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
4383 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4384 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4385 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4386 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4387 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4390 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
4391 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
4392 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
4393 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
4394 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
4398 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4399 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4400 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4401 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4402 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4403 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4404 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4405 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4406 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4407 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4408 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4409 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4410 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4411 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4412 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4415 StatusCode sc1 = m_nt1->write();
4416 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4421 phi1 = (*iter_trk)->getFFi0();
4422 r1 = (*iter_trk)->getFDr();
4423 z1 = (*iter_trk)->getFDz();
4424 kap1 = (*iter_trk)->getFCpa();
4425 tanl1 = (*iter_trk)->getFTanl();
4426 charge1 = kap1/fabs(kap1);
4429 p1 = sqrt(1+tanl1*tanl1)/kap1;
4430 the1 =
M_PI/2-atan(tanl1);
4433 pz1= tanl1/fabs(kap1);
4435 }
else if(jj == 2) {
4436 phi2 = (*iter_trk)->getFFi0();
4437 r2 = (*iter_trk)->getFDr();
4438 z2 = (*iter_trk)->getFDz();
4439 kap2 = (*iter_trk)->getFCpa();
4440 tanl2 = (*iter_trk)->getFTanl();
4441 charge2 = kap2/fabs(kap2);
4444 p2 = sqrt(1+tanl2*tanl2)/kap1;
4445 the2 =
M_PI/2-atan(tanl2);
4448 pz2= tanl2/fabs(kap2);
4456 m_delthe = the1 + the2;
4459 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
4460 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
4461 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
4463 StatusCode sc2 = m_nt2->write();
4464 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4473 cout <<
"Kalfitting finished " << std::endl;
4480 MsgStream log(
msgSvc(), name());
4481 double Pt_threshold(0.3);
4482 Hep3Vector IP(0,0,0);
4489 if ( !&whMgr )
return;
4492 int ntrk = mdcMgr->size();
4493 double* rPt =
new double[ntrk];
4494 int* rOM =
new int[ntrk];
4495 unsigned int* order =
new unsigned int[ntrk];
4496 unsigned int* rCont =
new unsigned int[ntrk];
4497 unsigned int* rGen =
new unsigned int[ntrk];
4500 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
4501 end = mdcMgr->end(); it != end; it++) {
4505 rPt[index] = 1 / fabs(it->helix[2]);
4506 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4507 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4509 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4510 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4512 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4513 ii !=pt.begin()-1; ii--) {
4514 int lyr((*ii)->geo->Lyr()->Id());
4515 if (outermost < lyr) outermost = lyr;
4516 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4518 rOM[index] = outermost;
4519 order[index] = index;
4524 for (
int j, k = ntrk - 1; k >= 0; k = j){
4526 for(
int i = 1; i <= k; i++)
4527 if(rPt[i - 1] < rPt[i]){
4529 std::swap(order[i], order[j]);
4530 std::swap(rPt[i], rPt[j]);
4531 std::swap(rOM[i], rOM[j]);
4532 std::swap(rCont[i], rCont[j]);
4533 std::swap(rGen[i], rGen[j]);
4540 DataObject *aReconEvent;
4541 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4545 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4546 if(sc!=StatusCode::SUCCESS) {
4547 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4555 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4558 for(
int l = 0; l < ntrk; l++) {
4560 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
4564 int trasqual = TrasanTRK_add.
quality;
4565 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4566 if (trasqual)
continue;
4570 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4574 if ((TrasanTRK_add.
decision & 32) == 32 ||
4575 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4580 TrasanTRK.
pivot[2]);
4583 for(
int i = 0; i < 5; i++)
4584 a[i] = TrasanTRK.
helix[i];
4587 for(
int i = 0, k = 0; i < 5; i++) {
4588 for(
int j = 0; j <= i; j++) {
4590 ea[j][i] = ea[i][j];
4596 double fiTerm = TrasanTRK.
fiTerm;
4602 int inlyr(999), outlyr(-1);
4603 int* rStat =
new int[43];
4604 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4605 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4607 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4610 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4611 ii != pt.begin()-1; ii--) {
4612 Num[(*ii)->geo->Lyr()->Id()]++;
4616 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4617 ii != pt.begin()-1; ii--) {
4620 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4622 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4623 <<
" hits in the layer "
4624 << (*ii)->geo->Lyr()->Id() << std::endl;
4651 double dist[2] = {rechit.
ddl, rechit.
ddr};
4652 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4657 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4659 else if (rechit.
lr==1) lr_decision=1;
4662 int ind(geo->
Lyr()->
Id());
4664 lr_decision, rechit.
tdc,
4669 if (inlyr>ind) inlyr = ind;
4670 if (outlyr<ind) outlyr = ind;
4673 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4676 int empty_between(0), empty(0);
4677 for (
int i= inlyr; i <= outlyr; i++)
4678 if (!rStat[i]) empty_between++;
4679 empty = empty_between+inlyr+(42-outlyr);
4684 track_lead.
type(type);
4685 unsigned int nhit = track_lead.
HitsMdc().size();
4686 if (!nhit &&
debug_ == 4) {
4687 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4692 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4694 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4697 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4699 track_lead.
pivot(outer_pivot);
4704 HepSymMatrix Eakal(5,0);
4708 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4718 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4724 std::cout <<
" dr = " << work.
a()[0]
4725 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4726 std::cout <<
" phi0 = " << work.
a()[1]
4727 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4728 std::cout <<
" PT = " << 1/work.
a()[2]
4729 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4730 std::cout <<
" dz = " << work.
a()[3]
4731 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4732 std::cout <<
" tanl = " << work.
a()[4]
4733 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4741 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4746 cout <<
" dr = " << work.
a()[0]
4747 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4748 cout <<
" phi0 = " << work.
a()[1]
4749 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4750 cout <<
" PT = " << 1/work.
a()[2]
4751 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4752 cout <<
" dz = " << work.
a()[3]
4753 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4754 cout <<
" tanl = " << work.
a()[4]
4755 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4762 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4768 DataObject *aRecKalEvent;
4769 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4770 if(aRecKalEvent!=NULL) {
4772 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4773 if(kalsc != StatusCode::SUCCESS) {
4774 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4779 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4780 if( kalsc.isFailure()) {
4781 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4784 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4790 DataObject *aRecKalSegEvent;
4791 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4792 if(aRecKalSegEvent!=NULL) {
4794 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4795 if(segsc != StatusCode::SUCCESS) {
4796 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4801 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4802 if( segsc.isFailure() ) {
4803 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4806 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4809 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.);
4810 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4813 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4815 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4818 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4819 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4820 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4821 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4822 <<
"Track Id: " << (*iter_trk)->getTrackId()
4823 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4824 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4825 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4826 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4827 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4828 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4829 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4830 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4835 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4838 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4839 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4841 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4842 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4843 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4844 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4845 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4846 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4849 for(
int i = 0; i<43; i++) {
4850 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4851 << (*iter_trk)->getPathl(i) <<endreq;
4855 m_trackid = (*iter_trk)->getTrackId();
4856 for(
int jj =0, iii=0; jj<5; jj++) {
4857 m_length[jj] = (*iter_trk)->getLength(jj);
4858 m_tof[jj] = (*iter_trk)->getTof(jj);
4859 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4860 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4861 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4862 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4863 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4864 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4865 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4866 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4867 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4868 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4869 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4870 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4871 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4872 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4873 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4874 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4876 for(
int kk=0; kk<=jj; kk++,iii++) {
4877 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4878 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4879 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4880 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4881 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4882 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4883 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4884 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4885 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4886 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4887 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4888 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4889 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4890 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4891 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4932 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4933 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4934 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4935 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4936 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4937 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4938 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4939 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4940 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4941 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4943 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4944 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4945 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4946 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4947 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4948 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4949 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4950 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4951 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4952 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4954 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4955 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4956 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4957 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4958 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4959 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4960 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4961 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4962 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4963 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4965 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4966 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4967 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4968 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4969 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4971 m_zpt = 1/m_zhelix[2];
4972 m_zpte = 1/m_zhelixe[2];
4973 m_zptmu = 1/m_zhelixmu[2];
4974 m_zptk = 1/m_zhelixk[2];
4975 m_zptp = 1/m_zhelixp[2];
4977 m_fpt = 1/m_fhelix[2];
4978 m_fpte = 1/m_fhelixe[2];
4979 m_fptmu = 1/m_fhelixmu[2];
4980 m_fptk = 1/m_fhelixk[2];
4981 m_fptp = 1/m_fhelixp[2];
4983 m_lpt = 1/m_lhelix[2];
4984 m_lpte = 1/m_lhelixe[2];
4985 m_lptmu = 1/m_lhelixmu[2];
4986 m_lptk = 1/m_lhelixk[2];
4987 m_lptp = 1/m_lhelixp[2];
4989 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4990 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4991 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4992 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4993 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4995 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4996 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4997 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4998 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4999 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5001 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5002 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5003 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5004 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5005 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5006 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5007 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5008 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5009 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5010 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5011 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5012 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5013 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5014 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5015 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5018 StatusCode sc1 = m_nt1->write();
5019 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5024 phi1 = (*iter_trk)->getFFi0();
5025 r1 = (*iter_trk)->getFDr();
5026 z1 = (*iter_trk)->getFDz();
5027 kap1 = (*iter_trk)->getFCpa();
5028 tanl1 = (*iter_trk)->getFTanl();
5031 p1 = sqrt(1+tanl1*tanl1)/kap1;
5032 the1 =
M_PI/2-atan(tanl1);
5033 }
else if(jj == 2) {
5034 phi2 = (*iter_trk)->getFFi0();
5035 r2 = (*iter_trk)->getFDr();
5036 z2 = (*iter_trk)->getFDz();
5037 kap2 = (*iter_trk)->getFCpa();
5038 tanl2 = (*iter_trk)->getFTanl();
5041 p2 = sqrt(1+tanl2*tanl2)/kap1;
5042 the2 =
M_PI/2-atan(tanl2);
5050 m_delthe = the1 + the2;
5053 StatusCode sc2 = m_nt2->write();
5054 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5062 cout <<
"Kalfitting finished " << std::endl;
5069 MsgStream log(
msgSvc(), name());
5070 double Pt_threshold(0.3);
5071 Hep3Vector IP(0,0,0);
5078 if ( !&whMgr )
return;
5081 int ntrk = mdcMgr->size();
5084 int nhits = whMgr->size();
5089 DataObject *aReconEvent;
5090 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5094 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5095 if(sc!=StatusCode::SUCCESS) {
5096 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5104 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5116 if ((TrasanTRK_add.
decision & 32) == 32 ||
5117 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5122 TrasanTRK.
pivot[2]);
5125 for(
int i = 0; i < 5; i++)
5126 a[i] = TrasanTRK.
helix[i];
5129 for(
int i = 0, k = 0; i < 5; i++) {
5130 for(
int j = 0; j <= i; j++) {
5132 ea[j][i] = ea[i][j];
5138 double fiTerm = TrasanTRK.
fiTerm;
5146 int trasqual = TrasanTRK_add.
quality;
5147 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5148 if (trasqual)
return;
5150 int inlyr(999), outlyr(-1);
5151 int* rStat =
new int[43];
5152 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5153 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
5155 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5158 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5159 ii != pt.begin()-1; ii--) {
5160 Num[(*ii)->geo->Lyr()->Id()]++;
5163 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5164 ii != pt.begin()-1; ii--) {
5167 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5169 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5170 <<
" hits in the layer "
5171 << (*ii)->geo->Lyr()->Id() << std::endl;
5177 double dist[2] = {rechit.
ddl, rechit.
ddr};
5178 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5183 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5185 else if (rechit.
lr==1) lr_decision=1;
5188 int ind(geo->
Lyr()->
Id());
5190 lr_decision, rechit.
tdc,
5195 if (inlyr>ind) inlyr = ind;
5196 if (outlyr<ind) outlyr = ind;
5199 int empty_between(0), empty(0);
5200 for (
int i= inlyr; i <= outlyr; i++)
5201 if (!rStat[i]) empty_between++;
5202 empty = empty_between+inlyr+(42-outlyr);
5206 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5211 track_lead.
type(type);
5212 unsigned int nhit = track_lead.
HitsMdc().size();
5214 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5219 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5221 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5224 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5226 track_lead.
pivot(outer_pivot);
5231 HepSymMatrix Eakal(5,0);
5235 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5245 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5251 std::cout <<
" dr = " << work.
a()[0]
5252 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5253 std::cout <<
" phi0 = " << work.
a()[1]
5254 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5255 std::cout <<
" PT = " << 1/work.
a()[2]
5256 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5257 std::cout <<
" dz = " << work.
a()[3]
5258 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5259 std::cout <<
" tanl = " << work.
a()[4]
5260 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5267 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5272 cout <<
" dr = " << work1.
a()[0]
5273 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5274 cout <<
" phi0 = " << work1.
a()[1]
5275 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5276 cout <<
" PT = " << 1/work1.
a()[2]
5277 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5278 cout <<
" dz = " << work1.
a()[3]
5279 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5280 cout <<
" tanl = " << work1.
a()[4]
5281 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5288 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5293 DataObject *aRecKalEvent;
5294 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5295 if(aRecKalEvent!=NULL) {
5297 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5298 if(kalsc != StatusCode::SUCCESS) {
5299 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5304 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5305 if( kalsc.isFailure()) {
5306 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5309 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5315 DataObject *aRecKalSegEvent;
5316 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5317 if(aRecKalSegEvent!=NULL) {
5319 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5320 if(segsc != StatusCode::SUCCESS) {
5321 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5326 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5327 if( segsc.isFailure() ) {
5328 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5331 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5334 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.);
5335 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5338 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5340 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5343 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5344 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5345 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5346 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5347 <<
"Track Id: " << (*iter_trk)->getTrackId()
5348 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5349 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5350 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5351 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5352 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5353 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5354 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5355 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5356 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5361 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5364 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5365 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5367 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5368 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5369 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5370 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5371 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5372 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5375 for(
int i = 0; i<43; i++) {
5376 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5377 << (*iter_trk)->getPathl(i) <<endreq;
5381 m_trackid = (*iter_trk)->getTrackId();
5382 for(
int jj =0, iii=0; jj<5; jj++) {
5383 m_length[jj] = (*iter_trk)->getLength(jj);
5384 m_tof[jj] = (*iter_trk)->getTof(jj);
5385 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5386 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5387 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5388 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5389 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5390 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5391 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5392 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5393 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5394 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5395 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5396 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5397 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5398 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5399 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5400 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5402 for(
int kk=0; kk<=jj; kk++,iii++) {
5403 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5404 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5405 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5406 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5407 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5408 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5409 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5410 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5411 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5412 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5413 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5414 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5415 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5416 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5417 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5424 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5425 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5426 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5427 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5428 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5429 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5430 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5431 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5432 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5433 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5435 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5436 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5437 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5438 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5439 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5440 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5441 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5442 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5443 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5444 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5446 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5447 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5448 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5449 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5450 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5451 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5452 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5453 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5454 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5455 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5457 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5458 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5459 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5460 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5461 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5463 m_zpt = 1/m_zhelix[2];
5464 m_zpte = 1/m_zhelixe[2];
5465 m_zptmu = 1/m_zhelixmu[2];
5466 m_zptk = 1/m_zhelixk[2];
5467 m_zptp = 1/m_zhelixp[2];
5469 m_fpt = 1/m_fhelix[2];
5470 m_fpte = 1/m_fhelixe[2];
5471 m_fptmu = 1/m_fhelixmu[2];
5472 m_fptk = 1/m_fhelixk[2];
5473 m_fptp = 1/m_fhelixp[2];
5475 m_lpt = 1/m_lhelix[2];
5476 m_lpte = 1/m_lhelixe[2];
5477 m_lptmu = 1/m_lhelixmu[2];
5478 m_lptk = 1/m_lhelixk[2];
5479 m_lptp = 1/m_lhelixp[2];
5481 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5482 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5483 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5484 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5485 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5487 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5488 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5489 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5490 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5491 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5493 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5494 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5495 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5496 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5497 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5498 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5499 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5500 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5501 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5502 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5503 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5504 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5505 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5506 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5507 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5510 StatusCode sc1 = m_nt1->write();
5511 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5516 phi1 = (*iter_trk)->getFFi0();
5517 r1 = (*iter_trk)->getFDr();
5518 z1 = (*iter_trk)->getFDz();
5519 kap1 = (*iter_trk)->getFCpa();
5520 tanl1 = (*iter_trk)->getFTanl();
5523 p1 = sqrt(1+tanl1*tanl1)/kap1;
5524 the1 =
M_PI/2-atan(tanl1);
5525 }
else if(jj == 2) {
5526 phi2 = (*iter_trk)->getFFi0();
5527 r2 = (*iter_trk)->getFDr();
5528 z2 = (*iter_trk)->getFDz();
5529 kap2 = (*iter_trk)->getFCpa();
5530 tanl2 = (*iter_trk)->getFTanl();
5533 p2 = sqrt(1+tanl2*tanl2)/kap1;
5534 the2 =
M_PI/2-atan(tanl2);
5542 m_delthe = the1 + the2;
5545 StatusCode sc2 = m_nt2->write();
5546 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5550 cout <<
"Kalfitting finished " << std::endl;
5558 MsgStream log(
msgSvc(), name());
5559 double Pt_threshold(0.3);
5560 Hep3Vector IP(0,0,0);
5567 if ( !&whMgr )
return;
5570 int ntrk = mdcMgr->size();
5573 int nhits = whMgr->size();
5577 double* rY =
new double[ntrk];
5578 double* rfiTerm =
new double[ntrk];
5579 double* rPt =
new double[ntrk];
5580 int* rOM =
new int[ntrk];
5581 unsigned int* order =
new unsigned int[ntrk];
5582 unsigned int* rCont =
new unsigned int[ntrk];
5583 unsigned int* rGen =
new unsigned int[ntrk];
5586 Hep3Vector csmp3[2];
5587 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
5588 end = mdcMgr->end(); it != end; it++) {
5590 rfiTerm[index]=it->fiTerm;
5595 rPt[index] = 1 / fabs(it->helix[2]);
5596 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
5597 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
5599 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
5600 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
5602 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5603 ii !=pt.begin()-1; ii--) {
5604 int lyr((*ii)->geo->Lyr()->Id());
5605 if (outermost < lyr) {
5607 rY[index] = (*ii)->geo->Forward().y();
5609 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
5611 rOM[index] = outermost;
5612 order[index] = index;
5617 for (
int j, k = ntrk - 1; k >= 0; k = j){
5619 for(
int i = 1; i <= k; i++)
5620 if(rY[i - 1] < rY[i]){
5622 std::swap(order[i], order[j]);
5623 std::swap(rY[i], rY[j]);
5624 std::swap(rOM[i], rOM[j]);
5625 std::swap(rCont[i], rCont[j]);
5626 std::swap(rGen[i], rGen[j]);
5635 DataObject *aReconEvent;
5636 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5640 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5641 if(sc!=StatusCode::SUCCESS) {
5642 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5650 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5657 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
5666 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
5670 if ((TrasanTRK_add.
decision & 32) == 32 ||
5671 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5676 TrasanTRK.
pivot[2]);
5679 for(
int i = 0; i < 5; i++)
5680 a[i] = TrasanTRK.
helix[i];
5683 for(
int i = 0, k = 0; i < 5; i++) {
5684 for(
int j = 0; j <= i; j++) {
5686 ea[j][i] = ea[i][j];
5692 double fiTerm = TrasanTRK.
fiTerm;
5699 for(
int l = 0; l < ntrk; l++) {
5700 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
5701 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
5703 int trasqual = TrasanTRK_add1.
quality;
5704 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5705 if (trasqual)
continue;
5707 int inlyr(999), outlyr(-1);
5708 int* rStat =
new int[43];
5709 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5710 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
5712 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5715 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5716 ii != pt.begin()-1; ii--) {
5717 Num[(*ii)->geo->Lyr()->Id()]++;
5720 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5721 ii != pt.begin()-1; ii--) {
5724 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5726 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5727 <<
" hits in the layer "
5728 << (*ii)->geo->Lyr()->Id() << std::endl;
5734 double dist[2] = {rechit.
ddl, rechit.
ddr};
5735 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5740 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5742 else if (rechit.
lr==1) lr_decision=1;
5745 int ind(geo->
Lyr()->
Id());
5747 lr_decision, rechit.
tdc,
5752 if (inlyr>ind) inlyr = ind;
5753 if (outlyr<ind) outlyr = ind;
5756 int empty_between(0), empty(0);
5757 for (
int i= inlyr; i <= outlyr; i++)
5758 if (!rStat[i]) empty_between++;
5759 empty = empty_between+inlyr+(42-outlyr);
5763 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5768 track_lead.
type(type);
5769 unsigned int nhit = track_lead.
HitsMdc().size();
5771 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5776 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5778 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5781 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5783 track_lead.
pivot(outer_pivot);
5788 HepSymMatrix Eakal(5,0);
5792 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5802 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5808 std::cout <<
" dr = " << work.
a()[0]
5809 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5810 std::cout <<
" phi0 = " << work.
a()[1]
5811 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5812 std::cout <<
" PT = " << 1/work.
a()[2]
5813 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5814 std::cout <<
" dz = " << work.
a()[3]
5815 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5816 std::cout <<
" tanl = " << work.
a()[4]
5817 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5824 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5829 cout <<
" dr = " << work1.
a()[0]
5830 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5831 cout <<
" phi0 = " << work1.
a()[1]
5832 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5833 cout <<
" PT = " << 1/work1.
a()[2]
5834 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5835 cout <<
" dz = " << work1.
a()[3]
5836 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5837 cout <<
" tanl = " << work1.
a()[4]
5838 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5845 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5851 DataObject *aRecKalEvent;
5852 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5853 if(aRecKalEvent!=NULL) {
5855 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5856 if(kalsc != StatusCode::SUCCESS) {
5857 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5862 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5863 if( kalsc.isFailure()) {
5864 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5867 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5873 DataObject *aRecKalSegEvent;
5874 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5875 if(aRecKalSegEvent!=NULL) {
5877 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5878 if(segsc != StatusCode::SUCCESS) {
5879 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5884 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5885 if( segsc.isFailure() ) {
5886 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5889 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5892 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.);
5893 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5896 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5898 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5901 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5902 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5903 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5904 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5905 <<
"Track Id: " << (*iter_trk)->getTrackId()
5906 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5907 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5908 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5909 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5910 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5911 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5912 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5913 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5914 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5919 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5922 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5923 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5925 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5926 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5927 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5928 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5929 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5930 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5933 for(
int i = 0; i<43; i++) {
5934 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5935 << (*iter_trk)->getPathl(i) <<endreq;
5939 m_trackid = (*iter_trk)->getTrackId();
5940 for(
int jj =0, iii=0; jj<5; jj++) {
5941 m_length[jj] = (*iter_trk)->getLength(jj);
5942 m_tof[jj] = (*iter_trk)->getTof(jj);
5943 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5944 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5945 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5946 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5947 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5948 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5949 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5950 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5951 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5952 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5953 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5954 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5955 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5956 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5957 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5958 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5960 for(
int kk=0; kk<=jj; kk++,iii++) {
5961 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5962 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5963 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5964 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5965 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5966 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5967 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5968 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5969 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5970 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5971 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5972 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5973 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5974 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5975 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5982 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5983 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5984 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5985 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5986 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5987 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5988 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5989 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5990 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5991 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5993 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5994 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5995 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5996 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5997 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5998 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5999 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
6000 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
6001 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
6002 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
6004 m_stat[0][0] = (*iter_trk)->getStat(0,0);
6005 m_stat[1][0] = (*iter_trk)->getStat(0,1);
6006 m_stat[2][0] = (*iter_trk)->getStat(0,2);
6007 m_stat[3][0] = (*iter_trk)->getStat(0,3);
6008 m_stat[4][0] = (*iter_trk)->getStat(0,4);
6009 m_stat[0][1] = (*iter_trk)->getStat(1,0);
6010 m_stat[1][1] = (*iter_trk)->getStat(1,1);
6011 m_stat[2][1] = (*iter_trk)->getStat(1,2);
6012 m_stat[3][1] = (*iter_trk)->getStat(1,3);
6013 m_stat[4][1] = (*iter_trk)->getStat(1,4);
6015 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
6016 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
6017 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
6018 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
6019 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
6021 m_zpt = 1/m_zhelix[2];
6022 m_zpte = 1/m_zhelixe[2];
6023 m_zptmu = 1/m_zhelixmu[2];
6024 m_zptk = 1/m_zhelixk[2];
6025 m_zptp = 1/m_zhelixp[2];
6027 m_fpt = 1/m_fhelix[2];
6028 m_fpte = 1/m_fhelixe[2];
6029 m_fptmu = 1/m_fhelixmu[2];
6030 m_fptk = 1/m_fhelixk[2];
6031 m_fptp = 1/m_fhelixp[2];
6033 m_lpt = 1/m_lhelix[2];
6034 m_lpte = 1/m_lhelixe[2];
6035 m_lptmu = 1/m_lhelixmu[2];
6036 m_lptk = 1/m_lhelixk[2];
6037 m_lptp = 1/m_lhelixp[2];
6039 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
6040 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
6041 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
6042 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
6043 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
6045 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
6046 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
6047 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
6048 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
6049 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
6051 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
6052 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
6053 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
6054 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
6055 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
6056 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
6057 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
6058 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
6059 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
6060 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
6061 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
6062 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
6063 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
6064 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
6065 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
6068 StatusCode sc1 = m_nt1->write();
6069 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
6074 phi1 = (*iter_trk)->getFFi0();
6075 r1 = (*iter_trk)->getFDr();
6076 z1 = (*iter_trk)->getFDz();
6077 kap1 = (*iter_trk)->getFCpa();
6078 tanl1 = (*iter_trk)->getFTanl();
6081 p1 = sqrt(1+tanl1*tanl1)/kap1;
6082 the1 =
M_PI/2-atan(tanl1);
6083 }
else if(jj == 2) {
6084 phi2 = (*iter_trk)->getFFi0();
6085 r2 = (*iter_trk)->getFDr();
6086 z2 = (*iter_trk)->getFDz();
6087 kap2 = (*iter_trk)->getFCpa();
6088 tanl2 = (*iter_trk)->getFTanl();
6091 p2 = sqrt(1+tanl2*tanl2)/kap1;
6092 the2 =
M_PI/2-atan(tanl2);
6100 m_delthe = the1 + the2;
6103 StatusCode sc2 = m_nt2->write();
6104 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
6112 cout <<
"Kalfitting finished " << std::endl;
6125 MsgStream log(
msgSvc(), name());
6143 double pp = track_lead.
momentum().mag();
6147 for(
int l_mass = 0, flg = 1; l_mass < nmass;
6148 l_mass++, flg <<= 1) {
6150 if (!(
mhyp_ & flg))
continue;
6151 if (l_mass ==
lead_)
continue;
6160 if(
debug_ == 4) cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
6168 TrasanTRK.
pivot[2]);
6169 HepVector a_trasan(5);
6170 for(
int i = 0; i < 5; i++)
6171 a_trasan[i] = TrasanTRK.
helix[i];
6173 HepSymMatrix ea_trasan(5);
6174 for(
int i = 0, k = 0; i < 5; i++) {
6175 for(
int j = 0; j <= i; j++) {
6177 ea_trasan[j][i] = ea_trasan[i][j];
6181 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
6184 double fiTerm = TrasanTRK.
fiTerm;
6185 track.
pivot(track.
x(fiTerm));
6186 HepSymMatrix Eakal(5,0);
6188 double costheta = track.
a()[4] / sqrt(1.0 + track.
a()[4]*track.
a()[4]);
6207 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
6209 fillTds(TrasanTRK, track_first, kaltrk, l_mass);
6217 if(
ntuple_&256) StatusCode sc7 = m_nt7 ->write();
6223 HepSymMatrix ea_first(kaltrk->
getFError());
6226 for(
int i = 0; i < 5; i++)
6227 for(
int j = 0; j <= i; j++)
6228 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
6229 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
6236 cout <<
" Backward fitting flag:" <<
back_<< endl;
6237 cout <<
"track_back pivot " << track_back.
pivot()
6238 <<
" track_lead kappa " << track_lead.
a()[2]
6242 if (
back_ && track_lead.
a()[2] != 0 &&
6243 1/fabs(track_lead.
a()[2]) >
pT_) {
6247 double p_kaon(0), p_proton(0);
6248 if (!(kaltrk->
getStat(0,3))) {
6251 track_back.
p_kaon(p_kaon);
6253 p_kaon = 1 / fabs(track_back.
a()[2]) *
6254 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
6255 track_back.
p_kaon(p_kaon);
6257 if (!(kaltrk->
getStat(0,4))) {
6258 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
6262 p_proton = 1 / fabs(track_back.
a()[2]) *
6263 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
6270 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
6284 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass, segcol, 1);
6311 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
6313 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
6314 <<
"Length of the track: "<< kaltrk->
getLength(2)
6315 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
6316 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
6317 <<
" "<< kaltrk->
getChisq(1,2) << endreq
6318 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
6319 <<
" "<< kaltrk->
getNdf(1,2) << endreq
6324 kalcol->push_back(kaltrk);
6337 MsgStream log(
msgSvc(), name());
6342 cout <<
"track_first pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
6349 cout <<
"after inner wall, track_ip pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
6355 double pp = track_lead.
momentum().mag();
6360 for(
int l_mass = 0, flg = 1; l_mass < nmass;
6361 l_mass++, flg <<= 1) {
6363 if (!(
mhyp_ & flg))
continue;
6364 if (l_mass ==
lead_)
continue;
6375 cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
6385 TrasanTRK.
pivot[2]);
6387 HepVector a_trasan(5);
6388 for(
int i = 0; i < 5; i++){
6389 a_trasan[i] = TrasanTRK.
helix[i];
6392 HepSymMatrix ea_trasan(5);
6393 for(
int i = 0, k = 0; i < 5; i++) {
6394 for(
int j = 0; j <= i; j++) {
6396 ea_trasan[j][i] = ea_trasan[i][j];
6400 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
6403 double fiTerm = TrasanTRK.
fiTerm;
6404 track.
pivot(track.
x(fiTerm));
6406 HepSymMatrix Eakal(5,0);
6408 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
6420 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
6422 fillTds(TrasanTRK, track, kaltrk, l_mass);
6431 HepSymMatrix ea_first(kaltrk->
getFError());
6435 for(
int i = 0; i < 5; i++)
6436 for(
int j = 0; j <= i; j++)
6437 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
6438 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
6454 cout <<
" Backward fitting flag:" <<
back_<< endl;
6455 cout <<
"track_back pivot " << track_back.
pivot()
6456 <<
" track_lead kappa " << track_lead.
a()[2]
6460 if (
back_ && track_lead.
a()[2] != 0 &&
6461 1/fabs(track_lead.
a()[2]) >
pT_) {
6466 double p_kaon(0), p_proton(0);
6468 if (!(kaltrk->
getStat(0,3))) {
6471 track_back.
p_kaon(p_kaon);
6473 p_kaon = 1 / fabs(track_back.
a()[2]) *
6474 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
6475 track_back.
p_kaon(p_kaon);
6477 if (!(kaltrk->
getStat(0,4))) {
6478 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
6482 p_proton = 1 / fabs(track_back.
a()[2]) *
6483 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
6491 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
6496 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass,segcol);
6519 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
6521 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
6522 <<
"Length of the track: "<< kaltrk->
getLength(2)
6523 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
6524 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
6525 <<
" "<< kaltrk->
getChisq(1,2) << endreq
6526 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
6527 <<
" "<< kaltrk->
getNdf(1,2) << endreq
6531 kalcol->push_back(kaltrk);
6540 for (
int i=0; i<5; i++) {
6541 for(
int j = 1; j<i+2;j++) {
6543 Eakal(j,i+1) = Eakal(i+1,j);
6547 if (
debug_ == 4) cout<<
"initialised Ea.. "<<Eakal<<endl;
6555 for (
int i=0; i<5; i++) {
6556 for(
int j = 1; j<i+2;j++) {
6558 Eakal(j,i+1) = Eakal(i+1,j);
6560 Eakal(1,1) = Eakal(1,1)*
gain1_;
6561 Eakal(2,2) = Eakal(2,2)*
gain2_;
6562 Eakal(3,3) = Eakal(3,3)*
gain3_;
6563 Eakal(4,4) = Eakal(4,4)*
gain4_;
6564 Eakal(5,5) = Eakal(5,5)*
gain5_;
6630 if (
debug_ == 4) cout<<
"initialised Eakal.. "<<Eakal<<endl;
6639 cout <<
"start_seed begin... " << std::endl;
6641 Hep3Vector x_init(track.
pivot());
6642 HepSymMatrix Ea_init(5,0);
6643 Ea_init = track.
Ea();
6644 HepVector a_init(5);
6648 unsigned int nhit_included(10);
6660 unsigned int nhit = track.
HitsMdc().size();
6663 for (
int ktrial = 0; ktrial < 8; ktrial++) {
6666 track.
pivot(x_init);
6676 HepSymMatrix Eakal(5,0);
6680 for(
unsigned i=0 ; i < nhit; i++ ){
6683 if (i<3) lr_decision = LR[ktrial][i];
6684 HitMdc.
LR(lr_decision);
6685 if (i<nhit_included)
6698 cout <<
"---- Result for " << ktrial <<
" case : chi2 = " << track.
chiSq()
6699 <<
", nhits included = " << track.
nchits() <<
" for nhits available = "
6700 << nhit << std::endl;
6702 if (track.
chiSq() < chi2_min &&
6703 (track.
nchits() == nhit_included || track.
nchits() == nhit)){
6704 chi2_min = track.
chiSq();
6711 cout <<
"*** i_min = " << i_min <<
" with a chi2 = " << chi2_min << std::endl;
6713 for(
unsigned i=0 ; i < nhit; i++ ){
6716 if (i_min >= 0 && i < 3)
6717 lr_decision = LR[i_min][i];
6718 HitMdc.
LR(lr_decision);
6722 track.
pivot(x_init);
6733 for(
unsigned i=0 ; i < 3; i++ ){
6735 cout <<
" LR(" << i <<
") = " << HitMdc.
LR()
6747 if(
debug_ == 4) cout<<
"Begining to clear Tables ...."<<endl;
6751 vector<MdcRec_trk>::iterator tit=mdcMgr->begin();
6752 for( ; tit != mdcMgr->end(); tit++) {
6753 vector<MdcRec_wirhit*>::iterator vit= tit->hitcol.begin() ;
6754 for(; vit != tit->hitcol.end(); vit++) {
6760 mdc_addMgr->clear();
6772bool KalFitAlg::order_rechits(
const SmartRef<RecMdcHit>& m1,
const SmartRef<RecMdcHit>& m2) {
6853// --- print out GEM hits
6854if(giveCout) cout<<endl<<"------- GEM hits: --------"<<endl;
6855//for(int i=0; i<useNCGem_; i++)
6856for(int i=0; i<4; i++)
6858 if(giveCout) cout<<"### layer "<<i<<endl;
6859 //int nhits = myGemHitCol[i].size();
6860 KalFitGemHitCol::const_iterator it_GemHit = (myGemHitCol[i]).begin();
6861 for(int j=0; it_GemHit!=(myGemHitCol[i]).end(); it_GemHit++, j++)
6863 if(giveCout) cout<<" hit "<<(j+1)<<": phi, z, r = "<<it_GemHit->getPhi()<<", "<<it_GemHit->getZ()<<", "<<it_GemHit->getR()<<endl;
6864 if(ntuple_&1) m_rGem[i] = it_GemHit->getR();
6874 double rMaxInnerWall=6.42005+0.00495;
6878 rMaxInnerWall=m_innerWall[0].radius();
6881 rMaxInnerWall=_BesKalmanFitWalls[0].radius();
6884 HepPoint3D posAtInnerWall = track.
x(dPhiToInnerWall);
6885 double pathToInnerWall;
6886 track.
pivot_numf(posAtInnerWall, pathToInnerWall);
6889 if(pathToInnerWall>0)
6891 if(
muls_) track.
ms(pathToInnerWall, _BesKalmanFitMaterials[0], way);
6892 if(
loss_) track.
eloss(pathToInnerWall, _BesKalmanFitMaterials[0], way);
6894 if(useNCGem_>0) m_innerWall[0].updateTrack(track, way);
7256 for(
int i=0; i<3; i++) myGemHitCol[i].clear();
7278 SmartDataPtr<RecCgemClusterCol> cgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
7281 RecCgemClusterCol::const_iterator it = cgemClusterCol->begin();
7282 for(;it!=cgemClusterCol->end();it++){
7284 if((*it)->getflag()==2){
7285 layer = (*it)->getlayerid();
7286 sheet = (*it)->getsheetid();
7287 phi = (*it)->getrecphi();
7288 v = (*it)->getrecv();
7290 myGemHitCol[layer].push_back(aKalFitGemHit);
7294 cout <<
"Gem information" << endl;
7295 cout <<
"first layer " << myGemHitCol[0].size() << endl;
7296 cout <<
"second layer" << myGemHitCol[1].size() << endl;
7297 cout <<
"third layer" << myGemHitCol[2].size() << endl;
7344 bool prtGemHits=
false;
7346 if(prtGemHits) cout<<endl<<
"------- GEM hits: --------"<<endl;
7347 for(
int i=0; i<3; i++)
7349 if(prtGemHits) cout<<
"### layer "<<i<<endl;
7351 KalFitGemHitCol::const_iterator it_GemHit = (myGemHitCol[i]).begin();
7352 for(
int j=0; it_GemHit!=(myGemHitCol[i]).end(); it_GemHit++, j++)
7354 if(prtGemHits) cout<<
" hit "<<(j+1)<<
": phi, z, r = "<<it_GemHit->getPhi()<<
", "<<it_GemHit->getZ()<<
", "<<it_GemHit->getR()<<endl;
7371 std::map<double,double> esti1_r;
7372 std::map<double,double> esti2_r;
7373 std::map<double,double> esti1_phi;
7374 std::map<double,double> esti2_phi;
7375 std::map<double,double> esti1_z;
7376 std::map<double,double> esti2_z;
7377 std::map<double,double> diff1_phi;
7378 std::map<double,double> diff1_z;
7379 std::map<double,double> diff2_phi;
7380 std::map<double,double> diff2_z;
7381 std::map<double,HepPoint3D> posnew;
7382 std::map<double,HepSymMatrix> eanew;
7383 std::map<double,HepVector> anew;
7384 std::map<double,double> meas_r;
7385 std::map<double,double> meas_phi;
7386 std::map<double,double> meas_z;
7387 std::map<double,double> meas_phierr;
7388 std::map<double,double> esti_phierr;
7389 std::map<double,double> meas_zerr;
7390 std::map<double,double> esti_zerr;
7392 if(
debug_==4){ cout <<
"wall size" << _BesKalmanFitWalls.size() << endl;}
7393 for(
int iLayer=2; iLayer>=0; iLayer--)
7411 meas_phierr.clear();
7413 esti_phierr.clear();
7416 int nHits=myGemHitCol[iLayer].size();
7419 if (nHits <= 0)
continue;
7422 vector<KalFitGemHit>::iterator
iter;
7426 innerwall(track,hypo,way,64+iLayer*(-32),88+iLayer*(-32));
7429 for(
iter=myGemHitCol[iLayer].begin();
iter!=myGemHitCol[iLayer].end();
iter++ )
7432 track.
set(track_start.
pivot(),track_start.
a(),track_start.
Ea());
7434 double rGem=(*iter).getR();
7435 double phiGem=(*iter).getPhi();
7436 double zGem=(*iter).getZ();
7437 HepVector v_measu=(*iter).getVecPos();
7450 double x0=posEstiAtGem.x();
7451 double y0=posEstiAtGem.y();
7459 if(
muls_) track.
ms(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7460 if(
loss_) track.
eloss(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7463 double rGemEsti=posEstiAtGem.perp();
7464 double phiGemEsti=posEstiAtGem.phi();
7465 double zGemEsti=posEstiAtGem.z();
7466 HepVector v_estim(2,0);
7467 v_estim(1)=phiGemEsti;
7468 v_estim(2)=zGemEsti;
7475 const HepSymMatrix& Ea = track.
Ea();
7476 HepVector v_a = track.
a();
7479 double kappa=v_a(3);
7484 HepMatrix
H(2, 5, 0);
7487 H(1,1)=-y0*
cos(phi0)/(y0*y0+x0*x0)+x0*
sin(phi0)/(x0*x0+y0*y0);
7490 H(1,2)=-y0/(y0*y0+x0*x0)*((-1)*drho*
sin(phi0)+
m_alpha/kappa*(
sin(phi0+dPhi)-
sin(phi0)))+x0/(x0*x0+y0*y0)*(drho*
cos(phi0)+
m_alpha/kappa*(
cos(phi0)-
cos(phi0+dPhi)));
7491 H(1,3)=-y0/(y0*y0+x0*x0)*(-1)*
m_alpha/(kappa*kappa)*(
cos(phi0)-
cos(phi0+dPhi))+x0/(x0*x0+y0*y0)*(-1)*
m_alpha/(kappa*kappa)*(
sin(phi0)-
sin(phi0+dPhi));
7492 H(2,3)=
m_alpha/(kappa*kappa)*tl*dPhi;
7496 HepMatrix H_T=
H.T();
7501 HepSymMatrix V=(*iter).getErrMat();
7504 HepMatrix HEaH_T=HEa*H_T;
7506 HepMatrix Vinv=(V+HEaH_T).inverse(ierr);
7507 if(ierr!=0) cout<<
"error in HEaH_T.inverse()!"<<endl;
7509 HepMatrix K=Ea*H_T*Vinv;
7511 HepSymMatrix EaNew(5,0);
7512 EaNew.assign(Ea-K*
H*Ea);
7514 HepVector v_diff=v_measu-v_estim;
7516 double delPhi=v_diff(1);
7517 if(fabs(v_diff(1))>6.28) v_measu(1)=-1*v_measu(1);
7524 HepVector aNew=v_a+K*v_diff;
7530 HepPoint3D posEstiAtGem_new = track_test.
x(dPhiToGem_New);
7531 double phiGemEsti_new=posEstiAtGem_new.phi();
7532 double zGemEsti_new=posEstiAtGem_new.z();
7533 double x0_new = posEstiAtGem_new.x();
7534 double y0_new = posEstiAtGem_new.y();
7535 HepVector v_estim_new(2,0);
7536 v_estim_new(1)=phiGemEsti_new;
7537 v_estim_new(2)=zGemEsti_new;
7541 v_diff=v_measu-v_estim_new;
7549 track_test.
pivot_numf(posEstiAtGem_new, pathInGem);
7550 HepVector v_a_new = track_test.
a();
7552 HepSymMatrix Ea_new = track_test.
Ea();
7554 double phi0_new = v_a_new(2);
7555 double kappa_new=v_a_new(3);
7556 HepMatrix H_new(2, 5, 0);
7565 HepMatrix H_new_T=H_new.T();
7572 HepSymMatrix R(2,0);
7574 R.assign(V-
H*EaNew*H_T);
7582 HepVector v_dChi2=v_diff.T()*R.inverse(ierr)*v_diff;
7583 if(ierr!=0) cout<<
"error in R.inverse()!"<<endl;
7586 double dChi2=v_dChi2(1);
7590 cout<<
"pivot: "<<posEstiAtGem<<endl;
7591 cout<<
"new pivot: "<<posEstiAtGem<<endl;
7592 cout<<
"a: "<<v_a<<endl;
7593 cout<<
"new a: "<<aNew<<endl;
7594 cout<<
"Ea: "<<Ea<<endl;
7595 cout<<
"new Ea: "<<EaNew<<endl;
7596 cout<<
"dchi2= "<<dChi2<<endl;
7631 esti1_r[dChi2]=rGemEsti;
7632 esti1_phi[dChi2]=phiGemEsti;
7633 esti1_z[dChi2]=zGemEsti;
7634 esti2_r[dChi2]=rGem;
7635 esti2_phi[dChi2]=phiGemEsti_new;
7636 esti2_z[dChi2]=zGemEsti_new;
7637 diff1_phi[dChi2]=phiGem-phiGemEsti;
7638 diff1_z[dChi2]=zGem-zGemEsti;
7639 diff2_phi[dChi2]=phiGem-phiGemEsti_new;
7640 diff2_z[dChi2]=zGem-zGemEsti_new;
7641 posnew[dChi2]=posEstiAtGem;
7645 meas_phi[dChi2]=phiGem;
7648 meas_phierr[dChi2]=sqrt(V[0][0]);
7650 meas_zerr[dChi2]=sqrt(V[1][1]);
7652 esti_phierr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[0][0]);
7654 esti_zerr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[1][1]);
7663 track.
set((posnew.begin())->second,(anew.begin())->second,(eanew.begin()->second));
7667 HepPoint3D new_posEstiAtGem = track.
x(new_dPhiToGem);
7668 double new_pathInGem;
7671 track.
pivot_numf(new_posEstiAtGem, new_pathInGem);
7676 if(
muls_) track.
ms(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7677 if(
loss_) track.
eloss(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7681 innerwall(track,hypo,way,90+(-32)*iLayer,95+(-32)*iLayer);}
7683 innerwall(track,hypo,way,90+(-32)*iLayer,94+(-32)*iLayer);}
7691 m_esti1_r[ii]=(esti1_r.begin())->second;
7693 m_esti1_phi[ii]=(esti1_phi.begin())->second;
7694 m_esti1_z[ii]=(esti1_z.begin())->second;
7695 m_esti2_r[ii]=(esti2_r.begin())->second;
7696 m_esti2_phi[ii]=(esti2_phi.begin())->second;
7697 m_esti2_z[ii]=(esti2_z.begin())->second;
7698 m_diff1_phi[ii]=(diff1_phi.begin())->second;
7699 m_diff1_z[ii]=(diff1_z.begin())->second;
7700 m_diff2_phi[ii]=(diff2_phi.begin())->second;
7701 m_diff2_z[ii]=(diff2_z.begin())->second;
7702 m_Gchi2[ii]=(esti1_r.begin())->first;
7703 m_meas_r[ii]=(meas_r.begin())->second;
7704 m_meas_phi[ii]=(meas_phi.begin())->second;
7706 m_meas_z[ii]=(meas_z.begin())->second;
7707 m_meas_phierr[ii]=(meas_phierr.begin())->second;
7708 m_meas_zerr[ii]=(meas_zerr.begin())->second;
7709 m_esti_phierr[ii]=(esti_phierr.begin())->second;
7710 m_esti_zerr[ii]=(esti_zerr.begin())->second;
7711 m_GemLayer[ii]=iLayer;
7723 if(m_meas_phi[0]>0 && m_meas_phi[0]<1.5707963) m_qua=0;
7724 if(m_meas_phi[0]>1.5707963 && m_meas_phi[0]<3.1415926) m_qua=1;
7725 if(m_meas_phi[0]>-3.1415926&& m_meas_phi[0]<-1.5707963) m_qua=2;
7726 if(m_meas_phi[0]>-1.5707963&& m_meas_phi[0]<0) m_qua=3;
7749 meas_phierr.clear();
7751 esti_phierr.clear();
7757 SmartDataPtr<Event::EventHeader> evtHeader(eventSvc(),
"/Event/EventHeader");
7759 m_evt3 = evtHeader->eventNumber();
7771 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
7773 cout <<
"Could not find Event Header" << endl;
7775 int eventNo = eventHeader->eventNumber();
7776 int runNo = eventHeader->runNumber();
7777 SmartDataPtr<RecMdcTrackCol> recMdcTrack(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
7779 cout <<
"Could not find RecMdcTrackCol" << endl;
7781 RecMdcTrackCol::iterator itrk = recMdcTrack->begin();
7783 for(; itrk != recMdcTrack->end(); itrk++){
7784 if((*itrk)->trackId() == track.
trasan_id() && (*itrk)->getNcluster() > 0) {
7785 clusterRefVec = (*itrk)->getVecClusters();
7793 if(
debug_==4){ cout<<
"wall size"<<_BesKalmanFitWalls.size()<<endl;}
7796 ClusterRefVec::iterator itCluster;
7797 ClusterRefVec::iterator itStartFlag;
7798 ClusterRefVec::iterator itStopFlag;
7802 double R_o_cgem=_BesKalmanFitWalls[0].radius();
7804 itStartFlag = clusterRefVec.begin();
7805 itStopFlag = clusterRefVec.end();
7810 itStartFlag = clusterRefVec.end()-1;
7811 itStopFlag = clusterRefVec.begin()-1;
7821 for(itCluster = itStartFlag; itCluster !=itStopFlag; itCluster = itCluster + step){
7822 int layer = (*itCluster)->getlayerid();
7825 innerwall(track,hypo,way,lastR,dmR,index);
7835 m_dchi2_hit_e[m_hit_no] = dchi2;
7836 m_residest_hit_e[m_hit_no] = 0;
7837 m_residnew_hit_e[m_hit_no] = 0;
7838 m_layer_hit_e[m_hit_no] = layer;
7840 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
7841 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
7842 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
7843 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
7844 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
7845 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
7846 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
7847 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
7848 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
7849 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
7854 m_dchi2_hit_mu[m_hit_no] = dchi2;
7855 m_residest_hit_mu[m_hit_no] = 0;
7856 m_residnew_hit_mu[m_hit_no] = 0;
7857 m_layer_hit_mu[m_hit_no] = layer;
7859 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
7860 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
7861 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
7862 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
7863 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
7864 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
7865 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
7866 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
7867 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
7868 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
7874 m_dchi2_hit_pi[m_hit_no] = dchi2;
7875 m_residest_hit_pi[m_hit_no] = 0;
7876 m_residnew_hit_pi[m_hit_no] = 0;
7877 m_layer_hit_pi[m_hit_no] = layer;
7879 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
7880 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
7881 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
7882 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
7883 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
7884 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
7885 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
7886 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
7887 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
7888 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
7893 m_dchi2_hit_k[m_hit_no] = dchi2;
7894 m_residest_hit_k[m_hit_no] = 0;
7895 m_residnew_hit_k[m_hit_no] = 0;
7896 m_layer_hit_k[m_hit_no] = layer;
7898 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
7899 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
7900 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
7901 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
7902 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
7903 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
7904 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
7905 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
7906 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
7907 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
7912 m_dchi2_hit_p[m_hit_no] = dchi2;
7913 m_residest_hit_p[m_hit_no] = 0;
7914 m_residnew_hit_p[m_hit_no] = 0;
7915 m_layer_hit_p[m_hit_no] = layer;
7917 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
7918 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
7919 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
7920 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
7921 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
7922 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
7923 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
7924 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
7925 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
7926 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
7934 innerwall(track,hypo,way,lastR,R_o_cgem,index);
7938 innerwall(track,hypo,way,lastR,0,index);
8085 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
8087 cout <<
"Could not find Event Header" << endl;
8089 int eventNo = eventHeader->eventNumber();
8090 int runNo = eventHeader->runNumber();
8094 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
8096 cout <<
"Could not find McParticle" << endl;
8099 HepVector mc_a(5,0);
8100 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
8101 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
8102 if(!(*i_mcTrk)->primaryParticle())
continue;
8103 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
8104 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
8106 double charge = 0.0;
8107 int pid = (*i_mcTrk)->particleProperty();
8109 charge = m_particleTable->particle( pid )->charge();
8110 }
else if ( pid <0 ) {
8111 charge = m_particleTable->particle( -pid )->charge();
8116 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
8121 for(
int j =0; j<5; j++) {
8123 mc_a[j] = mchelix.
a()[j];
8131 m_diff_dr = track.
a()[0]-MChelix.
a()[0];
8132 m_diff_phi0 = track.
a()[1]-MChelix.
a()[1];
8133 m_diff_kappa = track.
a()[2]-MChelix.
a()[2];
8134 m_diff_dz = track.
a()[3]-MChelix.
a()[3];
8135 m_diff_tanl = track.
a()[4]-MChelix.
a()[4];
8136 if(track.
a()[2] != 0 && MChelix.
a()[2] != 0)
8137 m_diff_p = sqrt(1+track.
a()[4]*track.
a()[4])/track.
a()[2] - sqrt(1+MChelix.
a()[4]*MChelix.
a()[4])/MChelix.
a()[2];
8140 StatusCode sc12 = m_nt12->write();
8141 if( sc12.isFailure() ) cout<<
"Ntuple12 filling failed!"<<endl;
8145 bool printwalls =
false;
8147 double material_quantity(0);
8148 int size = _BesKalmanFitWalls.size();
8149 for(
int i=0;i<size;i++)
8151 double ro = _BesKalmanFitWalls[i].radius();
8152 double ri = _BesKalmanFitWalls[i].rmin();
8153 double thick = ro-ri;
8155 double X0 = material.
X0();
8157 double Z = material.
Z();
8159 material_quantity += 100*thick/
X0;
8160 cout<<
"----------------------------------------------------------------------------------------------------"<<endl;
8161 cout<<
"wall"<<i+1<<
" Z:"<<setw(12)<<Z<<endl;
8162 cout<<
"Ri(cm):"<<setw(12)<<ri<<
"Ro(cm):"<<setw(12)<<ro<<
"T(cm):"<<setw(12)<<thick<<
"X0(cm):"<<setw(12)<<
X0<<
"T/X0(%)="<<setw(12)<<100*thick/
X0<<
"accumulation(%):"<<setw(12)<<material_quantity<<endl;
8163 cout<<
"----------------------------------------------------------------------------------------------------"<<endl;
8169 cout<<
"============================================================"<<endl;
8170 cout<<
"total material quantity="<<material_quantity<<endl;
8171 cout<<
"============================================================"<<endl;
double sin(const BesAngle a)
double cos(const BesAngle a)
*******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
HepGeom::Point3D< double > HepPoint3D
const double dor_layer[3]
const double dmr_layer[3]
HepGeom::Point3D< double > HepPoint3D
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
ObjectVector< RecMdcKalHelixSeg > RecMdcKalHelixSegCol
SmartRefVector< RecMdcKalHelixSeg > HelixSegRefVec
ObjectVector< RecMdcKalTrack > RecMdcKalTrackCol
SmartRefVector< RecCgemCluster > ClusterRefVec
SmartRefVector< RecMdcHit > HitRefVec
double getMiddleROfGapD() const
double getOuterROfGapD() const
double getInnerROfGapD() const
CgemGeoLayer * getCgemLayer(int i) const
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)
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 residual(KalFitTrack &track)
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 fromFHitToInnerWall(KalFitTrack &track, int way)
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 fitGemHits(KalFitTrack &track, int hypo, int way)
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.
KalFitTrack Cgem_filter_anal(KalFitTrack &track, int hypo, int way)
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
double get_density(void) const
double X0(void) const
Extractor.
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
int trasan_id(void) const
const HepPoint3D & pivot_numf(const HepPoint3D &newPivot)
Sets pivot position in a given mag field.
void ms(double path, const KalFitMaterial &m, int index)
static int Tof_correc_
Flag for TOF correction.
double intersect_cylinder(double r) const
Intersection with different geometry.
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 set(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
sets helix pivot position, parameters, and error matrix.
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)