MarlinTrk  02.08
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DiagnosticsController.cc
Go to the documentation of this file.
1 
2 #include "MarlinTrk/MarlinTrkDiagnostics.h"
3 
4 #ifdef MARLINTRK_DIAGNOSTICS_ON
5 
6 #include "MarlinTrk/DiagnosticsController.h"
7 
8 #define MarlinTrkNtuple_cxx
9 #include "MarlinTrk/MarlinTrkNtuple.h"
10 
11 #include "MarlinTrk/MarlinKalTestTrack.h"
12 
13 #include "streamlog/streamlog.h"
14 
15 #include "kaldet/ILDVTrackHit.h"
16 #include "kaltest/TKalTrackSite.h"
17 #include "kaltest/TKalTrack.h"
18 
19 #include "TFile.h"
20 #include "TTree.h"
21 
22 //#include "MarlinTrk/MarlinTrkNtuple.h"
23 
24 #include "MarlinTrk/HelixTrack.h"
25 
26 #include "EVENT/MCParticle.h"
27 #include <UTIL/BitField64.h>
28 #include "UTIL/LCTrackerConf.h"
29 #include <UTIL/ILDConf.h>
30 
31 
32 
33 
34 namespace MarlinTrk{
35 
36 
37  DiagnosticsController::DiagnosticsController() {
38 
39  _initialised = false;
40  _recording_on = false; // recording is set to off by default so that processors do not have to perform any action if they are not interested in diagnostics, i.e. no need to call init
41  _current_track = 0;
42  _currentMCP = 0;
43  _mcpInfoStored=false;
44  _skip_track = false;
45  _ntracks_skipped = 0;
46  _ntracks_written = 0;
47  _track_record = 0;
48  _tree = 0;
49  _root_file = 0;
50 
51  } // constructor
52 
53  DiagnosticsController::~DiagnosticsController(){
54  if(_track_record) { delete _track_record; _track_record = 0; }
55  }
56 
57  void DiagnosticsController::init(std::string root_file_name, std::string root_tree_name, bool recording_on){
58 
59  streamlog_out(DEBUG4) << " DiagnosticsController::init called " << "\n"
60  << "\t root file name = " << root_file_name << "\n"
61  << "\t root tree name = " << root_tree_name << "\n"
62  << "\t recording on = " << recording_on << "\n"
63  << std::endl;
64 
65  _recording_on = recording_on;
66 
67  if ( _recording_on == false ) {
68  _initialised = true;
69  return;
70  }
71 
72  _root_file_name = root_file_name+".root";
73  _root_tree_name = root_tree_name;
74 
75  _root_file = new TFile(_root_file_name.c_str(),"RECREATE");
76 
77  _tree = new TTree( _root_tree_name.c_str(), _root_tree_name.c_str());
78 
79  _track_record = new MarlinTrkNtuple();
80 
81  _track_record->CreateBranches(_tree);
82 
83  _initialised = true;
84 
85  }
86 
87 
88  void DiagnosticsController::skip_current_track(){
89 
90  streamlog_out(DEBUG) << " DiagnosticsController::skip_current_track called " << std::endl;
91 
92  if ( _recording_on == false ) {
93  return;
94  }
95 
96  _skip_track = true ;
97  this->clear_track_record();
98  }
99 
100  void DiagnosticsController::clear_track_record(){
101 
102  streamlog_out(DEBUG) << " DiagnosticsController::clear_track_record called " << std::endl;
103 
104  if ( _recording_on == false ) {
105  return;
106  }
107 
108  _track_record->error_code = 0 ;
109 
110  _track_record->nsites = 0 ;
111 
112  _track_record->nsites_vxd = 0 ;
113  _track_record->nsites_sit = 0 ;
114  _track_record->nsites_ftd = 0 ;
115  _track_record->nsites_tpc = 0 ;
116  _track_record->nsites_set = 0 ;
117 
118  _track_record->x_mcp = 0 ;
119  _track_record->y_mcp = 0 ;
120  _track_record->z_mcp = 0 ;
121 
122  _track_record->px_mcp = 0 ;
123  _track_record->py_mcp = 0 ;
124  _track_record->pz_mcp = 0 ;
125  _track_record->p_mcp = 0 ;
126  _track_record->theta_mcp = 0 ;
127  _track_record->phi_mcp = 0 ;
128  _track_record->pdg_mcp = 0 ;
129 
130  _track_record->d0_mcp = 0 ;
131  _track_record->phi0_mcp = 0 ;
132  _track_record->omega_mcp = 0 ;
133  _track_record->z0_mcp = 0 ;
134  _track_record->tanL_mcp = 0 ;
135 
136  _track_record->ndf = 0 ;
137  _track_record->chi2 = 0 ;
138  _track_record->prob = 0 ;
139 
140  _track_record->d0_seed = 0 ;
141  _track_record->phi0_seed = 0 ;
142  _track_record->omega_seed = 0 ;
143  _track_record->z0_seed = 0 ;
144  _track_record->tanL_seed = 0 ;
145 
146  _track_record->seed_ref_point_x = 0 ;
147  _track_record->seed_ref_point_y = 0 ;
148  _track_record->seed_ref_point_z = 0 ;
149 
150  _track_record->cov_seed_d0d0 = 0 ;
151  _track_record->cov_seed_phi0d0 = 0 ;
152  _track_record->cov_seed_phi0phi0 = 0 ;
153  _track_record->cov_seed_kappad0 = 0 ;
154  _track_record->cov_seed_kappaphi0 = 0 ;
155  _track_record->cov_seed_kappakappa = 0 ;
156  _track_record->cov_seed_z0d0 = 0 ;
157  _track_record->cov_seed_z0phi0 = 0 ;
158  _track_record->cov_seed_z0kappa = 0 ;
159  _track_record->cov_seed_z0z0 = 0 ;
160  _track_record->cov_seed_tanLd0 = 0 ;
161  _track_record->cov_seed_tanLphi0 = 0 ;
162  _track_record->cov_seed_tanLkappa = 0 ;
163  _track_record->cov_seed_tanLz0 = 0 ;
164  _track_record->cov_seed_tanLtanL = 0 ;
165 
166 
167  _track_record->d0_ip = 0 ;
168  _track_record->phi0_ip = 0 ;
169  _track_record->omega_ip = 0 ;
170  _track_record->z0_ip = 0 ;
171  _track_record->tanL_ip = 0 ;
172 
173  _track_record->cov_ip_d0d0 = 0 ;
174  _track_record->cov_ip_phi0d0 = 0 ;
175  _track_record->cov_ip_phi0phi0 = 0 ;
176  _track_record->cov_ip_omegad0 = 0 ;
177  _track_record->cov_ip_omegaphi0 = 0 ;
178  _track_record->cov_ip_omegaomega = 0 ;
179  _track_record->cov_ip_z0d0 = 0 ;
180  _track_record->cov_ip_z0phi0 = 0 ;
181  _track_record->cov_ip_z0omega = 0 ;
182  _track_record->cov_ip_z0z0 = 0 ;
183  _track_record->cov_ip_tanLd0 = 0 ;
184  _track_record->cov_ip_tanLphi0 = 0 ;
185  _track_record->cov_ip_tanLomega = 0 ;
186  _track_record->cov_ip_tanLz0 = 0 ;
187  _track_record->cov_ip_tanLtanL = 0 ;
188 
189  for ( int i = 0 ; i<MAX_SITES; ++i) {
190 
191  _track_record->CellID0[i] = 0 ;
192  _track_record->rejected[i] = 0 ;
193 
194  _track_record->site_x[i] = 0 ;
195  _track_record->site_y[i] = 0 ;
196  _track_record->site_z[i] = 0 ;
197 
198  _track_record->ref_point_x[i] = 0 ;
199  _track_record->ref_point_y[i] = 0 ;
200  _track_record->ref_point_z[i] = 0 ;
201 
202  _track_record->d0_mc[i] = 0 ;
203  _track_record->phi0_mc[i] = 0 ;
204  _track_record->omega_mc[i] = 0 ;
205  _track_record->z0_mc[i] = 0 ;
206  _track_record->tanL_mc[i] = 0 ;
207 
208  _track_record->d0_predicted[i] = 0 ;
209  _track_record->phi0_predicted[i] = 0 ;
210  _track_record->omega_predicted[i] = 0 ;
211  _track_record->z0_predicted[i] = 0 ;
212  _track_record->tanL_predicted[i] = 0 ;
213 
214  _track_record->d0_filtered[i] = 0 ;
215  _track_record->phi0_filtered[i] = 0 ;
216  _track_record->omega_filtered[i] = 0 ;
217  _track_record->z0_filtered[i] = 0 ;
218  _track_record->tanL_filtered[i] = 0 ;
219 
220  _track_record->d0_smoothed[i] = 0 ;
221  _track_record->phi0_smoothed[i] = 0 ;
222  _track_record->omega_smoothed[i] = 0 ;
223  _track_record->z0_smoothed[i] = 0 ;
224  _track_record->tanL_smoothed[i] = 0 ;
225 
226 
227  _track_record->chi2_inc_filtered[i] = 0 ;
228  _track_record->chi2_inc_smoothed[i] = 0 ;
229  _track_record->dim[i] = 0 ;
230 
231  _track_record->cov_predicted_d0d0[i] = 0 ;
232  _track_record->cov_predicted_phi0d0[i] = 0 ;
233  _track_record->cov_predicted_phi0phi0[i] = 0 ;
234  _track_record->cov_predicted_omegad0[i] = 0 ;
235  _track_record->cov_predicted_omegaphi0[i] = 0 ;
236  _track_record->cov_predicted_omegaomega[i] = 0 ;
237  _track_record->cov_predicted_z0d0[i] = 0 ;
238  _track_record->cov_predicted_z0phi0[i] = 0 ;
239  _track_record->cov_predicted_z0omega[i] = 0 ;
240  _track_record->cov_predicted_z0z0[i] = 0 ;
241  _track_record->cov_predicted_tanLd0[i] = 0 ;
242  _track_record->cov_predicted_tanLphi0[i] = 0 ;
243  _track_record->cov_predicted_tanLomega[i] = 0 ;
244  _track_record->cov_predicted_tanLz0[i] = 0 ;
245  _track_record->cov_predicted_tanLtanL[i] = 0 ;
246 
247  _track_record->cov_filtered_d0d0[i] = 0 ;
248  _track_record->cov_filtered_phi0d0[i] = 0 ;
249  _track_record->cov_filtered_phi0phi0[i] = 0 ;
250  _track_record->cov_filtered_omegad0[i] = 0 ;
251  _track_record->cov_filtered_omegaphi0[i] = 0 ;
252  _track_record->cov_filtered_omegaomega[i] = 0 ;
253  _track_record->cov_filtered_z0d0[i] = 0 ;
254  _track_record->cov_filtered_z0phi0[i] = 0 ;
255  _track_record->cov_filtered_z0omega[i] = 0 ;
256  _track_record->cov_filtered_z0z0[i] = 0 ;
257  _track_record->cov_filtered_tanLd0[i] = 0 ;
258  _track_record->cov_filtered_tanLphi0[i] = 0 ;
259  _track_record->cov_filtered_tanLomega[i] = 0 ;
260  _track_record->cov_filtered_tanLz0[i] = 0 ;
261  _track_record->cov_filtered_tanLtanL[i] = 0 ;
262 
263  _track_record->cov_smoothed_d0d0[i] = 0 ;
264  _track_record->cov_smoothed_phi0d0[i] = 0 ;
265  _track_record->cov_smoothed_phi0phi0[i] = 0 ;
266  _track_record->cov_smoothed_omegad0[i] = 0 ;
267  _track_record->cov_smoothed_omegaphi0[i] = 0 ;
268  _track_record->cov_smoothed_omegaomega[i] = 0 ;
269  _track_record->cov_smoothed_z0d0[i] = 0 ;
270  _track_record->cov_smoothed_z0phi0[i] = 0 ;
271  _track_record->cov_smoothed_z0omega[i] = 0 ;
272  _track_record->cov_smoothed_z0z0[i] = 0 ;
273  _track_record->cov_smoothed_tanLd0[i] = 0 ;
274  _track_record->cov_smoothed_tanLphi0[i] = 0 ;
275  _track_record->cov_smoothed_tanLomega[i] = 0 ;
276  _track_record->cov_smoothed_tanLz0[i] = 0 ;
277  _track_record->cov_smoothed_tanLtanL[i] = 0 ;
278 
279 
280 
281  }
282 
283  }
284 
285  void DiagnosticsController::new_track(MarlinKalTestTrack* trk){
286 
287  if ( _recording_on == false ) {
288  return;
289  }
290 
291  if ( _initialised == false ){
292  streamlog_out(ERROR) << "DiagnosticsController::new_track: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
293  << __FILE__
294  << " line "
295  << __LINE__
296  << std::endl;
297 
298  exit(1);
299  }
300 
301  streamlog_out(DEBUG3) << " DiagnosticsController::new_track called " << std::endl;
302 
303  this->clear_track_record();
304 
305  _current_track = trk;
306 
307  _skip_track = false;
308  _currentMCP = 0;
309  _mcpInfoStored=false;
310 
311  }
312 
313 
314 
315  void DiagnosticsController::set_intial_track_parameters(double d0, double phi0, double omega, double z0, double tanL, double pivot_x, double pivot_y, double pivot_z, TKalMatrix& cov){
316 
317  if ( _recording_on == false ) {
318  return;
319  }
320 
321  if ( _initialised == false ){
322  streamlog_out(ERROR) << "DiagnosticsController::set_intial_track_parameters: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
323  << __FILE__
324  << " line "
325  << __LINE__
326  << std::endl;
327 
328  exit(1);
329  }
330 
331  streamlog_out(DEBUG3) << " DiagnosticsController::set_intial_track_parameters called " << std::endl;
332 
333  _track_record->d0_seed= d0;
334  _track_record->phi0_seed= phi0;
335  _track_record->omega_seed= omega;
336  _track_record->z0_seed= z0;
337  _track_record->tanL_seed= tanL;
338 
339  _track_record->seed_ref_point_x = pivot_x ;
340  _track_record->seed_ref_point_y = pivot_y ;
341  _track_record->seed_ref_point_z = pivot_z ;
342 
343 
344  _track_record->cov_seed_d0d0 = cov( 0 , 0 ) ;
345  _track_record->cov_seed_phi0d0 = cov( 1 , 0 ) ;
346  _track_record->cov_seed_phi0phi0 = cov( 1 , 1 ) ;
347  _track_record->cov_seed_kappad0 = cov( 2 , 0 ) ;
348  _track_record->cov_seed_kappaphi0 = cov( 2 , 1 ) ;
349  _track_record->cov_seed_kappakappa = cov( 2 , 2 ) ;
350  _track_record->cov_seed_z0d0 = cov( 3 , 0 ) ;
351  _track_record->cov_seed_z0phi0 = cov( 3 , 1 ) ;
352  _track_record->cov_seed_z0kappa = cov( 3 , 2 ) ;
353  _track_record->cov_seed_z0z0 = cov( 3 , 3 ) ;
354  _track_record->cov_seed_tanLd0 = cov( 4 , 0 ) ;
355  _track_record->cov_seed_tanLphi0 = cov( 4 , 1 ) ;
356  _track_record->cov_seed_tanLkappa = cov( 4 , 2 ) ;
357  _track_record->cov_seed_tanLz0 = cov( 4 , 3 ) ;
358  _track_record->cov_seed_tanLtanL = cov( 4 , 4 ) ;
359 
360  // cov.Print();
361 
362 
363  streamlog_out(DEBUG3) << " $#$#$# Initial Track Parameters: "
364  << "d0 = " << d0 << " " << "[+/-" << sqrt( _track_record->cov_seed_d0d0 ) << "] "
365  << "phi0 = " << phi0 << " " << "[+/-" << sqrt( _track_record->cov_seed_phi0phi0 ) << "] "
366  << "omega = " << omega << " " << "[+/-" << sqrt( _track_record->cov_seed_kappakappa ) << "] "
367  << "z0 = " << z0 << " " << "[+/-" << sqrt( _track_record->cov_seed_z0z0 ) << "] "
368  << "tanL = " << tanL << " " << "[+/-" << sqrt( _track_record->cov_seed_tanLtanL ) << "] "
369  << "ref = " << pivot_x << " " << pivot_y << " " << pivot_z << " "
370  << std::endl;
371 
372  }
373 
374 
375  void DiagnosticsController::record_rejected_site(ILDVTrackHit* hit, TKalTrackSite* site){
376 
377  if ( _recording_on == false ) {
378  return;
379  }
380 
381  if ( _initialised == false ){
382  streamlog_out(ERROR) << "DiagnosticsController::record_rejected_site: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
383  << __FILE__
384  << " line "
385  << __LINE__
386  << std::endl;
387 
388  exit(1);
389  }
390 
391 
392  if(_skip_track) return;
393 
394  _track_record->rejected[_track_record->nsites] = 1;
395  _track_record->CellID0[_track_record->nsites] = hit->getLCIOTrackerHit()->getCellID0();
396 
397 
398  EVENT::MCParticle* mcp = hit->getLCIOTrackerHit()->ext<MarlinTrk::MCTruth4HitExt>()->simhit->getMCParticle();
399 
400  ++_track_record->nsites;
401  streamlog_out(DEBUG2) << "record_rejected_site _track_record->nsites = " << _track_record->nsites << " MCParticle of simhit = " << mcp << " pdg " << mcp->getPDG() << " energy = " << mcp->getEnergy() << " id = " << mcp->id() << std::endl;
402 
403  }
404 
405 
406  void DiagnosticsController::record_site(ILDVTrackHit* hit, TKalTrackSite* site){
407 
408  if ( _recording_on == false ) {
409  return;
410  }
411 
412  if ( _initialised == false ){
413  streamlog_out(ERROR) << "DiagnosticsController::record_site: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
414  << __FILE__
415  << " line "
416  << __LINE__
417  << std::endl;
418 
419  exit(1);
420  }
421 
422  streamlog_out(DEBUG2) << "DiagnosticsController::record_site called " << std::endl;
423 
424  if(_skip_track) return;
425 
426  EVENT::TrackerHit* trkhit = hit->getLCIOTrackerHit();
427 
428  MarlinTrk::MCTruth4HitExtStruct* ext = trkhit->ext<MarlinTrk::MCTruth4HitExt>();
429 
430  if ( ! ext ) {
431 
432  streamlog_out(ERROR) << "MCTruth4HitExt not attached to TrackerHit use processor SetTrackerHitExtensions to set them: exit(1) called from file "
433  << __FILE__
434  << " line "
435  << __LINE__
436  << std::endl;
437 
438  exit(1);
439 
440  }
441 
442  EVENT::SimTrackerHit* simhit = trkhit->ext<MarlinTrk::MCTruth4HitExt>()->simhit;
443 
444  if ( ! simhit ) {
445 
446  streamlog_out(ERROR) << "SimTrackerHit not attached to TrackerHit using MCTruth4HitExt use processor SetTrackerHitExtensions to set them: exit(1) called from file "
447  << __FILE__
448  << " line "
449  << __LINE__
450  << std::endl;
451 
452  exit(1);
453 
454  }
455 
456  EVENT::MCParticle* mcp = simhit->getMCParticle() ;
457 
458 
459  if( _track_record->nsites >= 0 ){
460 
461  if ( _mcpInfoStored == false ) {
462 
463  streamlog_out(DEBUG2) << "DiagnosticsController::record_site store MCParticle parameters for " << mcp << " pdg " << mcp->getPDG() << " energy = " << mcp->getEnergy() << " id = " << mcp->id() << std::endl;
464 
465  _currentMCP = mcp;
466 
467  _track_record->x_mcp = mcp->getVertex()[0];
468  _track_record->y_mcp = mcp->getVertex()[1];
469  _track_record->z_mcp = mcp->getVertex()[2];
470 
471  _track_record->px_mcp = mcp->getMomentum()[0];
472  _track_record->py_mcp = mcp->getMomentum()[1];
473  _track_record->pz_mcp = mcp->getMomentum()[2];
474 
475  double pt = sqrt(_track_record->px_mcp*_track_record->px_mcp +
476  _track_record->py_mcp*_track_record->py_mcp ) ;
477 
478  _track_record->p_mcp = sqrt( pt*pt + _track_record->pz_mcp*_track_record->pz_mcp );
479 
480 
481  _track_record->theta_mcp = atan2( pt, _track_record->pz_mcp );
482  _track_record->phi_mcp = atan2( _track_record->py_mcp, _track_record->px_mcp );
483 
484  _track_record->pdg_mcp = mcp->getPDG();
485 
486  // HelixTrack helixMC( sim_pos, sim_mom, mcp->getCharge(), ml.GetBz() ) ;
487  HelixTrack helixMC( mcp->getVertex(), mcp->getMomentum(), mcp->getCharge(), site->GetBfield() ) ;
488 
489  helixMC.moveRefPoint(0.0, 0.0, 0.0);
490 
491  _track_record->d0_mcp= helixMC.getD0();
492  _track_record->phi0_mcp = helixMC.getPhi0();
493  _track_record->omega_mcp = helixMC.getOmega();
494  _track_record->z0_mcp = helixMC.getZ0();
495  _track_record->tanL_mcp = helixMC.getTanLambda();
496 
497  _mcpInfoStored = true;
498 
499  }
500  else if( _currentMCP != mcp ) {
501  _skip_track = true ; // do not store tracks formed from more than one MCParticle
502  streamlog_out(WARNING) << "DiagnosticsController::record_site: Track skipped. Not storing tracks formed from more than one MCParticle " << std::endl;
503  return ;
504  }
505 
506  double sim_pos[3];
507  sim_pos[0] = simhit->getPosition()[0];
508  sim_pos[1] = simhit->getPosition()[1];
509  sim_pos[2] = simhit->getPosition()[2];
510 
511  double sim_mom[3];
512  sim_mom[0] = simhit->getMomentum()[0];
513  sim_mom[1] = simhit->getMomentum()[1];
514  sim_mom[2] = simhit->getMomentum()[2];
515 
516  if ( fabs(sim_mom[0]) < 1.e-09 && fabs(sim_mom[1]) < 1.e-09 && fabs(sim_mom[2]) < 1.e-09 ) {
517  // then the momentum is sub eV and therefore the momentum was certainly not set
518  streamlog_out(ERROR) << "Momentum not set in SimTrackerHit exit(1) called from file "
519  << __FILE__
520  << " line "
521  << __LINE__
522  << std::endl;
523 
524  exit(1);
525 
526  }
527 
528 
529  _track_record->CellID0[_track_record->nsites] = trkhit->getCellID0() ;
530 
531  UTIL::BitField64 encoder(lcio::LCTrackerCellID::encoding_string());
532  encoder.setValue( trkhit->getCellID0() );
533 
534  if (encoder[lcio::LCTrackerCellID::subdet()] == lcio::ILDDetID::VXD) {
535  ++_track_record->nsites_vxd;
536  }
537  else if (encoder[lcio::LCTrackerCellID::subdet()] == lcio::ILDDetID::SIT) {
538  ++_track_record->nsites_sit;
539  }
540  else if (encoder[lcio::LCTrackerCellID::subdet()] == lcio::ILDDetID::FTD) {
541  ++_track_record->nsites_ftd;
542  }
543  else if (encoder[lcio::LCTrackerCellID::subdet()] == lcio::ILDDetID::TPC) {
544  ++_track_record->nsites_tpc;
545  }
546  else if (encoder[lcio::LCTrackerCellID::subdet()] == lcio::ILDDetID::SET) {
547  ++_track_record->nsites_set;
548  }
549 
550  _track_record->chi2_inc_filtered[_track_record->nsites] = site->GetDeltaChi2();
551  _track_record->dim[_track_record->nsites] = site->GetDimension();
552 
553  // create helix from MCTruth at sim hit
554  // HelixTrack helixMC( sim_pos, sim_mom, mcp->getCharge(), ml.GetBz() ) ;
555  HelixTrack helixMC( sim_pos, sim_mom, mcp->getCharge(), site->GetBfield() ) ;
556 
557  // here perhaps we should move the helix to the hit to calculate the deltas or though this could still be done in the analysis code, as both sim and rec hit positions are stored ?
558  // helixMC.moveRefPoint(0.0, 0.0, 0.0);
559 
560  streamlog_out(DEBUG2) << " $#$#$# SimHit Track Parameters: "
561  << "x = " << sim_pos[0] << " "
562  << "y = " << sim_pos[1] << " "
563  << "z = " << sim_pos[2] << " "
564  << "px = " << sim_mom[0] << " "
565  << "py = " << sim_mom[1] << " "
566  << "pz = " << sim_mom[2] << " "
567  << "p = " << sqrt(sim_mom[0]*sim_mom[0] + sim_mom[1]*sim_mom[1] + sim_mom[2]*sim_mom[2]) << " "
568  << "d0 = " << helixMC.getD0() << " "
569  << "phi0 = " << helixMC.getPhi0() << " "
570  << "omega = " << helixMC.getOmega() << " "
571  << "z0 = " << helixMC.getZ0() << " "
572  << "tanL = " << helixMC.getTanLambda() << " "
573  << " mcp ID = " << simhit->getMCParticle()->id()
574  << std::endl;
575 
576  _track_record->d0_mc[_track_record->nsites] = helixMC.getD0();
577  _track_record->phi0_mc[_track_record->nsites] = helixMC.getPhi0();
578  _track_record->omega_mc[_track_record->nsites] = helixMC.getOmega();
579  _track_record->z0_mc[_track_record->nsites] = helixMC.getZ0();
580  _track_record->tanL_mc[_track_record->nsites] = helixMC.getTanLambda();
581 
582  double rec_x = trkhit->getPosition()[0];
583  double rec_y = trkhit->getPosition()[1];
584  double rec_z = trkhit->getPosition()[2];
585 
586  _track_record->site_x[_track_record->nsites] = rec_x;
587  _track_record->site_y[_track_record->nsites] = rec_y;
588  _track_record->site_z[_track_record->nsites] = rec_z;
589 
590 
591 // // move track state to the sim hit for comparison
592  const TVector3 tpoint( sim_pos[0], sim_pos[1], sim_pos[2] ) ;
593 
594  // move track state to the IP to make all comparisons straight forward
595  // const TVector3 tpoint( 0.0, 0.0, 0.0 ) ;
596 
597 
598  IMPL::TrackStateImpl ts;
599  int ndf;
600  double chi2;
601  double dPhi ;
602 
603 
604  TKalTrackState& trkState_predicted = (TKalTrackState&) site->GetState(TVKalSite::kPredicted); // this segfaults if no hits are present
605 
606  THelicalTrack helix_predicted = trkState_predicted.GetHelix() ;
607  TMatrixD c0_predicted(trkState_predicted.GetCovMat());
608 
609  Int_t sdim = trkState_predicted.GetDimension(); // dimensions of the track state, it will be 5 or 6
610 
611  // now move to the point
612  TKalMatrix DF(sdim,sdim);
613  DF.UnitMatrix();
614  helix_predicted.MoveTo( tpoint , dPhi , &DF , &c0_predicted) ; // move helix to desired point, and get propagator matrix
615 
616  streamlog_out(DEBUG2) << "DiagnosticsController::record_site get PREDICTED trackstate " << std::endl;
617 
618  _current_track->ToLCIOTrackState( helix_predicted, c0_predicted, ts, chi2, ndf );
619 
620  _track_record->d0_predicted[_track_record->nsites] = ts.getD0() ;
621  _track_record->phi0_predicted[_track_record->nsites] = ts.getPhi() ;
622  _track_record->omega_predicted[_track_record->nsites] = ts.getOmega() ;
623  _track_record->z0_predicted[_track_record->nsites] = ts.getZ0() ;
624  _track_record->tanL_predicted[_track_record->nsites] = ts.getTanLambda() ;
625 
626 
627  _track_record->cov_predicted_d0d0[_track_record->nsites] = ts.getCovMatrix()[0];
628  _track_record->cov_predicted_phi0d0[_track_record->nsites] = ts.getCovMatrix()[1];
629  _track_record->cov_predicted_phi0phi0[_track_record->nsites] = ts.getCovMatrix()[2];
630  _track_record->cov_predicted_omegad0[_track_record->nsites] = ts.getCovMatrix()[3];
631  _track_record->cov_predicted_omegaphi0[_track_record->nsites] = ts.getCovMatrix()[4];
632  _track_record->cov_predicted_omegaomega[_track_record->nsites] = ts.getCovMatrix()[5];
633  _track_record->cov_predicted_z0d0[_track_record->nsites] = ts.getCovMatrix()[6];
634  _track_record->cov_predicted_z0phi0[_track_record->nsites] = ts.getCovMatrix()[7];
635  _track_record->cov_predicted_z0omega[_track_record->nsites] = ts.getCovMatrix()[8];
636  _track_record->cov_predicted_z0z0[_track_record->nsites] = ts.getCovMatrix()[9];
637  _track_record->cov_predicted_tanLd0[_track_record->nsites] = ts.getCovMatrix()[10];
638  _track_record->cov_predicted_tanLphi0[_track_record->nsites] = ts.getCovMatrix()[11];
639  _track_record->cov_predicted_tanLomega[_track_record->nsites] = ts.getCovMatrix()[12];
640  _track_record->cov_predicted_tanLz0[_track_record->nsites] = ts.getCovMatrix()[13];
641  _track_record->cov_predicted_tanLtanL[_track_record->nsites] = ts.getCovMatrix()[14];
642 
643 
644  streamlog_out(DEBUG2) << "DiagnosticsController::record_site get FILTERED trackstate " << std::endl;
645 
646  TKalTrackState& trkState_filtered = (TKalTrackState&) site->GetState(TVKalSite::kFiltered); // this segfaults if no hits are present
647 
648  THelicalTrack helix_filtered = trkState_filtered.GetHelix() ;
649  TMatrixD c0_filtered(trkState_filtered.GetCovMat());
650 
651  DF.UnitMatrix();
652  helix_filtered.MoveTo( tpoint , dPhi , &DF , &c0_filtered) ; // move helix to desired point, and get propagator matrix
653 
654  IMPL::TrackStateImpl ts_f;
655 
656  _current_track->ToLCIOTrackState( helix_filtered, c0_filtered, ts_f, chi2, ndf );
657 
658  _track_record->d0_filtered[_track_record->nsites] = ts_f.getD0() ;
659  _track_record->phi0_filtered[_track_record->nsites] = ts_f.getPhi() ;
660  _track_record->omega_filtered[_track_record->nsites] = ts_f.getOmega() ;
661  _track_record->z0_filtered[_track_record->nsites] = ts_f.getZ0() ;
662  _track_record->tanL_filtered[_track_record->nsites] = ts_f.getTanLambda() ;
663 
664 
665  _track_record->cov_filtered_d0d0[_track_record->nsites] = ts_f.getCovMatrix()[0];
666  _track_record->cov_filtered_phi0d0[_track_record->nsites] = ts_f.getCovMatrix()[1];
667  _track_record->cov_filtered_phi0phi0[_track_record->nsites] = ts_f.getCovMatrix()[2];
668  _track_record->cov_filtered_omegad0[_track_record->nsites] = ts_f.getCovMatrix()[3];
669  _track_record->cov_filtered_omegaphi0[_track_record->nsites] = ts_f.getCovMatrix()[4];
670  _track_record->cov_filtered_omegaomega[_track_record->nsites] = ts_f.getCovMatrix()[5];
671  _track_record->cov_filtered_z0d0[_track_record->nsites] = ts_f.getCovMatrix()[6];
672  _track_record->cov_filtered_z0phi0[_track_record->nsites] = ts_f.getCovMatrix()[7];
673  _track_record->cov_filtered_z0omega[_track_record->nsites] = ts_f.getCovMatrix()[8];
674  _track_record->cov_filtered_z0z0[_track_record->nsites] = ts_f.getCovMatrix()[9];
675  _track_record->cov_filtered_tanLd0[_track_record->nsites] = ts_f.getCovMatrix()[10];
676  _track_record->cov_filtered_tanLphi0[_track_record->nsites] = ts_f.getCovMatrix()[11];
677  _track_record->cov_filtered_tanLomega[_track_record->nsites] = ts_f.getCovMatrix()[12];
678  _track_record->cov_filtered_tanLz0[_track_record->nsites] = ts_f.getCovMatrix()[13];
679  _track_record->cov_filtered_tanLtanL[_track_record->nsites] = ts_f.getCovMatrix()[14];
680 
681 
682 
683  _track_record->ref_point_x[_track_record->nsites] = tpoint.X();
684  _track_record->ref_point_y[_track_record->nsites] = tpoint.Y();
685  _track_record->ref_point_z[_track_record->nsites] = tpoint.Z();
686 
687 
688  }
689 
690  ++_track_record->nsites;
691 
692  if (_track_record->nsites > MAX_SITES) {
693  _skip_track = true ;
694  this->clear_track_record();
695  streamlog_out(DEBUG4) << " DiagnosticsController::end_track: Number of site too large. Track Skipped. nsites = " << _track_record->nsites << " : maximum number of site allowed = " << MAX_SITES << std::endl;
696  return;
697  }
698 
699 
700  streamlog_out(DEBUG2) << "_track_record->nsites = " << _track_record->nsites << std::endl;
701  }
702 
703  void DiagnosticsController::end_track(){
704 
705 
706 
707  if ( _recording_on == false ) {
708  return;
709  }
710 
711  if ( _initialised == false ){
712  streamlog_out(ERROR) << "DiagnosticsController::end_track: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
713  << __FILE__
714  << " line "
715  << __LINE__
716  << std::endl;
717 
718  exit(1);
719  }
720 
721  streamlog_out(DEBUG3) << " DiagnosticsController::end_track called " << std::endl;
722 
723  if ( _skip_track ) { // just clear the track buffers and return.
724  ++_ntracks_skipped;
725  this->clear_track_record();
726  return;
727  } else {
728 
729  _track_record->chi2 = _current_track->_kaltrack->GetChi2();
730  _track_record->ndf = _current_track->_kaltrack->GetNDF();
731  _track_record->prob = TMath::Prob(_track_record->chi2, _track_record->ndf);
732 
733  TIter it(_current_track->_kaltrack,kIterForward);
734 
735  Int_t nsites = _current_track->_kaltrack->GetEntries();
736 
737  if (nsites > MAX_SITES) {
738  _skip_track = true ;
739  ++_ntracks_skipped;
740  this->clear_track_record();
741  streamlog_out(DEBUG4) << " DiagnosticsController::end_track: Number of site too large. Track Skipped. nsites = " << nsites << " : maximum number of site allowed = " << MAX_SITES << std::endl;
742  return;
743  }
744 
745 
746  if( _current_track->_smoothed ){
747 
748  for (Int_t isite=1; isite<nsites; isite++) {
749 
750  UTIL::BitField64 encoder(lcio::LCTrackerCellID::encoding_string());
751  encoder.setValue( _track_record->CellID0[isite] );
752 
753 
754  TVKalSite* site = static_cast<TVKalSite *>( _current_track->_kaltrack->At(isite));
755 
756  if ( _track_record->rejected[isite] == 0 && encoder[lcio::LCTrackerCellID::subdet()] != 0 ) {
757 
758 
759  _track_record->chi2_inc_smoothed[isite] = site->GetDeltaChi2();
760 
761  TKalTrackState& trkState_smoothed = (TKalTrackState&) site->GetState(TVKalSite::kSmoothed); // this segfaults if no hits are present
762 
763  THelicalTrack helix_smoothed = trkState_smoothed.GetHelix() ;
764  TMatrixD c0_smoothed(trkState_smoothed.GetCovMat());
765 
766  Int_t sdim = trkState_smoothed.GetDimension(); // dimensions of the track state, it will be 5 or 6
767 
768  // move track state to the sim hit for comparison
769  const TVector3 tpoint( _track_record->ref_point_x[isite], _track_record->ref_point_y[isite], _track_record->ref_point_z[isite] ) ;
770 
771  // now move to the point
772  TKalMatrix DF(sdim,sdim);
773  DF.UnitMatrix();
774 
775  double dPhi=0;
776 
777  helix_smoothed.MoveTo( tpoint , dPhi , &DF , &c0_smoothed) ; // move helix to desired point, and get propagator matrix
778 
779  IMPL::TrackStateImpl ts;
780  int ndf;
781  double chi2;
782 
783  _current_track->ToLCIOTrackState( helix_smoothed, c0_smoothed, ts, chi2, ndf );
784 
785 
786  _track_record->d0_smoothed[isite] = ts.getD0() ;
787  _track_record->phi0_smoothed[isite] = ts.getPhi() ;
788  _track_record->omega_smoothed[isite] = ts.getOmega() ;
789  _track_record->z0_smoothed[isite] = ts.getZ0() ;
790  _track_record->tanL_smoothed[isite] = ts.getTanLambda() ;
791 
792 
793  _track_record->cov_smoothed_d0d0[isite] = ts.getCovMatrix()[0];
794  _track_record->cov_smoothed_phi0d0[isite] = ts.getCovMatrix()[1];
795  _track_record->cov_smoothed_phi0phi0[isite] = ts.getCovMatrix()[2];
796  _track_record->cov_smoothed_omegad0[isite] = ts.getCovMatrix()[3];
797  _track_record->cov_smoothed_omegaphi0[isite] = ts.getCovMatrix()[4];
798  _track_record->cov_smoothed_omegaomega[isite] = ts.getCovMatrix()[5];
799  _track_record->cov_smoothed_z0d0[isite] = ts.getCovMatrix()[6];
800  _track_record->cov_smoothed_z0phi0[isite] = ts.getCovMatrix()[7];
801  _track_record->cov_smoothed_z0omega[isite] = ts.getCovMatrix()[8];
802  _track_record->cov_smoothed_z0z0[isite] = ts.getCovMatrix()[9];
803  _track_record->cov_smoothed_tanLd0[isite] = ts.getCovMatrix()[10];
804  _track_record->cov_smoothed_tanLphi0[isite] = ts.getCovMatrix()[11];
805  _track_record->cov_smoothed_tanLomega[isite] = ts.getCovMatrix()[12];
806  _track_record->cov_smoothed_tanLz0[isite] = ts.getCovMatrix()[13];
807  _track_record->cov_smoothed_tanLtanL[isite] = ts.getCovMatrix()[14];
808 
809  }
810 
811  }
812  }
813 
814  IMPL::TrackStateImpl ts_at_ip;
815  int ndf;
816  double chi2;
817 
818  Vector3D point(0.0,0.0,0.0);
819  _current_track->propagate( point, ts_at_ip, chi2, ndf );
820 
821  _track_record->d0_ip = ts_at_ip.getD0() ;
822  _track_record->phi0_ip = ts_at_ip.getPhi() ;
823  _track_record->omega_ip = ts_at_ip.getOmega() ;
824  _track_record->z0_ip = ts_at_ip.getZ0() ;
825  _track_record->tanL_ip = ts_at_ip.getTanLambda() ;
826 
827 
828  _track_record->cov_ip_d0d0 = ts_at_ip.getCovMatrix()[0];
829  _track_record->cov_ip_phi0d0 = ts_at_ip.getCovMatrix()[1];
830  _track_record->cov_ip_phi0phi0 = ts_at_ip.getCovMatrix()[2];
831  _track_record->cov_ip_omegad0 = ts_at_ip.getCovMatrix()[3];
832  _track_record->cov_ip_omegaphi0 = ts_at_ip.getCovMatrix()[4];
833  _track_record->cov_ip_omegaomega = ts_at_ip.getCovMatrix()[5];
834  _track_record->cov_ip_z0d0 = ts_at_ip.getCovMatrix()[6];
835  _track_record->cov_ip_z0phi0 = ts_at_ip.getCovMatrix()[7];
836  _track_record->cov_ip_z0omega = ts_at_ip.getCovMatrix()[8];
837  _track_record->cov_ip_z0z0 = ts_at_ip.getCovMatrix()[9];
838  _track_record->cov_ip_tanLd0 = ts_at_ip.getCovMatrix()[10];
839  _track_record->cov_ip_tanLphi0 = ts_at_ip.getCovMatrix()[11];
840  _track_record->cov_ip_tanLomega = ts_at_ip.getCovMatrix()[12];
841  _track_record->cov_ip_tanLz0 = ts_at_ip.getCovMatrix()[13];
842  _track_record->cov_ip_tanLtanL = ts_at_ip.getCovMatrix()[14];
843 
844  _tree->Fill();
845 
846  ++_ntracks_written;
847 
848  }
849 
850  }
851 
852  void DiagnosticsController::end(){
853 
854  if ( _recording_on == false ) {
855  return;
856  }
857 
858  if ( _initialised == false ){
859  streamlog_out(ERROR) << "DiagnosticsController::end: Diagnostics not initialised call init(std::string root_file_name, std::string root_tree_name, bool recording_off) first : exit(1) called from file "
860  << __FILE__
861  << " line "
862  << __LINE__
863  << std::endl;
864 
865  exit(1);
866  }
867 
868  streamlog_out(DEBUG4) << " DiagnosticsController::end() called \n"
869  << "\t number of tracks written = " << _ntracks_written
870  << "\t number of tracks skipped = " << _ntracks_skipped
871  << std::endl;
872 
873 
874 
875  // _tree->Print();
876  _root_file->Write();
877  _root_file->Close();
878  delete _root_file; _root_file = 0;
879 
880  }
881 
882 
883 }
884 
885 
886 
887 #endif
EVENT::SimTrackerHit::getPosition
virtual const double * getPosition() const =0
IMPL::TrackStateImpl::getZ0
virtual float getZ0() const
EVENT::TrackerHit::getPosition
virtual const double * getPosition() const =0
IMPL::TrackStateImpl::getOmega
virtual float getOmega() const
std::endl
T endl(T...args)
EVENT::MCParticle::getPDG
virtual int getPDG() const =0
EVENT::SimTrackerHit::getMomentum
virtual const float * getMomentum() const =0
MarlinTrk::Vector3D
dd4hep::rec::Vector3D Vector3D
the Vector3D used for the tracking interface
Definition: IMarlinTrack.h:26
std::atan2
T atan2(T...args)
std::string
STL class.
EVENT::MCParticle::getEnergy
virtual double getEnergy() const =0
std::exit
T exit(T...args)
EVENT::SimTrackerHit::getMCParticle
virtual MCParticle * getMCParticle() const =0
EVENT::MCParticle::getVertex
virtual const double * getVertex() const =0
IMPL::TrackStateImpl::getTanLambda
virtual float getTanLambda() const
IMPL::TrackStateImpl::getD0
virtual float getD0() const
std::fabs
T fabs(T...args)
EVENT::TrackerHit::getCellID0
virtual int getCellID0() const =0
IMPL::TrackStateImpl::getCovMatrix
virtual const EVENT::FloatVec & getCovMatrix() const
EVENT::MCParticle::getCharge
virtual float getCharge() const =0
MAX_SITES
#define MAX_SITES
Definition: MarlinTrkNtuple.h:15
HelixTrack::moveRefPoint
double moveRefPoint(double x, double y, double z)
Definition: HelixTrack.cc:66
std::sqrt
T sqrt(T...args)
EVENT::MCParticle::getMomentum
virtual const double * getMomentum() const =0
IMPL::TrackStateImpl::getPhi
virtual float getPhi() const
Generated on Sun Sep 26 2021 03:45:50 for MarlinTrk by   doxygen 1.8.5