BOSS 6.6.4.p03
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 44 of file BeamParams.cxx.

44 :
45 Algorithm (name, pSvcLocator)
46{
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}

Member Function Documentation

◆ execute()

StatusCode BeamParams::execute ( )

Definition at line 262 of file BeamParams.cxx.

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

◆ finalize()

StatusCode BeamParams::finalize ( )

Definition at line 603 of file BeamParams.cxx.

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

◆ initialize()

StatusCode BeamParams::initialize ( )

Definition at line 88 of file BeamParams.cxx.

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

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