BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
BeamParams Class Reference

#include <BeamParams.h>

Inheritance diagram for BeamParams:

Public Member Functions

 BeamParams (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Detailed Description

Definition at line 11 of file BeamParams.h.

Constructor & Destructor Documentation

◆ BeamParams()

BeamParams::BeamParams ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 45 of file BeamParams.cxx.

46 : Algorithm( name, pSvcLocator ) {
47 // Declare the properties
48 // declareProperty("RunNo", m_runNo = 9999);
49
50 declareProperty( "ParVer", m_parVer = 1 );
51 declareProperty( "TrackNumberCut", m_trackNumberCut = 1 );
52 declareProperty( "Vz0Cut", m_vz0Cut = 50.0 );
53 declareProperty( "CosThetaCut", m_cosThetaCut = 0.93 );
54
55 // particle identification
56 declareProperty( "OptionPID", m_pid = 0 );
57 declareProperty( "PidUseDedx", m_useDedx = true );
58 declareProperty( "PidUseTof1", m_useTof1 = true );
59 declareProperty( "PidUseTof2", m_useTof2 = true );
60 declareProperty( "PidUseTofE", m_useTofE = false );
61 declareProperty( "PidUseTofQ", m_useTofQ = false );
62 declareProperty( "PidUseEmc", m_useEmc = false );
63 declareProperty( "PidUseMuc", m_useMuc = false );
64
65 // vertex finding
66 declareProperty( "OptionVertexFind", m_vertexFind = 0 );
67 declareProperty( "MinDistance", m_minDistance = 100.0 );
68 declareProperty( "MinPointX", m_minPointX = 0.1 );
69 declareProperty( "MinPointY", m_minPointY = 0.1 );
70 declareProperty( "MeanPointX", m_meanPointX = 0.1 );
71 declareProperty( "MeanPointY", m_meanPointY = -0.25 );
72
73 // Kalman vertex fitting
74 declareProperty( "ChisqCut", m_chisqCut = 200.0 );
75 declareProperty( "TrackIteration", m_trackIteration = 5 );
76 declareProperty( "VertexIteration", m_vertexIteration = 5 );
77 declareProperty( "Chi2CutforTrkIter", m_chi2CutforTrkIter = 0.1 );
78 declareProperty( "Chi2CutforSmooth", m_chi2CutforSmooth = 200.0 );
79
80 // output file name
81 declareProperty( "HadronFile", m_hadronFile = 0 );
82 declareProperty( "DstFileName", m_fileNameDst = std::string( "dst.txt" ) );
83 declareProperty( "HadronFileName", m_fileNameHadron = std::string( "hadron.txt" ) );
84 declareProperty( "FigsName", m_figsName = std::string( "figs.ps" ) );
85}

Referenced by BeamParams().

Member Function Documentation

◆ execute()

StatusCode BeamParams::execute ( )

Definition at line 294 of file BeamParams.cxx.

294 {
295 MsgStream log( msgSvc(), name() );
296 log << MSG::INFO << "in execute()" << endmsg;
297
298 ///////////////////
299 // Read REC data
300 //////////////////
301 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
302 SmartDataPtr<EvtRecEvent> recEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
303 SmartDataPtr<EvtRecTrackCol> recTrackCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
304 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber() << " "
305 << eventHeader->eventNumber() << endmsg;
306 log << MSG::DEBUG << "ncharg, nneu, tottks = " << recEvent->totalCharged() << " , "
307 << recEvent->totalNeutral() << " , " << recEvent->totalTracks() << endmsg;
308
309 m_runNo = eventHeader->runNumber();
310 // FIXME : cout
311 if ( eventHeader->eventNumber() % 5000 == 0 )
312 {
313 cout << "Event number = " << eventHeader->eventNumber() << " run number = " << m_runNo
314 << endl;
315 }
316 m_sel_number[0]++;
317
318 // Cut 1 : Number of tracks. For anything sample, at least 3 charged tracks
319 if ( recEvent->totalCharged() < m_trackNumberCut ) return StatusCode::SUCCESS;
320 m_sel_number[1]++;
321
322 ///////////////////////////////////////////////////
323 // Select good charged tracks in MDC ( option for PID )
324 //////////////////////////////////////////////////
325 Vint icp, icm, iGood, jGood;
326 icp.clear();
327 icm.clear();
328 iGood.clear();
329 jGood.clear();
330
331 map<int, int> ParticleType; // 0:e, 1:mu, 2:pi, 3:K, 4:p
332 ParticleID* pid = ParticleID::instance();
333
334 for ( unsigned int i = 0; i < recEvent->totalCharged(); i++ )
335 {
336 jGood.push_back( 0 );
337 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
338 if ( !( *itTrk )->isMdcTrackValid() ) continue;
339 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
340 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
341 if ( fabs( cos( mdcTrk->theta() ) ) >= m_cosThetaCut ) continue;
342 if ( fabs( mdcTrk->z() ) >= m_vz0Cut ) continue;
343 // iGood.push_back((*itTrk)->trackId());
344 iGood.push_back( i );
345
346 if ( m_pid == 0 )
347 {
348 if ( mdcTrk->charge() > 0 )
349 {
350 // icp.push_back((*itTrk)->trackId());
351 icp.push_back( i );
352 }
353 else if ( mdcTrk->charge() < 0 )
354 {
355 // icm.push_back((*itTrk)->trackId());
356 icm.push_back( i );
357 }
358 }
359 else if ( m_pid == 1 )
360 {
361 ParticleType[i] = 2; // FIXME default value is Pion ?
362 // pid info
363 pid->init();
364 pid->setChiMinCut( 4 );
365 pid->setRecTrack( *itTrk );
366 pid->setMethod( pid->methodProbability() );
367 // pid->setMethod(pid->methodLikelihood());
368 // use PID sub-system
369 if ( m_useDedx ) pid->usePidSys( pid->useDedx() );
370 if ( m_useTof1 ) pid->usePidSys( pid->useTof1() );
371 if ( m_useTof2 ) pid->usePidSys( pid->useTof2() );
372 if ( m_useTofE ) pid->usePidSys( pid->useTofE() );
373 if ( m_useTofQ ) pid->usePidSys( pid->useTofQ() );
374 if ( m_useEmc ) pid->usePidSys( pid->useEmc() );
375 if ( m_useMuc ) pid->usePidSys( pid->useMuc() );
376 pid->identify( pid->onlyElectron() );
377 pid->identify( pid->onlyMuon() );
378 pid->identify( pid->onlyPion() );
379 pid->identify( pid->onlyKaon() );
380 pid->identify( pid->onlyProton() );
381 pid->calculate();
382 if ( !pid->IsPidInfoValid() ) continue;
383 double prob_value[5];
384 prob_value[0] = pid->probElectron();
385 prob_value[1] = pid->probMuon();
386 prob_value[2] = pid->probPion();
387 prob_value[3] = pid->probKaon();
388 prob_value[4] = pid->probProton();
389
390 int max = 0;
391 for ( int i = 1; i < 5; i++ )
392 {
393 if ( prob_value[i] > prob_value[max] ) { max = i; }
394 }
395
396 // cout << "PID : n = " << n << endl;
397
398 if ( max == 0 ) ParticleType[i] = 0; // m_sel_number[3]++;}
399 if ( max == 1 ) ParticleType[i] = 1; // m_sel_number[4]++;}
400 if ( max == 2 ) ParticleType[i] = 2; // m_sel_number[5]++;}
401 if ( max == 3 ) ParticleType[i] = 3; // m_sel_number[6]++;}
402 if ( max == 4 ) ParticleType[i] = 4; // m_sel_number[7]++;}
403 }
404
405 // fill z position of track parameters
406 RecMdcKalTrack* kalTrk = ( *itTrk )->mdcKalTrack();
408 m_mdcTrk_x = mdcTrk->x();
409 m_mdcTrk_y = mdcTrk->y();
410 m_mdcTrk_z = mdcTrk->helix( 3 );
411 m_mdcTrk_r = mdcTrk->r();
412 m_mdcTrk_dr = mdcTrk->helix( 0 );
413 m_mdcKalTrk_z = kalTrk->helix()[3];
414 m_tuple6->write();
415 }
416
417 // Cut 2 : for anything sample, at least 3 good charged tracks
418 if ( iGood.size() < m_trackNumberCut ) return StatusCode::SUCCESS;
419 m_sel_number[2]++;
420
421 /////////////////////
422 // Vertex Finding
423 ////////////////////
424 if ( m_vertexFind == 1 )
425 {
426 int nGood = iGood.size();
427 for ( int i1 = 0; i1 < nGood; i1++ )
428 {
429 int id_trk1 = iGood[i1];
430 EvtRecTrackIterator trk1 = ( recTrackCol->begin() + id_trk1 );
431 RecMdcKalTrack* kalTrk1 = ( *trk1 )->mdcKalTrack();
432 // FIXME default set is ?
433 HTrackParameter htrk1, htrk2;
434 if ( m_pid == 0 )
435 {
437 htrk1 = HTrackParameter( kalTrk1->helix(), kalTrk1->err(), id_trk1, 2 );
438 }
439 else if ( m_pid == 1 )
440 {
441 RecMdcKalTrack::setPidType( (RecMdcKalTrack::PidType)ParticleType[id_trk1] );
442 htrk1 = HTrackParameter( kalTrk1->helix(), kalTrk1->err(), id_trk1,
443 ParticleType[id_trk1] );
444 }
445 for ( int i2 = i1 + 1; i2 < nGood; i2++ )
446 {
447 int id_trk2 = iGood[i2];
448 EvtRecTrackIterator trk2 = ( recTrackCol->begin() + id_trk2 );
449 RecMdcKalTrack* kalTrk2 = ( *trk2 )->mdcKalTrack();
450 // FIXME default set is ?
451 if ( m_pid == 0 )
452 {
454 htrk2 = HTrackParameter( kalTrk2->helix(), kalTrk2->err(), id_trk2, 2 );
455 }
456 else if ( m_pid == 1 )
457 {
458 RecMdcKalTrack::setPidType( (RecMdcKalTrack::PidType)ParticleType[id_trk2] );
459 htrk2 = HTrackParameter( kalTrk2->helix(), kalTrk2->err(), id_trk2,
460 ParticleType[id_trk2] );
461 }
462 HepPoint3D cross( 0., 0., 0. );
463 double dist = htrk1.minDistanceTwoHelix( htrk2, cross );
464 m_mind = dist;
465 m_xc = cross.x();
466 m_yc = cross.y();
467 m_zc = cross.z();
468 m_tuple1->write();
469
470 if ( sqrt( cross.x() * cross.x() + cross.y() * cross.y() ) > 2 ) continue;
471
472 if ( fabs( dist ) > m_minDistance ) continue; // Cut condition 2
473 // double cross_cut = (cross.x() - m_meanPointX) * (cross.x() - m_meanPointX)
474 // /m_minPointX/m_minPointX
475 // + (cross.y() - m_meanPointY) * (cross.y() - m_meanPointY)
476 // /m_minPointY/m_minPointY;
477 // //FIXME sigma value from Gaussian fit //0.071/0.059
478 double cross_cut = 0.;
479 if ( cross_cut < 3.5 * 3.5 )
480 { // Cut condition 3
481 jGood[id_trk1] = 1;
482 jGood[id_trk2] = 1;
483 } // end if
484 }
485 }
486
487 iGood.clear();
488 for ( int i = 0; i < jGood.size(); i++ )
489 {
490 if ( jGood[i] == 1 ) iGood.push_back( i );
491 }
492
493 } // end if vertex finding
494
495 if ( iGood.size() >= 2 ) m_sel_number[3]++;
496 if ( iGood.size() >= 3 ) m_sel_number[4]++;
497
498 ///////////////////////////////
499 // Vertex Fit
500 //////////////////////////////
501 HepPoint3D vx( 0., 0., 0. );
502 HepSymMatrix Evx( 3, 0 );
503 double bx = 1E+6;
504 double by = 1E+6;
505 double bz = 1E+6;
506 Evx[0][0] = bx * bx;
507 Evx[1][1] = by * by;
508 Evx[2][2] = bz * bz;
509 VertexParameter vxpar;
510 vxpar.setVx( vx );
511 vxpar.setEvx( Evx );
512
513 VertexFit* vtxfit = VertexFit::instance();
514 // Step 1: using MdcKalTrk parameters
515 vtxfit->init();
516 Vint tlis;
517 tlis.clear();
518 for ( int i = 0; i < iGood.size(); i++ )
519 {
520 int trk_id = iGood[i];
521 EvtRecTrackIterator itTrk = recTrackCol->begin() + trk_id;
522 RecMdcKalTrack* kalTrk = ( *itTrk )->mdcKalTrack();
523 if ( m_pid == 0 )
524 {
526 WTrackParameter wtrk( xmass[2], kalTrk->helix(), kalTrk->err() );
527 vtxfit->AddTrack( i, wtrk );
528 }
529 else if ( m_pid == 1 )
530 {
532 WTrackParameter wtrk( xmass[ParticleType[trk_id]], kalTrk->helix(), kalTrk->err() );
533 vtxfit->AddTrack( i, wtrk );
534 }
535 tlis.push_back( i );
536 }
537 // if(iGood.size() >= m_trackNumberCut) { //at least N tracks
538 vtxfit->AddVertex( 0, vxpar, tlis );
539 if ( !vtxfit->Fit( 0 ) ) return StatusCode::SUCCESS;
540 m_chig = vtxfit->chisq( 0 );
541 m_ndofg = 2 * iGood.size() - 3;
542 m_probg = TMath::Prob( m_chig, m_ndofg );
543 HepVector glvtx = vtxfit->Vx( 0 );
544 m_gvx = glvtx[0];
545 m_gvy = glvtx[1];
546 m_gvz = glvtx[2];
547 m_tuple4->write();
548
549 if ( !( m_vertex_x->Fill( glvtx[0] ) ) )
550 {
551 log << MSG::FATAL << "Couldn't Fill x of vertex" << endmsg;
552 exit( 1 );
553 }
554 if ( !( m_vertex_y->Fill( glvtx[1] ) ) )
555 {
556 log << MSG::FATAL << "Couldn't Fill y of vertex" << endmsg;
557 exit( 1 );
558 }
559 if ( !( m_vertex_z->Fill( glvtx[2] ) ) )
560 {
561 log << MSG::FATAL << "Couldn't Fill z of vertex" << endmsg;
562 exit( 1 );
563 }
564 //}
565 for ( int j = 0; j < iGood.size(); j++ )
566 {
567 HepVector pull( 5, 0 );
568 bool success = vtxfit->pull( 0, j, pull );
569 if ( success )
570 {
571 m_gpull_drho = pull[0];
572 m_gpull_phi = pull[1];
573 m_gpull_kapha = pull[2];
574 m_gpull_dz = pull[3];
575 m_gpull_lamb = pull[4];
576 m_tuple7->write();
577 }
578 }
579
580 // Step 2 : Kalman Vertex Fitting
581 KalmanVertexFit* kalvtx = KalmanVertexFit::instance();
582 kalvtx->init();
583 kalvtx->initVertex( vxpar );
584 kalvtx->setChisqCut( m_chisqCut );
585 kalvtx->setTrackIteration( m_trackIteration );
586 kalvtx->setVertexIteration( m_vertexIteration );
587 kalvtx->setTrackIterationCut( m_chi2CutforTrkIter );
588 for ( int i = 0; i < iGood.size(); i++ )
589 {
590 int trk_id = iGood[i];
591 EvtRecTrackIterator itTrk = recTrackCol->begin() + trk_id;
592 RecMdcKalTrack* kalTrk = ( *itTrk )->mdcKalTrack();
593 if ( m_pid == 0 )
594 {
596 HTrackParameter htrk( kalTrk->helix(), kalTrk->err(), trk_id, 2 );
597 kalvtx->addTrack( htrk );
598 }
599 else if ( m_pid == 1 )
600 {
602 HTrackParameter htrk( kalTrk->helix(), kalTrk->err(), trk_id, ParticleType[trk_id] );
603 kalvtx->addTrack( htrk );
604 }
605 }
606 kalvtx->filter();
607 // int freedomCut = -3 + 2*m_trackNumberCut;
608 int freedomCut = 1;
609 if ( kalvtx->ndof() >= freedomCut )
610 { // Cut condition 5 //add 2 when add every track
611 // for(int i = 0; i < kalvtx->numTrack(); i++) {
612 for ( int i = 0; i < iGood.size(); i++ )
613 {
614 m_chif = kalvtx->chiF( i );
615 m_chis = kalvtx->chiS( i );
616 m_probs = TMath::Prob( m_chis, 2 );
617 m_probf = TMath::Prob( m_chif, 2 );
618 m_tuple2->write();
619
620 if ( kalvtx->chiS( i ) > m_chi2CutforSmooth ) kalvtx->remove( i );
621 }
622 }
623 if ( kalvtx->ndof() >= freedomCut )
624 { // FIXME to Rhopi is 0 , others are 1
625 for ( int i = 0; i < iGood.size(); i++ )
626 {
627 kalvtx->smooth( i );
628
629 HepVector Pull( 5, 0 );
630 Pull = kalvtx->pull( i );
631 m_pull_drho = Pull[0];
632 m_pull_phi = Pull[1];
633 m_pull_kapha = Pull[2];
634 m_pull_dz = Pull[3];
635 m_pull_lamb = Pull[4];
636 m_pull_momentum = kalvtx->pullmomentum( i );
637 m_tuple5->write();
638 }
639
640 m_chik = kalvtx->chisq();
641 m_ndofk = kalvtx->ndof();
642 m_probk = TMath::Prob( m_chik, m_ndofk );
643 HepVector xp( 3, 0 );
644 xp = kalvtx->x();
645 m_kvx = xp[0];
646 m_kvy = xp[1];
647 m_kvz = xp[2];
648 m_tuple3->write();
649
650 if ( !( m_vertex_x_kal->Fill( xp[0] ) ) )
651 {
652 log << MSG::FATAL << "Couldn't Fill kal x of vertex" << endmsg;
653 exit( 1 );
654 }
655 if ( !( m_vertex_y_kal->Fill( xp[1] ) ) )
656 {
657 log << MSG::FATAL << "Couldn't Fill kal y of vertex" << endmsg;
658 exit( 1 );
659 }
660 if ( !( m_vertex_z_kal->Fill( xp[2] ) ) )
661 {
662 log << MSG::FATAL << "Couldn't Fill kal z of vertex" << endmsg;
663 exit( 1 );
664 }
665
666 m_sel_number[3]++;
667 }
668 return StatusCode::SUCCESS;
669}
HepGeom::Point3D< double > HepPoint3D
#define max(a, b)
EvtRecTrackCol::iterator EvtRecTrackIterator
EvtVector3R cross(const EvtVector3R &p1, const EvtVector3R &p2)
const double xmass[5]
Definition Gam4pikp.cxx:35
std::vector< int > Vint
Definition Gam4pikp.cxx:37
IMessageSvc * msgSvc()
const HepVector helix() const
......
double minDistanceTwoHelix(const HTrackParameter G, HepPoint3D &pos)
double pullmomentum(const int k)
void addTrack(const HTrackParameter)
void smooth(const int k)
int filter(const int k)
double chiS(const int k) const
void initVertex(const VertexParameter vtx)
HepVector pull(const int k)
static KalmanVertexFit * instance()
void remove(const int k)
static ParticleID * instance()
bool IsPidInfoValid() const
void calculate()
void init()
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
void init()
Definition VertexFit.cxx:27
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
bool pull(int n, int itk, HepVector &p)
static VertexFit * instance()
Definition VertexFit.cxx:15
bool Fit()

◆ finalize()

StatusCode BeamParams::finalize ( )

Definition at line 673 of file BeamParams.cxx.

673 {
674 MsgStream log( msgSvc(), name() );
675 log << MSG::INFO << "in finalize()" << endmsg;
676
677 /*TF1 *func = new TF1("func", "gaus", -0.6, 0.6);
678 m_vertex_x->Fit("func", "NR+");
679 Double_t MeanX = func->GetParameter(1);
680 Double_t SigmaX = func->GetParameter(2);
681 Double_t ErrMeanX = func->GetParError(1);
682 Double_t ErrSigmaX = func->GetParError(2);
683
684 TF1 *funcY = new TF1("funcY", "gaus", -0.6, 0.1);
685 m_vertex_y->Fit("funcY", "NR+");
686 Double_t MeanY = funcY->GetParameter(1);
687 Double_t SigmaY = funcY->GetParameter(2);
688 Double_t ErrMeanY = funcY->GetParError(1);
689 Double_t ErrSigmaY = funcY->GetParError(2);
690
691 TF1 *funcZ = new TF1("funcZ", "gaus", -6, 6);
692 m_vertex_z->Fit("funcZ", "NR+");
693 Double_t MeanZ = funcZ->GetParameter(1);
694 Double_t SigmaZ = funcZ->GetParameter(2);
695 Double_t ErrMeanZ = funcZ->GetParError(1);
696 Double_t ErrSigmaZ = funcZ->GetParError(2);
697
698 TCanvas *myC = new TCanvas("myC", "myC", 1200, 400);
699 TPad *background = (TPad*)gPad;
700 background->SetFillColor(10);
701 myC->Divide(3,1);
702 myC->cd(1);
703
704 m_vertex_x_kal->Fit("func", "R+");
705 Double_t MeanXKal = func->GetParameter(1);
706 Double_t SigmaXKal = func->GetParameter(2);
707 Double_t ErrMeanXKal = func->GetParError(1);
708 Double_t ErrSigmaXKal = func->GetParError(2);
709
710 myC->cd(2);
711 m_vertex_y_kal->Fit("funcY", "R+");
712 Double_t MeanYKal = funcY->GetParameter(1);
713 Double_t SigmaYKal = funcY->GetParameter(2);
714 Double_t ErrMeanYKal = funcY->GetParError(1);
715 Double_t ErrSigmaYKal = funcY->GetParError(2);
716
717 myC->cd(3);
718 m_vertex_z_kal->Fit("funcZ", "R+");
719 Double_t MeanZKal = funcZ->GetParameter(1);
720 Double_t SigmaZKal = funcZ->GetParameter(2);
721 Double_t ErrMeanZKal = funcZ->GetParError(1);
722 Double_t ErrSigmaZKal = funcZ->GetParError(2);
723
724 cout << "Kal: Mean X, Y, Z = "<<MeanXKal << " " << MeanYKal << " " << MeanZKal << endl;
725 cout << "Kal: Sigma X, Y, Z = "<<SigmaXKal<<" " <<SigmaYKal <<" " << SigmaZKal << endl;
726 cout << "Gl: Mean X, Y, Z = " << MeanX << " " << MeanY << " " << MeanZ << endl;
727 cout << "Gl: Sigma X, Y, Z = " << SigmaX << " " << SigmaY << " " << SigmaZ << endl;
728
729 //////////////////////////////////////////////////////////////////
730 // Output a TXT file and a .ps figure for storing the fit results
731 //////////////////////////////////////////////////////////////////
732 const char *file_name, *figs_name;
733 if (m_hadronFile == 0) {
734 file_name = m_fileNameDst.c_str();
735 } else if (m_hadronFile == 1) {
736 file_name = m_fileNameHadron.c_str();
737 }
738
739 figs_name = m_figsName.c_str();
740 myC->SaveAs(figs_name);
741
742 ofstream file(file_name);
743
744 file << getenv("BES_RELEASE") << " " << m_parVer << " " << m_runNo << endl;
745 file << MeanX << " " << MeanY << " " << MeanZ << " "<< SigmaX << " "<< SigmaY << " " <<
746 SigmaZ << endl; file << ErrMeanX << " " << ErrMeanY<< " " << ErrMeanZ << " " << ErrSigmaX <<
747 " " << ErrSigmaY << " " << ErrSigmaZ << endl; file << MeanXKal << " "<< MeanYKal << " "<<
748 MeanZKal << " "<< SigmaXKal << " " << SigmaYKal << " " << SigmaZKal << endl; file <<
749 ErrMeanXKal << " " << ErrMeanYKal<< " " << ErrMeanZKal << " " << ErrSigmaXKal << " " <<
750 ErrSigmaYKal << " " << ErrSigmaZKal << endl;*/
751
752 /*delete m_vertex_x;
753 delete m_vertex_y;
754 delete m_vertex_z;
755 delete m_vertex_x_kal;
756 delete m_vertex_y_kal;
757 delete m_vertex_z_kal;
758 */
759 ////////////////////////////////////////////////
760 log << MSG::ALWAYS << "==================================" << endmsg;
761 log << MSG::ALWAYS << "survived event :";
762 for ( int i = 0; i < 10; i++ ) { log << MSG::ALWAYS << m_sel_number[i] << " "; }
763 log << MSG::ALWAYS << endmsg;
764 log << MSG::ALWAYS << "==================================" << endmsg;
765 return StatusCode::SUCCESS;
766}

◆ initialize()

StatusCode BeamParams::initialize ( )

Definition at line 88 of file BeamParams.cxx.

88 {
89 MsgStream log( msgSvc(), name() );
90 log << MSG::INFO << "in initialize()" << endmsg;
91
92 StatusCode sc = service( "THistSvc", m_thsvc );
93 if ( sc.isFailure() )
94 {
95 log << MSG::FATAL << "Couldn't get THistSvc" << endmsg;
96 exit( 1 );
97 }
98
99 for ( int i = 0; i < 15; i++ ) { m_sel_number[i] = 0; }
100
101 //////////////////////
102 // Book histogram
103 /////////////////////
104 m_vertex_x = new TH1D( "x_of_vertex", "x of vertex", 200, -5, 5 );
105 m_vertex_y = new TH1D( "y_of_vertex", "y of vertex", 200, -5, 5 );
106 m_vertex_z = new TH1D( "z_of_vertex", "z of vertex", 200, -10, 10 );
107 m_vertex_x_kal = new TH1D( "x_of_vertex_in_kal", "x of vertex in kal", 200, -5, 5 );
108 m_vertex_y_kal = new TH1D( "y_of_vertex_in_kal", "y of vertex in kal", 200, -5, 5 );
109 m_vertex_z_kal = new TH1D( "z_of_vertex_in_kal", "z of vertex in kal", 200, -10, 10 );
110 if ( m_thsvc->regHist( "/DQAHist/zhsVER/x_of_vertex", m_vertex_x ).isFailure() )
111 {
112 log << MSG::FATAL << "Couldn't register x of vertex " << endmsg;
113 exit( 1 );
114 }
115 if ( m_thsvc->regHist( "/DQAHist/zhsVER/y_of_vertex", m_vertex_y ).isFailure() )
116 {
117 log << MSG::FATAL << "Couldn't register y of vertex" << endmsg;
118 exit( 1 );
119 }
120 if ( m_thsvc->regHist( "/DQAHist/zhsVER/z_of_vertex", m_vertex_z ).isFailure() )
121 {
122 log << MSG::FATAL << "Couldn't register z of vertex" << endmsg;
123 exit( 1 );
124 }
125 if ( m_thsvc->regHist( "/DQAHist/zhsVER/x_of_vertex_in_kal", m_vertex_x_kal ).isFailure() )
126 {
127 log << MSG::FATAL << "Couldn't register x of vertex in kal" << endmsg;
128 exit( 1 );
129 }
130 if ( m_thsvc->regHist( "/DQAHist/zhsVER/y_of_vertex_in_kal", m_vertex_y_kal ).isFailure() )
131 {
132 log << MSG::FATAL << "Couldn't register y of vertex in kal" << endmsg;
133 exit( 1 );
134 }
135 if ( m_thsvc->regHist( "/DQAHist/zhsVER/z_of_vertex_in_kal", m_vertex_z_kal ).isFailure() )
136 {
137 log << MSG::FATAL << "Couldn't register z of vertex in kal" << endmsg;
138 exit( 1 );
139 }
140 // end book
141
142 StatusCode status;
143
144 NTuplePtr nt1( ntupleSvc(), "FILE1/minid" ); // minimal distance
145 if ( nt1 ) m_tuple1 = nt1;
146 else
147 {
148 m_tuple1 = ntupleSvc()->book( "FILE1/minid", CLID_ColumnWiseTuple, "minimal distance" );
149 if ( m_tuple1 )
150 {
151 status = m_tuple1->addItem( "xc", m_xc );
152 status = m_tuple1->addItem( "yc", m_yc );
153 status = m_tuple1->addItem( "zc", m_zc );
154 status = m_tuple1->addItem( "mind", m_mind );
155 }
156 else
157 {
158 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
159 return StatusCode::FAILURE;
160 }
161 }
162
163 NTuplePtr nt2( ntupleSvc(), "FILE1/chisq" ); // chisq of Kalman filter
164 if ( nt2 ) m_tuple2 = nt2;
165 else
166 {
167 m_tuple2 = ntupleSvc()->book( "FILE1/chisq", CLID_ColumnWiseTuple, "chi-square of " );
168 if ( m_tuple2 )
169 {
170 status = m_tuple2->addItem( "chis", m_chis );
171 status = m_tuple2->addItem( "probs", m_probs );
172 status = m_tuple2->addItem( "chif", m_chif );
173 status = m_tuple2->addItem( "probf", m_probf );
174 }
175 else
176 {
177 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
178 return StatusCode::FAILURE;
179 }
180 }
181
182 NTuplePtr nt3( ntupleSvc(), "FILE1/kalvtx" );
183 if ( nt3 ) m_tuple3 = nt3;
184 else
185 {
186 m_tuple3 = ntupleSvc()->book( "FILE1/kalvtx", CLID_ColumnWiseTuple, "kalman vertex" );
187 if ( m_tuple3 )
188 {
189 status = m_tuple3->addItem( "kvx", m_kvx );
190 status = m_tuple3->addItem( "kvy", m_kvy );
191 status = m_tuple3->addItem( "kvz", m_kvz );
192 status = m_tuple3->addItem( "chik", m_chik );
193 status = m_tuple3->addItem( "ndofk", m_ndofk );
194 status = m_tuple3->addItem( "probk", m_probk );
195 }
196 else
197 {
198 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple3 ) << endmsg;
199 return StatusCode::FAILURE;
200 }
201 }
202
203 NTuplePtr nt4( ntupleSvc(), "FILE1/glbvtx" );
204 if ( nt4 ) m_tuple4 = nt4;
205 else
206 {
207 m_tuple4 = ntupleSvc()->book( "FILE1/glbvtx", CLID_ColumnWiseTuple, "global vertex" );
208 if ( m_tuple4 )
209 {
210 status = m_tuple4->addItem( "gvx", m_gvx );
211 status = m_tuple4->addItem( "gvy", m_gvy );
212 status = m_tuple4->addItem( "gvz", m_gvz );
213 status = m_tuple4->addItem( "chig", m_chig );
214 status = m_tuple4->addItem( "ndofg", m_ndofg );
215 status = m_tuple4->addItem( "probg", m_probg );
216 }
217 else
218 {
219 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
220 return StatusCode::FAILURE;
221 }
222 }
223
224 NTuplePtr nt5( ntupleSvc(), "FILE1/Pull" );
225 if ( nt5 ) m_tuple5 = nt5;
226 else
227 {
228 m_tuple5 = ntupleSvc()->book( "FILE1/Pull", CLID_ColumnWiseTuple, "Pull" );
229 if ( m_tuple5 )
230 {
231 status = m_tuple5->addItem( "pull_drho", m_pull_drho );
232 status = m_tuple5->addItem( "pull_phi", m_pull_phi );
233 status = m_tuple5->addItem( "pull_kapha", m_pull_kapha );
234 status = m_tuple5->addItem( "pull_dz", m_pull_dz );
235 status = m_tuple5->addItem( "pull_lamb", m_pull_lamb );
236 status = m_tuple5->addItem( "pull_momentum", m_pull_momentum );
237 }
238 else
239 {
240 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple5 ) << endmsg;
241 return StatusCode::FAILURE;
242 }
243 }
244
245 NTuplePtr nt6( ntupleSvc(), "FILE1/MdcTrack" );
246 if ( nt6 ) m_tuple6 = nt6;
247 else
248 {
249 m_tuple6 = ntupleSvc()->book( "FILE1/MdcTrack", CLID_ColumnWiseTuple, "MdcTrack" );
250 if ( m_tuple6 )
251 {
252 status = m_tuple6->addItem( "MdcTrkX", m_mdcTrk_x );
253 status = m_tuple6->addItem( "MdcTrkY", m_mdcTrk_y );
254 status = m_tuple6->addItem( "MdcTrkZ", m_mdcTrk_z );
255 status = m_tuple6->addItem( "MdcTrkR", m_mdcTrk_r );
256 status = m_tuple6->addItem( "MdcTrkDrho", m_mdcTrk_dr );
257 status = m_tuple6->addItem( "Rxy", m_rxy );
258 status = m_tuple6->addItem( "MdcKalTrkZ", m_mdcKalTrk_z );
259 }
260 else
261 {
262 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple6 ) << endmsg;
263 return StatusCode::FAILURE;
264 }
265 }
266
267 NTuplePtr nt7( ntupleSvc(), "FILE1/PullG" );
268 if ( nt7 ) m_tuple7 = nt7;
269 else
270 {
271 m_tuple7 = ntupleSvc()->book( "FILE1/PullG", CLID_ColumnWiseTuple, "Pull" );
272 if ( m_tuple7 )
273 {
274 status = m_tuple7->addItem( "gpull_drho", m_gpull_drho );
275 status = m_tuple7->addItem( "gpull_phi", m_gpull_phi );
276 status = m_tuple7->addItem( "gpull_kapha", m_gpull_kapha );
277 status = m_tuple7->addItem( "gpull_dz", m_gpull_dz );
278 status = m_tuple7->addItem( "gpull_lamb", m_gpull_lamb );
279 }
280 else
281 {
282 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple7 ) << endmsg;
283 return StatusCode::FAILURE;
284 }
285 }
286
287 log << MSG::INFO << "successfully return from initialize()" << endmsg;
288 return StatusCode::SUCCESS;
289}
INTupleSvc * ntupleSvc()

The documentation for this class was generated from the following files: