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"
43#include "KalFitAlg/helix/Helix.h"
45#include "KalFitAlg/lpav/Lpav.h"
47#include "KalFitAlg/coil/Bfield.h"
48#include "KalFitAlg/KalFitTrack.h"
49#include "KalFitAlg/KalFitHitMdc.h"
50#include "KalFitAlg/KalFitHelixSeg.h"
51#include "KalFitAlg/KalFitElement.h"
52#include "KalFitAlg/KalFitAlg.h"
53#include "McTruth/McParticle.h"
54#include "McTruth/MdcMcHit.h"
55#include "EventModel/EventHeader.h"
56#include "EvTimeEvent/RecEsTime.h"
57#include "ReconEvent/ReconEvent.h"
58#include "MdcRawEvent/MdcDigi.h"
59#include "MdcRecEvent/RecMdcHit.h"
60#include "MdcRecEvent/RecMdcTrack.h"
61#include "MdcRecEvent/RecMdcKalHelixSeg.h"
62#include "MdcRecEvent/RecMdcKalTrack.h"
63#include "MdcGeomSvc/MdcGeomSvc.h"
64#include "MagneticField/IMagneticFieldSvc.h"
65#include "MagneticField/MagneticFieldSvc.h"
67#include "VertexFit/IVertexDbSvc.h"
69#include "Identifier/Identifier.h"
70#include "Identifier/MdcID.h"
71#include "GaudiKernel/IPartPropSvc.h"
72#include "GaudiKernel/INTupleSvc.h"
73#include "RawEvent/RawDataUtil.h"
74#include "McTruth/CgemMcHit.h"
75#include "CgemRecEvent/RecCgemCluster.h"
76#include "CgemGeomSvc/ICgemGeomSvc.h"
77#include "CgemGeomSvc/CgemGeomSvc.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;
329 return StatusCode::SUCCESS;
335 MsgStream log(
msgSvc(), name());
336 log << MSG::INFO <<
"in beginRun()" << endreq;
337 log << MSG::INFO <<
"Present Parameters: muls: " <<
muls_ <<
" loss: "
346 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
347 if(sc!=StatusCode::SUCCESS) {
348 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
357 std::cout<<
" beginRun, referField from MagneticFieldSvc:"<< (IMFSvc->
getReferField())*10000 <<std::endl;
362 return StatusCode::SUCCESS;
372 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101");
374 if ( nt1 ) m_nt1 = nt1;
376 m_nt1=
ntupleSvc()->book(
"FILE104/n101",CLID_ColumnWiseTuple,
"KalFit");
379 status = m_nt1->addItem(
"trackid",m_trackid);
380 status = m_nt1->addItem(
"stat",5,2,m_stat);
381 status = m_nt1->addItem(
"ndf",5,2,m_ndf);
382 status = m_nt1->addItem(
"chisq",5,2,m_chisq);
383 status = m_nt1->addItem(
"length",5,m_length);
384 status = m_nt1->addItem(
"tof",5,m_tof);
385 status = m_nt1->addItem(
"nhits",5,m_nhits);
386 status = m_nt1->addItem(
"zhelix",5,m_zhelix);
387 status = m_nt1->addItem(
"zhelixe",5,m_zhelixe);
388 status = m_nt1->addItem(
"zhelixmu",5,m_zhelixmu);
389 status = m_nt1->addItem(
"zhelixk",5,m_zhelixk);
390 status = m_nt1->addItem(
"zhelixp",5,m_zhelixp);
391 status = m_nt1->addItem(
"zptot",m_zptot);
392 status = m_nt1->addItem(
"zptote",m_zptote);
393 status = m_nt1->addItem(
"zptotmu",m_zptotmu);
394 status = m_nt1->addItem(
"zptotk",m_zptotk);
395 status = m_nt1->addItem(
"zptotp",m_zptotp);
397 status = m_nt1->addItem(
"zpt",m_zpt);
398 status = m_nt1->addItem(
"zpte",m_zpte);
399 status = m_nt1->addItem(
"zptmu",m_zptmu);
400 status = m_nt1->addItem(
"zptk",m_zptk);
401 status = m_nt1->addItem(
"zptp",m_zptp);
403 status = m_nt1->addItem(
"fptot",m_fptot);
404 status = m_nt1->addItem(
"fptote",m_fptote);
405 status = m_nt1->addItem(
"fptotmu",m_fptotmu);
406 status = m_nt1->addItem(
"fptotk",m_fptotk);
407 status = m_nt1->addItem(
"fptotp",m_fptotp);
408 status = m_nt1->addItem(
"fpt",m_fpt);
409 status = m_nt1->addItem(
"fpte",m_fpte);
410 status = m_nt1->addItem(
"fptmu",m_fptmu);
411 status = m_nt1->addItem(
"fptk",m_fptk);
412 status = m_nt1->addItem(
"fptp",m_fptp);
414 status = m_nt1->addItem(
"lptot",m_lptot);
415 status = m_nt1->addItem(
"lptote",m_lptote);
416 status = m_nt1->addItem(
"lptotmu",m_lptotmu);
417 status = m_nt1->addItem(
"lptotk",m_lptotk);
418 status = m_nt1->addItem(
"lptotp",m_lptotp);
419 status = m_nt1->addItem(
"lpt",m_lpt);
420 status = m_nt1->addItem(
"lpte",m_lpte);
421 status = m_nt1->addItem(
"lptmu",m_lptmu);
422 status = m_nt1->addItem(
"lptk",m_lptk);
423 status = m_nt1->addItem(
"lptp",m_lptp);
425 status = m_nt1->addItem(
"zsigp",m_zsigp);
426 status = m_nt1->addItem(
"zsigpe",m_zsigpe);
427 status = m_nt1->addItem(
"zsigpmu",m_zsigpmu);
428 status = m_nt1->addItem(
"zsigpk",m_zsigpk);
429 status = m_nt1->addItem(
"zsigpp",m_zsigpp);
430 status = m_nt1->addItem(
"fhelix",5,m_fhelix);
431 status = m_nt1->addItem(
"fhelixe",5,m_fhelixe);
432 status = m_nt1->addItem(
"fhelixmu",5,m_fhelixmu);
433 status = m_nt1->addItem(
"fhelixk",5,m_fhelixk);
434 status = m_nt1->addItem(
"fhelixp",5,m_fhelixp);
435 status = m_nt1->addItem(
"lhelix",5,m_lhelix);
436 status = m_nt1->addItem(
"lhelixe",5,m_lhelixe);
437 status = m_nt1->addItem(
"lhelixmu",5,m_lhelixmu);
438 status = m_nt1->addItem(
"lhelixk",5,m_lhelixk);
439 status = m_nt1->addItem(
"lhelixp",5,m_lhelixp);
440 status = m_nt1->addItem(
"rGem",4,m_rGem);
441 status = m_nt1->addItem(
"chi2Gem",4,m_chi2Gem);
442 status = m_nt1->addItem(
"phiGemExp",4,m_phiGemExp);
443 status = m_nt1->addItem(
"phiGemHit",4,m_phiGemHit);
444 status = m_nt1->addItem(
"zGemExp",4,m_zGemExp);
445 status = m_nt1->addItem(
"zGemHit",4,m_zGemHit);
447 status = m_nt1->addItem(
"zerror",15,m_zerror);
448 status = m_nt1->addItem(
"zerrore",15,m_zerrore);
449 status = m_nt1->addItem(
"zerrormu",15,m_zerrormu);
450 status = m_nt1->addItem(
"zerrork",15,m_zerrork);
451 status = m_nt1->addItem(
"zerrorp",15,m_zerrorp);
452 status = m_nt1->addItem(
"ferror",15,m_ferror);
453 status = m_nt1->addItem(
"ferrore",15,m_ferrore);
454 status = m_nt1->addItem(
"ferrormu",15,m_ferrormu);
455 status = m_nt1->addItem(
"ferrork",15,m_ferrork);
456 status = m_nt1->addItem(
"ferrorp",15,m_ferrorp);
457 status = m_nt1->addItem(
"lerror",15,m_lerror);
458 status = m_nt1->addItem(
"lerrore",15,m_lerrore);
459 status = m_nt1->addItem(
"lerrormu",15,m_lerrormu);
460 status = m_nt1->addItem(
"lerrork",15,m_lerrork);
461 status = m_nt1->addItem(
"lerrorp",15,m_lerrorp);
464 status = m_nt1->addItem(
"evtid",m_evtid);
465 status = m_nt1->addItem(
"mchelix",5,m_mchelix);
466 status = m_nt1->addItem(
"mcptot",m_mcptot);
467 status = m_nt1->addItem(
"mcpid",m_mcpid);
469 if( status.isFailure() ) cout<<
"Ntuple1 add item failed!"<<endl;
475 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102");
477 if ( nt2 ) m_nt2 = nt2;
479 m_nt2=
ntupleSvc()->book(
"FILE104/n102",CLID_ColumnWiseTuple,
"KalFitComp");
481 status2 = m_nt2->addItem(
"delx",m_delx);
482 status2 = m_nt2->addItem(
"dely",m_dely);
483 status2 = m_nt2->addItem(
"delz",m_delz);
484 status2 = m_nt2->addItem(
"delthe",m_delthe);
485 status2 = m_nt2->addItem(
"delphi",m_delphi);
486 status2 = m_nt2->addItem(
"delp",m_delp);
487 status2 = m_nt2->addItem(
"delpx",m_delpx);
488 status2 = m_nt2->addItem(
"delpy",m_delpy);
489 status2 = m_nt2->addItem(
"delpz",m_delpz);
491 if( status2.isFailure() ) cout<<
"Ntuple2 add item failed!"<<endl;
497 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103");
499 if ( nt3 ) m_nt3 = nt3;
501 m_nt3=
ntupleSvc()->book(
"FILE104/n103",CLID_ColumnWiseTuple,
"PatRec");
503 status3 = m_nt3->addItem(
"trkhelix",5,m_trkhelix);
504 status3 = m_nt3->addItem(
"trkptot",m_trkptot);
506 status3 = m_nt3->addItem(
"trkerror",15,m_trkerror);
507 status3 = m_nt3->addItem(
"trksigp",m_trksigp);
509 status3 = m_nt3->addItem(
"trkndf",m_trkndf);
510 status3 = m_nt3->addItem(
"trkchisq",m_trkchisq);
511 if( status3.isFailure() ) cout<<
"Ntuple3 add item failed!"<<endl;
517 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104");
519 if ( nt4 ) m_nt4 = nt4;
521 m_nt4=
ntupleSvc()->book(
"FILE104/n104",CLID_ColumnWiseTuple,
"PatRecComp");
523 status4 = m_nt4->addItem(
"trkdelx",m_trkdelx);
524 status4 = m_nt4->addItem(
"trkdely",m_trkdely);
525 status4 = m_nt4->addItem(
"trkdelz",m_trkdelz);
526 status4 = m_nt4->addItem(
"trkdelthe",m_trkdelthe);
527 status4 = m_nt4->addItem(
"trkdelphi",m_trkdelphi);
528 status4 = m_nt4->addItem(
"trkdelp",m_trkdelp);
529 if( status4.isFailure() ) cout<<
"Ntuple4 add item failed!"<<endl;
534 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105");
536 if ( nt5 ) m_nt5 = nt5;
538 m_nt5=
ntupleSvc()->book(
"FILE104/n105",CLID_ColumnWiseTuple,
"KalFitdChisq");
540 status5 = m_nt5->addItem(
"dchi2",m_dchi2);
541 status5 = m_nt5->addItem(
"masshyp",m_masshyp);
542 status5 = m_nt5->addItem(
"residual_estim",m_residest);
543 status5 = m_nt5->addItem(
"residual",m_residnew);
544 status5 = m_nt5->addItem(
"layer",m_layer);
545 status5 = m_nt5->addItem(
"kaldr",m_anal_dr);
546 status5 = m_nt5->addItem(
"kalphi0",m_anal_phi0);
547 status5 = m_nt5->addItem(
"kalkappa",m_anal_kappa);
548 status5 = m_nt5->addItem(
"kaldz",m_anal_dz);
549 status5 = m_nt5->addItem(
"kaltanl",m_anal_tanl);
550 status5 = m_nt5->addItem(
"dr_ea",m_anal_ea_dr);
551 status5 = m_nt5->addItem(
"phi0_ea",m_anal_ea_phi0);
552 status5 = m_nt5->addItem(
"kappa_ea",m_anal_ea_kappa);
553 status5 = m_nt5->addItem(
"dz_ea",m_anal_ea_dz);
554 status5 = m_nt5->addItem(
"tanl_ea",m_anal_ea_tanl);
555 if( status5.isFailure() ) cout<<
"Ntuple5 add item failed!"<<endl;
560 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
562 if ( nt6 ) m_nt6 = nt6;
564 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
566 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
567 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
568 status6 = m_nt6->addItem(
"tdr",m_tdrift);
569 status6 = m_nt6->addItem(
"layerid", m_layerid);
570 status6 = m_nt6->addItem(
"event", m_eventNo);
571 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
572 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
573 status6 = m_nt6->addItem(
"lr",m_lr);
574 status6 = m_nt6->addItem(
"dd",m_dd);
575 status6 = m_nt6->addItem(
"y",m_yposition);
577 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
582 NTuplePtr nt10(
ntupleSvc(),
"FILE104/n110");
584 if ( nt10 ) m_nt10 = nt10;
586 m_nt10=
ntupleSvc()->book(
"FILE104/n110",CLID_ColumnWiseTuple,
"test");
588 status10 = m_nt10->addItem(
"evt",m_evt3);
589 status10 = m_nt10->addItem(
"qua",m_qua);
590 status10 = m_nt10->addItem(
"nhit",m_nGemHits,0,4000000);
591 status10 = m_nt10->addIndexedItem(
"meas_r",m_nGemHits,m_meas_r);
592 status10 = m_nt10->addIndexedItem(
"meas_phi",m_nGemHits,m_meas_phi);
593 status10 = m_nt10->addIndexedItem(
"meas_z",m_nGemHits,m_meas_z);
594 status10 = m_nt10->addIndexedItem(
"esti1_r",m_nGemHits,m_esti1_r);
595 status10 = m_nt10->addIndexedItem(
"esti1_phi",m_nGemHits,m_esti1_phi);
596 status10 = m_nt10->addIndexedItem(
"esti1_z",m_nGemHits,m_esti1_z);
597 status10 = m_nt10->addIndexedItem(
"esti2_r",m_nGemHits,m_esti2_r);
598 status10 = m_nt10->addIndexedItem(
"esti2_phi",m_nGemHits,m_esti2_phi);
599 status10 = m_nt10->addIndexedItem(
"esti2_z",m_nGemHits,m_esti2_z);
600 status10 = m_nt10->addIndexedItem(
"diff1_phi",m_nGemHits,m_diff1_phi);
601 status10 = m_nt10->addIndexedItem(
"diff1_z",m_nGemHits,m_diff1_z);
602 status10 = m_nt10->addIndexedItem(
"diff2_phi",m_nGemHits,m_diff2_phi);
603 status10 = m_nt10->addIndexedItem(
"diff2_z",m_nGemHits,m_diff2_z);
604 status10 = m_nt10->addIndexedItem(
"layer",m_nGemHits,m_GemLayer);
605 status10 = m_nt10->addIndexedItem(
"mass",m_nGemHits,m_mass);
606 status10 = m_nt10->addIndexedItem(
"dchi2",m_nGemHits,m_Gchi2);
607 status10 = m_nt10->addIndexedItem(
"meas_phierr",m_nGemHits,m_meas_phierr);
608 status10 = m_nt10->addIndexedItem(
"meas_zerr",m_nGemHits,m_meas_zerr);
609 status10 = m_nt10->addIndexedItem(
"esti_phierr",m_nGemHits,m_esti_phierr);
610 status10 = m_nt10->addIndexedItem(
"esti_zerr",m_nGemHits,m_esti_zerr);
611 if( status10.isFailure() ) cout<<
"Ntuple10 add item failed!"<<endl;
617 NTuplePtr nt11(
ntupleSvc(),
"FILE104/n111");
619 if ( nt11 ) m_nt11 = nt11;
621 m_nt11=
ntupleSvc()->book(
"FILE104/n111",CLID_ColumnWiseTuple,
"truth");
623 status11 = m_nt11->addItem(
"evt",m_evt4);
624 status11 = m_nt11->addItem(
"ntruth",m_ntruth,0,400000);
625 status11 = m_nt11->addIndexedItem(
"phi",m_ntruth,m_dtphi);
626 status11 = m_nt11->addIndexedItem(
"z",m_ntruth,m_dtv);
627 status11 = m_nt11->addIndexedItem(
"postphi",m_ntruth,m_dtpostphi);
628 status11 = m_nt11->addIndexedItem(
"postz",m_ntruth,m_dtpostz);
629 status11 = m_nt11->addIndexedItem(
"layer",m_ntruth,m_tlayer);
630 if( status11.isFailure() ) cout<<
"Ntuple11 add item failed!"<<endl;
636 NTuplePtr nt12(
ntupleSvc(),
"FILE104/n112");
638 if ( nt12 ) m_nt12 = nt12;
640 m_nt12=
ntupleSvc()->book(
"FILE104/n112",CLID_ColumnWiseTuple,
"PatRecComp");
642 status12 = m_nt12->addItem(
"evt",m_evt);
643 status12 = m_nt12->addItem(
"track",m_track);
644 status12 = m_nt12->addItem(
"diff_dr",m_diff_dr);
645 status12 = m_nt12->addItem(
"diff_phi0",m_diff_phi0);
646 status12 = m_nt12->addItem(
"diff_kappa",m_diff_kappa);
647 status12 = m_nt12->addItem(
"diff_dz",m_diff_dz);
648 status12 = m_nt12->addItem(
"diff_tanl",m_diff_tanl);
649 status12 = m_nt12->addItem(
"diff_p",m_diff_p);
650 if( status12.isFailure() ) cout<<
"Ntuple12 add item failed!"<<endl;
663 for(layid = 0; layid<2; layid++) {
668 for(layid = 4; layid<12; layid++) {
671 for(layid = 12; layid<20; layid++) {
674 for(layid = 20; layid<43; layid++) {
680 for(layid = 0; layid<2; layid++) {
686 for(layid = 4; layid<12; layid++) {
689 for(layid = 12; layid<20; layid++) {
692 for(layid = 20; layid<43; layid++) {
703 for(layid = 0; layid<2; layid++) {
710 for(layid = 4; layid<12; layid++) {
714 for(layid = 12; layid<20; layid++) {
718 for(layid = 20; layid<43; layid++) {
724 for(layid = 0; layid<43; layid++) {
731 for(layid = 0; layid<2; layid++) {
738 for(layid = 4; layid<12; layid++) {
742 for(layid = 12; layid<20; layid++) {
746 for(layid = 20; layid<43; layid++) {
752 for(layid = 0; layid<43; layid++) {
763 MsgStream log(
msgSvc(), name());
764 log << MSG::INFO <<
"in execute()" << endreq;
765 if(
testOutput) std::cout<<
"begin to deal with EVENT ..."<<(++
eventno)<<std::endl;
811 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
812 if(sc!=StatusCode::SUCCESS) {
813 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
822 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
829 IDataProviderSvc* evtSvc = NULL;
830 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
832 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
834 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
835 return StatusCode::SUCCESS;
839 IDataManagerSvc *dataManSvc;
840 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
841 DataObject *aKalTrackCol;
842 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
843 if(aKalTrackCol != NULL) {
844 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
845 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
848 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
849 if( kalsc.isFailure() ) {
850 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
851 return StatusCode::SUCCESS;
853 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
855 DataObject *aKalHelixSegCol;
856 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
857 if(aKalHelixSegCol != NULL){
858 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
859 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
862 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
863 if( kalsc.isFailure()){
864 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
865 return StatusCode::SUCCESS;
867 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
877 ISvcLocator* svcLocator = Gaudi::svcLocator();
879 StatusCode Cgem_sc=svcLocator->service(
"CgemGeomSvc", ISvc);
881 if (!Cgem_sc.isSuccess()) log<< MSG::INFO <<
"KalFitAlg::execute(): Could not open CGEM geometry file" << endreq;
887 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
890 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
892 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
893 return StatusCode::FAILURE;
895 int eventNo = eventHeader->eventNumber();
896 int runNo = eventHeader->runNumber();
900 if(
testOutput) cout<<endl<<
"$$$$$$$$$$$ run="<<
runNo<<
", evt="<<
eventNo<<
" $$$$$$$$$$$$$$$$$"<<endl<<endl;
903 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
904 if (estimeCol && estimeCol->size()) {
905 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
906 t0 = (*iter_evt)->getTest();
909 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
910 return StatusCode::SUCCESS;
915 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
919 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
920 if (sc!=StatusCode::SUCCESS) {
921 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
922 return StatusCode::SUCCESS;
933 m_evtid = eventHeader->eventNumber();
936 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
938 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
943 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
944 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
945 if(!(*i_mcTrk)->primaryParticle())
continue;
946 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
947 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
948 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
949 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
950 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
953 int pid = (*i_mcTrk)->particleProperty();
954 int mother_id = ((*i_mcTrk)->mother()).particleProperty();
956 charge = m_particleTable->particle( pid )->charge();
957 }
else if ( pid <0 ) {
958 charge = m_particleTable->particle( -pid )->charge();
961 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
964 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
966 Helix mchelix(pos2, mom2, charge);
967 log << MSG::DEBUG <<
"charge of the track " << charge << endreq;
968 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
970 for(
int j =0; j<5; j++) {
971 m_mchelix[j] = mchelix.
a()[j];
974 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
977 cout<<
"MC pid, mother_id = "<<pid<<
", "<<mother_id<<endl;
978 cout<<
" p4 = "<<(*i_mcTrk)->initialFourMomentum()<<endl;
979 cout<<
" start from "<<(*i_mcTrk)->initialPosition()<<endl;
980 cout<<
" end at "<<(*i_mcTrk)->finalPosition()<<endl;
981 cout<<
" Helix: "<<mchelix.
a()<<endl;
982 cout<<
"mc ptot, theta, phi, R = "<<m_mcptot<<
", "<<mom2.theta()/acos(-1.)*180<<
", "<<mom2.phi()/acos(-1.)*180<<
", "<<mchelix.
radius()<<endl;
991 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
993 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
994 return( StatusCode::SUCCESS);
996 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
1001 mtrkadd_mgr->clear();
1005 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
1006 Hep3Vector csmp3[2];
1009 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
1010 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
1012 csmp3[kj-1]=(*iter_trk)->p3();
1013 csmphi[kj-1] = (*iter_trk)->phi();
1017 for(
int j = 0, ij = 0; j<5; j++) {
1018 m_trkhelix[j] = (*iter_trk)->helix()[j];
1020 for(
int k=0; k<=j; k++,ij++) {
1021 m_trkerror[ij] = (*iter_trk)->err()[j][k];
1025 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
1027 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
1028 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
1029 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
1031 m_trkndf = (*iter_trk)->ndof();
1032 m_trkchisq = (*iter_trk)->chi2();
1034 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
1036 StatusCode sc3 = m_nt3->write();
1037 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
1067 log << MSG::DEBUG <<
"retrieved MDC tracks:"
1068 <<
" Nhits " <<(*iter_trk)->getNhits()
1069 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
1071 HitRefVec gothits = (*iter_trk)->getVecHits();
1075 rectrk->
id = (*iter_trk)->trackId();
1076 rectrk->
chiSq = (*iter_trk)->chi2();
1077 rectrk->
ndf = (*iter_trk)->ndof();
1078 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
1079 rectrk->
nhits = (*iter_trk)->getNhits();
1080 rectrk->
nster = (*iter_trk)->nster();
1082 rectrk->
stat = (*iter_trk)->stat();
1083 status_temp = (*iter_trk)->stat();
1085 trkadd->
id = (*iter_trk)->trackId();
1089 trkadd->
body = rectrk;
1090 rectrk->
add = trkadd;
1092 for (
int i=0; i<5; i++) {
1093 rectrk->
helix[i] = (*iter_trk)->helix()[i];
1094 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
1095 for(
int j = 1; j<i+2;j++) {
1096 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
1099 std::sort(gothits.begin(), gothits.end(), order_rechits);
1100 HitRefVec::iterator it_gothit = gothits.begin();
1101 for( ; it_gothit != gothits.end(); it_gothit++) {
1103 if( (*it_gothit)->getStat() != 1 ) {
1105 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
1110 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
1111 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
1112 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
1113 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
1114 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
1115 <<
" id of hit "<<(*it_gothit)->getId()
1116 <<
" track id of hit "<<(*it_gothit)->getTrkId()
1117 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
1120 whit->
id = (*it_gothit)->getId();
1121 whit->
ddl = (*it_gothit)->getDriftDistLeft();
1122 whit->
ddr = (*it_gothit)->getDriftDistRight();
1123 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
1124 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
1125 whit->
pChiSq = (*it_gothit)->getChisqAdd();
1126 whit->
lr = (*it_gothit)->getFlagLR();
1127 whit->
stat = (*it_gothit)->getStat();
1128 mdcid = (*it_gothit)->getMdcId();
1132 int wid = w0id + localwid;
1134 <<
"lr from PR: "<<whit->
lr
1135 <<
" layerId = " << layid
1136 <<
" wireId = " << localwid
1146 whit->
tdc = (*it_gothit)->getTdc();
1147 whit->
adc= (*it_gothit)->getAdc();
1148 rectrk->
hitcol.push_back(whit);
1149 mhit_mgr->push_back(*whit);
1151 mtrk_mgr->push_back(*rectrk);
1152 mtrkadd_mgr->push_back(*trkadd);
1160 m_trkdelx = trkx1 - trkx2;
1161 m_trkdely = trky1 - trky2;
1162 m_trkdelz = trkz1 - trkz2;
1163 m_trkdelthe = trkthe1 + trkthe2;
1164 m_trkdelphi = trkphi1- trkphi2;
1165 m_trkdelp = trkp1 - trkp2;
1166 StatusCode sc4 = m_nt4->write();
1167 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1170 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1171 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1182 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1183 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1188 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1195 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1200 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1202 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++){
1204 for(
int hypo=0; hypo<5; hypo++)
1206 if((*KalTrk)->getStat(0,hypo)==1) nFailedTrks[hypo]++;
1209 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1211 int nhitofthistrk=0;
1212 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1219 iter_hit=gothelixsegs.begin();
1267 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
1269 log << MSG::WARNING <<
"Could not retrieve Cgem MC truth" << endreq;
1270 return StatusCode::FAILURE;
1272 m_evt4=eventHeader->eventNumber();
1273 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
1275 double dprex,dprey,dprez,dpostx,dposty,dpostz;
1277 for(; iter_truth != cgemMcHitCol->end(); ++iter_truth){
1278 layer = (*iter_truth)->GetLayerID();
1279 dprex = (*iter_truth)->GetPositionXOfPrePoint();
1280 dprey = (*iter_truth)->GetPositionYOfPrePoint();
1281 dprez = (*iter_truth)->GetPositionZOfPrePoint();
1282 dpostx = (*iter_truth)->GetPositionXOfPostPoint();
1283 dposty = (*iter_truth)->GetPositionYOfPostPoint();
1284 dpostz = (*iter_truth)->GetPositionZOfPostPoint();
1292 if((dprex >=0&&dprey>=0)||(dprex <0&&dprey>=0)){
1294 m_dtphi[jj] = acos(dprex/diR);
1297 m_dtpostphi[jj] = acos(dpostx/doR);
1298 m_dtpostz[jj]= dpostz;
1300 if((dprex <0&&dprey<0)||(dprex >=0&&dprey<0)){
1302 m_dtphi[jj] = -acos(dprex/diR);
1305 m_dtpostphi[jj] = -acos(dpostx/doR);
1306 m_dtpostz[jj]= dpostz;
1308 midx=(dprex+dpostx)/2;
1309 midy=(dprey+dposty)/2;
1310 if((midx>=0&&midy>=0)||(midx<0&&midy>=0)){
1312 m_dtphi[jj]=acos(midx/dmR);
1314 if((midx<0&&midy<0)||(midx>=0&&midy<0)){
1316 m_dtphi[jj]=-acos(midx/dmR);
1318 m_dtv[jj] = (m_dtv[jj]+m_dtpostz[jj])/20;
1319 m_tlayer[jj] = layer;
1329 return StatusCode::SUCCESS;
1340 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1341 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1343 TrasanTRK.
helix[2]*track.
a()[2]<0)
1347 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1348 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1349 <<
" "<<TrasanTRK.
helix[2]<<
" "<<track.
a()[2]<<
" "<<iqual<<endl;
1350 cout<<
"FillTds> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1351 <<
" nster "<<track.
nster()<<
" iqual "<<iqual<<
" track.Ea "<< track.
Ea()<<endl;
1353 cout<<
"fillTds>.....track.Ea[2][2] "<<track.
Ea()[2][2]<<endl;
1354 cout <<
" TRASAN stereo = " << trasster
1355 <<
" and KalFitTrack = " << trakster << std::endl;
1356 cout <<
" TRASAN axial = " << trasax
1357 <<
" and KalFitTrack = " << trakax << std::endl;
1360 cout <<
"...there is a problem during fit !! " << std::endl;
1361 if (trasster-trakster>5)
1362 cout <<
" because stereo " << trasster-trakster << std::endl;
1363 if (trasax-trakax >5)
1364 cout <<
" because axial " << std::endl;
1365 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1366 cout <<
" because kappa sign " << std::endl;
1372 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1373 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1374 track.
Ea()[4][4] > 0 && iqual) {
1375 if(
debug_ == 4) cout<<
"fillTds>.....going on "<<endl;
1389 if(
debug_) cout<<
"ALARM: FillTds Not refit with KalFilter!!!"<<endl;
1398 double a_trasan[5], ea_trasan[15];
1399 for(
int i =0 ; i <5; i++){
1400 a_trasan[i] = TrasanTRK.
helix[i];
1402 for(
int j =0 ; j <15; j++){
1403 ea_trasan[j] = TrasanTRK.
error[j];
1418 int trasster = TrasanTRK.
nster, trakster = track.
nster(),
1419 trasax(TrasanTRK.
nhits-trasster), trakax(track.
nchits()-trakster);
1421 TrasanTRK.
helix[2]*track.
a()[2]<0)
1425 cout<<
"Nhit from PR "<<TrasanTRK.
nhits<<
" nhit "<<track.
nchits()<<endl;
1426 cout<<
"trasster trakster trasax trakax TrasK trackK iqual"<<endl
1427 <<trasster<<
" "<<trakster<<
" "<<trasax<<
" "<<trakax
1429 cout<<
"FillTds_lead> track.chiSq..."<<track.
chiSq()<<
" nchits "<<track.
nchits()
1430 <<
" nster "<<track.
nster()<<
" iqual_front_[l_mass] "<<
iqual_front_[l_mass]<<
" track.Ea "<<track.
Ea()<<endl;
1432 cout <<
" TRASAN stereo = " << trasster
1433 <<
" and KalFitTrack = " << trakster << std::endl;
1434 cout <<
" TRASAN axial = " << trasax
1435 <<
" and KalFitTrack = " << trakax << std::endl;
1438 cout <<
"...there is a problem during fit !! " << std::endl;
1439 if (trasster-trakster>5)
1440 cout <<
" because stereo " << trasster-trakster << std::endl;
1441 if (trasax-trakax >5)
1442 cout <<
" because axial " << std::endl;
1443 if (TrasanTRK.
helix[2]*track.
a()[2]<0)
1444 cout <<
" because kappa sign " << std::endl;
1450 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1451 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1462 if (
debug_ == 4) cout<<
" trasan id...1 "<<TrasanTRK.
id<<endl;
1471 if(
debug_) cout<<
"ALARM: FillTds_forMdc Not refit with KalFilter!!!"<<endl;
1482 if (
debug_ ==4) cout<<
" trasan id...2 "<<TrasanTRK.
id<<endl;
1486 double a_trasan[5], ea_trasan[15];
1487 for(
int i =0 ; i <5; i++){
1488 a_trasan[i] = TrasanTRK.
helix[i];
1490 for(
int j =0 ; j <15; j++){
1491 ea_trasan[j] = TrasanTRK.
error[j];
1510 cout <<
"fillTds_IP>......"<<endl;
1511 cout <<
" dr = " << track.
a()[0]
1512 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1513 cout <<
" phi0 = " << track.
a()[1]
1514 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1515 cout <<
" PT = " << 1/track.
a()[2]
1516 <<
", Er_kappa =" << sqrt(track.
Ea()[2][2]) << std::endl;
1517 cout <<
" dz = " << track.
a()[3]
1518 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1519 cout <<
" tanl = " << track.
a()[4]
1520 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1524 TrasanTRK.
helix[2]*track.
a()[2]<0)
1530 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1531 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1535 double dr = track.
a()[0];
1536 double phi0 = track.
a()[1];
1537 double kappa = track.
a()[2];
1538 double dz = track.
a()[3];
1539 double tanl = track.
a()[4];
1542 double vx = dr*
cos(phi0);
1543 double vy = dr*
sin(phi0);
1548 if(0==kappa) kappa = 10e-10;
1549 double px = -
sin(phi0)/fabs(kappa);
1550 double py =
cos(phi0)/fabs(kappa);
1551 double pz = tanl/fabs(kappa);
1553 trk->
setX(vx, l_mass);
1554 trk->
setY(vy, l_mass);
1555 trk->
setZ(vz, l_mass);
1556 trk->
setPx(px, l_mass);
1557 trk->
setPy(py, l_mass);
1558 trk->
setPz(pz, l_mass);
1569 if (kappa > 0.0000000001)
1571 else if (kappa < -0.0000000001)
1576 double ptot = sqrt(px*px+py*py+pz*pz);
1577 trk->
setTheta(acos(pz/ptot),l_mass);
1584 double dr = TrasanTRK.
helix[0];
1585 double phi0 = TrasanTRK.
helix[1];
1586 double kappa = TrasanTRK.
helix[2];
1587 double dz = TrasanTRK.
helix[3];
1588 double tanl = TrasanTRK.
helix[4];
1590 double vx = dr*
cos(phi0);
1591 double vy = dr*
sin(phi0);
1594 if(0==kappa) kappa = 10e-10;
1595 double px = -
sin(phi0)/fabs(kappa);
1596 double py =
cos(phi0)/fabs(kappa);
1597 double pz = tanl/fabs(kappa);
1599 trk->
setX(vx, l_mass);
1600 trk->
setY(vy, l_mass);
1601 trk->
setZ(vz, l_mass);
1603 trk->
setPx(px, l_mass);
1604 trk->
setPy(py, l_mass);
1605 trk->
setPz(pz, l_mass);
1612 double a_trasan[5], ea_trasan[15];
1613 for(
int i =0 ; i <5; i++){
1614 a_trasan[i] = TrasanTRK.
helix[i];
1616 for(
int j =0 ; j <15; j++){
1617 ea_trasan[j] = TrasanTRK.
error[j];
1624 if (kappa > 0.0000000001)
1626 else if (kappa < -0.0000000001)
1631 double ptot = sqrt(px*px+py*py+pz*pz);
1632 trk->
setTheta(acos(pz/ptot),l_mass);
1637 std::cout<<
"px: "<<trk->
px()<<
" py: "<<trk->
py()<<
" pz: "<<trk->
pz()<<std::endl;
1638 std::cout<<
"vx: "<<trk->
x()<<
" vy: "<<trk->
y()<<
" vz: "<<trk->
z()<<std::endl;
1655 if(
debug_ == 4) cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0] "<<trk->
getNdf(0,2)<<endl;
1659 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1660 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1661 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
1672 if(
debug_ == 4) cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
1682 for (
int i = 0; i<43; i++) {
1696 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
1697 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1698 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1699 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1703 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
1708 TrasanTRK.
pivot[2]);
1711 for(
int i = 0; i < 5; i++)
1712 a[i] = TrasanTRK.
helix[i];
1715 for(
int i = 0, k = 0; i < 5; i++) {
1716 for(
int j = 0; j <= i; j++) {
1718 ea[j][i] = ea[i][j];
1724 double fiTerm = TrasanTRK.
fiTerm;
1726 double fi0 = track_rep.
phi0();
1732 if( fabs(
x ) > 1.0e-10 ){
1733 phi_x = atan2( y,
x );
1734 if( phi_x < 0 ) phi_x += 2*
M_PI;
1736 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
1738 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
1739 double dphi = fabs( fiTerm + fi0 - phi_x );
1740 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
1741 double tanl = track_rep.
tanl();
1742 double cosl_inv = sqrt( tanl*tanl + 1.0 );
1744 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
1745 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
1747 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
1748 double light_speed( 29.9792458 );
1749 double pt( 1.0 / track_rep.
kappa() );
1750 double p( pt * sqrt( 1.0 + tanl*tanl ) );
1756 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<endl;
1757 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<endl;
1762 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
1763 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
1765 track_rep.
pivot(IP);
1781 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1782 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1783 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1785 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
1786 cout <<
" dr = " << track.
a()[0]
1787 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
1788 cout<<
" phi0 = " << track.
a()[1]
1789 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
1790 cout <<
" PT = " << 1/track.
a()[2]
1791 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
1792 cout <<
" dz = " << track.
a()[3]
1793 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
1794 cout <<
" tanl = " << track.
a()[4]
1795 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
1796 cout <<
" Ea = " << track.
Ea() <<endl;
1823 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
1824 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
1825 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;
1829 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
1830 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
1831 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
1840 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
1851 helixseg->
setResIncl(it->residual_include());
1852 helixseg->
setResExcl(it->residual_exclude());
1854 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
1856 helixseg->
setDrIncl(it->a_include()[0]);
1859 helixseg->
setDzIncl(it->a_include()[3]);
1863 helixseg->
setDrExcl(it->a_exclude()[0]);
1866 helixseg->
setDzExcl(it->a_exclude()[3]);
1881 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
1882 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
1883 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
1884 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
1886 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
1887 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
1888 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
1890 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
1894 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
1895 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
1896 helixseg->
setFlagLR(it->HitMdc()->LR());
1897 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
1898 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
1899 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
1900 helixseg->
setTof(it->tof());
1903 helixseg->
setDD(it->dd());
1904 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
1905 helixseg->
setDT(it->dt());
1906 segcol->push_back(helixseg);
1907 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
1908 helixsegrefvec.push_back(refhelixseg);
1911 m_docaInc = helixseg -> getDocaIncl();
1912 m_docaExc = helixseg -> getDocaExcl();
1913 m_residualInc = helixseg -> getResIncl();
1914 m_residualExc = helixseg -> getResExcl();
1915 m_dd = helixseg -> getDD();
1917 m_tdrift = helixseg -> getDT();
1918 m_layerid = helixseg -> getLayerId();
1919 m_yposition= it->HitMdc()->wire().fwd().y();
1921 StatusCode sc6 = m_nt6->write();
1922 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
1929 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
1937 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
1941 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
1950 for (
int i = 0; i<43; i++) {
1961 double a_trasan[5], ea_trasan[15];
1962 for(
int i =0 ; i <5; i++){
1963 a_trasan[i] = TrasanTRK.
helix[i];
1965 for(
int j =0 ; j <15; j++){
1966 ea_trasan[j] = TrasanTRK.
helix[j];
1972 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
1973 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
1974 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
1975 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
1980 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
1986 TrasanTRK.
pivot[2]);
1989 for(
int i = 0; i < 5; i++)
1990 a[i] = TrasanTRK.
helix[i];
1993 for(
int i = 0, k = 0; i < 5; i++) {
1994 for(
int j = 0; j <= i; j++) {
1996 ea[j][i] = ea[i][j];
2002 double fiTerm = TrasanTRK.
fiTerm;
2004 double fi0 = track_rep.
phi0();
2010 if( fabs(
x ) > 1.0e-10 ){
2011 phi_x = atan2( y,
x );
2012 if( phi_x < 0 ) phi_x += 2*
M_PI;
2014 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
2016 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
2017 double dphi = fabs( fiTerm + fi0 - phi_x );
2018 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
2019 double tanl = track_rep.
tanl();
2020 double cosl_inv = sqrt( tanl*tanl + 1.0 );
2022 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
2023 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
2025 double track_len(fabs( track_rep.
radius() * dphi * cosl_inv ));
2026 double light_speed( 29.9792458 );
2027 double pt( 1.0 / track_rep.
kappa() );
2028 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2036 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
2037 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
2042 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2043 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2045 track_rep.
pivot(IP);
2063 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2064 cout <<
" dr = " << track.
a()[0]
2065 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2066 cout<<
" phi0 = " << track.
a()[1]
2067 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2068 cout <<
" PT = " << 1/track.
a()[2]
2069 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2070 cout <<
" dz = " << track.
a()[3]
2071 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2072 cout <<
" tanl = " << track.
a()[4]
2073 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2074 cout <<
" Ea = " << track.
Ea() <<endl;
2099 std::cout<<
"fillTds_back> mass "<<trk->
getMass(2)<<
" ndf[0][l_mass] "<<trk->
getNdf(0,l_mass)<<endl;
2100 std::cout<<
"ndf_back "<< track.
ndf_back() <<
" chi2_back " << track.
chiSq_back()<<endl;
2105 track.
Ea()[0][0] > 0 && track.
Ea()[1][1] > 0 &&
2106 track.
Ea()[2][2] > 0 && track.
Ea()[3][3] > 0 &&
2107 track.
Ea()[4][4] > 0 && fabs(track.
a()[0]) <
DBL_MAX &&
2116 for (vector<KalFitHelixSeg>::iterator it = track.
HelixSegs().begin(); it!=track.
HelixSegs().end();it ++)
2127 helixseg->
setResIncl(it->residual_include());
2128 helixseg->
setResExcl(it->residual_exclude());
2130 std::cout<<
"helixseg->Res_inc ..."<<helixseg->
getResIncl()<<std::endl;
2157 std::cout<<
"KalFitHelixSeg track id .."<<it->HitMdc()->rechitptr()->getTrkId()<<std::endl;
2158 std::cout<<
"helixseg a: "<<it->a()<<std::endl;
2159 std::cout<<
"helixseg a_excl: "<<helixseg->
getHelixExcl()<<std::endl;
2160 std::cout<<
"helixseg a_incl: "<<helixseg->
getHelixIncl()<<std::endl;
2162 std::cout<<
"helixseg Ea: "<<it->Ea()<<std::endl;
2163 std::cout<<
"helixseg Ea_excl: "<<helixseg->
getErrorExcl()<<std::endl;
2164 std::cout<<
"helixseg Ea_incl: "<<helixseg->
getErrorIncl()<<std::endl;
2166 std::cout<<
"helixseg layer: "<<it->layer()<<std::endl;
2170 helixseg->
setTrackId(it->HitMdc()->rechitptr()->getTrkId());
2171 helixseg->
setMdcId(it->HitMdc()->rechitptr()->getMdcId());
2172 helixseg->
setFlagLR(it->HitMdc()->LR());
2173 helixseg->
setTdc(it->HitMdc()->rechitptr()->getTdc());
2174 helixseg->
setAdc(it->HitMdc()->rechitptr()->getAdc());
2175 helixseg->
setZhit(it->HitMdc()->rechitptr()->getZhit());
2176 helixseg->
setTof(it->tof());
2179 helixseg->
setDD(it->dd());
2180 helixseg->
setEntra(it->HitMdc()->rechitptr()->getEntra());
2181 helixseg->
setDT(it->dt());
2183 segcol->push_back(helixseg);
2184 SmartRef<RecMdcKalHelixSeg> refhelixseg(helixseg);
2185 helixsegrefvec.push_back(refhelixseg);
2187 m_docaInc = helixseg -> getDocaIncl();
2188 m_docaExc = helixseg -> getDocaExcl();
2189 m_residualInc = helixseg -> getResIncl();
2190 m_residualExc = helixseg -> getResExcl();
2191 m_dd = helixseg -> getDD();
2193 m_tdrift = helixseg -> getDT();
2194 m_layerid = helixseg -> getLayerId();
2195 m_yposition= it->HitMdc()->wire().fwd().y();
2197 StatusCode sc6 = m_nt6->write();
2198 if( sc6.isFailure() ) cout<<
"Ntuple6 helixseg filling failed!"<<endl;
2206 std::cout<<
"trk->getVecHelixSegs size..."<<(trk->
getVecHelixSegs()).size()<<std::endl;
2214 std::cout<<
"THEY ARE NOT EQUALL!!!"<<std::endl;
2218 std::cout<<
"l_mass "<<l_mass<<
" path set as "<<track.
pathip()<<endl;
2227 for (
int i = 0; i<43; i++) {
2238 double a_trasan[5], ea_trasan[15];
2239 for(
int i =0 ; i <5; i++){
2240 a_trasan[i] = TrasanTRK.
helix[i];
2242 for(
int j =0 ; j <15; j++){
2243 ea_trasan[j] = TrasanTRK.
helix[j];
2249 std::cout<<
" last pivot: "<< trk->
getLPivot(0)<<std::endl;
2250 std::cout<<
" pathl in SM: "<< trk->
getPathSM(0)<<std::endl;
2251 std::cout<<
" fiTerm: "<< trk->
getFiTerm(0)<<std::endl;
2252 std::cout<<
" last point: "<< trk->
getLPoint(0)<<std::endl;
2257 if(
debug_) cout<<
"ALARM: FillTds_back Not refit with KalFilter!!!"<<endl;
2263 TrasanTRK.
pivot[2]);
2266 for(
int i = 0; i < 5; i++)
2267 a[i] = TrasanTRK.
helix[i];
2270 for(
int i = 0, k = 0; i < 5; i++) {
2271 for(
int j = 0; j <= i; j++) {
2273 ea[j][i] = ea[i][j];
2279 double fiTerm = TrasanTRK.
fiTerm;
2281 double fi0 = track_rep.
phi0();
2287 if( fabs(
x ) > 1.0e-10 ){
2288 phi_x = atan2( y,
x );
2289 if( phi_x < 0 ) phi_x += 2*
M_PI;
2291 phi_x = ( y > 0 ) ? M_PI_4: 3.0*M_PI_4;
2293 if(
debug_ == 4) cout<<
"fiterm "<<fiTerm<<
" fi0 "<<fi0<<
" phi_x "<<phi_x<<endl;
2294 double dphi = fabs( fiTerm + fi0 - phi_x );
2295 if( dphi >= 2*
M_PI ) dphi -= 2*
M_PI;
2296 double tanl = track_rep.
tanl();
2297 double cosl_inv = sqrt( tanl*tanl + 1.0 );
2299 cout<<
"tanl= "<<tanl<<
" radius "<<track_rep.
radius()<<
" dphi "<<dphi<<endl;
2300 cout<<
" cosl_inv "<<cosl_inv<<
" radius_numf "<<track_rep.
radius_numf()<<endl;
2303 double track_len(fabs( track_rep.
radius() * fiTerm * cosl_inv ));
2305 double light_speed( 29.9792458 );
2306 double pt( 1.0 / track_rep.
kappa() );
2307 double p( pt * sqrt( 1.0 + tanl*tanl ) );
2315 std::cout<<
".....fillTds_back...chiSq..."<< TrasanTRK.
chiSq<<std::endl;
2316 std::cout<<
"...track_len..."<<track_len<<
" ndf[1] "<< trk->
getNdf(0,l_mass)<<std::endl;
2321 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2322 trk->
setTof(track_len / ( light_speed * beta ), l_mass) ;
2326 track_rep.
pivot(LPiovt);
2413 cout<<
"Now let us see results after smoothering at IP:........."<<endl;
2414 cout <<
" dr = " << track.
a()[0]
2415 <<
", Er_dr = " << sqrt(track.
Ea()[0][0]) << std::endl;
2416 cout<<
" phi0 = " << track.
a()[1]
2417 <<
", Er_phi0 = " << sqrt(track.
Ea()[1][1]) << std::endl;
2418 cout <<
" PT = " << 1/track.
a()[2]
2419 <<
", Er_kappa = " << sqrt(track.
Ea()[2][2]) << std::endl;
2420 cout <<
" dz = " << track.
a()[3]
2421 <<
", Er_dz = " << sqrt(track.
Ea()[3][3]) << std::endl;
2422 cout <<
" tanl = " << track.
a()[4]
2423 <<
", Er_tanl = " << sqrt(track.
Ea()[4][4]) << std::endl;
2424 cout <<
" Ea = " << track.
Ea() <<endl;
2440 for(
int jj = 0; jj<2; jj++) {
2454 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2455 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2458 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2463 Hep3Vector xorigin(0,0,0);
2465 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
2469 xorigin.setX(dbv[0]);
2470 xorigin.setY(dbv[1]);
2471 xorigin.setZ(dbv[2]);
2485 double tanl = track.
tanl();
2486 double phi_old = work.
phi0();
2487 double phi = track.
phi0();
2489 if (fabs(phi - phi_old) >
M_PI) {
2490 if (phi > phi_old) phi -= 2 *
M_PI;
2491 else phi_old -= 2 *
M_PI;
2494 double path_zero = fabs(track.
radius() * (phi_old-phi)* sqrt(1 + tanl * tanl));
2498 HepSymMatrix Eakal(5,0);
2506 unsigned int nhit = track.
HitsMdc().size();
2507 int layer_prev = -1;
2509 HepVector pos_old(3,0);
2510 double r0kal_prec(0);
2512 for(
unsigned i=0 ; i < nhit; i++ ) {
2513 int ihit = (nhit-1)-i;
2517 int wireid = Wire.
geoID();
2521 if (
pathl_ && layer != layer_prev) {
2523 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
2524 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
2532 Hep3Vector wire = (CLHEP::Hep3Vector)fwd - (CLHEP::Hep3Vector)bck;
2535 work.
pivot((fwd + bck) * .5);
2536 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
2538 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2542 double tension = geowire->
Tension();
2545 double zinit(x0kal.z()), lzx(Wire.
lzx());
2548 double A = 47.35E-6/tension;
2549 double Zp = (zinit - bck.z())*lzx/wire.z();
2552 std::cout<<
" sag in smoother_anal: "<<std::endl;
2553 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2554 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2555 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2556 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2557 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2558 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2561 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2562 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2563 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2565 wire.setX(wire.x()/wire.z());
2566 wire.setY(result.z());
2569 x0kal.setX(result.x());
2570 x0kal.setY(result.y());
2573 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2576 double r0kal = x0kal.perp();
2578 cout<<
"wire direction "<<wire<<endl;
2579 cout<<
"x0kal "<<x0kal<<endl;
2580 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
2597 double pmag( sqrt( 1.0 + track.
a()[4]*track.
a()[4]) / track.
a()[2]);
2598 double mass_over_p( track.
mass()/ pmag );
2599 double beta( 1.0 / sqrt( 1.0 + mass_over_p * mass_over_p ) );
2600 double tofest = pathl / ( 29.9792458 * beta );
2610 if(!(way<0&&fabs(track.
kappa())>1000.0)) {
2620 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2622 HepSymMatrix Ma(5,0);
2625 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2639 point.setX(x0kal.x() + track.
a()[0]*
cos(track.
a()[1]));
2640 point.setY(x0kal.y() + track.
a()[0]*
sin(track.
a()[1]));
2641 point.setZ(x0kal.z() + track.
a()[3]);
2645 double phi_old = track.
a()[1];
2648 double phi_new = temp.
a()[1];
2649 double fi = phi_new - phi_old;
2653 if(fabs(fi) >= CLHEP::twopi) fi = fmod(fi+2*CLHEP::twopi,CLHEP::twopi);
2658 if (
debug_) cout<<
"track----7-----"<<track.
a()<<endl;
2670 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2671 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2674 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2677 HepSymMatrix Eakal(5,0);
2683 std::cout<<
"the initial error matrix in smoothing is "<<Eakal<<std::endl;
2688 unsigned int nseg = track.
HelixSegs().size();
2689 int layer_prev = -1;
2691 HepVector pos_old(3,0);
2692 double r0kal_prec(0);
2694 for(
unsigned i=0 ; i < nseg; i++ ) {
2697 int iseg = (nseg-1)-i;
2701 int wireid = Wire.
geoID();
2705 if (
pathl_ && layer != layer_prev) {
2707 if (
debug_ == 4) cout<<
"in smoother,layerid "<<layer<<
" layer_prev "
2708 <<layer_prev <<
" pathl_ "<<
pathl_<<endl;
2716 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
2719 work.
pivot((fwd + bck) * .5);
2720 HepPoint3D x0kal = (work.
x(0).z() - bck.z()) / wire.z() * wire + bck;
2723 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2729 double tension = geowire->
Tension();
2732 double zinit(x0kal.z()), lzx(Wire.
lzx());
2736 double A = 47.35E-6/tension;
2737 double Zp = (zinit - bck.z())*lzx/wire.z();
2741 std::cout<<
" sag in smoother_calib: "<<std::endl;
2742 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2743 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2744 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2745 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2746 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2747 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2750 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2751 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2752 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2754 wire.setX(wire.x()/wire.z());
2755 wire.setY(result.z());
2758 x0kal.setX(result.x());
2759 x0kal.setY(result.y());
2762 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2766 double r0kal = x0kal.perp();
2768 cout<<
"wire direction "<<wire<<endl;
2769 cout<<
"x0kal "<<x0kal<<endl;
2770 cout<<
"smoother::r0kal "<<r0kal<<
" r0kal_prec "<<r0kal_prec <<endl;
2777 if (
debug_ == 4) cout<<
"track----6-----"<<track.
a()<<
" ...path..."<<pathl
2778 <<
"momentum"<<track.
momentum(0)<<endl;
2783 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2785 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2792 if (
debug_) cout<<
"track----7-----"<<track.
a()<<endl;
2811 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
2812 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
2815 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
2818 Hep3Vector x0inner(track.
pivot());
2819 HepVector pos_old(3,0);
2820 double r0kal_prec(0);
2822 int nhit = track.
HitsMdc().size();
2823 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
2824 for(
int i=0 ; i < nhit; i++ ) {
2827 if (HitMdc.
chi2()<0)
continue;
2833 int wireid = Wire.
geoID();
2837 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
2840 work.
pivot((fwd + bck) * .5);
2847 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
2849 if(4 ==
debug_) std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
2893 double tension = geowire->
Tension();
2895 double zinit(x0kal.z()), lzx(Wire.
lzx());
2897 double A = 47.35E-6/tension;
2898 double Zp = (zinit - bck.z())*lzx/wire.z();
2901 std::cout<<
" sag in filter_fwd_anal: "<<std::endl;
2902 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
2903 std::cout<<
"zinit: "<<zinit<<
" bck.z(): "<<bck.z()<<std::endl;
2904 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
2905 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
2906 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
2907 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
2908 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
2911 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
2912 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
2913 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
2915 wire.setX(wire.x()/wire.z());
2916 wire.setY(result.z());
2918 x0kal.setX(result.x());
2919 x0kal.setY(result.y());
2922 if(4 ==
debug_) std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
2925 double r0kal = x0kal.perp();
2932 if (nhits_read == 1) {
2939 double dtracknew = 0.;
2943 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
2944 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
2945 double diff_chi2 = track.
chiSq();
2946 Hep3Vector IP(0,0,0);
2952 for(
unsigned k=i+1 ; k < nhit; k++ )
2957 double dchi2 = -1.0;
3008 std::cout<<
" ... ERROR OF dchi2... "<<std::endl;
3014 m_residest = dtrack-dtdc;
3015 m_residnew = dtracknew -dtdc;
3018 m_anal_dr = worktemp.
a()[0];
3019 m_anal_phi0 = worktemp.
a()[1];
3020 m_anal_kappa = worktemp.
a()[2];
3021 m_anal_dz = worktemp.
a()[3];
3022 m_anal_tanl = worktemp.
a()[4];
3023 m_anal_ea_dr = worktemp.
Ea()[0][0];
3024 m_anal_ea_phi0 = worktemp.
Ea()[1][1];
3025 m_anal_ea_kappa = worktemp.
Ea()[2][2];
3026 m_anal_ea_dz = worktemp.
Ea()[3][3];
3027 m_anal_ea_tanl = worktemp.
Ea()[4][4];
3028 StatusCode sc5 = m_nt5->write();
3029 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3034 diff_chi2 = chi2 - diff_chi2;
3035 HitMdc.
chi2(diff_chi2);
3049 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3050 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3053 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3057 std::cout<<
"filter_fwd_calib starting ...the inital error Matirx is "<<track.
Ea()<<std::endl;
3059 Hep3Vector x0inner(track.
pivot());
3060 HepVector pos_old(3,0);
3061 double r0kal_prec(0);
3064 unsigned int nhit = track.
HitsMdc().size();
3065 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
3066 for(
unsigned i=0 ; i < nhit; i++ ) {
3070 cout<<
"filter_fwd...........222 chi2="<<HitMdc.
chi2()<<endl;
3072 if (HitMdc.
chi2()<0)
continue;
3075 int wireid = Wire.
geoID();
3081 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3084 work.
pivot((fwd + bck) * .5);
3085 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
3088 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3093 double tension = geowire->
Tension();
3095 double zinit(x0kal.z()), lzx(Wire.
lzx());
3097 double A = 47.35E-6/tension;
3098 double Zp = (zinit - bck.z())*lzx/wire.z();
3102 std::cout<<
" sag in filter_fwd_calib: "<<std::endl;
3103 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3104 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3105 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3106 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3107 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3108 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3111 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3112 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3113 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3114 wire.setX(wire.x()/wire.z());
3115 wire.setY(result.z());
3117 x0kal.setX(result.x());
3118 x0kal.setY(result.y());
3122 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3125 double r0kal = x0kal.perp();
3134 cout <<
"*** move from " << track.
pivot() <<
" ( Kappa = "
3135 << track.
a()[2] <<
")" << endl;
3143 cout <<
"*** to " << track.
pivot() <<
" ( Kappa = "
3144 << track.
a()[2] <<
")" << std::endl;
3146 if (nhits_read == 1) {
3148 if(
debug_==4) cout <<
"filter_fwd::Ea = " << track.
Ea()<<endl;
3161 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3162 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3164 double diff_chi2 = track.
chiSq();
3166 Hep3Vector IP(0,0,0);
3172 for(
unsigned k=i+1 ; k < nhit; k++ )
3177 double dchi2 = -1.0;
3180 std::cout<<
"the value of x0kal is "<<x0kal<<std::endl;
3181 std::cout<<
"the value of track.pivot() is "<<track.
pivot()<<std::endl;
3185 HepSymMatrix Ma(5,0);
3189 std::cout<<
"HelixSeg.Ea_pre_fwd() ..."<<HelixSeg.
Ea_pre_fwd()<<std::endl;
3190 std::cout<<
"HelixSeg.a_pre_fwd() ..."<<HelixSeg.
a_pre_fwd()<<std::endl;
3191 std::cout<<
"HelixSeg.Ea_filt_fwd() ..."<<HelixSeg.
Ea_filt_fwd()<<std::endl;
3204 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;
3206 if(
debug_ == 4) cout<<
"****inext***"<<inext<<
" !!!! dchi2= "<< dchi2
3207 <<
" chisq= "<< chi2<< endl;
3212 StatusCode sc5 = m_nt5->write();
3213 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3219 diff_chi2 = chi2 - diff_chi2;
3220 HitMdc.
chi2(diff_chi2);
3224 cout <<
" -> after include meas = " << meas
3225 <<
" at R = " << track.
pivot().perp() << std::endl;
3226 cout <<
" chi2 = " << chi2 <<
", diff_chi2 = "
3227 << diff_chi2 <<
", LR = "
3228 << HitMdc.
LR() <<
", stereo = " << HitMdc.
wire().
stereo()
3230 cout<<
"filter_fwd..........HitMdc.chi2... "<<HitMdc.
chi2()<<endl;
3233 if(dchi2>0.0 && (way!=999)) {
3244 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3245 for(
int i = 0; i < _BesKalmanFitWalls.size(); i++) {
3246 _BesKalmanFitWalls[i].updateTrack(track, way);
3247 if (
debug_ == 4) cout<<
"Wall "<<i<<
" update the track!(filter)"<<endl;
3252 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3253 for(
int i = begin; i <= end; i++) {
3254 _BesKalmanFitWalls[i].updateTrack(track, way);
3255 if(
debug_ == 4){cout <<
"wall" << i << endl;
3256 cout <<
"track" <<
"pivot=" << track.
pivot() <<
"a=" << track.
a() <<
"ea=" <<track.
Ea() << endl;}
3257 if (
debug_ == 4) cout<<
"Wall "<<i<<
" update the track!(filter)"<<endl;
3262 if (
debug_ ==4) cout<<
"innerwall....."<<endl;
3263 double ri=(r1<r2)?r1:r2;
3264 double ro=(r1<r2)?r2:r1;
3267 vector<KalFitCylinder>::iterator it;
3268 vector<KalFitCylinder>::iterator itStart;
3269 vector<KalFitCylinder>::iterator itStop;
3272 itStart=_BesKalmanFitWalls.end()-index-1;
3273 itStop =_BesKalmanFitWalls.begin()-1;
3277 itStart=_BesKalmanFitWalls.begin()+index;
3278 itStop =_BesKalmanFitWalls.end();
3280 for(it=itStart; it != itStop; it=it+step) {
3281 if(ri < it->rmin() && it->radius() < ro) {
3282 it->updateTrack(track,way);
3287 else if(it->rmin() <= ri && ri <= it->radius() && it->radius() < ro){
3288 it->updateTrack(track,way,ri,it->
radius());
3297 else if(ri < it->rmin() && it->rmin() <= ro && ro <= it->radius()){
3298 it->updateTrack(track,way,it->rmin(),ro);
3307 else if(it->radius()<ri || ro < it->rmin()){
3311 it->updateTrack(track,way,ri,ro);
3323 index = _BesKalmanFitWalls.end() - it-1;
3330 index = it - _BesKalmanFitWalls.begin();
3345 MsgStream log(
msgSvc(), name());
3346 double Pt_threshold(0.3);
3347 Hep3Vector IP(0,0,0);
3353 if ( !&whMgr )
return;
3356 int ntrk = mdcMgr->size();
3357 double* rPt =
new double[ntrk];
3358 int* rOM =
new int[ntrk];
3359 unsigned int* order =
new unsigned int[ntrk];
3360 unsigned int* rCont =
new unsigned int[ntrk];
3361 unsigned int* rGen =
new unsigned int[ntrk];
3363 if(
testOutput) cout<<
"nMdcTrk = "<<ntrk<<endl;
3366 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3367 end = mdcMgr->end(); it != end; it++) {
3371 rPt[index] = 1 / fabs(it->helix[2]);
3372 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3373 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3375 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3376 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3378 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3379 ii !=pt.begin()-1; ii--) {
3380 int lyr((*ii)->geo->Lyr()->Id());
3381 if (outermost < lyr) outermost = lyr;
3382 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3384 rOM[index] = outermost;
3385 order[index] = index;
3390 for (
int j, k = ntrk - 1; k >= 0; k = j){
3392 for(
int i = 1; i <= k; i++)
3393 if(rPt[i - 1] < rPt[i]){
3395 std::swap(order[i], order[j]);
3396 std::swap(rPt[i], rPt[j]);
3397 std::swap(rOM[i], rOM[j]);
3398 std::swap(rCont[i], rCont[j]);
3399 std::swap(rGen[i], rGen[j]);
3407 DataObject *aReconEvent;
3408 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3412 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3413 if(sc!=StatusCode::SUCCESS) {
3414 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3422 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3424 for(
int l = 0; l < ntrk; l++) {
3433 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3437 int trasqual = TrasanTRK_add.
quality;
3438 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3439 if (trasqual)
continue;
3443 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3447 if ((TrasanTRK_add.
decision & 32) == 32 ||
3448 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3453 TrasanTRK.
pivot[2]);
3456 for(
int i = 0; i < 5; i++)
3457 a[i] = TrasanTRK.
helix[i];
3460 for(
int i = 0, k = 0; i < 5; i++) {
3461 for(
int j = 0; j <= i; j++) {
3463 ea[j][i] = ea[i][j];
3466 double fiTerm = TrasanTRK.
fiTerm;
3474 int inlyr(999), outlyr(-1);
3475 int* rStat =
new int[43];
3476 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3477 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3479 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3482 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3483 ii != pt.begin()-1; ii--) {
3484 Num[(*ii)->geo->Lyr()->Id()]++;
3487 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3488 ii != pt.begin()-1; ii--) {
3490 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
3492 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3493 <<
" hits in the layer "
3494 << (*ii)->geo->Lyr()->Id() << std::endl;
3512 double dist[2] = {rechit.
ddl, rechit.
ddr};
3513 double erdist[2] = {rechit.
erddl, rechit.
erddr};
3518 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
3520 else if (rechit.
lr==1) lr_decision=1;
3523 int ind(geo->
Lyr()->
Id());
3527 lr_decision, rechit.
tdc,
3532 if (inlyr>ind) inlyr = ind;
3533 if (outlyr<ind) outlyr = ind;
3536 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3539 int empty_between(0), empty(0);
3540 for (
int i= inlyr; i <= outlyr; i++)
3541 if (!rStat[i]) empty_between++;
3542 empty = empty_between+inlyr+(42-outlyr);
3549 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
3555 track_lead.
type(type);
3557 unsigned int nhit = track_lead.
HitsMdc().size();
3558 if (!nhit &&
debug_ == 4) {
3559 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3564 double KalFitst(0), KalFitax(0), KalFitschi2(0);
3567 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
3569 track_lead.
pivot(outer_pivot);
3575 HepSymMatrix Eakal(5,0);
3578 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
3588 cout <<
"from Mdc Pattern Recognition: " << std::endl;
3594 cout <<
" dr = " << work.
a()[0]
3595 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3596 cout <<
" phi0 = " << work.
a()[1]
3597 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3598 cout <<
" PT = " << 1/work.
a()[2]
3599 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3600 cout <<
" dz = " << work.
a()[3]
3601 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3602 cout <<
" tanl = " << work.
a()[4]
3603 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3607 cout<<
"----------------------"<<endl;
3608 cout<<
"track "<<l<<
", pid = "<<
lead_<<
": "<<endl;
3625 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3630 cout <<
" dr = " << work.
a()[0]
3631 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3632 cout <<
" phi0 = " << work.
a()[1]
3633 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3634 cout <<
" PT = " << 1/work.
a()[2]
3635 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3636 cout <<
" dz = " << work.
a()[3]
3637 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3638 cout <<
" tanl = " << work.
a()[4]
3639 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3645 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
3649 IDataProviderSvc* eventSvc = NULL;
3650 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
3652 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
3654 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
3658 IDataManagerSvc *dataManSvc;
3659 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
3660 DataObject *aKalTrackCol;
3661 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
3662 if(aKalTrackCol != NULL) {
3663 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
3664 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
3666 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
3667 if( kalsc.isFailure() ) {
3668 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
3671 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
3674 DataObject *aRecKalSegEvent;
3675 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
3676 if(aRecKalSegEvent!=NULL) {
3678 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
3679 if(segsc != StatusCode::SUCCESS) {
3680 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
3684 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
3685 if( segsc.isFailure() ) {
3686 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
3689 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
3691 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.);
3692 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
3693 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
3696 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
3698 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
3701 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
3702 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
3703 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
3704 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
3705 <<
"Track Id: " << (*iter_trk)->getTrackId()
3706 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
3707 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
3708 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
3709 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
3710 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
3711 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
3712 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
3713 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
3715 for(
int i = 0; i<43; i++) {
3716 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
3717 << (*iter_trk)->getPathl(i) <<endreq;
3720 m_trackid = (*iter_trk)->getTrackId();
3722 for(
int jj =0, iii=0; jj<5; jj++) {
3724 m_length[jj] = (*iter_trk)->getLength(jj);
3725 m_tof[jj] = (*iter_trk)->getTof(jj);
3726 m_nhits[jj] = (*iter_trk)->getNhits(jj);
3727 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
3728 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
3729 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
3730 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
3731 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
3732 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
3733 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
3734 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
3735 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
3736 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
3737 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
3738 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
3739 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
3740 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
3741 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
3744 for(
int kk=0; kk<=jj; kk++,iii++) {
3745 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
3746 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
3747 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
3748 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
3749 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
3750 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
3751 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
3752 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
3753 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
3754 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
3755 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
3756 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
3757 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
3758 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
3759 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
3799 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
3800 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
3801 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
3802 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
3803 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
3804 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
3805 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
3806 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
3807 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
3808 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
3810 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
3811 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
3812 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
3813 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
3814 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
3815 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
3816 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
3817 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
3818 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
3819 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
3821 m_stat[0][0] = (*iter_trk)->getStat(0,0);
3822 m_stat[1][0] = (*iter_trk)->getStat(0,1);
3823 m_stat[2][0] = (*iter_trk)->getStat(0,2);
3824 m_stat[3][0] = (*iter_trk)->getStat(0,3);
3825 m_stat[4][0] = (*iter_trk)->getStat(0,4);
3826 m_stat[0][1] = (*iter_trk)->getStat(1,0);
3827 m_stat[1][1] = (*iter_trk)->getStat(1,1);
3828 m_stat[2][1] = (*iter_trk)->getStat(1,2);
3829 m_stat[3][1] = (*iter_trk)->getStat(1,3);
3830 m_stat[4][1] = (*iter_trk)->getStat(1,4);
3832 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
3833 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
3834 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
3835 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
3836 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
3838 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
3839 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
3840 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
3841 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
3842 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
3844 m_lpt = 1/m_lhelix[2];
3845 m_lpte = 1/m_lhelixe[2];
3846 m_lptmu = 1/m_lhelixmu[2];
3847 m_lptk = 1/m_lhelixk[2];
3848 m_lptp = 1/m_lhelixp[2];
3850 m_fpt = 1/m_fhelix[2];
3851 m_fpte = 1/m_fhelixe[2];
3852 m_fptmu = 1/m_fhelixmu[2];
3853 m_fptk = 1/m_fhelixk[2];
3854 m_fptp = 1/m_fhelixp[2];
3857 std::cout<<
" "<<std::endl;
3858 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
3859 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
3860 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
3861 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
3862 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
3865 m_zpt = 1/m_zhelix[2];
3866 m_zpte = 1/m_zhelixe[2];
3867 m_zptmu = 1/m_zhelixmu[2];
3868 m_zptk = 1/m_zhelixk[2];
3869 m_zptp = 1/m_zhelixp[2];
3872 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
3873 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
3874 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
3875 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
3876 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
3878 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
3879 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
3880 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
3881 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
3882 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
3885 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
3886 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
3887 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
3888 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
3889 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
3893 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
3894 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
3895 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
3896 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
3897 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
3898 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
3899 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
3900 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
3901 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
3902 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
3903 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
3904 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
3905 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
3906 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
3907 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
3910 StatusCode sc1 = m_nt1->write();
3911 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
3916 phi1 = (*iter_trk)->getFFi0();
3917 r1 = (*iter_trk)->getFDr();
3918 z1 = (*iter_trk)->getFDz();
3919 kap1 = (*iter_trk)->getFCpa();
3920 tanl1 = (*iter_trk)->getFTanl();
3921 charge1 = kap1/fabs(kap1);
3924 p1 = sqrt(1+tanl1*tanl1)/kap1;
3925 the1 =
M_PI/2-atan(tanl1);
3928 pz1= tanl1/fabs(kap1);
3930 }
else if(jj == 2) {
3931 phi2 = (*iter_trk)->getFFi0();
3932 r2 = (*iter_trk)->getFDr();
3933 z2 = (*iter_trk)->getFDz();
3934 kap2 = (*iter_trk)->getFCpa();
3935 tanl2 = (*iter_trk)->getFTanl();
3936 charge2 = kap2/fabs(kap2);
3939 p2 = sqrt(1+tanl2*tanl2)/kap1;
3940 the2 =
M_PI/2-atan(tanl2);
3943 pz2= tanl2/fabs(kap2);
3951 m_delthe = the1 + the2;
3954 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
3955 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
3956 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
3958 StatusCode sc2 = m_nt2->write();
3959 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
3968 cout <<
"Kalfitting finished " << std::endl;
3975 MsgStream log(
msgSvc(), name());
3976 double Pt_threshold(0.3);
3977 Hep3Vector IP(0,0,0);
3984 if ( !&whMgr )
return;
3987 int ntrk = mdcMgr->size();
3988 double* rPt =
new double[ntrk];
3989 int* rOM =
new int[ntrk];
3990 unsigned int* order =
new unsigned int[ntrk];
3991 unsigned int* rCont =
new unsigned int[ntrk];
3992 unsigned int* rGen =
new unsigned int[ntrk];
3995 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3996 end = mdcMgr->end(); it != end; it++) {
4000 rPt[index] = 1 / fabs(it->helix[2]);
4001 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4002 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4004 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4005 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4007 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4008 ii !=pt.begin()-1; ii--) {
4009 int lyr((*ii)->geo->Lyr()->Id());
4010 if (outermost < lyr) outermost = lyr;
4011 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4013 rOM[index] = outermost;
4014 order[index] = index;
4019 for (
int j, k = ntrk - 1; k >= 0; k = j){
4021 for(
int i = 1; i <= k; i++)
4022 if(rPt[i - 1] < rPt[i]){
4024 std::swap(order[i], order[j]);
4025 std::swap(rPt[i], rPt[j]);
4026 std::swap(rOM[i], rOM[j]);
4027 std::swap(rCont[i], rCont[j]);
4028 std::swap(rGen[i], rGen[j]);
4035 DataObject *aReconEvent;
4036 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4040 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4041 if(sc!=StatusCode::SUCCESS) {
4042 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4050 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4053 for(
int l = 0; l < ntrk; l++) {
4055 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
4059 int trasqual = TrasanTRK_add.
quality;
4060 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4061 if (trasqual)
continue;
4065 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4069 if ((TrasanTRK_add.
decision & 32) == 32 ||
4070 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4075 TrasanTRK.
pivot[2]);
4078 for(
int i = 0; i < 5; i++)
4079 a[i] = TrasanTRK.
helix[i];
4082 for(
int i = 0, k = 0; i < 5; i++) {
4083 for(
int j = 0; j <= i; j++) {
4085 ea[j][i] = ea[i][j];
4091 double fiTerm = TrasanTRK.
fiTerm;
4097 int inlyr(999), outlyr(-1);
4098 int* rStat =
new int[43];
4099 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4100 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4102 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4105 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4106 ii != pt.begin()-1; ii--) {
4107 Num[(*ii)->geo->Lyr()->Id()]++;
4111 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4112 ii != pt.begin()-1; ii--) {
4115 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4117 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4118 <<
" hits in the layer "
4119 << (*ii)->geo->Lyr()->Id() << std::endl;
4146 double dist[2] = {rechit.
ddl, rechit.
ddr};
4147 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4152 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4154 else if (rechit.
lr==1) lr_decision=1;
4157 int ind(geo->
Lyr()->
Id());
4159 lr_decision, rechit.
tdc,
4164 if (inlyr>ind) inlyr = ind;
4165 if (outlyr<ind) outlyr = ind;
4168 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4171 int empty_between(0), empty(0);
4172 for (
int i= inlyr; i <= outlyr; i++)
4173 if (!rStat[i]) empty_between++;
4174 empty = empty_between+inlyr+(42-outlyr);
4179 track_lead.
type(type);
4180 unsigned int nhit = track_lead.
HitsMdc().size();
4181 if (!nhit &&
debug_ == 4) {
4182 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4187 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4189 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4192 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4194 track_lead.
pivot(outer_pivot);
4199 HepSymMatrix Eakal(5,0);
4203 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4213 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4219 std::cout <<
" dr = " << work.
a()[0]
4220 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4221 std::cout <<
" phi0 = " << work.
a()[1]
4222 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4223 std::cout <<
" PT = " << 1/work.
a()[2]
4224 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4225 std::cout <<
" dz = " << work.
a()[3]
4226 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4227 std::cout <<
" tanl = " << work.
a()[4]
4228 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4236 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4241 cout <<
" dr = " << work.
a()[0]
4242 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4243 cout <<
" phi0 = " << work.
a()[1]
4244 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4245 cout <<
" PT = " << 1/work.
a()[2]
4246 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4247 cout <<
" dz = " << work.
a()[3]
4248 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4249 cout <<
" tanl = " << work.
a()[4]
4250 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4257 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4263 DataObject *aRecKalEvent;
4264 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4265 if(aRecKalEvent!=NULL) {
4267 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4268 if(kalsc != StatusCode::SUCCESS) {
4269 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4274 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4275 if( kalsc.isFailure()) {
4276 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4279 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4285 DataObject *aRecKalSegEvent;
4286 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4287 if(aRecKalSegEvent!=NULL) {
4289 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4290 if(segsc != StatusCode::SUCCESS) {
4291 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4296 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4297 if( segsc.isFailure() ) {
4298 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4301 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4304 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.);
4305 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4308 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4310 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4313 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4314 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4315 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4316 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4317 <<
"Track Id: " << (*iter_trk)->getTrackId()
4318 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4319 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4320 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4321 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4322 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4323 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4324 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4325 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4330 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4333 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4334 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4336 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4337 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4338 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4339 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4340 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4341 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4344 for(
int i = 0; i<43; i++) {
4345 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4346 << (*iter_trk)->getPathl(i) <<endreq;
4350 m_trackid = (*iter_trk)->getTrackId();
4351 for(
int jj =0, iii=0; jj<5; jj++) {
4352 m_length[jj] = (*iter_trk)->getLength(jj);
4353 m_tof[jj] = (*iter_trk)->getTof(jj);
4354 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4355 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4356 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4357 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4358 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4359 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4360 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4361 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4362 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4363 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4364 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4365 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4366 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4367 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4368 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4369 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4371 for(
int kk=0; kk<=jj; kk++,iii++) {
4372 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4373 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4374 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4375 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4376 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4377 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4378 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4379 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4380 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4381 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4382 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4383 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4384 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4385 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4386 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4427 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4428 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4429 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4430 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4431 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4432 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4433 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4434 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4435 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4436 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4438 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4439 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4440 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4441 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4442 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4443 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4444 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4445 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4446 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4447 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4449 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4450 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4451 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4452 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4453 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4454 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4455 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4456 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4457 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4458 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4460 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4461 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4462 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4463 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4464 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4466 m_zpt = 1/m_zhelix[2];
4467 m_zpte = 1/m_zhelixe[2];
4468 m_zptmu = 1/m_zhelixmu[2];
4469 m_zptk = 1/m_zhelixk[2];
4470 m_zptp = 1/m_zhelixp[2];
4472 m_fpt = 1/m_fhelix[2];
4473 m_fpte = 1/m_fhelixe[2];
4474 m_fptmu = 1/m_fhelixmu[2];
4475 m_fptk = 1/m_fhelixk[2];
4476 m_fptp = 1/m_fhelixp[2];
4478 m_lpt = 1/m_lhelix[2];
4479 m_lpte = 1/m_lhelixe[2];
4480 m_lptmu = 1/m_lhelixmu[2];
4481 m_lptk = 1/m_lhelixk[2];
4482 m_lptp = 1/m_lhelixp[2];
4484 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4485 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4486 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4487 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4488 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4490 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4491 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4492 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4493 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4494 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4496 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4497 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4498 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4499 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4500 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4501 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4502 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4503 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4504 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4505 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4506 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4507 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4508 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4509 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4510 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4513 StatusCode sc1 = m_nt1->write();
4514 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4519 phi1 = (*iter_trk)->getFFi0();
4520 r1 = (*iter_trk)->getFDr();
4521 z1 = (*iter_trk)->getFDz();
4522 kap1 = (*iter_trk)->getFCpa();
4523 tanl1 = (*iter_trk)->getFTanl();
4526 p1 = sqrt(1+tanl1*tanl1)/kap1;
4527 the1 =
M_PI/2-atan(tanl1);
4528 }
else if(jj == 2) {
4529 phi2 = (*iter_trk)->getFFi0();
4530 r2 = (*iter_trk)->getFDr();
4531 z2 = (*iter_trk)->getFDz();
4532 kap2 = (*iter_trk)->getFCpa();
4533 tanl2 = (*iter_trk)->getFTanl();
4536 p2 = sqrt(1+tanl2*tanl2)/kap1;
4537 the2 =
M_PI/2-atan(tanl2);
4545 m_delthe = the1 + the2;
4548 StatusCode sc2 = m_nt2->write();
4549 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4557 cout <<
"Kalfitting finished " << std::endl;
4564 MsgStream log(
msgSvc(), name());
4565 double Pt_threshold(0.3);
4566 Hep3Vector IP(0,0,0);
4573 if ( !&whMgr )
return;
4576 int ntrk = mdcMgr->size();
4579 int nhits = whMgr->size();
4584 DataObject *aReconEvent;
4585 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4589 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4590 if(sc!=StatusCode::SUCCESS) {
4591 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4599 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4611 if ((TrasanTRK_add.
decision & 32) == 32 ||
4612 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4617 TrasanTRK.
pivot[2]);
4620 for(
int i = 0; i < 5; i++)
4621 a[i] = TrasanTRK.
helix[i];
4624 for(
int i = 0, k = 0; i < 5; i++) {
4625 for(
int j = 0; j <= i; j++) {
4627 ea[j][i] = ea[i][j];
4633 double fiTerm = TrasanTRK.
fiTerm;
4641 int trasqual = TrasanTRK_add.
quality;
4642 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4643 if (trasqual)
return;
4645 int inlyr(999), outlyr(-1);
4646 int* rStat =
new int[43];
4647 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4648 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4650 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4653 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4654 ii != pt.begin()-1; ii--) {
4655 Num[(*ii)->geo->Lyr()->Id()]++;
4658 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4659 ii != pt.begin()-1; ii--) {
4662 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4664 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4665 <<
" hits in the layer "
4666 << (*ii)->geo->Lyr()->Id() << std::endl;
4672 double dist[2] = {rechit.
ddl, rechit.
ddr};
4673 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4678 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4680 else if (rechit.
lr==1) lr_decision=1;
4683 int ind(geo->
Lyr()->
Id());
4685 lr_decision, rechit.
tdc,
4690 if (inlyr>ind) inlyr = ind;
4691 if (outlyr<ind) outlyr = ind;
4694 int empty_between(0), empty(0);
4695 for (
int i= inlyr; i <= outlyr; i++)
4696 if (!rStat[i]) empty_between++;
4697 empty = empty_between+inlyr+(42-outlyr);
4701 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4706 track_lead.
type(type);
4707 unsigned int nhit = track_lead.
HitsMdc().size();
4709 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4714 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4716 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4719 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4721 track_lead.
pivot(outer_pivot);
4726 HepSymMatrix Eakal(5,0);
4730 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4740 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4746 std::cout <<
" dr = " << work.
a()[0]
4747 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4748 std::cout <<
" phi0 = " << work.
a()[1]
4749 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4750 std::cout <<
" PT = " << 1/work.
a()[2]
4751 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4752 std::cout <<
" dz = " << work.
a()[3]
4753 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4754 std::cout <<
" tanl = " << work.
a()[4]
4755 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4762 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4767 cout <<
" dr = " << work1.
a()[0]
4768 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
4769 cout <<
" phi0 = " << work1.
a()[1]
4770 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
4771 cout <<
" PT = " << 1/work1.
a()[2]
4772 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
4773 cout <<
" dz = " << work1.
a()[3]
4774 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
4775 cout <<
" tanl = " << work1.
a()[4]
4776 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
4783 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4788 DataObject *aRecKalEvent;
4789 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4790 if(aRecKalEvent!=NULL) {
4792 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4793 if(kalsc != StatusCode::SUCCESS) {
4794 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4799 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4800 if( kalsc.isFailure()) {
4801 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4804 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4810 DataObject *aRecKalSegEvent;
4811 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4812 if(aRecKalSegEvent!=NULL) {
4814 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4815 if(segsc != StatusCode::SUCCESS) {
4816 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4821 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4822 if( segsc.isFailure() ) {
4823 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4826 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4829 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.);
4830 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4833 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4835 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4838 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4839 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4840 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4841 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4842 <<
"Track Id: " << (*iter_trk)->getTrackId()
4843 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4844 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4845 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4846 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4847 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4848 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4849 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4850 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4851 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
4856 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4859 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4860 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4862 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4863 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4864 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4865 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4866 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4867 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4870 for(
int i = 0; i<43; i++) {
4871 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4872 << (*iter_trk)->getPathl(i) <<endreq;
4876 m_trackid = (*iter_trk)->getTrackId();
4877 for(
int jj =0, iii=0; jj<5; jj++) {
4878 m_length[jj] = (*iter_trk)->getLength(jj);
4879 m_tof[jj] = (*iter_trk)->getTof(jj);
4880 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4881 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4882 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4883 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4884 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4885 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4886 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4887 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4888 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4889 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4890 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4891 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4892 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4893 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4894 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4895 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4897 for(
int kk=0; kk<=jj; kk++,iii++) {
4898 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4899 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4900 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4901 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4902 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4903 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4904 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4905 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4906 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4907 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4908 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4909 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4910 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4911 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4912 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4919 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4920 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4921 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4922 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4923 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4924 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4925 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4926 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4927 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4928 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4930 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4931 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4932 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4933 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4934 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4935 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4936 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4937 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4938 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4939 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4941 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4942 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4943 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4944 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4945 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4946 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4947 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4948 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4949 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4950 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4952 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4953 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4954 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4955 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4956 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4958 m_zpt = 1/m_zhelix[2];
4959 m_zpte = 1/m_zhelixe[2];
4960 m_zptmu = 1/m_zhelixmu[2];
4961 m_zptk = 1/m_zhelixk[2];
4962 m_zptp = 1/m_zhelixp[2];
4964 m_fpt = 1/m_fhelix[2];
4965 m_fpte = 1/m_fhelixe[2];
4966 m_fptmu = 1/m_fhelixmu[2];
4967 m_fptk = 1/m_fhelixk[2];
4968 m_fptp = 1/m_fhelixp[2];
4970 m_lpt = 1/m_lhelix[2];
4971 m_lpte = 1/m_lhelixe[2];
4972 m_lptmu = 1/m_lhelixmu[2];
4973 m_lptk = 1/m_lhelixk[2];
4974 m_lptp = 1/m_lhelixp[2];
4976 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4977 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4978 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4979 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4980 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4982 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4983 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4984 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4985 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4986 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4988 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4989 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4990 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4991 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4992 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4993 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4994 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4995 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4996 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4997 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4998 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4999 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5000 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5001 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5002 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5005 StatusCode sc1 = m_nt1->write();
5006 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5011 phi1 = (*iter_trk)->getFFi0();
5012 r1 = (*iter_trk)->getFDr();
5013 z1 = (*iter_trk)->getFDz();
5014 kap1 = (*iter_trk)->getFCpa();
5015 tanl1 = (*iter_trk)->getFTanl();
5018 p1 = sqrt(1+tanl1*tanl1)/kap1;
5019 the1 =
M_PI/2-atan(tanl1);
5020 }
else if(jj == 2) {
5021 phi2 = (*iter_trk)->getFFi0();
5022 r2 = (*iter_trk)->getFDr();
5023 z2 = (*iter_trk)->getFDz();
5024 kap2 = (*iter_trk)->getFCpa();
5025 tanl2 = (*iter_trk)->getFTanl();
5028 p2 = sqrt(1+tanl2*tanl2)/kap1;
5029 the2 =
M_PI/2-atan(tanl2);
5037 m_delthe = the1 + the2;
5040 StatusCode sc2 = m_nt2->write();
5041 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5045 cout <<
"Kalfitting finished " << std::endl;
5053 MsgStream log(
msgSvc(), name());
5054 double Pt_threshold(0.3);
5055 Hep3Vector IP(0,0,0);
5062 if ( !&whMgr )
return;
5065 int ntrk = mdcMgr->size();
5068 int nhits = whMgr->size();
5072 double* rY =
new double[ntrk];
5073 double* rfiTerm =
new double[ntrk];
5074 double* rPt =
new double[ntrk];
5075 int* rOM =
new int[ntrk];
5076 unsigned int* order =
new unsigned int[ntrk];
5077 unsigned int* rCont =
new unsigned int[ntrk];
5078 unsigned int* rGen =
new unsigned int[ntrk];
5081 Hep3Vector csmp3[2];
5082 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
5083 end = mdcMgr->end(); it != end; it++) {
5085 rfiTerm[index]=it->fiTerm;
5090 rPt[index] = 1 / fabs(it->helix[2]);
5091 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
5092 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
5094 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
5095 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
5097 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5098 ii !=pt.begin()-1; ii--) {
5099 int lyr((*ii)->geo->Lyr()->Id());
5100 if (outermost < lyr) {
5102 rY[index] = (*ii)->geo->Forward().y();
5104 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
5106 rOM[index] = outermost;
5107 order[index] = index;
5112 for (
int j, k = ntrk - 1; k >= 0; k = j){
5114 for(
int i = 1; i <= k; i++)
5115 if(rY[i - 1] < rY[i]){
5117 std::swap(order[i], order[j]);
5118 std::swap(rY[i], rY[j]);
5119 std::swap(rOM[i], rOM[j]);
5120 std::swap(rCont[i], rCont[j]);
5121 std::swap(rGen[i], rGen[j]);
5130 DataObject *aReconEvent;
5131 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5135 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5136 if(sc!=StatusCode::SUCCESS) {
5137 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5145 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5152 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
5161 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
5165 if ((TrasanTRK_add.
decision & 32) == 32 ||
5166 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5171 TrasanTRK.
pivot[2]);
5174 for(
int i = 0; i < 5; i++)
5175 a[i] = TrasanTRK.
helix[i];
5178 for(
int i = 0, k = 0; i < 5; i++) {
5179 for(
int j = 0; j <= i; j++) {
5181 ea[j][i] = ea[i][j];
5187 double fiTerm = TrasanTRK.
fiTerm;
5194 for(
int l = 0; l < ntrk; l++) {
5195 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
5196 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
5198 int trasqual = TrasanTRK_add1.
quality;
5199 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5200 if (trasqual)
continue;
5202 int inlyr(999), outlyr(-1);
5203 int* rStat =
new int[43];
5204 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5205 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
5207 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5210 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5211 ii != pt.begin()-1; ii--) {
5212 Num[(*ii)->geo->Lyr()->Id()]++;
5215 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5216 ii != pt.begin()-1; ii--) {
5219 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5221 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5222 <<
" hits in the layer "
5223 << (*ii)->geo->Lyr()->Id() << std::endl;
5229 double dist[2] = {rechit.
ddl, rechit.
ddr};
5230 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5235 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5237 else if (rechit.
lr==1) lr_decision=1;
5240 int ind(geo->
Lyr()->
Id());
5242 lr_decision, rechit.
tdc,
5247 if (inlyr>ind) inlyr = ind;
5248 if (outlyr<ind) outlyr = ind;
5251 int empty_between(0), empty(0);
5252 for (
int i= inlyr; i <= outlyr; i++)
5253 if (!rStat[i]) empty_between++;
5254 empty = empty_between+inlyr+(42-outlyr);
5258 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5263 track_lead.
type(type);
5264 unsigned int nhit = track_lead.
HitsMdc().size();
5266 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5271 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5273 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5276 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5278 track_lead.
pivot(outer_pivot);
5283 HepSymMatrix Eakal(5,0);
5287 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5297 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5303 std::cout <<
" dr = " << work.
a()[0]
5304 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5305 std::cout <<
" phi0 = " << work.
a()[1]
5306 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5307 std::cout <<
" PT = " << 1/work.
a()[2]
5308 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5309 std::cout <<
" dz = " << work.
a()[3]
5310 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5311 std::cout <<
" tanl = " << work.
a()[4]
5312 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5319 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5324 cout <<
" dr = " << work1.
a()[0]
5325 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5326 cout <<
" phi0 = " << work1.
a()[1]
5327 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5328 cout <<
" PT = " << 1/work1.
a()[2]
5329 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5330 cout <<
" dz = " << work1.
a()[3]
5331 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5332 cout <<
" tanl = " << work1.
a()[4]
5333 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5340 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5346 DataObject *aRecKalEvent;
5347 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5348 if(aRecKalEvent!=NULL) {
5350 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5351 if(kalsc != StatusCode::SUCCESS) {
5352 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5357 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5358 if( kalsc.isFailure()) {
5359 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5362 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5368 DataObject *aRecKalSegEvent;
5369 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5370 if(aRecKalSegEvent!=NULL) {
5372 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5373 if(segsc != StatusCode::SUCCESS) {
5374 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5379 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5380 if( segsc.isFailure() ) {
5381 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5384 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5387 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.);
5388 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5391 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5393 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5396 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5397 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5398 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5399 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5400 <<
"Track Id: " << (*iter_trk)->getTrackId()
5401 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5402 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5403 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5404 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5405 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5406 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5407 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5408 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5409 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5414 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5417 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5418 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5420 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5421 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5422 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5423 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5424 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5425 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5428 for(
int i = 0; i<43; i++) {
5429 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5430 << (*iter_trk)->getPathl(i) <<endreq;
5434 m_trackid = (*iter_trk)->getTrackId();
5435 for(
int jj =0, iii=0; jj<5; jj++) {
5436 m_length[jj] = (*iter_trk)->getLength(jj);
5437 m_tof[jj] = (*iter_trk)->getTof(jj);
5438 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5439 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5440 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5441 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5442 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5443 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5444 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5445 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5446 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5447 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5448 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5449 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5450 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5451 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5452 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5453 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5455 for(
int kk=0; kk<=jj; kk++,iii++) {
5456 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5457 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5458 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5459 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5460 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5461 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5462 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5463 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5464 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5465 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5466 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5467 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5468 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5469 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5470 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5477 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5478 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5479 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5480 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5481 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5482 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5483 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5484 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5485 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5486 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5488 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5489 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5490 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5491 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5492 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5493 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5494 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5495 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5496 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5497 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5499 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5500 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5501 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5502 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5503 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5504 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5505 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5506 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5507 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5508 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5510 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5511 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5512 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5513 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5514 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5516 m_zpt = 1/m_zhelix[2];
5517 m_zpte = 1/m_zhelixe[2];
5518 m_zptmu = 1/m_zhelixmu[2];
5519 m_zptk = 1/m_zhelixk[2];
5520 m_zptp = 1/m_zhelixp[2];
5522 m_fpt = 1/m_fhelix[2];
5523 m_fpte = 1/m_fhelixe[2];
5524 m_fptmu = 1/m_fhelixmu[2];
5525 m_fptk = 1/m_fhelixk[2];
5526 m_fptp = 1/m_fhelixp[2];
5528 m_lpt = 1/m_lhelix[2];
5529 m_lpte = 1/m_lhelixe[2];
5530 m_lptmu = 1/m_lhelixmu[2];
5531 m_lptk = 1/m_lhelixk[2];
5532 m_lptp = 1/m_lhelixp[2];
5534 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5535 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5536 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5537 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5538 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5540 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5541 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5542 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5543 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5544 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5546 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5547 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5548 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5549 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5550 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5551 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5552 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5553 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5554 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5555 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5556 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5557 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5558 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5559 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5560 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5563 StatusCode sc1 = m_nt1->write();
5564 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5569 phi1 = (*iter_trk)->getFFi0();
5570 r1 = (*iter_trk)->getFDr();
5571 z1 = (*iter_trk)->getFDz();
5572 kap1 = (*iter_trk)->getFCpa();
5573 tanl1 = (*iter_trk)->getFTanl();
5576 p1 = sqrt(1+tanl1*tanl1)/kap1;
5577 the1 =
M_PI/2-atan(tanl1);
5578 }
else if(jj == 2) {
5579 phi2 = (*iter_trk)->getFFi0();
5580 r2 = (*iter_trk)->getFDr();
5581 z2 = (*iter_trk)->getFDz();
5582 kap2 = (*iter_trk)->getFCpa();
5583 tanl2 = (*iter_trk)->getFTanl();
5586 p2 = sqrt(1+tanl2*tanl2)/kap1;
5587 the2 =
M_PI/2-atan(tanl2);
5595 m_delthe = the1 + the2;
5598 StatusCode sc2 = m_nt2->write();
5599 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5607 cout <<
"Kalfitting finished " << std::endl;
5620 MsgStream log(
msgSvc(), name());
5637 double pp = track_lead.
momentum().mag();
5641 for(
int l_mass = 0, flg = 1; l_mass < nmass;
5642 l_mass++, flg <<= 1) {
5644 if (!(
mhyp_ & flg))
continue;
5645 if (l_mass ==
lead_)
continue;
5654 if(
debug_ == 4) cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
5662 TrasanTRK.
pivot[2]);
5663 HepVector a_trasan(5);
5664 for(
int i = 0; i < 5; i++)
5665 a_trasan[i] = TrasanTRK.
helix[i];
5667 HepSymMatrix ea_trasan(5);
5668 for(
int i = 0, k = 0; i < 5; i++) {
5669 for(
int j = 0; j <= i; j++) {
5671 ea_trasan[j][i] = ea_trasan[i][j];
5675 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
5678 double fiTerm = TrasanTRK.
fiTerm;
5679 track.
pivot(track.
x(fiTerm));
5680 HepSymMatrix Eakal(5,0);
5682 double costheta = track.
a()[4] / sqrt(1.0 + track.
a()[4]*track.
a()[4]);
5696 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
5698 fillTds(TrasanTRK, track, kaltrk, l_mass);
5711 HepSymMatrix ea_first(kaltrk->
getFError());
5714 for(
int i = 0; i < 5; i++)
5715 for(
int j = 0; j <= i; j++)
5716 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
5717 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
5723 cout <<
" Backward fitting flag:" <<
back_<< endl;
5724 cout <<
"track_back pivot " << track_back.
pivot()
5725 <<
" track_lead kappa " << track_lead.
a()[2]
5729 if (
back_ && track_lead.
a()[2] != 0 &&
5730 1/fabs(track_lead.
a()[2]) >
pT_) {
5734 double p_kaon(0), p_proton(0);
5735 if (!(kaltrk->
getStat(0,3))) {
5738 track_back.
p_kaon(p_kaon);
5740 p_kaon = 1 / fabs(track_back.
a()[2]) *
5741 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5742 track_back.
p_kaon(p_kaon);
5744 if (!(kaltrk->
getStat(0,4))) {
5745 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
5749 p_proton = 1 / fabs(track_back.
a()[2]) *
5750 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5757 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
5771 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass, segcol, 1);
5786 for(
int pid = 0; pid < nmass;
5788 if (pid ==
lead_)
continue;
5796 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
5798 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
5799 <<
"Length of the track: "<< kaltrk->
getLength(2)
5800 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
5801 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
5802 <<
" "<< kaltrk->
getChisq(1,2) << endreq
5803 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
5804 <<
" "<< kaltrk->
getNdf(1,2) << endreq
5808 kalcol->push_back(kaltrk);
5821 MsgStream log(
msgSvc(), name());
5826 cout <<
"track_first pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
5833 cout <<
"after inner wall, track_ip pivot "<<track_first.
pivot()<<
" helix "<<track_first.
a()<<endl;
5839 double pp = track_lead.
momentum().mag();
5844 for(
int l_mass = 0, flg = 1; l_mass < nmass;
5845 l_mass++, flg <<= 1) {
5847 if (!(
mhyp_ & flg))
continue;
5848 if (l_mass ==
lead_)
continue;
5859 cout<<
"complete_track..REFIT ASSUMPION " << l_mass << endl;
5869 TrasanTRK.
pivot[2]);
5871 HepVector a_trasan(5);
5872 for(
int i = 0; i < 5; i++){
5873 a_trasan[i] = TrasanTRK.
helix[i];
5876 HepSymMatrix ea_trasan(5);
5877 for(
int i = 0, k = 0; i < 5; i++) {
5878 for(
int j = 0; j <= i; j++) {
5880 ea_trasan[j][i] = ea_trasan[i][j];
5884 KalFitTrack track(x_trasan,a_trasan, ea_trasan, l_mass, chiSq, nhits);
5887 double fiTerm = TrasanTRK.
fiTerm;
5888 track.
pivot(track.
x(fiTerm));
5890 HepSymMatrix Eakal(5,0);
5892 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5904 fillTds_ip(TrasanTRK, track_z, kaltrk, l_mass);
5906 fillTds(TrasanTRK, track, kaltrk, l_mass);
5915 HepSymMatrix ea_first(kaltrk->
getFError());
5919 for(
int i = 0; i < 5; i++)
5920 for(
int j = 0; j <= i; j++)
5921 ea_first[i][j] = fac[i]*fac[j]*ea_first[i][j];
5922 KalFitTrack track(x_first, a_first, ea_first, 2, 0, 0);
5938 cout <<
" Backward fitting flag:" <<
back_<< endl;
5939 cout <<
"track_back pivot " << track_back.
pivot()
5940 <<
" track_lead kappa " << track_lead.
a()[2]
5944 if (
back_ && track_lead.
a()[2] != 0 &&
5945 1/fabs(track_lead.
a()[2]) >
pT_) {
5950 double p_kaon(0), p_proton(0);
5952 if (!(kaltrk->
getStat(0,3))) {
5955 track_back.
p_kaon(p_kaon);
5957 p_kaon = 1 / fabs(track_back.
a()[2]) *
5958 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5959 track_back.
p_kaon(p_kaon);
5961 if (!(kaltrk->
getStat(0,4))) {
5962 p_proton = 1 / fabs(kaltrk->
getZHelixP()[2]) *
5966 p_proton = 1 / fabs(track_back.
a()[2]) *
5967 sqrt(1 + track_back.
a()[4]*track_back.
a()[4]);
5975 for(
int l_mass = 0; l_mass < nmass; l_mass++) {
5980 fillTds_back(track_seed, kaltrk, TrasanTRK, l_mass,segcol);
6003 log << MSG::DEBUG <<
"registered MDC Kalmantrack:"
6005 <<
" Mass of the fit: "<< kaltrk->
getMass(2)<< endreq
6006 <<
"Length of the track: "<< kaltrk->
getLength(2)
6007 <<
" Tof of the track: "<< kaltrk->
getTof(2) << endreq
6008 <<
"Chisq of the fit: "<< kaltrk->
getChisq(0,2)
6009 <<
" "<< kaltrk->
getChisq(1,2) << endreq
6010 <<
"Ndf of the fit: "<< kaltrk->
getNdf(0,2)
6011 <<
" "<< kaltrk->
getNdf(1,2) << endreq
6015 kalcol->push_back(kaltrk);
6024 for (
int i=0; i<5; i++) {
6025 for(
int j = 1; j<i+2;j++) {
6027 Eakal(j,i+1) = Eakal(i+1,j);
6031 if (
debug_ == 4) cout<<
"initialised Ea.. "<<Eakal<<endl;
6039 for (
int i=0; i<5; i++) {
6040 for(
int j = 1; j<i+2;j++) {
6042 Eakal(j,i+1) = Eakal(i+1,j);
6044 Eakal(1,1) = Eakal(1,1)*
gain1_;
6045 Eakal(2,2) = Eakal(2,2)*
gain2_;
6046 Eakal(3,3) = Eakal(3,3)*
gain3_;
6047 Eakal(4,4) = Eakal(4,4)*
gain4_;
6048 Eakal(5,5) = Eakal(5,5)*
gain5_;
6114 if (
debug_ == 4) cout<<
"initialised Eakal.. "<<Eakal<<endl;
6123 cout <<
"start_seed begin... " << std::endl;
6125 Hep3Vector x_init(track.
pivot());
6126 HepSymMatrix Ea_init(5,0);
6127 Ea_init = track.
Ea();
6128 HepVector a_init(5);
6132 unsigned int nhit_included(10);
6144 unsigned int nhit = track.
HitsMdc().size();
6147 for (
int ktrial = 0; ktrial < 8; ktrial++) {
6150 track.
pivot(x_init);
6160 HepSymMatrix Eakal(5,0);
6164 for(
unsigned i=0 ; i < nhit; i++ ){
6167 if (i<3) lr_decision = LR[ktrial][i];
6168 HitMdc.
LR(lr_decision);
6169 if (i<nhit_included)
6182 cout <<
"---- Result for " << ktrial <<
" case : chi2 = " << track.
chiSq()
6183 <<
", nhits included = " << track.
nchits() <<
" for nhits available = "
6184 << nhit << std::endl;
6186 if (track.
chiSq() < chi2_min &&
6187 (track.
nchits() == nhit_included || track.
nchits() == nhit)){
6188 chi2_min = track.
chiSq();
6195 cout <<
"*** i_min = " << i_min <<
" with a chi2 = " << chi2_min << std::endl;
6197 for(
unsigned i=0 ; i < nhit; i++ ){
6200 if (i_min >= 0 && i < 3)
6201 lr_decision = LR[i_min][i];
6202 HitMdc.
LR(lr_decision);
6206 track.
pivot(x_init);
6217 for(
unsigned i=0 ; i < 3; i++ ){
6219 cout <<
" LR(" << i <<
") = " << HitMdc.
LR()
6231 if(
debug_ == 4) cout<<
"Begining to clear Tables ...."<<endl;
6235 vector<MdcRec_trk>::iterator tit=mdcMgr->begin();
6236 for( ; tit != mdcMgr->end(); tit++) {
6237 vector<MdcRec_wirhit*>::iterator vit= tit->hitcol.begin() ;
6238 for(; vit != tit->hitcol.end(); vit++) {
6244 mdc_addMgr->clear();
6256bool KalFitAlg::order_rechits(
const SmartRef<RecMdcHit>& m1,
const SmartRef<RecMdcHit>& m2) {
6337// --- print out GEM hits
6338if(giveCout) cout<<endl<<"------- GEM hits: --------"<<endl;
6339//for(int i=0; i<useNCGem_; i++)
6340for(int i=0; i<4; i++)
6342 if(giveCout) cout<<"### layer "<<i<<endl;
6343 //int nhits = myGemHitCol[i].size();
6344 KalFitGemHitCol::const_iterator it_GemHit = (myGemHitCol[i]).begin();
6345 for(int j=0; it_GemHit!=(myGemHitCol[i]).end(); it_GemHit++, j++)
6347 if(giveCout) cout<<" hit "<<(j+1)<<": phi, z, r = "<<it_GemHit->getPhi()<<", "<<it_GemHit->getZ()<<", "<<it_GemHit->getR()<<endl;
6348 if(ntuple_&1) m_rGem[i] = it_GemHit->getR();
6358 double rMaxInnerWall=6.42005+0.00495;
6362 rMaxInnerWall=m_innerWall[0].radius();
6365 rMaxInnerWall=_BesKalmanFitWalls[0].radius();
6368 HepPoint3D posAtInnerWall = track.
x(dPhiToInnerWall);
6369 double pathToInnerWall;
6370 track.
pivot_numf(posAtInnerWall, pathToInnerWall);
6373 if(pathToInnerWall>0)
6375 if(
muls_) track.
ms(pathToInnerWall, _BesKalmanFitMaterials[0], way);
6376 if(
loss_) track.
eloss(pathToInnerWall, _BesKalmanFitMaterials[0], way);
6378 if(useNCGem_>0) m_innerWall[0].updateTrack(track, way);
6740 for(
int i=0; i<3; i++) myGemHitCol[i].clear();
6762 SmartDataPtr<RecCgemClusterCol> cgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
6765 RecCgemClusterCol::const_iterator it = cgemClusterCol->begin();
6766 for(;it!=cgemClusterCol->end();it++){
6768 if((*it)->getflag()==2){
6769 layer = (*it)->getlayerid();
6770 sheet = (*it)->getsheetid();
6771 phi = (*it)->getrecphi();
6772 v = (*it)->getrecv();
6774 myGemHitCol[layer].push_back(aKalFitGemHit);
6778 cout <<
"Gem information" << endl;
6779 cout <<
"first layer " << myGemHitCol[0].size() << endl;
6780 cout <<
"second layer" << myGemHitCol[1].size() << endl;
6781 cout <<
"third layer" << myGemHitCol[2].size() << endl;
6828 bool prtGemHits=
false;
6830 if(prtGemHits) cout<<endl<<
"------- GEM hits: --------"<<endl;
6831 for(
int i=0; i<3; i++)
6833 if(prtGemHits) cout<<
"### layer "<<i<<endl;
6835 KalFitGemHitCol::const_iterator it_GemHit = (myGemHitCol[i]).begin();
6836 for(
int j=0; it_GemHit!=(myGemHitCol[i]).end(); it_GemHit++, j++)
6838 if(prtGemHits) cout<<
" hit "<<(j+1)<<
": phi, z, r = "<<it_GemHit->getPhi()<<
", "<<it_GemHit->getZ()<<
", "<<it_GemHit->getR()<<endl;
6855 std::map<double,double> esti1_r;
6856 std::map<double,double> esti2_r;
6857 std::map<double,double> esti1_phi;
6858 std::map<double,double> esti2_phi;
6859 std::map<double,double> esti1_z;
6860 std::map<double,double> esti2_z;
6861 std::map<double,double> diff1_phi;
6862 std::map<double,double> diff1_z;
6863 std::map<double,double> diff2_phi;
6864 std::map<double,double> diff2_z;
6865 std::map<double,HepPoint3D> posnew;
6866 std::map<double,HepSymMatrix> eanew;
6867 std::map<double,HepVector> anew;
6868 std::map<double,double> meas_r;
6869 std::map<double,double> meas_phi;
6870 std::map<double,double> meas_z;
6871 std::map<double,double> meas_phierr;
6872 std::map<double,double> esti_phierr;
6873 std::map<double,double> meas_zerr;
6874 std::map<double,double> esti_zerr;
6876 if(
debug_==4){ cout <<
"wall size" << _BesKalmanFitWalls.size() << endl;}
6877 for(
int iLayer=2; iLayer>=0; iLayer--)
6895 meas_phierr.clear();
6897 esti_phierr.clear();
6900 int nHits=myGemHitCol[iLayer].size();
6903 if (nHits <= 0)
continue;
6906 vector<KalFitGemHit>::iterator
iter;
6910 innerwall(track,hypo,way,64+iLayer*(-32),88+iLayer*(-32));
6913 for(
iter=myGemHitCol[iLayer].begin();
iter!=myGemHitCol[iLayer].end();
iter++ )
6916 track.
set(track_start.
pivot(),track_start.
a(),track_start.
Ea());
6918 double rGem=(*iter).getR();
6919 double phiGem=(*iter).getPhi();
6920 double zGem=(*iter).getZ();
6921 HepVector v_measu=(*iter).getVecPos();
6934 double x0=posEstiAtGem.x();
6935 double y0=posEstiAtGem.y();
6943 if(
muls_) track.
ms(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
6944 if(
loss_) track.
eloss(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
6947 double rGemEsti=posEstiAtGem.perp();
6948 double phiGemEsti=posEstiAtGem.phi();
6949 double zGemEsti=posEstiAtGem.z();
6950 HepVector v_estim(2,0);
6951 v_estim(1)=phiGemEsti;
6952 v_estim(2)=zGemEsti;
6959 const HepSymMatrix& Ea = track.
Ea();
6960 HepVector v_a = track.
a();
6963 double kappa=v_a(3);
6968 HepMatrix
H(2, 5, 0);
6971 H(1,1)=-y0*
cos(phi0)/(y0*y0+x0*x0)+x0*
sin(phi0)/(x0*x0+y0*y0);
6974 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)));
6975 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));
6976 H(2,3)=
m_alpha/(kappa*kappa)*tl*dPhi;
6980 HepMatrix H_T=
H.T();
6985 HepSymMatrix V=(*iter).getErrMat();
6988 HepMatrix HEaH_T=HEa*H_T;
6990 HepMatrix Vinv=(V+HEaH_T).inverse(ierr);
6991 if(ierr!=0) cout<<
"error in HEaH_T.inverse()!"<<endl;
6993 HepMatrix K=Ea*H_T*Vinv;
6995 HepSymMatrix EaNew(5,0);
6996 EaNew.assign(Ea-K*
H*Ea);
6998 HepVector v_diff=v_measu-v_estim;
7000 double delPhi=v_diff(1);
7001 if(fabs(v_diff(1))>6.28) v_measu(1)=-1*v_measu(1);
7008 HepVector aNew=v_a+K*v_diff;
7014 HepPoint3D posEstiAtGem_new = track_test.
x(dPhiToGem_New);
7015 double phiGemEsti_new=posEstiAtGem_new.phi();
7016 double zGemEsti_new=posEstiAtGem_new.z();
7017 double x0_new = posEstiAtGem_new.x();
7018 double y0_new = posEstiAtGem_new.y();
7019 HepVector v_estim_new(2,0);
7020 v_estim_new(1)=phiGemEsti_new;
7021 v_estim_new(2)=zGemEsti_new;
7025 v_diff=v_measu-v_estim_new;
7033 track_test.
pivot_numf(posEstiAtGem_new, pathInGem);
7034 HepVector v_a_new = track_test.
a();
7036 HepSymMatrix Ea_new = track_test.
Ea();
7038 double phi0_new = v_a_new(2);
7039 double kappa_new=v_a_new(3);
7040 HepMatrix H_new(2, 5, 0);
7049 HepMatrix H_new_T=H_new.T();
7056 HepSymMatrix R(2,0);
7058 R.assign(V-
H*EaNew*H_T);
7066 HepVector v_dChi2=v_diff.T()*R.inverse(ierr)*v_diff;
7067 if(ierr!=0) cout<<
"error in R.inverse()!"<<endl;
7070 double dChi2=v_dChi2(1);
7074 cout<<
"pivot: "<<posEstiAtGem<<endl;
7075 cout<<
"new pivot: "<<posEstiAtGem<<endl;
7076 cout<<
"a: "<<v_a<<endl;
7077 cout<<
"new a: "<<aNew<<endl;
7078 cout<<
"Ea: "<<Ea<<endl;
7079 cout<<
"new Ea: "<<EaNew<<endl;
7080 cout<<
"dchi2= "<<dChi2<<endl;
7115 esti1_r[dChi2]=rGemEsti;
7116 esti1_phi[dChi2]=phiGemEsti;
7117 esti1_z[dChi2]=zGemEsti;
7118 esti2_r[dChi2]=rGem;
7119 esti2_phi[dChi2]=phiGemEsti_new;
7120 esti2_z[dChi2]=zGemEsti_new;
7121 diff1_phi[dChi2]=phiGem-phiGemEsti;
7122 diff1_z[dChi2]=zGem-zGemEsti;
7123 diff2_phi[dChi2]=phiGem-phiGemEsti_new;
7124 diff2_z[dChi2]=zGem-zGemEsti_new;
7125 posnew[dChi2]=posEstiAtGem;
7129 meas_phi[dChi2]=phiGem;
7132 meas_phierr[dChi2]=sqrt(V[0][0]);
7134 meas_zerr[dChi2]=sqrt(V[1][1]);
7136 esti_phierr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[0][0]);
7138 esti_zerr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[1][1]);
7147 track.
set((posnew.begin())->second,(anew.begin())->second,(eanew.begin()->second));
7151 HepPoint3D new_posEstiAtGem = track.
x(new_dPhiToGem);
7152 double new_pathInGem;
7155 track.
pivot_numf(new_posEstiAtGem, new_pathInGem);
7160 if(
muls_) track.
ms(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7161 if(
loss_) track.
eloss(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7165 innerwall(track,hypo,way,90+(-32)*iLayer,95+(-32)*iLayer);}
7167 innerwall(track,hypo,way,90+(-32)*iLayer,94+(-32)*iLayer);}
7175 m_esti1_r[ii]=(esti1_r.begin())->second;
7177 m_esti1_phi[ii]=(esti1_phi.begin())->second;
7178 m_esti1_z[ii]=(esti1_z.begin())->second;
7179 m_esti2_r[ii]=(esti2_r.begin())->second;
7180 m_esti2_phi[ii]=(esti2_phi.begin())->second;
7181 m_esti2_z[ii]=(esti2_z.begin())->second;
7182 m_diff1_phi[ii]=(diff1_phi.begin())->second;
7183 m_diff1_z[ii]=(diff1_z.begin())->second;
7184 m_diff2_phi[ii]=(diff2_phi.begin())->second;
7185 m_diff2_z[ii]=(diff2_z.begin())->second;
7186 m_Gchi2[ii]=(esti1_r.begin())->first;
7187 m_meas_r[ii]=(meas_r.begin())->second;
7188 m_meas_phi[ii]=(meas_phi.begin())->second;
7190 m_meas_z[ii]=(meas_z.begin())->second;
7191 m_meas_phierr[ii]=(meas_phierr.begin())->second;
7192 m_meas_zerr[ii]=(meas_zerr.begin())->second;
7193 m_esti_phierr[ii]=(esti_phierr.begin())->second;
7194 m_esti_zerr[ii]=(esti_zerr.begin())->second;
7195 m_GemLayer[ii]=iLayer;
7207 if(m_meas_phi[0]>0 && m_meas_phi[0]<1.5707963) m_qua=0;
7208 if(m_meas_phi[0]>1.5707963 && m_meas_phi[0]<3.1415926) m_qua=1;
7209 if(m_meas_phi[0]>-3.1415926&& m_meas_phi[0]<-1.5707963) m_qua=2;
7210 if(m_meas_phi[0]>-1.5707963&& m_meas_phi[0]<0) m_qua=3;
7233 meas_phierr.clear();
7235 esti_phierr.clear();
7241 SmartDataPtr<Event::EventHeader> evtHeader(eventSvc(),
"/Event/EventHeader");
7243 m_evt3 = evtHeader->eventNumber();
7255 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
7257 cout <<
"Could not find Event Header" << endl;
7259 int eventNo = eventHeader->eventNumber();
7260 int runNo = eventHeader->runNumber();
7261 SmartDataPtr<RecMdcTrackCol> recMdcTrack(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
7263 cout <<
"Could not find RecMdcTrackCol" << endl;
7265 RecMdcTrackCol::iterator itrk = recMdcTrack->begin();
7267 for(; itrk != recMdcTrack->end(); itrk++){
7268 if((*itrk)->trackId() == track.
trasan_id() && (*itrk)->getNcluster() > 0) {
7269 clusterRefVec = (*itrk)->getVecClusters();
7277 if(
debug_==4){ cout<<
"wall size"<<_BesKalmanFitWalls.size()<<endl;}
7280 ClusterRefVec::iterator itCluster;
7281 ClusterRefVec::iterator itStartFlag;
7282 ClusterRefVec::iterator itStopFlag;
7286 double R_o_cgem=_BesKalmanFitWalls[0].radius();
7288 itStartFlag = clusterRefVec.begin();
7289 itStopFlag = clusterRefVec.end();
7294 itStartFlag = clusterRefVec.end()-1;
7295 itStopFlag = clusterRefVec.begin()-1;
7305 for(itCluster = itStartFlag; itCluster !=itStopFlag; itCluster = itCluster + step){
7306 int layer = (*itCluster)->getlayerid();
7309 innerwall(track,hypo,way,lastR,dmR,index);
7315 innerwall(track,hypo,way,lastR,R_o_cgem,index);
7319 innerwall(track,hypo,way,lastR,0,index);
7466 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
7468 cout <<
"Could not find Event Header" << endl;
7470 int eventNo = eventHeader->eventNumber();
7471 int runNo = eventHeader->runNumber();
7475 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
7477 cout <<
"Could not find McParticle" << endl;
7480 HepVector mc_a(5,0);
7481 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
7482 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
7483 if(!(*i_mcTrk)->primaryParticle())
continue;
7484 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
7485 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
7487 double charge = 0.0;
7488 int pid = (*i_mcTrk)->particleProperty();
7490 charge = m_particleTable->particle( pid )->charge();
7491 }
else if ( pid <0 ) {
7492 charge = m_particleTable->particle( -pid )->charge();
7497 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
7499 Helix mchelix(pos2, mom2, charge);
7502 for(
int j =0; j<5; j++) {
7504 mc_a[j] = mchelix.
a()[j];
7510 Helix MChelix(IP,mc_a);
7512 m_diff_dr = track.
a()[0]-MChelix.
a()[0];
7513 m_diff_phi0 = track.
a()[1]-MChelix.
a()[1];
7514 m_diff_kappa = track.
a()[2]-MChelix.
a()[2];
7515 m_diff_dz = track.
a()[3]-MChelix.
a()[3];
7516 m_diff_tanl = track.
a()[4]-MChelix.
a()[4];
7517 if(track.
a()[2] != 0 && MChelix.
a()[2] != 0)
7518 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];
7521 StatusCode sc12 = m_nt12->write();
7522 if( sc12.isFailure() ) cout<<
"Ntuple12 filling failed!"<<endl;
7526 bool printwalls =
false;
7528 double material_quantity(0);
7529 int size = _BesKalmanFitWalls.size();
7530 for(
int i=0;i<size;i++)
7532 double ro = _BesKalmanFitWalls[i].radius();
7533 double ri = _BesKalmanFitWalls[i].rmin();
7534 double thick = ro-ri;
7536 double X0 = material.
X0();
7538 double Z = material.
Z();
7540 material_quantity += 100*thick/
X0;
7541 cout<<
"----------------------------------------------------------------------------------------------------"<<endl;
7542 cout<<
"wall"<<i+1<<
" Z:"<<setw(12)<<Z<<endl;
7543 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;
7544 cout<<
"----------------------------------------------------------------------------------------------------"<<endl;
7550 cout<<
"============================================================"<<endl;
7551 cout<<
"total material quantity="<<material_quantity<<endl;
7552 cout<<
"============================================================"<<endl;
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
double sin(const BesAngle a)
double cos(const BesAngle a)
ObjectVector< RecMdcKalHelixSeg > RecMdcKalHelixSegCol
SmartRefVector< RecMdcKalHelixSeg > HelixSegRefVec
ObjectVector< RecMdcKalTrack > RecMdcKalTrackCol
SmartRefVector< RecCgemCluster > ClusterRefVec
SmartRefVector< RecMdcHit > HitRefVec
HepGeom::Point3D< double > HepPoint3D
const double dor_layer[3]
const double dmr_layer[3]
**********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
HepGeom::Point3D< double > HepPoint3D
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)