CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemMdcFitAlg.cxx
Go to the documentation of this file.
1#include "CgemMdcFitAlg/CgemMdcFitAlg.h"
2#include "GaudiKernel/MsgStream.h"
3
4#include "GaudiKernel/IPartPropSvc.h"
5
6//#include "GaudiKernel/MsgStream.h"
7//#include "GaudiKernel/AlgFactory.h"
8//#include "GaudiKernel/ISvcLocator.h"
9//#include "GaudiKernel/SmartDataPtr.h"
10//#include "GaudiKernel/IDataProviderSvc.h"
11//#include "GaudiKernel/PropertyMgr.h"
12//#include "GaudiKernel/IService.h"
13//#include "MdcRawEvent/MdcDigi.h"
14//#include "Identifier/Identifier.h"
15//#include "Identifier/MdcID.h"
16//#include "CLHEP/Units/PhysicalConstants.h"
17
18#include "HepPDT/ParticleDataTable.hh"
19#include "Identifier/MdcID.h"
20#include "BField/BField.h"
21#include "MagneticField/IMagneticFieldSvc.h"
22#include "TrkFitter/TrkContextEv.h"
23#include "EventModel/EventHeader.h"
24#include "EvTimeEvent/RecEsTime.h"
25#include "McTruth/McParticle.h"
26#include "TrackUtil/Helix.h"
27#include "MdcRawEvent/MdcDigi.h"
28#include "MdcData/MdcHit.h"
29#include "MdcData/MdcRecoHitOnTrack.h"
30#include "CgemData/CgemHitOnTrack.h"
31#include "TrkFitter/TrkHelixMaker.h"
32#include "TrkFitter/TrkCircleMaker.h"
33#include "TrkFitter/TrkLineMaker.h"
34#include "TrkBase/TrkHitList.h"
35#include "TrkBase/TrkExchangePar.h"
36#include "TrkBase/TrkFit.h"
37#include "TrkBase/TrkRecoTrk.h"
38#include "MdcTrkRecon/MdcSeg.h"
39#include "MdcTrkRecon/MdcTrackParams.h"
40
41#include <iomanip>
42//#include <iostream>
43#include <fstream>
44
45using namespace std;
46using namespace Event;
47const TrkContextEv* CgemMdcFitAlg::m_context=NULL;
48const MdcCalibFunSvc* CgemMdcFitAlg::m_mdcCalibFunSvc=NULL;
49const CgemCalibFunSvc* CgemMdcFitAlg::m_cgemCalibFunSvc=NULL;
50const CgemGeomSvc* CgemMdcFitAlg::m_cgemGeomSvc=NULL;
51const MdcDetector* CgemMdcFitAlg::m_mdcDetector=NULL;
52
53//double CgemMdcFitAlg::m_qualityFactor=3;
54//double CgemMdcFitAlg::m_dropTrkDrCut=1; //no limit of dr dz
55//double CgemMdcFitAlg::m_dropTrkDzCut=10;
56//double CgemMdcFitAlg::m_dropTrkNhitCut=9;
57//double CgemMdcFitAlg::m_dropTrkChi2Cut=99999;
58//double CgemMdcFitAlg::m_dropTrkChi2NdfCut=99999;
59
60CgemMdcFitAlg::CgemMdcFitAlg(const std::string& name, ISvcLocator* pSvcLocator):
61 Algorithm(name, pSvcLocator)
62{
63 declareProperty("debug" ,m_debug = 0 );
64 declareProperty("tuple" ,m_tuple = 0 );
65 declareProperty("change" ,m_changeTDS = 2 );// 0:not change, 1:update, 2:replace
66 declareProperty("debugArbHit" ,m_debugArbHit= 0 );
67 declareProperty("fit2D" ,m_fit2D= 1 );
68 declareProperty("drCut" ,m_drCut= 1. );
69 declareProperty("dzCut" ,m_dzCut= 10. );
70 declareProperty("filter" ,m_filter= 0 );
71 declareProperty("eventFile" ,m_evtFile= "EventList" );
72 declareProperty("compareTrack" ,m_compareTrack= 0 );
73 declareProperty("removeBadTrack",m_removeBadTrack= 0 );
74}
75
76StatusCode CgemMdcFitAlg::initialize()
77{
78 MsgStream log(msgSvc(), name());
79 log << MSG::INFO << "In CgemMdcFitAlg initialize()" << endreq;
80
81 StatusCode sc;
82
83 IMagneticFieldSvc* pIMF;
84 sc = service ("MagneticFieldSvc",pIMF);
85 if(sc!=StatusCode::SUCCESS) {
86 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
87 }
88 BField* bfield = new BField(pIMF);
89 log << MSG::INFO << "field z = "<<bfield->bFieldNominal()<< endreq;
90 m_context = new TrkContextEv(bfield);
91
92 IMdcCalibFunSvc* imdcCalibSvc;
93 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
94 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
95 if ( sc.isFailure() ){
96 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
97 return StatusCode::FAILURE;
98 }
99
100 ICgemCalibFunSvc* icgemCalibSvc;
101 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
102 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
103 if ( sc.isFailure() ){
104 log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
105 return StatusCode::FAILURE;
106 }
107
108 ICgemGeomSvc* icgemGeomSvc;
109 sc = service ("CgemGeomSvc",icgemGeomSvc);
110 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*>(icgemGeomSvc);
111 if ( sc.isFailure() ){
112 log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
113 return StatusCode::FAILURE;
114 }
115
116 m_mdcDetector = new MdcDetector();
117 if(m_tuple) sc = bookTuple();
118
119 return StatusCode::SUCCESS;
120}
121//cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
122
123StatusCode CgemMdcFitAlg::execute()
124{
125//cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
126
127 MsgStream log(msgSvc(), name());
128 log << MSG::INFO << "In CgemMdcFitAlg execute()" << endreq;
129
130 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
131 if (!eventHeader) {
132 log << MSG::FATAL << "Could not find Event Header" << endreq;
133 return StatusCode::FAILURE;
134 }
135 int run = eventHeader->runNumber();
136 int evt = eventHeader->eventNumber();
137 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
138 if (!recMdcTrackCol){
139 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
140 return StatusCode::SUCCESS;
141 }
142
143 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
144 if (!evTimeCol) {
145 log << MSG::WARNING<< "Could not retrieve RecEsTimeCol" << endreq;
146 return StatusCode::SUCCESS;
147 }else{
148 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
149 if (iter_evt != evTimeCol->end()){
150 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
151 }
152 }
153
154 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(),"/Event/Digi/MdcDigiCol");
155 if (!mdcDigiCol) {
156 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
157 return StatusCode::SUCCESS;
158 }
159
160 if(m_filter){
161 ifstream lowPt_Evt;
162 lowPt_Evt.open(m_evtFile.c_str());
163 vector<int> evtlist;
164 int evtNum;
165 while( lowPt_Evt >> evtNum) {
166 evtlist.push_back(evtNum);
167 }
168 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),evt);
169 if( iter_evt == evtlist.end() ) { setFilterPassed(false); return StatusCode::SUCCESS; }
170 else setFilterPassed(true);
171 }
172
173 if(m_debug)cout<<endl;
174 if(m_debug)cout<<"run:"<<run<<" event:"<<evt<<endl;
175 if(m_tuple)getMcTruth();
176 if(m_debug)cout<<"tracks before fitting: "<<recMdcTrackCol->size()<<endl;
177 //if(recMdcTrackCol->size()==0)return StatusCode::SUCCESS;
178 if(m_tuple){
179 m_run = run;
180 m_event = evt;
181 m_nTrkRec = recMdcTrackCol->size();
182 }
183
184 for(RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();itTrk != recMdcTrackCol->end();itTrk++)
185 {
186 HitRefVec vechits = (*itTrk)->getVecHits();
187 HitRefVec::iterator itHit = vechits.begin();
188 for( ; itHit != vechits.end(); itHit++)
189 {
190 double tdc = (*itHit)->getTdc();
191 double adc = (*itHit)->getAdc();
192 const Identifier mdcid = (*itHit)->getMdcId();
193 int layer = MdcID::layer(mdcid);
194 int wire = MdcID::wire(mdcid);
195 //cout<<"layer: "<<layer<<" wire: "<<wire<<" tdc: "<<tdc<<" adc: "<<adc<<endl;
196 MdcDigiCol::iterator iterDigi = mdcDigiCol->begin();
197 for(; iterDigi != mdcDigiCol->end(); iterDigi++ ) {
198 if((*iterDigi)->identify() == mdcid) break;
199 }
200 if(tdc==0)tdc = (*iterDigi)->getTimeChannel();
201 //cout<<" tdc= "<<(*iterDigi)->getTimeChannel()<<" adc= "<<(*iterDigi)->getChargeChannel()<<endl;
202
203 bool used = false;
204 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
205 for(;imdcHit != mdcHitCol.end();imdcHit++){
206 if((*imdcHit)->mdcId() == mdcid){
207 used = true;
208 //cout<<"layer: "<<layer<<" view: "<<mdcHit->layer()->view()<<" wire: "<<wire<<endl;
209 break;
210 }
211 }
212 if(used)continue;
213 //cout<<"layer: "<<layer<<" wire: "<<wire<<" tdc: "<<tdc<<" adc: "<<adc<<endl;
214 const MdcDigi *mdcDigi = new MdcDigi(mdcid,tdc,adc);
215 MdcHit *mdcHit = new MdcHit(mdcDigi,m_mdcDetector);
216 mdcHitCol.push_back(mdcHit);
217 }
218 }
219 //cout<<mdcHitCol.size()<<endl;
220 MdcTrackParams m_par;
221 if(m_debugArbHit>0 ) m_par.lPrint=8;
222 m_par.lRemoveInActive=1;
223 //m_par.lUseQualCuts=0;
224 //m_par.maxGapLength=99;
225 //m_par.maxChisq=99;
226 //m_par.minHits=99;
227 MdcTrackList mdcTrackList(m_par);
228
229 RecMdcTrackCol* trackList_tds;
230 RecMdcHitCol* hitList_tds;
231 //vector<MdcTrack*> vec_MdcTrack;
232 RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
233 for(int itrk=0;itTrk != recMdcTrackCol->end(); itrk++,itTrk++)
234 {
235/*
236 int trkId = (*itTrk)->trackId();
237 double dr = (*itTrk)->helix(0);
238 double phi0 = (*itTrk)->helix(1);
239 double kappa = (*itTrk)->helix(2);
240 double dz = (*itTrk)->helix(3);
241 double tanl = (*itTrk)->helix(4);
242 int chargeRec = (*itTrk)->charge() ;
243 int statRec = (*itTrk)->stat() ;
244 int nhitsRec = (*itTrk)->getNhits() ;
245 int nclusterRec = (*itTrk)->getNcluster() ;
246 int nsterRec = (*itTrk)->nster() ;
247 int ndofRec = (*itTrk)->ndof() ;
248 double chi2Rec = (*itTrk)->chi2() ;
249 double pxRec = (*itTrk)->px() ;
250 double pyRec = (*itTrk)->py() ;
251 double pzRec = (*itTrk)->pz() ;
252 double pxyRec = (*itTrk)->pxy() ;
253 double pRec = (*itTrk)->p() ;
254 double xRec = (*itTrk)->x() ;
255 double yRec = (*itTrk)->y() ;
256 double zRec = (*itTrk)->z() ;
257 double rRec = (*itTrk)->r() ;
258 double thetaRec = (*itTrk)->theta() ;
259 double phiRec = (*itTrk)->phi() ;
260 double fiTermRec = (*itTrk)->getFiTerm() ;
261 if(m_tuple){
262 m_trkidRec[itrk] = trkId;
263 m_drRec[itrk] = dr ;
264 m_phi0Rec[itrk] = phi0 ;
265 m_kappaRec[itrk] = kappa;
266 m_dzRec[itrk] = dz ;
267 m_tanlRec[itrk] = tanl ;
268 m_chargeRec[itrk] = chargeRec ;
269 m_statRec[itrk] = statRec ;
270 m_nhitsRec[itrk] = nhitsRec ;
271 m_nclusterRec[itrk]= nclusterRec;
272 m_nsterRec[itrk] = nsterRec ;
273 m_ndofRec[itrk] = ndofRec ;
274 m_chi2Rec[itrk] = chi2Rec ;
275 m_pxRec[itrk] = pxRec ;
276 m_pyRec[itrk] = pyRec ;
277 m_pzRec[itrk] = pzRec ;
278 m_pxyRec[itrk] = pxyRec ;
279 m_pRec[itrk] = pRec ;
280 m_xRec[itrk] = xRec ;
281 m_yRec[itrk] = yRec ;
282 m_zRec[itrk] = zRec ;
283 m_rRec[itrk] = rRec ;
284 m_thetaRec[itrk] = thetaRec ;
285 m_phiRec[itrk] = phiRec ;
286 m_fiTermRec[itrk] = fiTermRec ;
287 }
288*/
289 if(m_tuple){
290 m_trkidRec[itrk] = (*itTrk)->trackId();
291 m_drRec[itrk] = (*itTrk)->helix(0);
292 m_phi0Rec[itrk] = (*itTrk)->helix(1);
293 m_kappaRec[itrk] = (*itTrk)->helix(2);
294 m_dzRec[itrk] = (*itTrk)->helix(3);
295 m_tanlRec[itrk] = (*itTrk)->helix(4);
296 m_chargeRec[itrk] = (*itTrk)->charge() ;
297 m_statRec[itrk] = (*itTrk)->stat() ;
298 m_nhitsRec[itrk] = (*itTrk)->getNhits() ;
299 m_nclusterRec[itrk]= (*itTrk)->getNcluster() ;
300 m_nsterRec[itrk] = (*itTrk)->nster() ;
301 m_ndofRec[itrk] = (*itTrk)->ndof() ;
302 m_chi2Rec[itrk] = (*itTrk)->chi2() ;
303 m_pxRec[itrk] = (*itTrk)->px() ;
304 m_pyRec[itrk] = (*itTrk)->py() ;
305 m_pzRec[itrk] = (*itTrk)->pz() ;
306 m_pxyRec[itrk] = (*itTrk)->pxy() ;
307 m_pRec[itrk] = (*itTrk)->p() ;
308 m_xRec[itrk] = (*itTrk)->x() ;
309 m_yRec[itrk] = (*itTrk)->y() ;
310 m_zRec[itrk] = (*itTrk)->z() ;
311 m_rRec[itrk] = (*itTrk)->r() ;
312 m_thetaRec[itrk] = (*itTrk)->theta() ;
313 m_phiRec[itrk] = (*itTrk)->phi() ;
314 m_fiTermRec[itrk] = (*itTrk)->getFiTerm() ;
315 }
316
317 if(m_debug)cout<<"fitting track: "<<itrk<<endl;
318
319 //---make tracks
320 float chisum =0.;
321 bool permitFlips = true;
322 bool lPickHits = true;
323 int trkId = (*itTrk)->trackId();
324 double dr = (*itTrk)->helix(0);
325 double phi0 = (*itTrk)->helix(1);
326 double kappa = (*itTrk)->helix(2);
327 double dz = (*itTrk)->helix(3);
328 double tanl = (*itTrk)->helix(4);
329 if(m_debug)cout<<"helix parameters before fitting: "<<setw(12)<<dr<<" "<<setw(12)<<phi0<<" "<<setw(12)<<kappa<<" "<<setw(12)<<dz<<" "<<setw(12)<<tanl<<endl;
330
331 //TrkLineMaker linefactory;
332 //---fitting circle
333 if(m_fit2D){
334 TrkExchangePar circle(-dr,phi0+M_PI/2,kappa/333.567,0,0);
335
336 trackType = circleType;
337 TrkCircleMaker circlefactory;
338 TrkRecoTrk* circleTrk = circlefactory.makeTrack(circle, chisum, *m_context, m_bunchT0);
339 circlefactory.setFlipAndDrop(*circleTrk, permitFlips, lPickHits);
340 TrkErrCode fit_circle = fit((*itTrk), circleTrk, trackType);
341 TrkErrCode check_circle = check(circleTrk,trackType);
342 if (m_debug && fit_circle.failure()){cout<<"fit: "<<fit_circle<<endl;}//FIXME
343 if (m_debug && check_circle.failure()){cout<<"check: "<<check_circle<<endl;}//FIXME
344 if (fit_circle.failure() || check_circle.failure()){continue;}//FIXME
345
346 //---update initial fitting parameters of helix dr, phi0, kappa
347 TrkExchangePar par = circleTrk->fitResult()->helix(0.);
348 dr = -par.d0();
349 phi0 = par.phi0()-M_PI/2;
350 if(phi0<0)phi0+=2*M_PI;
351 kappa = par.omega()*333.567;
352 if(m_debug)cout<<"helix parameters after 2D fitting: "<<setw(12)<<dr<<" "<<setw(12)<<phi0<<" "<<setw(12)<<kappa<<" "<<setw(12)<<dz<<" "<<setw(12)<<tanl<<endl;
353 }
354
355 //--- fitting helix
356 TrkExchangePar helix(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
357
358 trackType = helixType;
359 TrkHelixMaker helixfactory;
360 TrkRecoTrk* helixTrk = helixfactory.makeTrack(helix, chisum, *m_context, m_bunchT0);
361 //const TrkId id = helixTrk->id();
362 //cout<<"======================================================"<<long(id)<<endl;
363 helixfactory.setFlipAndDrop(*helixTrk, permitFlips, lPickHits);
364 TrkErrCode fit_helix = fit((*itTrk), helixTrk, trackType);
365 TrkErrCode check_helix = check(helixTrk,trackType);
366 if (m_debug && fit_helix.failure()){cout<<"fit: "<<fit_helix<<endl;}//FIXME
367 if (m_debug && check_helix.failure()){cout<<"check: "<<check_helix<<endl;}//FIXME
368 if (fit_helix.failure() || check_helix.failure()){continue;}//FIXME
369
370 TrkExchangePar par = helixTrk->fitResult()->helix(0.);
371 dr = -par.d0();
372 phi0 = par.phi0()-M_PI/2;
373 if(phi0<0)phi0+=2*M_PI;
374 kappa = par.omega()*333.567;
375 dz=par.z0();
376 tanl=par.tanDip();
377 if(m_debug)cout<<"helix parameters after fitting: "<<setw(12)<<dr<<" "<<setw(12)<<phi0<<" "<<setw(12)<<kappa<<" "<<setw(12)<<dz<<" "<<setw(12)<<tanl<<endl<<endl;
378
379 //RecMdcTrack* recMdcTrack = (*itTrk);
380 int trkStat = 5;//track fi by CgemMdcFitAlg set stat=5
381 if(m_changeTDS ==1)updateTracks(trkId,helixTrk,(*itTrk),trkStat);
382 //mdcTrack->storeTrack(trkId, trackList_tds, hitList_tds, trkStat);
383 //vec_MdcTrack.push_back(mdcTrack);
384 //trackList_tds->push_back(*itTrk);
385 //MdcTrack *mdcTrack ;
386 if(m_changeTDS ==2){
387 MdcTrack *mdcTrack = new MdcTrack(helixTrk);
388 //mdcTrack = new MdcTrack(helixTrk);
389 mdcTrackList.append(mdcTrack);
390 }
391
392 }
393
394 if(m_changeTDS ==2){
395 mdcTrackList.setNoInner(true);
396 if(m_removeBadTrack){
397 int nDeleted=mdcTrackList.arbitrateHits();//FIXME
398 if(m_debug>0 && nDeleted>0)cout<<nDeleted<<" tracks have been deleted by MdcTrackList..arbitrateHits()"<<endl;//FIXME
399 }
400 if(m_compareTrack>0)compareTracks(mdcTrackList, m_dzCut);
401 storeTracks(trackList_tds, hitList_tds,mdcTrackList);
402 }
403
404 if(m_tuple){
405 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
406 if (!recMdcTrackCol){
407 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
408 return StatusCode::SUCCESS;
409 }
410 m_nTrkFit = recMdcTrackCol->size();
411 RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
412 for(int itrk=0;itTrk != recMdcTrackCol->end(); itrk++,itTrk++)
413 {
414 //RecMdcTrack* recMdcTrack = (*itTrk);
415 m_trkidFit[itrk] = (*itTrk)->trackId();
416 m_drFit[itrk] = (*itTrk)->helix(0);
417 m_phi0Fit[itrk] = (*itTrk)->helix(1);
418 m_kappaFit[itrk] = (*itTrk)->helix(2);
419 m_dzFit[itrk] = (*itTrk)->helix(3);
420 m_tanlFit[itrk] = (*itTrk)->helix(4);
421 m_chargeFit[itrk] = (*itTrk)->charge() ;
422 m_statFit[itrk] = (*itTrk)->stat() ;
423 m_nhitsFit[itrk] = (*itTrk)->getNhits() ;
424 m_nclusterFit[itrk]= (*itTrk)->getNcluster() ;
425 m_nsterFit[itrk] = (*itTrk)->nster() ;
426 m_ndofFit[itrk] = (*itTrk)->ndof() ;
427 m_chi2Fit[itrk] = (*itTrk)->chi2() ;
428 m_pxFit[itrk] = (*itTrk)->px() ;
429 m_pyFit[itrk] = (*itTrk)->py() ;
430 m_pzFit[itrk] = (*itTrk)->pz() ;
431 m_pxyFit[itrk] = (*itTrk)->pxy() ;
432 m_pFit[itrk] = (*itTrk)->p() ;
433 m_xFit[itrk] = (*itTrk)->x() ;
434 m_yFit[itrk] = (*itTrk)->y() ;
435 m_zFit[itrk] = (*itTrk)->z() ;
436 m_rFit[itrk] = (*itTrk)->r() ;
437 m_thetaFit[itrk] = (*itTrk)->theta() ;
438 m_phiFit[itrk] = (*itTrk)->phi() ;
439 m_fiTermFit[itrk] = (*itTrk)->getFiTerm() ;
440 }
441 }
442
443//cout<<"--------------------------------------------------------------------------------"<<endl;
444 if(m_tuple)ntuple_evt->write();
445 //if(m_debug)cout<<endl;
446 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
447 for(;imdcHit != mdcHitCol.end();imdcHit++){
448 delete (*imdcHit);
449 }
450 mdcHitCol.clear();
451 return StatusCode::SUCCESS;
452}
453//cout<<"--------------------------------------------------------------------------------"<<endl;
454
455StatusCode CgemMdcFitAlg::finalize()
456{
457 MsgStream log(msgSvc(), name());
458 log << MSG::INFO<< "In CgemMdcFitAlg finalize()" << endreq;
459 //delete m_bfield ;
460 return StatusCode::SUCCESS;
461}
462
463TrkErrCode CgemMdcFitAlg::fit(RecMdcTrack* recMdcTrack, TrkRecoTrk* trkRecoTrk, TrackType trackType)
464//template <typename T>
465//TrkErrCode CgemMdcFitAlg::fit(double& dr,double& phi0, double& kappa, double& dz, double& tanl, RecMdcTrack* recMdcTrack, TrkRecoTrk* trkRecoTrk, T &trackfactory, TrackType trkType)
466{
467 //---make tracks
468 //double dr = recMdcTrack->helix(0);
469 //double phi0 = recMdcTrack->helix(1);
470 //double kappa = recMdcTrack->helix(2);
471 //double dz = recMdcTrack->helix(3);
472 //double tanl = recMdcTrack->helix(4);
473 //TrkHelixMaker trackfactory;
474 //TrkLineMaker trackfactory;
475 //TrkCircleMaker trackfactory;
476 //if(typeid(trackfactory)==typeid(TrkHelixMaker)){
477 //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
478 //}else if(typeid(trackfactory)==typeid(TrkCircleMaker)){
479 //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,0,0);
480 //}else if(typeid(trackfactory)==typeid(TrkLineMaker)){
481 //FIXME
482 //}else return TrkErrCode(TrkErrCode::fail);
483 //float chisum =0.;
484 //bool permitFlips = true;
485 //bool lPickHits = true;
486 //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
487 //trkRecoTrk = trackfactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
488 //trackfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
489 TrkHitList* trkHitList = trkRecoTrk->hits();
490
491 //--- for ODC hits
492 double ddCut=1;
493 HitRefVec vechits = recMdcTrack->getVecHits();
494 HitRefVec::iterator itHit = vechits.begin();
495 for( ; itHit != vechits.end(); itHit++)
496 {
497 double flight = (*itHit)->getFltLen();
498 double tdc = (*itHit)->getTdc();
499 double adc = (*itHit)->getAdc();
500 const Identifier mdcid = (*itHit)->getMdcId();
501 int layer = MdcID::layer(mdcid);
502 int wire = MdcID::wire(mdcid);
503 MdcHit* mdcHit;
504 if(trackType == circleType){
505 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
506 for(;imdcHit != mdcHitCol.end();imdcHit++){
507 if((*imdcHit)->mdcId() == mdcid){
508 //mdcHit = *imdcHit;
509 tdc = (*imdcHit)->tdcIndex();
510 break;
511 }
512 }
513 const MdcDigi *mdcDigi = new MdcDigi(mdcid,tdc,adc);
514 mdcHit = new MdcHit(mdcDigi,m_mdcDetector);
515 }else if(trackType == helixType){
516 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
517 for(;imdcHit != mdcHitCol.end();imdcHit++){
518 if((*imdcHit)->mdcId() == mdcid){
519 mdcHit = *imdcHit;
520 break;
521 }
522 }
523 }else{
524 TrkErrCode err;
525 err.setFailure(10,"fit failed,track type error!");
526 return err;
527 }
528 //cout<<"layer: "<<layer<<" wire: "<<wire<<" tdc= "<<tdc<<" adc: "<<adc<<endl;
529
530 mdcHit->setCalibSvc(m_mdcCalibFunSvc);
531 mdcHit->setCountPropTime(true);
532
533 int newAmbig = 0;
534 MdcRecoHitOnTrack temp( *mdcHit, newAmbig, m_bunchT0*1.e9);
535 MdcHitOnTrack* newhot = &temp;
536 newhot->setFltLen( flight );
537 if(m_debug>1)cout<<"layer: "<<layer<<" view: "<<mdcHit->layer()->view()<<" wire: "<<wire<<endl;
538 if(mdcHit->driftTime(m_bunchT0,0)>1000)continue;
539 if(trackType==circleType && mdcHit->layer()->view())continue;
540 if(trackType==circleType && mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;
541 trkHitList->appendHot(newhot);
542 //hitList_tds->push_back((*itHit));
543 }
544 //---for CGEM clusters
545 ClusterRefVec clusterRefVec = recMdcTrack->getVecClusters();
546 ClusterRefVec::iterator itCluster = clusterRefVec.begin();
547 for( ; itCluster != clusterRefVec.end(); itCluster++)
548 {
549 RecCgemCluster *recCgemCluster = *itCluster;
550 CgemHitOnTrack* newhot = new CgemHitOnTrack(*recCgemCluster,m_bunchT0*1.e9);
551 newhot->setCgemGeomSvc(m_cgemGeomSvc);
552 newhot->setCgemCalibFunSvc(m_cgemCalibFunSvc);
553 trkHitList->appendHot(newhot);
554 if(m_debug>1)cout<<"layer: "<<recCgemCluster->getlayerid()<<" flag: "<<recCgemCluster->getflag()<<" clusterId: "<<recCgemCluster->getclusterid()<<endl;
555 }
556
557 //---fitting
558 TrkHitList* qhits = trkRecoTrk->hits();
559 TrkErrCode err=qhits->fit();
560 return err;
561}
562
563TrkErrCode CgemMdcFitAlg::check(TrkRecoTrk* trkRecoTrk, TrackType trackType)
564{
565 m_qualityFactor=3;
566 m_dropTrkDrCut=1;
567 m_dropTrkDzCut=10;
568 if(trackType==circleType)m_dropTrkNhitCut=5;
569 if(trackType==helixType)m_dropTrkNhitCut=9;
570 m_dropTrkChi2Cut=99999;
571 m_dropTrkChi2NdfCut=99999;
572
573 TrkErrCode err;
574 int nActiveHit = trkRecoTrk->hots()->nActive();
575 TrkHitList* qhits = trkRecoTrk->hits();
576 const TrkFit* trkFit = trkRecoTrk->fitResult();
577 //if (!trkFit){cout<<"fit failed"<<endl;}
578 //else{cout<<"fit success"<<endl;}
579 //cout<<trkFit<<endl;
580 //TrkExchangePar par = trkFit->helix(0.);
581 int fit_stat=false;
582 double chi2=-999.;
583 if (!trkFit){
584 err.setFailure(12,"fit failed");
585 return err;
586 }else{
587 TrkExchangePar par = trkFit->helix(0.);
588 if( abs( 1/(par.omega()) )>0.03) {
589 fit_stat = 1;
590 chi2=trkFit->chisq();
591 }
592 else {
593 fit_stat = 0;
594 chi2=-999;
595 if(m_debug>0){
596 std::cout<<__FILE__<<" "<<__LINE__<<" drop track radius cut "<<abs(1/(par.omega()))<<" < "<<0.03 <<std::endl;
597 }
598 }
599
600 bool badQuality = false;
601 if(fabs(chi2)>m_qualityFactor*m_dropTrkChi2Cut ){
602 if(m_debug>0){
603 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2 "<<chi2<<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2Cut <<std::endl;
604 }
605 badQuality=1;
606 }
607 if( fabs(par.d0())>m_qualityFactor*m_dropTrkDrCut) {
608 if(m_debug>0){
609 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by d0 "<<par.d0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDrCut <<std::endl;
610 }
611 badQuality=1;
612 }
613 if( fabs(par.z0())>m_qualityFactor*m_dropTrkDzCut && trackType == helixType) {
614 if(m_debug>0){
615 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by z0 "<<par.z0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDzCut <<std::endl;
616 }
617 badQuality=1;
618 }
619 if( (fabs(chi2)/qhits->nHit()) > m_qualityFactor*m_dropTrkChi2NdfCut) {
620 if(m_debug>0){
621 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2/ndf "<<(fabs(chi2)/qhits->nHit()) <<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2NdfCut<<std::endl;
622 }
623 badQuality=1;
624 }
625 if( nActiveHit <= m_dropTrkNhitCut) {
626 if(m_debug>0){
627 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by nhit "<<nActiveHit <<" <= "<<m_dropTrkNhitCut<<std::endl;
628 }
629 badQuality=1;
630 }
631 //badQuality = false; //not delete track in this stage
632 if( badQuality) fit_stat=0;
633 }
634 //if( fit_stat!=1 ){
635 //delete trkRecoTrk;
636 //vectrk_for_clean.pop_back();
637 //}
638 if( fit_stat==1 ){
639
640 if(m_debug>1){
641 TrkExchangePar par = trkFit->helix(0.);
642 double _d0=par.d0();
643 double _phi0=par.phi0();
644 double _omega=par.omega();
645 double _tanl=0;
646 double _z0=0;
647 if(trackType == helixType){
648 _tanl=par.tanDip();
649 _z0=par.z0();
650 }
651 //double _tanl=par.tanDip();
652 //double _z0=par.z0();
653 int _charge=trkFit->charge();
654 cout << " global 3d fit success"<<endl;
655 //cout<<__FILE__<<__LINE__<<"AFTER fit 3d "<<endl;
656 trkRecoTrk->print(std::cout);
657 double _circleR=_charge/par.omega();
658 double _circleX= sin(par.phi0()) *(par.d0()+1/par.omega()) * -1.; //z axis verse,x_babar = -x_belle
659 double _circleY= -1.*cos(par.phi0()) *(par.d0()+1/par.omega()) * -1;//???
660 double pxy=fabs(_circleR/333.567);
661 double pz=pxy*_tanl;
662 double pxyz=_charge*sqrt(pxy*pxy+pz*pz);
663 double _pt=pxy;
664 double _pz=pz;
665 double _p=pxyz;
666 cout<<" circle after globle 3d: "<<"("<<_circleX<<","<<_circleY<<","<<_circleR<<")"<<endl;
667 cout<<"pt_rec: "<<_pt<<endl;
668 cout<<"pz_rec: "<<_pz<<endl;
669 cout<<"p_rec: "<<_p<<endl;
670
671 int nfit3d=0;
672 cout<<" hot list:"<<endl;
673 TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
674 int lay=((MdcHit*)(hotIter->hit()))->layernumber();
675 int wir=((MdcHit*)(hotIter->hit()))->wirenumber();
676 while (hotIter!=qhits->hotList().end()) {
677 if( typeid(*hotIter)==typeid(MdcRecoHitOnTrack) || typeid(*hotIter)==typeid(MdcHitOnTrack)){
678 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
679 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
680 <<":"<<hotIter->isActive()<<") ";
681 }else if(typeid(*hotIter)==typeid(CgemHitOnTrack)){
682 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hotIter));
683 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
684 cout<<"(" <<recCgemCluster->getlayerid()
685 <<","<<recCgemCluster->getsheetid()
686 <<":"<<hotIter->isActive()<<") ";
687 //<<"clusterID:"<<recCgemCluster->getclusterid()
688 //<<" flag:"<<recCgemCluster->getflag()<<endl;
689 //hotIter++;
690 //continue;
691 }
692 // cout << "nuse "<<hotIter->hit()->nUsedHits()<<endl;
693 if (hotIter->isActive()==1) nfit3d++;
694 hotIter++;
695 }
696 double _nfit=nfit3d;
697 cout<<" in 3D fit number of Active hits : "<<_nfit<<endl;
698 }
699 err.setSuccess(11,"good track");
700 return err;
701 }else {
702 err.setFailure(12,"drop by cut");
703 return err;
704 }
705 //_chi2_aver=chi2/_nfit;
706 // for(int i=0;i<t_hitCol.size();i++){
707 // delete t_hitCol[i];
708 // }
709 //return fit_stat;
710
711}
712
713StatusCode CgemMdcFitAlg::storeTracks(RecMdcTrackCol*& trackList_tds ,RecMdcHitCol*& hitList_tds, MdcTrackList& mdcTrackList){
714 MsgStream log(msgSvc(), name());
715 StatusCode sc;
716 //--- clear subtrees out of tds
717 DataObject *aTrackCol;
718 DataObject *aRecHitCol;
719 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
720 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
721 if(aTrackCol != NULL) {
722 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
723 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
724 }
725 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
726 if(aRecHitCol != NULL) {
727 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
728 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
729 }
730
731 //---register subtrees into tds
732 //RecMdcTrackCol* trackList_tds;
733 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
734 if (aTrackCol) {
735 trackList_tds = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
736 }else{
737 trackList_tds = new RecMdcTrackCol;
738 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList_tds);
739 if(!sc.isSuccess()) {
740 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
741 return StatusCode::FAILURE;
742 }
743 }
744 //RecMdcHitCol* hitList_tds;
745 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
746 if (aRecHitCol) {
747 hitList_tds = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
748 }else{
749 hitList_tds = new RecMdcHitCol;
750 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList_tds);
751 if(!sc.isSuccess()) {
752 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
753 return StatusCode::FAILURE;
754 }
755 }
756
757 //---store fitting result into tds
758 if(m_debug>0)cout<<"tracks after fitting: "<<mdcTrackList.length()<<endl;
759 for(int i =0;i<mdcTrackList.length();i++){
760 int trkStat = 5;//track fi by CgemMdcFitAlg set stat=5
761 mdcTrackList[i]->storeTrack(i, trackList_tds, hitList_tds, trkStat);
762 delete mdcTrackList[i];
763 }
764
765 return StatusCode::SUCCESS;
766}
767
768bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
769 return clusterA->getlayerid()<clusterB->getlayerid();
770}
771void CgemMdcFitAlg::updateTracks(int trackId, TrkRecoTrk* trkRecoTrk, RecMdcTrack* recMdcTrack, int tkStat){
772
773 //get the results of the fit to this track
774 const TrkFit* fitresult = trkRecoTrk->fitResult();
775 //check the fit worked
776 if (fitresult == 0) return;
777
778 //new a Rec. Track MdcTrack
779 //RecMdcTrack* recMdcTrack = new RecMdcTrack();
780 const TrkHitList* aList = trkRecoTrk->hits();
781 //some track info
782 const BField& theField = trkRecoTrk->bField();
783 double Bz = theField.bFieldZ();
784 //Fit info
785 int hitId = 0;
786 int nHits=0;
787 int nSt=0;
788 //int nAct=0; int nSeg=0;
789 //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
790 //****************************
791 //* active flag open this
792 //****************************
793 //* if (allHit>0){
794 //* nHits = aList->nHit();//add inActive hit also
795 //* }else{
796 //* nHits = fitresult->nMdc();// = nActive()
797 //* }
798 //****************************
799 //for 5.1.0 use all hits
800 nHits = aList->nHit();
801 //****************************
802 // use active only
803 // nHits = fitresult->nMdc();// = nActive()
804 //****************************
805
806 int q = fitresult->charge();
807 double chisq = fitresult->chisq();
808 int nDof = fitresult->nDof();
809 //track parameters at closest approach to beamline
810 double fltLenPoca = 0.0;
811 TrkExchangePar helix = fitresult->helix(fltLenPoca);
812 //std::cout<< __FILE__ << " phi0 " << helix.phi0() << " omega " <<helix.omega()<<std::endl;
813 double phi0 = helix.phi0();
814 double tanDip = helix.tanDip();
815 //double z0 = helix.z0();
816 double d0 = helix.d0();
817 //momenta and position
818 //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
819 HepPoint3D poca = fitresult->position(fltLenPoca);
820
821 double helixPar[5];
822 //dro =-d0
823 helixPar[0] = -d0;
824 //phi0 = phi0 - pi/2 [0,2pi)
825 double tphi0 = phi0 - Constants::pi/2.;
826 if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
827 helixPar[1] = tphi0;
828 //kappa = q/pxy
829 double pxy = fitresult->pt();
830 if (pxy == 0.) helixPar[2] = 9999.;
831 else helixPar[2] = q/fabs(pxy);
832 if(pxy>9999.) helixPar[2] = 0.00001;
833 //dz = z0
834 helixPar[3] = helix.z0();
835 //tanl
836 helixPar[4] = tanDip;
837 //error
838 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
839 HepSymMatrix mS(helix.params().num_row(),0);
840 mS[0][0]=-1.;//dr0=-d0
841 mS[1][1]=1.;
842 mS[2][2]=-333.567/Bz;//pxy -> cpa
843 mS[3][3]=1.;
844 mS[4][4]=1.;
845 HepSymMatrix mVy = helix.covariance().similarity(mS);
846 double errorMat[15];
847 int k = 0;
848 for (int ie = 0 ; ie < 5 ; ie ++){
849 for (int je = ie ; je < 5 ; je ++) {
850 errorMat[k] = mVy[ie][je];
851 k++;
852 }
853 }
854 double p,px,py,pz;
855 px = pxy * (-sin(helixPar[1]));
856 py = pxy * cos(helixPar[1]);
857 pz = pxy * helixPar[4];
858 p = sqrt(pxy*pxy + pz*pz);
859 //theta, phi
860 double theta = acos(pz/p);
861 double phi = atan2(py,px);
862 recMdcTrack->setTrackId(trackId);
863 recMdcTrack->setHelix(helixPar);
864 recMdcTrack->setCharge(q);
865 recMdcTrack->setPxy(pxy);
866 recMdcTrack->setPx(px);
867 recMdcTrack->setPy(py);
868 recMdcTrack->setPz(pz);
869 recMdcTrack->setP(p);
870 recMdcTrack->setTheta(theta);
871 recMdcTrack->setPhi(phi);
872 recMdcTrack->setPoca(poca);
873 recMdcTrack->setX(poca.x());//poca
874 recMdcTrack->setY(poca.y());
875 recMdcTrack->setZ(poca.z());
876 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
877 HepPoint3D pivot(0.,0.,0.);
878 recMdcTrack->setPivot(pivot);
879 recMdcTrack->setVX0(0.);//pivot
880 recMdcTrack->setVY0(0.);
881 recMdcTrack->setVZ0(0.);
882 recMdcTrack->setError(errorMat);
883 recMdcTrack->setError(mVy);
884 recMdcTrack->setChi2(chisq);
885 recMdcTrack->setNdof(nDof);
886 recMdcTrack->setStat(tkStat);
887 recMdcTrack->setNhits(nHits);
888 //-----------hit list-------------
889 HitRefVec hitRefVec;
890 ClusterRefVec clusterRefVec;
891 //vector<const RecCgemCluster*> clusterCol;
892 vector<int> vecHits;
893 //for fiterm
894 int maxLayerId = -1;
895 int minLayerId = 43;
896 double fiTerm = 999.;
897 const MdcRecoHitOnTrack* fiTermHot = NULL;
898 TrkHitList::hot_iterator hot = aList->begin();
899 for (;hot!=aList->end();hot++){
900 //cout<<__FILE__ <<" "<<__LINE__ <<" "<< typeid(*hot).name()<<endl;
901 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
902 const MdcRecoHitOnTrack* recoHot;
903 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
904 bool findhit = false;
905 RecMdcHit* recMdcHit;
906 Identifier mdcId = recoHot->mdcHit()->digi()->identify();
907 HitRefVec vechits = recMdcTrack->getVecHits();
908 HitRefVec::iterator itHit = vechits.begin();
909 for( ; itHit != vechits.end(); itHit++)
910 {
911 Identifier mdcid = (*itHit)->getMdcId();
912 if(mdcid == mdcId){
913 recMdcHit = (*itHit);
914 findhit = true;
915 break;
916 }
917 }
918 if(!findhit)continue;
919 //if (!recoHot->mdcHit()) return;
920 //id
921 recMdcHit->setId(hitId);
922 //---------BESIII----------
923 //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
924 //phiWire - phiHit >0 flagLR=1 right driftright>0;
925 //flag = 2 ambig;
926 //-----Babar wireAmbig----
927 //-1 -> right, 0 -> don't know, +1 -> left
928 //+1 phiWire-phiHit<0 Left,
929 //-1 phiWire-phiHit>0 Right,
930 //0 don't know
931 //ambig() w.r.t track direction
932 //wireAmbig() w.r.t. wire location
933 double hotWireAmbig = recoHot->wireAmbig();
934 double driftDist = fabs(recoHot->drift());
935 double sigma = recoHot->hitRms();
936 double doca = fabs(recoHot->dcaToWire());
937 int flagLR=2;
938 if ( hotWireAmbig == 1){
939 flagLR = 0; //left driftDist <0
940 doca *= -1.;//2012-07-19
941 }else if( hotWireAmbig == -1){
942 flagLR = 1; //right driftDist >0
943 }else if( hotWireAmbig == 0){
944 flagLR = 2; //don't know
945 }
946 recMdcHit->setFlagLR(flagLR);
947 recMdcHit->setDriftDistLeft((-1)*driftDist);
948 recMdcHit->setErrDriftDistLeft(sigma);
949 recMdcHit->setDriftDistRight(driftDist);
950 recMdcHit->setErrDriftDistRight(sigma);
951 //Mdc Id
952 recMdcHit->setMdcId(mdcId);
953 //corrected ADC
954
955 //contribution to chisquare
956 //fill chisq
957 double res=999.,rese=999.;
958 if (recoHot->resid(res,rese,false)){
959 }else{}
960 double deltaChi=0;
961 recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20
962 recMdcHit->setChisqAdd( deltaChi * deltaChi );
963 //set tracks containing this hit
964 recMdcHit->setTrkId(trackId);
965 //doca
966 recMdcHit->setDoca(doca);//doca sign left<0
967 //entrance angle
968
969 recMdcHit->setEntra(recoHot->entranceAngle());
970 //z of hit
971 HepPoint3D pos; Hep3Vector dir;
972 double fltLen = recoHot->fltLen();
973 recMdcHit->setFltLen(fltLen);
974 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
975 recMdcHit->setZhit(pos.z());
976 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
977 recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
978 recMdcHit->setAdc(recoHot->mdcHit()->adcIndex());
979 //drift time
980 recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
981 //for fiterm
982 int layerId = recoHot->mdcHit()->layernumber();
983 int wireId = recoHot->mdcHit()->wirenumber();
984 //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig "<<hotWireAmbig<<" y "<<pos.y()<<std::endl;
985 // <<recoHot->entranceAngle()<<std::endl;
986 if (layerId >= maxLayerId){
987 maxLayerId = layerId;
988 fiTermHot = recoHot;
989 }
990 if (layerId < minLayerId){
991 minLayerId = layerId;
992 }
993 // status flag
994 if (recoHot->isActive()) {
995 recMdcHit->setStat(1);
996 }else{
997 recMdcHit->setStat(0);
998 }
999 // for 5.1.0 use all hits
1000 if (recoHot->layer()->view()) {
1001 ++nSt;
1002 }
1003 //hitList->push_back(recMdcHit);
1004 SmartRef<RecMdcHit> refHit(recMdcHit);
1005 hitRefVec.push_back(refHit);
1006 vecHits.push_back(mdcId.get_value());
1007 ++hitId;
1008 }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
1009 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
1010 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
1011 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
1012 clusterRefVec.push_back(recCgemCluster);
1013 //clusterCol.push_back(recCgemCluster);
1014 }
1015 }
1016 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
1017 //std::sort(clusterCol.begin(),clusterCol.end(),sortCluster);
1018 //fi terminal angle
1019 if (fiTermHot!=NULL){
1020 fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
1021 }
1022 recMdcTrack->setFiTerm(fiTerm);
1023 // number of stereo hits contained
1024 recMdcTrack->setNster(nSt);
1025 recMdcTrack->setFirstLayer(maxLayerId);
1026 recMdcTrack->setLastLayer(minLayerId);
1027 recMdcTrack->setVecHits(hitRefVec);
1028 recMdcTrack->setVecClusters(clusterRefVec);
1029 recMdcTrack->setNcluster(clusterRefVec.size());
1030 //trackList->push_back(recMdcTrack);
1031
1032 //cout<<__FILE__ <<" "<<__LINE__ << endl;
1033 //for(int i=0;i<clusterRefVec.size();i++)
1034 //{
1035 //cout<<clusterRefVec[i]->getlayerid()<<endl;
1036 //}
1037 //hot = aList->begin();
1038 //for (;hot!=aList->end();hot++){
1039 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
1040 //const MdcRecoHitOnTrack* recoHot;
1041 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
1042 //int layerId = recoHot->mdcHit()->layernumber();
1043 //int wireId = recoHot->mdcHit()->wirenumber();
1044 //cout<<"hit ("<<layerId<<","<<wireId<<")"<<endl;
1045 //}else if(typeid(*hot)==typeid(CgemHitOnTrack)){
1046 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
1047 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
1048 //cout<<"hit ("<<recCgemCluster->getlayerid()<<","<<recCgemCluster->getclusterid()<<")"<<endl;
1049 //clusterRefVec.push_back(recCgemCluster);
1050 //}
1051 //}
1052}//end of storeTrack
1053
1054void CgemMdcFitAlg::getMcTruth(){
1055 MsgStream log(msgSvc(), name());
1056 StatusCode sc;
1057 t_t0Truth=-1;
1058 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
1059 if (!mcParticleCol) {
1060 log << MSG::ERROR<< "Could not find McParticle" << endreq;
1061 } else{
1062 int itrk= 0;
1063 t_nTrkMC=0;
1064 bool findJpsi=false;
1065 if(m_debug>1) cout<< "size_mcParticleCol "<<mcParticleCol->size()<< endl;
1066 McParticleCol::iterator iter_mc = mcParticleCol->begin();
1067 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
1068 if(!(*iter_mc)->primaryParticle()) continue;
1069 int pdg=(*iter_mc)->particleProperty();
1070 if(pdg==443) findJpsi=true;
1071 //if(!findJpsi) continue;
1072 //if(abs(pdg)!=211) continue;
1073 t_pidTruth[itrk] = (*iter_mc)->particleProperty();
1074
1075 if(m_debug>1) cout<< "tk "<<itrk<< " pid="<< t_pidTruth[itrk]<< endl;
1076 //if((m_pid!=0) && (abs(t_pidTruth[itrk]) != m_pid)){ continue; }
1077
1078 t_t0Truth=(*iter_mc)->initialPosition().t();
1079
1080 if(itrk>=10) break;
1081 //if(itrk>=10) continue;
1082
1083 int pid = t_pidTruth[itrk];
1084 if( pid == 0 ) {
1085 log << MSG::WARNING << "Wrong particle id" <<endreq;
1086 t_qMC[itrk]=0;
1087 }else{
1088 IPartPropSvc* p_PartPropSvc;
1089 static const bool CREATEIFNOTTHERE(true);
1090 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
1091 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
1092 cout<< " Could not initialize Particle Properties Service" << endl;
1093 }
1094 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
1095 std::string name;
1096 if( p_particleTable->particle(pid) ){
1097 name = p_particleTable->particle(pid)->name();
1098 t_qMC[itrk] = p_particleTable->particle(pid)->charge();
1099 }else if( p_particleTable->particle(-pid) ){
1100 name = "anti " + p_particleTable->particle(-pid)->name();
1101 t_qMC[itrk] = (-1.)*p_particleTable->particle(-pid)->charge();
1102 }
1103 }
1104
1105 t_vrMC[itrk] = sqrt((*iter_mc)->initialPosition().x()*(*iter_mc)->initialPosition().x() + (*iter_mc)->initialPosition().y()*(*iter_mc)->initialPosition().y()) ;
1106 t_vzMC[itrk] = (*iter_mc)->initialPosition().z();
1107 t_pxMC[itrk] = (*iter_mc)->initialFourMomentum().px();
1108 t_pyMC[itrk] = (*iter_mc)->initialFourMomentum().py();
1109 t_pzMC[itrk] = (*iter_mc)->initialFourMomentum().pz();
1110 t_ptMC[itrk] = sqrt(t_pxMC[itrk]*t_pxMC[itrk] + t_pyMC[itrk]*t_pyMC[itrk]);
1111 t_pMC[itrk] = sqrt(t_ptMC[itrk]*t_ptMC[itrk] + t_pzMC[itrk]*t_pzMC[itrk]);
1112 Hep3Vector p3(t_pxMC[itrk],t_pyMC[itrk],t_pzMC[itrk]);
1113
1114 t_costhetaMC[itrk] = t_pzMC[itrk]/t_pMC[itrk];
1115
1116 //double phi=9999.;
1117 t_phiMC[itrk] = p3.phi();
1118
1119 Helix mchelix = Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qMC[itrk]);//cm
1120 mchelix.pivot( HepPoint3D(0.,0.,0.) );
1121 int trkId = (*iter_mc)->trackIndex();
1122 double dr = mchelix.dr() ;
1123 double phi0 = mchelix.phi0() ;
1124 double kappa = mchelix.kappa();
1125 double dz = mchelix.dz() ;
1126 double tanl = mchelix.tanl() ;
1127 if(m_debug>0)cout<<"MCtruth helix parameters: dr, phi0, kappa, dz, tanl ---> "<<dr<<" "<<phi0<<" "<<kappa<<" "<<dz<<" "<<tanl<<endl;
1128 //cout<<"MC: "<<setw(12)<<dr<<" "<<setw(12)<<phi0<<" "<<setw(12)<<kappa<<" "<<setw(12)<<dz<<" "<<setw(12)<<tanl<<endl;
1129 m_trkidMC[itrk] = trkId;
1130 m_drMC[itrk] = dr ;
1131 m_phi0MC[itrk] = phi0 ;
1132 m_kappaMC[itrk] = kappa;
1133 m_dzMC[itrk] = dz ;
1134 m_tanlMC[itrk] = tanl ;
1135 m_qMC[itrk] = t_qMC[itrk];
1136 m_costhetaMC[itrk] = t_costhetaMC[itrk];
1137 m_phiMC[itrk] = t_phiMC[itrk];
1138 m_vzMC[itrk] = t_vzMC[itrk];
1139 m_vrMC[itrk] = t_vrMC[itrk];
1140 m_pxMC[itrk] = t_pxMC[itrk];
1141 m_pyMC[itrk] = t_pyMC[itrk];
1142 m_pzMC[itrk] = t_pzMC[itrk];
1143 m_ptMC[itrk] = t_ptMC[itrk];
1144 m_pMC[itrk] = t_pMC[itrk];
1145 itrk++;
1146 t_nTrkMC++;
1147 }
1148 m_nTrkMC = t_nTrkMC;
1149 }
1150
1151 // //mc hit
1152 // SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
1153 // if (!mdcMcHitCol) {
1154 // log << MSG::ERROR<< "Could not find MdcMcHitCol" << endreq;
1155 // } else{
1156 // int ihit= 0;
1157 // t_nMdcMcHit = mdcMcHitCol->size();
1158 // if(m_debug>1) cout<< "nMcHit="<<t_nMdcMcHit << endl;
1159 // Event::MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
1160 // for (;iter_mchit != mdcMcHitCol->end(); iter_mchit++,ihit++ ) {
1161 // const Identifier id= (*iter_mchit)->identify();
1162 // //int layer = MdcID::layer(id);
1163 // //int wire = MdcID::wire(id);
1164 // int iMcTk = (*iter_mchit)->getTrackIndex();
1165 // double mcX = (*iter_mchit)->getPositionX()/10.;
1166 // double mcY = (*iter_mchit)->getPositionY()/10.;
1167 // double mcZ = (*iter_mchit)->getPositionZ()/10.;
1168 // double mcLR = (*iter_mchit)->getPositionFlag();
1169 // double mcDrift = (*iter_mchit)->getDriftDistance()/10.;
1170 // double mcEnergy = (*iter_mchit)->getDepositEnergy()/10.;
1171 // if (mcLR == 0) mcLR = -1;
1172 // if(m_debug){
1173 // cout<<"Truth hit "<<ihit <<" tk "<< iMcTk
1174 // <<" xyz("<< mcX<<"," <<mcY<<","<< mcZ
1175 // <<") lr " <<mcLR
1176 // <<" drift " <<mcDrift
1177 // <<" energy " <<mcEnergy
1178 // <<endl;
1179 // }
1180 // }
1181 // }
1182
1183}
1184
1185void CgemMdcFitAlg::compareTracks(MdcTrackList& mdcTrackList, double dzCut){
1186 if(m_debug>1) cout<<"ntrack : "<<mdcTrackList.length()<<endl;
1187 for(int itrack=0; itrack<mdcTrackList.length(); itrack++){
1188 MdcTrack* track = mdcTrackList[itrack];
1189 TrkExchangePar par = track->track().fitResult()->helix(0.);
1190 int nhit = track->track().hots()->nHit();
1191 double pt = (1./par.omega())/333.567;
1192 if(m_debug>1) cout<<"i Track : "<<itrack<<" nHit: "<<nhit<<" pt: "<<pt<<endl;
1193 }
1194
1195 for(int itrack1=0; itrack1<mdcTrackList.length(); itrack1++){
1196 MdcTrack* track_1 = mdcTrackList[itrack1];
1197 TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
1198 int nhit1 = track_1->track().hots()->nHit();
1199 double d0_1 = par1.d0();
1200 double phi0_1 = par1.phi0();
1201 double omega_1 = par1.omega();
1202 double z0_1= par1.z0();
1203 double cr=fabs(1./par1.omega());
1204 double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
1205 double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
1206 if(m_debug>1) cout<<"circle 1 : "<<cr<<","<<cx<<","<<cy<<endl;
1207
1208 for(int itrack2=itrack1+1; itrack2<mdcTrackList.length(); itrack2++){
1209 MdcTrack* track_2 = mdcTrackList[itrack2];
1210 TrkExchangePar par2 = track_2->track().fitResult()->helix(0.);
1211 int nhit2 = track_2->track().hots()->nHit();
1212 double d0_2 = par2.d0();
1213 double phi0_2 = par2.phi0();
1214 double omega_2 = par2.omega();
1215 double z0_2= par2.z0();
1216 double cR=fabs(1./par2.omega());
1217 double cX= sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle
1218 double cY= -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//???
1219 if(m_debug>1) cout<<"circle 2 : "<<cR<<","<<cX<<","<<cY<<endl;
1220
1221 double bestDiff = 1.0e+20;
1222 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
1223 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
1224 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
1225 if(diff < bestDiff){
1226 bestDiff = diff;
1227 }
1228 }
1229 }
1230
1231 double zdiff = z0_1-z0_2;
1232 if(bestDiff != 1.0e20 && fabs(zdiff)<25.){
1233 if(m_debug>1) cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
1234 if( nhit1<nhit2 && (fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>dzCut) ){ //FIXME
1235 if(m_debug>0) cout<<"remove track1 "<<1./par1.omega()/333.567<<endl;
1236 //remove 1
1237 mdcTrackList.remove( track_1 );
1238 itrack1--;
1239 break;
1240 }
1241 else{
1242 //remove 2
1243 if(m_debug>0) cout<<"remove track2 "<<1./par2.omega()/333.567<<endl;
1244 //remove 1
1245 mdcTrackList.remove( track_2 );
1246 itrack2--;
1247 }
1248 }
1249 if(bestDiff == 1.0e20) { if(m_debug>1) cout<<" no same track in hough"<<endl; }
1250 }
1251 }
1252}
1253
1254StatusCode CgemMdcFitAlg::bookTuple(){
1255 MsgStream log(msgSvc(), name());
1256 NTuplePtr nt1(ntupleSvc(), "CgemMdcFitAlg/evt");
1257 if ( nt1 ){
1258 ntuple_evt = nt1;
1259 } else {
1260 ntuple_evt = ntupleSvc()->book("CgemMdcFitAlg/evt", CLID_ColumnWiseTuple, "evt");
1261 if(ntuple_evt){
1262 ntuple_evt->addItem("run" , m_run);
1263 ntuple_evt->addItem("event" , m_event);
1264 ntuple_evt->addItem("nTrkRec" , m_nTrkRec , 0,100);
1265 ntuple_evt->addItem("trkidRec", m_nTrkRec , m_trkidRec);
1266 ntuple_evt->addItem("drRec" , m_nTrkRec , m_drRec );
1267 ntuple_evt->addItem("phi0Rec" , m_nTrkRec , m_phi0Rec );
1268 ntuple_evt->addItem("kappaRec", m_nTrkRec , m_kappaRec);
1269 ntuple_evt->addItem("dzRec" , m_nTrkRec , m_dzRec );
1270 ntuple_evt->addItem("tanlRec" , m_nTrkRec , m_tanlRec );
1271 ntuple_evt->addItem("chargeRec" , m_nTrkRec , m_chargeRec );
1272 ntuple_evt->addItem("statRec" , m_nTrkRec , m_statRec );
1273 ntuple_evt->addItem("nhitsRec" , m_nTrkRec , m_nhitsRec );
1274 ntuple_evt->addItem("nclusterRec" , m_nTrkRec , m_nclusterRec );
1275 ntuple_evt->addItem("nsterRec" , m_nTrkRec , m_nsterRec );
1276 ntuple_evt->addItem("ndofRec" , m_nTrkRec , m_ndofRec );
1277 ntuple_evt->addItem("chi2Rec" , m_nTrkRec , m_chi2Rec );
1278 ntuple_evt->addItem("pxRec" , m_nTrkRec , m_pxRec );
1279 ntuple_evt->addItem("pyRec" , m_nTrkRec , m_pyRec );
1280 ntuple_evt->addItem("pzRec" , m_nTrkRec , m_pzRec );
1281 ntuple_evt->addItem("pxyRec" , m_nTrkRec , m_pxyRec );
1282 ntuple_evt->addItem("pRec" , m_nTrkRec , m_pRec );
1283 ntuple_evt->addItem("xRec" , m_nTrkRec , m_xRec );
1284 ntuple_evt->addItem("yRec" , m_nTrkRec , m_yRec );
1285 ntuple_evt->addItem("zRec" , m_nTrkRec , m_zRec );
1286 ntuple_evt->addItem("rRec" , m_nTrkRec , m_rRec );
1287 ntuple_evt->addItem("thetaRec" , m_nTrkRec , m_thetaRec );
1288 ntuple_evt->addItem("phiRec" , m_nTrkRec , m_phiRec );
1289 ntuple_evt->addItem("fiTermRec" , m_nTrkRec , m_fiTermRec );
1290
1291 ntuple_evt->addItem("nTrkFit" , m_nTrkFit , 0,100);
1292 ntuple_evt->addItem("trkidFit", m_nTrkFit , m_trkidFit);
1293 ntuple_evt->addItem("drFit" , m_nTrkFit , m_drFit );
1294 ntuple_evt->addItem("phi0Fit" , m_nTrkFit , m_phi0Fit );
1295 ntuple_evt->addItem("kappaFit", m_nTrkFit , m_kappaFit);
1296 ntuple_evt->addItem("dzFit" , m_nTrkFit , m_dzFit );
1297 ntuple_evt->addItem("tanlFit" , m_nTrkFit , m_tanlFit );
1298 ntuple_evt->addItem("chargeFit" , m_nTrkFit , m_chargeFit );
1299 ntuple_evt->addItem("statFit" , m_nTrkFit , m_statFit );
1300 ntuple_evt->addItem("nhitsFit" , m_nTrkFit , m_nhitsFit );
1301 ntuple_evt->addItem("nclusterFit" , m_nTrkFit , m_nclusterFit );
1302 ntuple_evt->addItem("nsterFit" , m_nTrkFit , m_nsterFit );
1303 ntuple_evt->addItem("ndofFit" , m_nTrkFit , m_ndofFit );
1304 ntuple_evt->addItem("chi2Fit" , m_nTrkFit , m_chi2Fit );
1305 ntuple_evt->addItem("pxFit" , m_nTrkFit , m_pxFit );
1306 ntuple_evt->addItem("pyFit" , m_nTrkFit , m_pyFit );
1307 ntuple_evt->addItem("pzFit" , m_nTrkFit , m_pzFit );
1308 ntuple_evt->addItem("pxyFit" , m_nTrkFit , m_pxyFit );
1309 ntuple_evt->addItem("pFit" , m_nTrkFit , m_pFit );
1310 ntuple_evt->addItem("xFit" , m_nTrkFit , m_xFit );
1311 ntuple_evt->addItem("yFit" , m_nTrkFit , m_yFit );
1312 ntuple_evt->addItem("zFit" , m_nTrkFit , m_zFit );
1313 ntuple_evt->addItem("rFit" , m_nTrkFit , m_rFit );
1314 ntuple_evt->addItem("thetaFit" , m_nTrkFit , m_thetaFit );
1315 ntuple_evt->addItem("phiFit" , m_nTrkFit , m_phiFit );
1316 ntuple_evt->addItem("fiTermFit" , m_nTrkFit , m_fiTermFit );
1317
1318 ntuple_evt->addItem("nTrkMC" , m_nTrkMC , 0,100);
1319 ntuple_evt->addItem("trkidMC", m_nTrkMC , m_trkidMC);
1320 ntuple_evt->addItem("drMC" , m_nTrkMC , m_drMC );
1321 ntuple_evt->addItem("phi0MC" , m_nTrkMC , m_phi0MC );
1322 ntuple_evt->addItem("kappaMC", m_nTrkMC , m_kappaMC);
1323 ntuple_evt->addItem("dzMC" , m_nTrkMC , m_dzMC );
1324 ntuple_evt->addItem("tanlMC" , m_nTrkMC , m_tanlMC );
1325 ntuple_evt->addItem("qMC" , m_nTrkMC , m_qMC );
1326 ntuple_evt->addItem("costhetaMC" , m_nTrkMC , m_costhetaMC );
1327 ntuple_evt->addItem("phiMC" , m_nTrkMC , m_phiMC );
1328 ntuple_evt->addItem("vzMC" , m_nTrkMC , m_vzMC );
1329 ntuple_evt->addItem("vrMC" , m_nTrkMC , m_vrMC );
1330 ntuple_evt->addItem("pxMC" , m_nTrkMC , m_pxMC );
1331 ntuple_evt->addItem("pyMC" , m_nTrkMC , m_pyMC );
1332 ntuple_evt->addItem("pzMC" , m_nTrkMC , m_pzMC );
1333 ntuple_evt->addItem("ptMC" , m_nTrkMC , m_ptMC );
1334 ntuple_evt->addItem("pMC" , m_nTrkMC , m_pMC );
1335 }
1336 }
1337 return StatusCode::SUCCESS;
1338}
1339
1340//cout<<__FILE__<<" "<<__LINE__<<endl;
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
sortCluster
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition CgemMdcFitAlg.cxx:768
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition RecMdcHit.h:99
ClusterRefVec
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition RecMdcTrack.h:28
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition RecMdcTrack.h:102
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition RecMdcTrack.h:26
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
M_PI
#define M_PI
Definition TConstant.h:4
CgemHitOnTrack::setCgemCalibFunSvc
void setCgemCalibFunSvc(const CgemCalibFunSvc *svc)
Definition CgemHitOnTrack.h:115
CgemHitOnTrack::baseHit
const RecCgemCluster * baseHit() const
Definition CgemHitOnTrack.h:119
CgemHitOnTrack::setCgemGeomSvc
void setCgemGeomSvc(const CgemGeomSvc *svc)
Definition CgemHitOnTrack.h:114
CgemMdcFitAlg::CgemMdcFitAlg
CgemMdcFitAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition CgemMdcFitAlg.cxx:60
CgemMdcFitAlg::execute
StatusCode execute()
Definition CgemMdcFitAlg.cxx:123
CgemMdcFitAlg::updateTracks
void updateTracks(int trackId, TrkRecoTrk *trkRecoTrk, RecMdcTrack *recMdcTrack, int trkStat)
Definition CgemMdcFitAlg.cxx:771
CgemMdcFitAlg::compareTracks
void compareTracks(MdcTrackList &mdcTrackList, double dzCut)
Definition CgemMdcFitAlg.cxx:1185
CgemMdcFitAlg::check
TrkErrCode check(TrkRecoTrk *trkRecoTrk, TrackType trackType)
Definition CgemMdcFitAlg.cxx:563
CgemMdcFitAlg::initialize
StatusCode initialize()
Definition CgemMdcFitAlg.cxx:76
CgemMdcFitAlg::finalize
StatusCode finalize()
Definition CgemMdcFitAlg.cxx:455
CgemMdcFitAlg::fit
TrkErrCode fit(RecMdcTrack *recMdcTrack, TrkRecoTrk *trkRecoTrk, TrackType trackType)
Definition CgemMdcFitAlg.cxx:463
Constants::pi
static const double pi
Definition Constants.h:38
DstMdcTrack::setFirstLayer
void setFirstLayer(const int id)
Definition DstMdcTrack.h:101
DstMdcTrack::setPxy
void setPxy(const double pxy)
Definition DstMdcTrack.h:86
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition DstMdcTrack.h:84
DstMdcTrack::setPy
void setPy(const double py)
Definition DstMdcTrack.h:88
DstMdcTrack::setZ
void setZ(const double z)
Definition DstMdcTrack.h:95
DstMdcTrack::setNster
void setNster(const int ns)
Definition DstMdcTrack.h:100
DstMdcTrack::setX
void setX(const double x)
Definition DstMdcTrack.h:93
DstMdcTrack::setError
void setError(double err[15])
Definition DstMdcTrack.cxx:118
DstMdcTrack::setNdof
void setNdof(const int ndof)
Definition DstMdcTrack.h:99
DstMdcTrack::setTheta
void setTheta(const double theta)
Definition DstMdcTrack.h:91
DstMdcTrack::setStat
void setStat(const int stat)
Definition DstMdcTrack.h:97
DstMdcTrack::setP
void setP(const double p)
Definition DstMdcTrack.h:90
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition DstMdcTrack.cxx:104
DstMdcTrack::setPoca
void setPoca(double poca[3])
Definition DstMdcTrack.cxx:111
DstMdcTrack::setR
void setR(const double r)
Definition DstMdcTrack.h:96
DstMdcTrack::setCharge
void setCharge(const int charge)
Definition DstMdcTrack.h:85
DstMdcTrack::setLastLayer
void setLastLayer(const int id)
Definition DstMdcTrack.h:102
DstMdcTrack::setY
void setY(const double y)
Definition DstMdcTrack.h:94
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition DstMdcTrack.h:98
DstMdcTrack::setPhi
void setPhi(const double phi)
Definition DstMdcTrack.h:92
DstMdcTrack::setPz
void setPz(const double pz)
Definition DstMdcTrack.h:89
DstMdcTrack::setPx
void setPx(const double px)
Definition DstMdcTrack.h:87
Helix::dr
double dr(void) const
returns an element of parameters.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:236
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:218
Identifier::get_value
value_type get_value() const
Definition Identifier.h:163
MdcHitOnTrack::wireAmbig
int wireAmbig() const
Definition MdcHitOnTrack.cxx:415
MdcHitOnTrack::dcaToWire
double dcaToWire() const
Definition MdcHitOnTrack.cxx:123
MdcHitOnTrack::drift
double drift() const
Definition MdcHitOnTrack.h:75
MdcHitOnTrack::entranceAngle
double entranceAngle() const
Definition MdcHitOnTrack.cxx:172
MdcHitOnTrack::layer
const MdcLayer * layer() const
Definition MdcHitOnTrack.cxx:354
MdcHit::digi
const MdcDigi * digi() const
Definition MdcHit.h:55
MdcHit::setCalibSvc
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition MdcHit.cxx:136
MdcHit::layernumber
unsigned layernumber() const
Definition MdcHit.h:61
MdcHit::wirenumber
unsigned wirenumber() const
Definition MdcHit.h:62
MdcHit::adcIndex
unsigned adcIndex() const
Definition MdcHit.h:64
MdcHit::driftTime
double driftTime(double tof, double z) const
Definition MdcHit.cxx:142
MdcHit::setCountPropTime
void setCountPropTime(const bool count)
Definition MdcHit.h:86
MdcHit::driftDist
double driftDist(double, int, double, double, double) const
Definition MdcHit.cxx:156
MdcHit::layer
const MdcLayer * layer() const
Definition MdcHit.h:56
MdcHit::tdcIndex
unsigned tdcIndex() const
Definition MdcHit.h:63
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:54
MdcLayer::view
int view(void) const
Definition MdcLayer.h:28
MdcRecoHitOnTrack::mdcHit
const MdcHit * mdcHit() const
Definition MdcRecoHitOnTrack.cxx:38
MdcTrackListBase::setNoInner
void setNoInner(bool noInnerFlag)
Definition MdcTrackListBase.h:63
MdcTrackListBase::remove
void remove(MdcTrack *atrack)
Definition MdcTrackListBase.cxx:617
MdcTrack::track
TrkRecoTrk & track()
Definition MdcTrack.h:27
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition RecCgemCluster.h:54
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition RecCgemCluster.h:56
RecCgemCluster::getflag
int getflag(void) const
Definition RecCgemCluster.h:58
RecCgemCluster::getsheetid
int getsheetid(void) const
Definition RecCgemCluster.h:57
RecMdcHit::setMdcId
void setMdcId(Identifier mdcid)
Definition RecMdcHit.h:67
RecMdcHit::setErrDriftDistRight
void setErrDriftDistRight(double erddr)
Definition RecMdcHit.h:63
RecMdcHit::setFltLen
void setFltLen(double fltLen)
Definition RecMdcHit.h:74
RecMdcHit::setErrDriftDistLeft
void setErrDriftDistLeft(double erddl)
Definition RecMdcHit.h:62
RecMdcHit::setDriftDistLeft
void setDriftDistLeft(double ddl)
Definition RecMdcHit.h:60
RecMdcHit::setDoca
void setDoca(double doca)
Definition RecMdcHit.h:71
RecMdcHit::setStat
void setStat(int stat)
Definition RecMdcHit.h:66
RecMdcHit::setTdc
void setTdc(double tdc)
Definition RecMdcHit.h:68
RecMdcHit::setAdc
void setAdc(double adc)
Definition RecMdcHit.h:69
RecMdcHit::setFlagLR
void setFlagLR(int lr)
Definition RecMdcHit.h:65
RecMdcHit::setChisqAdd
void setChisqAdd(double pChisq)
Definition RecMdcHit.h:64
RecMdcHit::setZhit
void setZhit(double zhit)
Definition RecMdcHit.h:73
RecMdcHit::setDriftT
void setDriftT(double driftT)
Definition RecMdcHit.h:70
RecMdcHit::setDriftDistRight
void setDriftDistRight(double ddr)
Definition RecMdcHit.h:61
RecMdcHit::setTrkId
void setTrkId(int trkid)
Definition RecMdcHit.h:59
RecMdcHit::setId
void setId(int id)
Definition RecMdcHit.h:58
RecMdcHit::setEntra
void setEntra(double entra)
Definition RecMdcHit.h:72
RecMdcTrack::getVecHits
const HitRefVec getVecHits(void) const
Definition RecMdcTrack.h:67
RecMdcTrack::setVecClusters
void setVecClusters(ClusterRefVec vecclusters)
Definition RecMdcTrack.cxx:86
RecMdcTrack::setPivot
void setPivot(const HepPoint3D &pivot)
Definition RecMdcTrack.h:79
RecMdcTrack::setVZ0
void setVZ0(double z0)
Definition RecMdcTrack.h:76
RecMdcTrack::setVY0
void setVY0(double y0)
Definition RecMdcTrack.h:75
RecMdcTrack::setVecHits
void setVecHits(HitRefVec vechits)
Definition RecMdcTrack.cxx:80
RecMdcTrack::setNhits
void setNhits(int nhits)
Definition RecMdcTrack.h:78
RecMdcTrack::setFiTerm
void setFiTerm(double fiterm)
Definition RecMdcTrack.h:77
RecMdcTrack::getVecClusters
const ClusterRefVec getVecClusters() const
Definition RecMdcTrack.h:70
RecMdcTrack::setNcluster
void setNcluster(int ncluster)
Definition RecMdcTrack.h:83
RecMdcTrack::setVX0
void setVX0(double x0)
Definition RecMdcTrack.h:74
Trajectory::getInfo
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &direction) const =0
TrkAbsFit::pt
virtual double pt(double fltL=0.) const =0
TrkAbsFit::chisq
virtual double chisq() const =0
TrkAbsFit::charge
virtual int charge() const =0
TrkAbsFit::nDof
virtual int nDof() const =0
TrkAbsFit::traj
virtual const TrkDifTraj & traj() const =0
TrkAbsFit::position
virtual HepPoint3D position(double fltL) const =0
TrkErrCode::failure
int failure() const
Definition TrkErrCode.h:61
TrkErrCode::setSuccess
void setSuccess(int i, const char *str=0)
Definition TrkErrCode.h:84
TrkErrCode::setFailure
void setFailure(int i, const char *str=0)
Definition TrkErrCode.h:79
TrkExchangePar::phi0
double phi0() const
Definition TrkExchangePar.h:43
TrkExchangePar::omega
double omega() const
Definition TrkExchangePar.h:44
TrkExchangePar::z0
double z0() const
Definition TrkExchangePar.h:45
TrkExchangePar::params
const HepVector & params() const
Definition TrkExchangePar.h:48
TrkExchangePar::d0
double d0() const
Definition TrkExchangePar.h:42
TrkExchangePar::tanDip
double tanDip() const
Definition TrkExchangePar.h:46
TrkExchangePar::covariance
const HepSymMatrix & covariance() const
Definition TrkExchangePar.h:50
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitList::fit
TrkErrCode fit()
Definition TrkHitList.cxx:59
TrkHitList::begin
hot_iterator begin() const
Definition TrkHitList.h:45
TrkHitList::nHit
unsigned nHit() const
Definition TrkHitList.h:44
TrkHitList::appendHot
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition TrkHitList.cxx:91
TrkHitList::hotList
const TrkHotList & hotList() const
Definition TrkHitList.cxx:51
TrkHitList::end
hot_iterator end() const
Definition TrkHitList.h:46
TrkHitOnTrk::fltLen
double fltLen() const
Definition TrkHitOnTrk.h:91
TrkHitOnTrk::resid
double resid(bool exclude=false) const
Definition TrkHitOnTrk.cxx:191
TrkHitOnTrk::hitRms
double hitRms() const
Definition TrkHitOnTrk.h:89
TrkHitOnTrk::setFltLen
void setFltLen(double f)
Definition TrkHitOnTrk.h:147
TrkHitOnTrk::getParentRep
const TrkRep * getParentRep() const
Definition TrkHitOnTrk.h:73
TrkHitOnTrk::isActive
bool isActive() const
Definition TrkHitOnTrk.h:200
TrkHitOnTrk::getFitStuff
TrkErrCode getFitStuff(HepVector &derivs, double &deltaChi) const
Definition TrkHitOnTrk.cxx:248
TrkHotList::end
hot_iterator end() const
Definition TrkHotList.h:45
TrkHotList::begin
hot_iterator begin() const
Definition TrkHotList.h:44
TrkHotList::nActive
virtual int nActive(TrkEnums::TrkViewInfo view=TrkEnums::bothView) const =0
TrkHotList::nHit
virtual int nHit(TrkEnums::TrkViewInfo view=TrkEnums::bothView) const =0
TrkRecoTrk::hits
TrkHitList * hits()
Definition TrkRecoTrk.h:107
TrkRecoTrk::bField
const BField & bField() const
Definition TrkRecoTrk.h:82
TrkRecoTrk::hots
TrkHotList * hots()
Definition TrkRecoTrk.h:113
TrkRecoTrk::print
virtual void print(std::ostream &) const
Definition TrkRecoTrk.cxx:166
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition TrkRecoTrk.cxx:387
TrkRep::arrivalTime
virtual double arrivalTime(double fltL) const
Definition TrkRep.cxx:192
TrkSimpleMaker::setFlipAndDrop
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkSimpleMaker::makeTrack
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const
EventModel::Recon::RecMdcTrackCol
_EXTERN_ std::string RecMdcTrackCol
Definition EventModel.h:92
EventModel::Recon::RecMdcHitCol
_EXTERN_ std::string RecMdcHitCol
Definition EventModel.h:91
Event
Definition Event.h:21