BOSS 7.1.2
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 12 of file BeamParams.h.

Constructor & Destructor Documentation

◆ BeamParams()

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

Definition at line 45 of file BeamParams.cxx.

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

Member Function Documentation

◆ execute()

StatusCode BeamParams::execute ( )

Definition at line 263 of file BeamParams.cxx.

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

◆ finalize()

StatusCode BeamParams::finalize ( )

Definition at line 604 of file BeamParams.cxx.

605{
606 MsgStream log(msgSvc(), name());
607 log << MSG::INFO << "in finalize()" << endmsg;
608
609 /*TF1 *func = new TF1("func", "gaus", -0.6, 0.6);
610 m_vertex_x->Fit("func", "NR+");
611 Double_t MeanX = func->GetParameter(1);
612 Double_t SigmaX = func->GetParameter(2);
613 Double_t ErrMeanX = func->GetParError(1);
614 Double_t ErrSigmaX = func->GetParError(2);
615
616 TF1 *funcY = new TF1("funcY", "gaus", -0.6, 0.1);
617 m_vertex_y->Fit("funcY", "NR+");
618 Double_t MeanY = funcY->GetParameter(1);
619 Double_t SigmaY = funcY->GetParameter(2);
620 Double_t ErrMeanY = funcY->GetParError(1);
621 Double_t ErrSigmaY = funcY->GetParError(2);
622
623 TF1 *funcZ = new TF1("funcZ", "gaus", -6, 6);
624 m_vertex_z->Fit("funcZ", "NR+");
625 Double_t MeanZ = funcZ->GetParameter(1);
626 Double_t SigmaZ = funcZ->GetParameter(2);
627 Double_t ErrMeanZ = funcZ->GetParError(1);
628 Double_t ErrSigmaZ = funcZ->GetParError(2);
629
630 TCanvas *myC = new TCanvas("myC", "myC", 1200, 400);
631 TPad *background = (TPad*)gPad;
632 background->SetFillColor(10);
633 myC->Divide(3,1);
634 myC->cd(1);
635
636 m_vertex_x_kal->Fit("func", "R+");
637 Double_t MeanXKal = func->GetParameter(1);
638 Double_t SigmaXKal = func->GetParameter(2);
639 Double_t ErrMeanXKal = func->GetParError(1);
640 Double_t ErrSigmaXKal = func->GetParError(2);
641
642 myC->cd(2);
643 m_vertex_y_kal->Fit("funcY", "R+");
644 Double_t MeanYKal = funcY->GetParameter(1);
645 Double_t SigmaYKal = funcY->GetParameter(2);
646 Double_t ErrMeanYKal = funcY->GetParError(1);
647 Double_t ErrSigmaYKal = funcY->GetParError(2);
648
649 myC->cd(3);
650 m_vertex_z_kal->Fit("funcZ", "R+");
651 Double_t MeanZKal = funcZ->GetParameter(1);
652 Double_t SigmaZKal = funcZ->GetParameter(2);
653 Double_t ErrMeanZKal = funcZ->GetParError(1);
654 Double_t ErrSigmaZKal = funcZ->GetParError(2);
655
656 cout << "Kal: Mean X, Y, Z = "<<MeanXKal << " " << MeanYKal << " " << MeanZKal << endl;
657 cout << "Kal: Sigma X, Y, Z = "<<SigmaXKal<<" " <<SigmaYKal <<" " << SigmaZKal << endl;
658 cout << "Gl: Mean X, Y, Z = " << MeanX << " " << MeanY << " " << MeanZ << endl;
659 cout << "Gl: Sigma X, Y, Z = " << SigmaX << " " << SigmaY << " " << SigmaZ << endl;
660
661 //////////////////////////////////////////////////////////////////
662 // Output a TXT file and a .ps figure for storing the fit results
663 //////////////////////////////////////////////////////////////////
664 const char *file_name, *figs_name;
665 if (m_hadronFile == 0) {
666 file_name = m_fileNameDst.c_str();
667 } else if (m_hadronFile == 1) {
668 file_name = m_fileNameHadron.c_str();
669 }
670
671 figs_name = m_figsName.c_str();
672 myC->SaveAs(figs_name);
673
674 ofstream file(file_name);
675
676 file << getenv("BES_RELEASE") << " " << m_parVer << " " << m_runNo << endl;
677 file << MeanX << " " << MeanY << " " << MeanZ << " "<< SigmaX << " "<< SigmaY << " " << SigmaZ << endl;
678 file << ErrMeanX << " " << ErrMeanY<< " " << ErrMeanZ << " " << ErrSigmaX << " " << ErrSigmaY << " " << ErrSigmaZ << endl;
679 file << MeanXKal << " "<< MeanYKal << " "<< MeanZKal << " "<< SigmaXKal << " " << SigmaYKal << " " << SigmaZKal << endl;
680 file << ErrMeanXKal << " " << ErrMeanYKal<< " " << ErrMeanZKal << " " << ErrSigmaXKal << " " << ErrSigmaYKal << " " << ErrSigmaZKal << endl;*/
681
682 /*delete m_vertex_x;
683 delete m_vertex_y;
684 delete m_vertex_z;
685 delete m_vertex_x_kal;
686 delete m_vertex_y_kal;
687 delete m_vertex_z_kal;
688 */
689 ////////////////////////////////////////////////
690 log << MSG::ALWAYS << "==================================" << endreq;
691 log << MSG::ALWAYS << "survived event :";
692 for(int i = 0; i < 10; i++){
693 log << MSG::ALWAYS << m_sel_number[i] << " ";
694 }
695 log << MSG::ALWAYS << endreq;
696 log << MSG::ALWAYS << "==================================" << endreq;
697 return StatusCode::SUCCESS;
698}

◆ initialize()

StatusCode BeamParams::initialize ( )

Definition at line 89 of file BeamParams.cxx.

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

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