BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
MdcDedxRecon Class Reference

#include <MdcDedxRecon.h>

+ Inheritance diagram for MdcDedxRecon:

Public Member Functions

 MdcDedxRecon (const std::string &name, ISvcLocator *pSvcLocator)
 
 ~MdcDedxRecon ()
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
StatusCode beginRun ()
 
const std::vector< MdcDedxTrk > & tracks (void) const
 
void clearTables ()
 
void mdctrackrec (RecMdcTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
 
void kaltrackrec (RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
 
void switchtomdctrack (int trkid, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
 

Detailed Description

Definition at line 29 of file MdcDedxRecon.h.

Constructor & Destructor Documentation

◆ MdcDedxRecon()

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

Definition at line 53 of file MdcDedxRecon.cxx.

53 : Algorithm(name, pSvcLocator)
54{
55 ex_calib=0; // get MdcDedxCorrection
56 calib_flag = 0x7F; //all calib on
57 landau = 1; //0: gauss distribution; 1:landau distribution;
58 ntpFlag = 1;
59 doNewGlobal = 0;
60 recalg = 0; //reconstruction method: 0:RecMdcTrack; 1:RecMdcKalTrack;
61 //2:use RecMdcTrack when no RecMdcKalTrack
62 ParticleFlag = -1; //e:0, mu:1, pi:2, K:3, p:4, cosmic:5
63 m_alg = 1; //calculation method of dE/dx of one track; 1: truncation and get average;
64 truncate = 0.70; // truncation rate
65
66 // Declare the properties
67 declareProperty("CalibFlag",calib_flag);
68 declareProperty("NTupleFlag",ntpFlag);
69 declareProperty("RecAlg",recalg);
70 declareProperty("ParticleType",ParticleFlag);
71 declareProperty("dEdxMethod",m_alg);
72 declareProperty("TruncRate",truncate);
73 declareProperty("RootFile",m_rootfile = std::string("no rootfile"));
74}

◆ ~MdcDedxRecon()

MdcDedxRecon::~MdcDedxRecon ( )
inline

Definition at line 33 of file MdcDedxRecon.h.

33{};

Member Function Documentation

◆ beginRun()

StatusCode MdcDedxRecon::beginRun ( )

Definition at line 239 of file MdcDedxRecon.cxx.

240{
241 MsgStream log(msgSvc(), name());
242 log << MSG::DEBUG <<"in MdcDedxRecon::beginrun()!!!"<< endreq;
243
244 StatusCode gesc = service("MdcGeomSvc", geosvc);
245 if (gesc == StatusCode::SUCCESS)
246 {
247 log << MSG::INFO <<"MdcGeomSvc is running"<<endreq;
248 }
249 else
250 {
251 log << MSG::ERROR <<"MdcGeomSvc is failed"<<endreq;
252 return StatusCode::SUCCESS;
253 }
254
255 return StatusCode::SUCCESS;
256}
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13

◆ clearTables()

void MdcDedxRecon::clearTables ( )

Definition at line 1119 of file MdcDedxRecon.cxx.

1119{}

◆ execute()

StatusCode MdcDedxRecon::execute ( )

Definition at line 288 of file MdcDedxRecon.cxx.

289{
290 MsgStream log(msgSvc(), name());
291 log << MSG::INFO << "in execute()" << endreq;
292
293 vector<double> Curve_Para, Sigma_Para;
294 int vFlag[3];//vFlag[0]:dedx curve version; vFlag[1]:dedx sigma version; vFlag[2]: 0:data; 1:MC
295
296 for(int i=0; i< dedxcursvc->getCurveSize(); i++) // get the dedx curve parameters
297 {
298 if(i==0) vFlag[0] = (int) dedxcursvc->getCurve(i); //first element is dedx curve version
299 else Curve_Para.push_back( dedxcursvc->getCurve(i) ); //dedx curve parameters
300 }
301 for(int i=0; i< dedxcursvc->getSigmaSize(); i++)
302 {
303 if(i==0) vFlag[1] = (int) dedxcursvc->getSigma(i); //dedx sigma version: 2: psip; 3:jpsi
304 else Sigma_Para.push_back( dedxcursvc->getSigma(i) ); //dedx sigma parameters
305 }
306
307 //check if size of parameters is right
308 if(vFlag[0] ==1 && Curve_Para.size() != 5) //version 1: 5 parameters 652a psip data
309 cout<<" size of dedx curve paramters for version 1 is not right!"<<endl;
310 if(vFlag[0] ==2 && Curve_Para.size() != 16) //version 2: increase from 5 parameters of 652 to 16
311 cout<<" size of dedx curve paramters for version 2 is not right!"<<endl;
312
313 if(vFlag[1] ==1 && Sigma_Para.size() != 14) //vesion 1: 14 parameters 652a psip data (old way)
314 cout<<" size of dedx sigma paramters for version 1 is not right!"<<endl;
315 if(vFlag[1] ==2 && Sigma_Para.size() != 21) //version 2: include t0 correction (for psip09 data)
316 cout<<" size of dedx sigma paramters for version 2 is not right!"<<endl;
317 if(vFlag[1] ==3 && Sigma_Para.size() != 18) //version 3: no t0 correction (for jpsi09 data and beyond)
318 cout<<" size of dedx sigma paramters for version 3 is not right!"<<endl;
319 if(vFlag[1] ==4 && Sigma_Para.size() != 19) //version 4: data with mdc track defaulted when kal track not ok(no t0 correction)
320 cout<<" size of dedx sigma paramters for version 4 is not right!"<<endl;
321 if(vFlag[1] ==5 && Sigma_Para.size() != 22) //version 5: data with mdc track defaulted when kal track not ok(with t0 correction)
322 cout<<" size of dedx sigma paramters for version 5 is not right!"<<endl;
323
324
325 //---------- register RecMdcDedxCol and RecMdcDedxHitCol to tds---------//
326 DataObject *aReconEvent;
327 eventSvc()->findObject("/Event/Recon",aReconEvent);
328 if(aReconEvent==NULL)
329 {
330 aReconEvent = new ReconEvent();
331 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
332 if(sc!=StatusCode::SUCCESS)
333 {
334 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
335 return( StatusCode::FAILURE);
336 }
337 }
338
339 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
340
341 DataObject *aDedxcol;
342 eventSvc()->findObject("/Event/Recon/RecMdcDedxCol",aDedxcol);
343 if(aDedxcol != NULL)
344 {
345 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxCol");
346 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxCol");
347 }
348 dedxList = new RecMdcDedxCol;
349 StatusCode dedxsc;
350 dedxsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxCol, dedxList);
351 if(!dedxsc.isSuccess())
352 {
353 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
354 return ( StatusCode::FAILURE);
355 }
356
357 DataObject *aDedxhitcol;
358 eventSvc()->findObject("/Event/Recon/RecMdcDedxHitCol",aDedxhitcol);
359 if(aDedxhitcol != NULL)
360 {
361 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxHitCol");
362 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxHitCol");
363 }
364 dedxhitList = new RecMdcDedxHitCol;
365 StatusCode dedxhitsc;
366 dedxhitsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxHitCol, dedxhitList);
367 if(!dedxhitsc.isSuccess())
368 {
369 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
370 return ( StatusCode::FAILURE);
371 }
372
373 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
374 if (!eventHeader)
375 {
376 log << MSG::INFO << "Could not find Event Header" << endreq;
377 return( StatusCode::FAILURE);
378 }
379 int eventNO = eventHeader->eventNumber();
380 int runNO = eventHeader->runNumber();
381 log << MSG::INFO << "now begin the event: runNO= "<<runNO<<" eventNO= "<< eventNO<< endreq;
382 RunNO2= runNO;
383 if(RunNO1==0) RunNO1=runNO;
384 if(RunNO2 != RunNO1)
385 {
386 cout<<"RunNO2 = "<<RunNO2 <<" RunNO1 = "<<RunNO1<<endl;
387 }
388 RunNO1 = runNO;
389 int runFlag; //data type flag
390 if(runNO<MdcDedxParam::RUN0) runFlag =0; //MC: less than 0
391 else if(runNO>=MdcDedxParam::RUN1 && runNO<MdcDedxParam::RUN2) runFlag =1;
392 else if(runNO>=MdcDedxParam::RUN2 && runNO<MdcDedxParam::RUN3) runFlag =2;
393 else if(runNO>=MdcDedxParam::RUN3 && runNO<MdcDedxParam::RUN4) runFlag =3;
394 else runFlag =4;
395
396 //vFlag[2] identify MC or data: 1:Mc; 0:data
397 if(runNO < 0) vFlag[2]=1;
398 else vFlag[2]=0;
399
400 double tes = -999.0;
401 int esTimeflag = -1;
402 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
403 if( ! estimeCol)
404 {
405 log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
406 }
407 else
408 {
409 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
410 for(; iter_evt!=estimeCol->end(); iter_evt++)
411 {
412 tes = (*iter_evt)->getTest();
413 esTimeflag = (*iter_evt)->getStat();
414 }
415 }
416 //cout<<"estime : "<<tes<<endl;
417
418
419 Identifier mdcid;
420 int ntrk;
421 ex_trks.clear(); // clear the vector of MdcDedxTrk,when read a new event
422 m_trkalgs.clear();
423 if( !doNewGlobal )
424 {
425 log << MSG::DEBUG << "recalg: "<<recalg<<endreq;
426
427 //---------using kal algorithm by default, switch to seek mdc track if no kaltrack hits //
428 if(recalg ==2 )
429 {
430 //retrieve MdcKalTrackCol from TDS
431 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
432 if (!newtrkCol)
433 {
434 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
435 return StatusCode::SUCCESS;
436 }
437 if(ntpFlag>0) eventNo++;
438 ntrk = newtrkCol->size();
439 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
440
441 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
442 //cout<<"MdcDedxRecon newtrkCol.size() "<<newtrkCol->size()<<endl;
443 for( ; trk != newtrkCol->end(); trk++)
444 {
445 if(ntpFlag>0) trackNO1++;
446
447 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
448 //if not set ParticleFlag, we will get the last succefully hypothesis;
449 if(gothits.size()==0)
450 {
451 m_trkalg=0;
452 int trkid=(*trk)->trackId();
453 switchtomdctrack(trkid, mdcid, tes, runNO, eventNO, runFlag, log);
454 }
455 else
456 {
457 m_trkalg=1;
458 if(gothits.size()<2) continue;
459 kaltrackrec(trk, mdcid, tes, runNO, eventNO, runFlag, log );
460 }
461 }//end of track loop
462 }//end of recalg==2
463
464 //------------------------ Use information of MDC track rec --------------------------//
465 else if(recalg ==0 )
466 {
467 m_trkalg=0;
468
469 //retrieve MdcTrackCol from TDS
470 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
471 if (!newtrkCol)
472 {
473 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
474 return StatusCode::SUCCESS;
475 }
476 if(ntpFlag>0) eventNo++;
477 ntrk = newtrkCol->size();
478 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
479
480 vector<double> phlist; //dE/dx only after hit level correction
481 vector<double> phlist_hit; //dE/dx after hit and track level correction
482 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
483 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
484 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
485 int Nhits=0;
486 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
487
488 RecMdcTrackCol::iterator trk = newtrkCol->begin();
489 for( ; trk != newtrkCol->end(); trk++)
490 {
491 if(ntpFlag>0) trackNO1++;
492
493 MdcDedxTrk trk_ex( **trk);
494 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
495 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
496
497 HepVector a = (*trk)->helix();
498 HepSymMatrix tkErrM = (*trk)->err();
499 m_d0E = tkErrM[0][0];
500 m_phi0E = tkErrM[1][1];
501 m_cpaE = tkErrM[2][2];
502 m_z0E = tkErrM[3][3];
503
504 HepPoint3D IP(0,0,0);
505 Dedx_Helix exhel(IP, a); //initialize class Dedx_Helix for DedxCorrecSvc
506 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
507 //cout<<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl;
508 phi0= a[1];
509 costheta = cos(M_PI_2-atan(a[4]));
510 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack costheta: "<<trk_ex.theta()<<endl;
511 sintheta = sin(M_PI_2-atan(a[4]));
512 m_Pt = 1.0/fabs( a[2] );
513 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
514 charge = ( a[2] > 0 )? 1 : -1;
515 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
516 dedxhitrefvec = new DedxHitRefVec;
517
518
519 HitRefVec gothits= (*trk)->getVecHits();
520 Nhits = gothits.size();
521 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
522 if(gothits.size()<2)
523 {
524 delete dedxhitrefvec;
525 continue;
526 }
527 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
528
529 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
530 HitRefVec::iterator it_gothit = gothits.begin();
531 for( ; it_gothit != gothits.end(); it_gothit++)
532 {
533 mdcid = (*it_gothit)->getMdcId();
534 layid = MdcID::layer(mdcid);
535 localwid = MdcID::wire(mdcid);
536 w0id = geosvc->Layer(layid)->Wirst();
537 wid = w0id + localwid;
538 adc_raw = (*it_gothit)->getAdc();
539 tdc_raw = (*it_gothit)->getTdc();
540 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
541 zhit = (*it_gothit)->getZhit();
542 lr = (*it_gothit)->getFlagLR();
543 if(lr == 2) cout<<"lr = "<<lr<<endl;
544 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
545 else driftD = (*it_gothit)->getDriftDistRight();
546 //cout<<"lr: "<<lr<<" driftD: "<<driftD<<endl;
547 driftd = abs(driftD);
548 dd = (*it_gothit)->getDoca();
549 if(lr == 0 || lr == 2 ) dd = -abs(dd);
550 if(lr == 1 ) dd = abs(dd);
551 driftT = (*it_gothit)->getDriftT();
552
553 eangle = (*it_gothit)->getEntra();
554 eangle = eangle/M_PI;
555 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
556 rphi_path=pathlength * sintheta;
557
558 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd: "<<dd<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
559 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
560 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
561 //if(pathlength == -1)
562 //cout<<"parameter0: "<<"eventNO: "<<eventNO<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
563
564 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
565 if(runNO<0)
566 {
567 if (layid<8)
568 {
569 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
570 }
571 else if(layid >= 8)
572 {
573 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
574 }
575 }
576 else if(runNO>=0)
577 {
578 if(layid <8)
579 {
580 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
581 }
582 else if(layid >= 8)
583 {
584 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
585 }
586 }
587 //cout<<"recAlg=0 para mdc: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
588
589 if (ph<0.0||ph==0) continue;
590 else
591 {
592 //-----------------------put data into TDS-----------------------------//
593 dedxhit = new RecMdcDedxHit;
594 dedxhit->setMdcHit(*it_gothit);
595 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg );
596 dedxhit->setMdcId(mdcid);
597 dedxhit->setFlagLR(lr);
598 dedxhit->setTrkId(trk_ex.trk_id());
599 dedxhit->setDedx(ph_hit);
600 dedxhit->setPathLength(pathlength);
601
602 //---------------------------Fill root file----------------------------//
603 if(m_rootfile!= "no rootfile")
604 {
605 m_hitNo_wire->Fill(wid);
606 }
607
608 //-------------------------Fill ntuple n102---------------------------//
609 if ( ntpFlag ==2 )
610 {
611 m_charge1 = (*trk)->charge();
612 m_t01 = tes;
613 m_driftT = driftT;
614 m_tdc_raw = tdc_raw;
615 m_phraw = adc_raw;
616 m_exraw = ph_hit;
617 m_localwid = localwid;
618 m_wire = wid;
619 m_runNO1 = runNO;
620 m_evtNO1 = eventNO;
621 m_nhit_hit = Nhits;
622 m_doca_in = dd;
623 m_doca_ex = dd;
624 m_driftdist = driftD;
625 m_eangle = eangle;
626 m_costheta1 = costheta;
627 m_pathL = pathlength;
628 m_layer = layid;
629 m_ptrk1 = m_P;
630 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
631 m_trackId1 = trk_ex.trk_id();
632 StatusCode status2 = m_nt2->write();
633 if ( status2.isFailure() )
634 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
635 }
636 if(layid>3)
637 {
638 phlist.push_back(ph);
639 phlist_hit.push_back(ph_hit);
640 }
641 }
642 dedxhitList->push_back( dedxhit );
643 dedxhitrefvec->push_back( dedxhit );
644 }//end of hit loop
645 trk_ex.set_phlist( phlist );
646 trk_ex.set_phlist_hit( phlist_hit );
647 trk_ex.setVecDedxHits( *dedxhitrefvec );
648 ex_trks.push_back(trk_ex );
649 m_trkalgs.push_back(m_trkalg);
650
651 delete dedxhitrefvec;
652 phlist.clear();
653 phlist_hit.clear();
654 if(ntpFlag>0) trackNO2++;
655 }//end of track loop
656 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
657 }//end of recalg==0
658
659 //------------------------ Use information of MDC KAL track rec -----------------------//
660 else if(recalg ==1 )
661 {
662 m_trkalg=1;
663
664 //retrieve MdcKalTrackCol from TDS
665 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
666 if (!newtrkCol)
667 {
668 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
669 return StatusCode::SUCCESS;
670 }
671 if(ntpFlag>0) eventNo++;
672 ntrk = newtrkCol->size();
673 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
674
675 vector<double> phlist; //dE/dx only after hit level correction
676 vector<double> phlist_hit; //dE/dx after hit and track level correction
677 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
678 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
679 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
680 int Nhits=0;
681 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
682
683
684 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
685 for( ; trk != newtrkCol->end(); trk++)
686 {
687 if(ntpFlag>0) trackNO1++;
688
689 MdcDedxTrk trk_ex( **trk, ParticleFlag);
690 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
691 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
692
693 HepVector a;
694 HepSymMatrix tkErrM;
695 if(ParticleFlag>-1&&ParticleFlag<5)
696 {
697 DstMdcKalTrack* dstTrk = *trk;
698 a = dstTrk->getFHelix(ParticleFlag);
699 tkErrM = dstTrk->getFError(ParticleFlag);
700 }
701 else
702 {
703 a = (*trk)->getFHelix();
704 tkErrM = (*trk)->getFError();
705 }
706 //cout<<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl; //getFHelix is first layer of MdcKalTrack;
707 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
708
709 m_d0E = tkErrM[0][0];
710 m_phi0E = tkErrM[1][1];
711 m_cpaE = tkErrM[2][2];
712 m_z0E = tkErrM[3][3];
713
714 phi0= a[1];
715 costheta = cos(M_PI_2-atan(a[4]));
716 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl; //theta in trk_ex is got by getFHelix();
717 sintheta = sin(M_PI_2-atan(a[4]));
718 m_Pt = 1.0/fabs( a[2] );
719 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
720 //cout<<"track ptrk: "<<m_P<<" extrack ptrk: "<<trk_ex.P()<<endl;
721 charge = ( a[2] > 0 )? 1 : -1;
722 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
723 dedxhitrefvec = new DedxHitRefVec;
724
725
726 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
727 //HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
728 //if not set ParticleFlag, we will get the last succefully hypothesis;
729 Nhits = gothits.size();
730 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
731 if(gothits.size()<2)
732 {
733 delete dedxhitrefvec;
734 continue;
735 }
736 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
737
738 RecMdcHit* mdcHit = 0;
739 HelixSegRefVec::iterator it_gothit = gothits.begin();
740 for( ; it_gothit != gothits.end(); it_gothit++)
741 {
742 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
743 //HepVector a_hit_ex = (*it_gothit)->getHelixExcl(); //same with getHelixIncl()
744 HepPoint3D IP(0,0,0);
745 Dedx_Helix exhel_hit_in(IP, a_hit_in);
746 //Dedx_Helix exhel_hit_ex(IP, a_hit_ex);
747 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
748 //cout<<"getHelixIncl 5 para: "<<a_hit_in[0]<<" "<<a_hit_in[1] <<" "<<a_hit_in[2]<<" "<<a_hit_in[3]<<" "<<a_hit_in[4]<<endl;
749
750 mdcid = (*it_gothit)->getMdcId();
751 layid = MdcID::layer(mdcid);
752 localwid = MdcID::wire(mdcid);
753 w0id = geosvc->Layer(layid)->Wirst();
754 wid = w0id + localwid;
755 adc_raw = (*it_gothit)->getAdc();
756 tdc_raw = (*it_gothit)->getTdc();
757 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
758 zhit = (*it_gothit)->getZhit();
759 lr = (*it_gothit)->getFlagLR();
760 if(lr == 2) cout<<"lr = "<<lr<<endl;
761 driftD = (*it_gothit)->getDD();
762 driftd = abs(driftD);
763 driftT = (*it_gothit)->getDT();
764 dd_in = (*it_gothit)->getDocaIncl(); //getDocaIncl() include fit unused hit
765 dd_ex = (*it_gothit)->getDocaExcl(); //getDocaExcl() exclude fit unused hit
766
767 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
768 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
769 //dd = dd/doca_norm[layid];
770 //cout<<"lr "<<lr<<" dd_in: "<<dd_in<<" dd_ex: "<<dd_ex<<endl;
771
772 eangle = (*it_gothit)->getEntra();
773 eangle = eangle/M_PI;
774 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
775 rphi_path=pathlength * sintheta;
776 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd_in: "<<dd_in<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
777 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
778 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
779 //if(pathlength == -1)
780 //cout<<"parameter1: "<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
781
782 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
783 if(runNO<0)
784 {
785 if (layid<8)
786 {
787 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
788 }
789 else if(layid >= 8)
790 {
791 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
792 }
793 }
794 else if(runNO>=0)
795 {
796 if(layid <8)
797 {
798 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
799 }
800 else if(layid >= 8)
801 {
802 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
803 }
804 }
805 //cout<<"recAlg=1 para kal: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd_in<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
806
807 if (ph<0.0||ph==0) continue;
808 else
809 {
810 //-----------------------put data into TDS-----------------------------//
811 dedxhit = new RecMdcDedxHit;
812 dedxhit->setMdcHit(mdcHit);
813 dedxhit->setMdcKalHelixSeg(*it_gothit);
814 dedxhit->setMdcId(mdcid);
815 dedxhit->setFlagLR(lr);
816 dedxhit->setTrkId(trk_ex.trk_id());
817 dedxhit->setDedx(ph_hit);
818 dedxhit->setPathLength(pathlength);
819
820 //---------------------------Fill root file----------------------------//
821 if(m_rootfile!= "no rootfile")
822 {
823 m_hitNo_wire->Fill(wid);
824 }
825
826 //-------------------------Fill ntuple n102---------------------------//
827 if ( ntpFlag ==2 )
828 {
829 m_charge1 = (*trk)->charge();
830 m_t01 = tes;
831 m_driftT = driftT;
832 m_tdc_raw = tdc_raw;
833 m_phraw = adc_raw;
834 m_exraw = ph_hit;
835 m_localwid = localwid;
836 m_wire = wid;
837 m_runNO1 = runNO;
838 m_evtNO1 = eventNO;
839 m_nhit_hit = Nhits;
840 m_doca_in = dd_in;
841 m_doca_ex = dd_ex;
842 m_driftdist = driftD;
843 m_eangle = eangle;
844 m_costheta1 = costheta;
845 m_pathL = pathlength;
846 m_layer = layid;
847 m_ptrk1 = m_P;
848 m_ptrk_hit = p_hit;
849 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
850 m_trackId1 = trk_ex.trk_id();
851 StatusCode status2 = m_nt2->write();
852 if ( status2.isFailure() )
853 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
854 }
855 if(layid>3)
856 {
857 phlist.push_back(ph);
858 phlist_hit.push_back(ph_hit);
859 }
860 }
861 dedxhitList->push_back( dedxhit );
862 dedxhitrefvec->push_back( dedxhit );
863 }//end of hit loop
864 trk_ex.set_phlist( phlist );
865 trk_ex.set_phlist_hit( phlist_hit );
866 trk_ex.setVecDedxHits( *dedxhitrefvec );
867 ex_trks.push_back(trk_ex );
868 m_trkalgs.push_back(m_trkalg);
869
870 delete dedxhitrefvec;
871 phlist.clear();
872 phlist_hit.clear();
873 if(ntpFlag>0) trackNO3++;
874 }//end of track loop
875 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
876 }//end of recalg==1
877 //--------------------------------- Hit level finished ---------------------------------//
878
879 //-------------------------------- Track level begin --------------------------------//
880 if( ntrk != ex_trks.size())
881 log << MSG::DEBUG <<"ntrkCol size: "<<ntrk<<" dedxtrk size: "<<ex_trks.size()<< endreq;
882
883 double poca_x=0,poca_y=0,poca_z=0;
884 float sintheta=0,costheta=0,ptrk=0,charge=0;
885 int Nhit=0,Nphlisthit=0;
886 int usedhit = 0;
887 double phtm=0,median=0,geometric=0,geometric_trunc=0,harmonic=0,harmonic_trunc=0,transform=0,logtrunc=0;
888
889 enum pid_dedx parType_temp(electron);
890 float probpid_temp=-0.01,expect_temp=-0.01,sigma_temp=-0.01,chidedx_temp=-0.01;
891
892 double dEdx_meas_hit=0;
893 double dEdx_meas=0,dEdx_meas_esat=0,dEdx_meas_norun=0;
894 int trk_recalg=-1;
895
896 int idedxid = 0;
897 for(std::vector<MdcDedxTrk>::iterator it = ex_trks.begin(); it != ex_trks.end(); it++, idedxid++ )
898 {
899 log << MSG::DEBUG <<"-------------------------------------------------------"<< endreq;
900 log << MSG::DEBUG <<" trk id ="<< it->trk_id()<<" : P ="<<it->P() <<" Pt ="<<it->Pt()<<" : theta ="<<it->theta() <<" : phi ="<<it->phi()<< " : charge = "<<it->charge()<<endreq;
901 log << MSG::DEBUG <<"all hits on track: "<<it->nsample()<<" phlist size: "<<it->get_phlist().size()<<endreq;
902
903 if(it->trk_ptr_kal()!=0)
904 {
905 poca_x = it->trk_ptr_kal()->x(); //get poca, default pid is pion; change pid using setPidType();
906 poca_y = it->trk_ptr_kal()->y();
907 poca_z = it->trk_ptr_kal()->z();
908 }
909 else if(it->trk_ptr()!=0)
910 {
911 poca_x = it->trk_ptr()->x();
912 poca_y = it->trk_ptr()->y();
913 poca_z = it->trk_ptr()->z();
914 }
915 //cout<<"poca_x: "<<poca_x<<" poca_y: "<<poca_y<<" poca_z: "<<poca_z<<endl;
916
917 sintheta = sin(it->theta());
918 costheta = cos(it->theta());
919 ptrk = it->P();
920 charge = it->charge();
921 Nhit = it->nsample(); //total hits on track used as sample;
922 Nphlisthit = it->get_phlist_hit().size(); //hits in phlist_hit, exclude first 4 layers;
923 //usedhit: hits after truncting phlist and used in cal dE/dx value;
924
925 //--------------------------extrk truncation--------------------------------//
926 phtm = it->cal_dedx( truncate );
927 //cout<<"phtm: "<<phtm<<endl;
928 //median = it->cal_dedx_median( truncate );
929 //geometric = it->cal_dedx_geometric( truncate );
930 //geometric_trunc = it->cal_dedx_geometric_trunc( truncate );
931 //harmonic = it->cal_dedx_harmonic( truncate );
932 //harmonic_trunc = it->cal_dedx_harmonic_trunc( truncate );
933 //transform = it->cal_dedx_transform( 0 );
934 //logtrunc = it->cal_dedx_log( 1.0, 0);
935
936 if(m_alg==1) dEdx_meas_hit = it->cal_dedx_bitrunc(truncate, 0, usedhit);
937 else if(m_alg==2) dEdx_meas_hit = it-> cal_dedx_transform(1);
938 else if(m_alg==3) dEdx_meas_hit = it->cal_dedx_log(1.0, 1);
939 else cout<<"-------------Truncate Algorithm Error!!!------------------------"<<endl;
940 if(m_alg==1 && usedhit==0)
941 {
942 cout<<"***************bad extrk with no hits!!!!!******************"<<endl;
943 continue;
944 }
945 // force to use the same definition of usedhit in TDS and what used in setDedx() function
946 usedhit = (int)(Nphlisthit*truncate);
947
948 //--------------------track level correction of extrk dE/dx Value---------------//
949 dEdx_meas = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
950 dEdx_meas_esat = exsvc->StandardTrackCorrec(1,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
951 dEdx_meas_norun = exsvc->StandardTrackCorrec(2,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
952 log << MSG::INFO << "===================== " << endreq << endreq;
953 log << MSG::DEBUG <<"dEdx_meas_hit: "<<dEdx_meas_hit<<" dEdx_meas: "<<dEdx_meas<<" dEdx_meas_esat: "<<dEdx_meas_esat<<" dEdx_meas_norun: "<<dEdx_meas_norun<<" ptrk: "<<it->P()<<endreq;
954 log << MSG::INFO << "ptrk " << ptrk << " costhe " << costheta << " runno " << runNO << " evtno " << eventNO << endreq;
955 //cout<<"dEdx_meas_hit: "<<dEdx_meas_hit<<" dEdx_meas: "<<dEdx_meas<<" dEdx_meas_esat: "<<dEdx_meas_esat<<" dEdx_meas_norun: "<<dEdx_meas_norun<<" ptrk: "<<it->P()<<endl;
956
957 trk_recalg = m_trkalgs[idedxid];
958
959 if(dEdx_meas<0 || dEdx_meas==0) continue;
960 it->set_dEdx( landau , dEdx_meas, trk_recalg, runFlag, vFlag , tes, Curve_Para, Sigma_Para, ex_calib); // calculate expect
961 parType_temp = electron;
962 probpid_temp=-0.01;
963 expect_temp=-0.01;
964 sigma_temp=-0.01;
965 chidedx_temp=-0.01;
966 for(int i=0;i<5;i++)
967 {
968 if(it->pprob()[i]>probpid_temp)
969 {
970 parType_temp = (enum pid_dedx) i; //e:0, mu:1, pi:2, K:3, p:4
971 probpid_temp = it->pprob()[i];
972 expect_temp = it->pexpect()[i];
973 sigma_temp = it->pexp_sigma()[i];
974 chidedx_temp = it->pchi_dedx()[i];
975 }
976 }
977 log<< MSG::INFO<<"expect dE/dx: type: "<<parType_temp<<" probpid: "<<probpid_temp<<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq;
978 //cout<<"dEdx_meas: "<<dEdx_meas<<endl;
979
980 //-----------------------put data into TDS-----------------------------//
981 resdedx = new RecMdcDedx;
982 resdedx->setDedxHit(dEdx_meas_hit);
983 resdedx->setDedxEsat(dEdx_meas_esat);
984 resdedx->setDedxNoRun(dEdx_meas_norun);
985 resdedx->setDedxMoment(it->P());
986 resdedx->setTrackId( it->trk_id() );
987 resdedx->setProbPH( dEdx_meas );
988 resdedx->setNormPH( dEdx_meas/550.0 );
989 resdedx->setDedxExpect( it->pexpect() );
990 resdedx->setSigmaDedx( it->pexp_sigma() );
991 resdedx->setPidProb( it->pprob() );
992 resdedx->setChi( it->pchi_dedx() );
993 //cout<<"recdedx: "<<" "<<resdedx->getPidProb(parType_temp)<<" "<<resdedx->getDedxExpect(parType_temp)<<" "<<resdedx->getSigmaDedx(parType_temp)<<" "<<resdedx->chi(parType_temp)<<endl;
994 resdedx->setNumTotalHits(it->nsample() ); //all hits on track;
995 resdedx->setNumGoodHits(usedhit); //hits after truncing the phlist
996 //phlist are all hits on track excluding first 4 layers;
997 //resdedx->setStatus( it->quality() );
998 resdedx->setStatus(trk_recalg );
999 //cout<<"trk_recalg: "<<trk_recalg<<" trk stat: "<<it->quality()<<endl;
1000 resdedx->setTruncAlg( m_alg );
1001 resdedx->setParticleId(parType_temp);
1002 //cout<<"ParticleType: "<<parType_temp<<" "<<resdedx->particleType()<<endl;
1003 resdedx->setVecDedxHits(it->getVecDedxHits());
1004 resdedx->setMdcTrack(it->trk_ptr());
1005 resdedx->setMdcKalTrack(it->trk_ptr_kal());
1006
1007 //-------------------------Fill ntuple n103---------------------------//
1008 if(ntpFlag ==2)
1009 {
1010 m_phtm=phtm;
1011 //m_median=median;
1012 //m_geometric=geometric;
1013 //m_geometric_trunc=geometric_trunc;
1014 //m_harmonic=harmonic;
1015 //m_harmonic_trunc=harmonic_trunc;
1016 //m_transform=transform;
1017 //m_logtrunc=logtrunc;
1018 m_dEdx_meas = dEdx_meas;
1019
1020 m_poca_x = poca_x;
1021 m_poca_y = poca_y;
1022 m_poca_z = poca_z;
1023 m_ptrk=ptrk;
1024 m_sintheta=sintheta;
1025 m_costheta=costheta;
1026 m_charge=charge;
1027 m_runNO = runNO;
1028 m_evtNO = eventNO;
1029
1030 m_t0 = tes;
1031 m_trackId = it->trk_id();
1032 m_recalg = trk_recalg;
1033
1034 m_nhit=Nhit;
1035 m_nphlisthit=Nphlisthit;
1036 m_usedhit=usedhit;
1037 for(int i=0; i<Nphlisthit; i++) m_dEdx_hit[i] = it->get_phlist_hit()[i];
1038
1039 m_parttype = parType_temp;
1040 m_prob = probpid_temp;
1041 m_expect = expect_temp;
1042 m_sigma = sigma_temp;
1043 m_chidedx = chidedx_temp;
1044 m_chidedxe=it->pchi_dedx()[0];
1045 m_chidedxmu=it->pchi_dedx()[1];
1046 m_chidedxpi=it->pchi_dedx()[2];
1047 m_chidedxk=it->pchi_dedx()[3];
1048 m_chidedxp=it->pchi_dedx()[4];
1049 for(int i=0;i<5;i++)
1050 {
1051 m_probpid[i]= it->pprob()[i];
1052 m_expectid[i]= it->pexpect()[i];
1053 m_sigmaid[i]= it->pexp_sigma()[i];
1054 }
1055
1056 StatusCode status = m_nt1->write();
1057 if ( status.isFailure() )
1058 {
1059 log << MSG::ERROR << "Cannot fill Ntuple n103" << endreq;
1060 }
1061 }
1062
1063 log<< MSG::INFO<<"-----------------Summary of this ExTrack----------------"<<endreq
1064 <<"dEdx_mean: "<<dEdx_meas<<" type: "<<parType_temp<<" probpid: "<<probpid_temp
1065 <<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq<<endreq;
1066
1067 dedxList->push_back( resdedx );
1068 }//ExTrk loop end
1069 } //doNewGlobal==0 END . . .
1070
1071
1072 // check MdcDedxCol registered
1073 SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
1074 if (!newexCol)
1075 {
1076 log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
1077 return( StatusCode::SUCCESS);
1078 }
1079 log << MSG::DEBUG << "----------------Begin to check RecMdcDedxCol-----------------"<<endreq;
1080 RecMdcDedxCol::iterator it_dedx = newexCol->begin();
1081 for( ; it_dedx != newexCol->end(); it_dedx++)
1082 {
1083 log << MSG::INFO << "retrieved MDC dE/dx:" << endreq
1084 << "dEdx Id: " << (*it_dedx)->trackId()
1085 << " part Id: " << (*it_dedx)->particleType()
1086 << " measured dEdx: " << (*it_dedx)->probPH()
1087 << " dEdx std: " << (*it_dedx)->normPH() << endreq
1088 << "hits on track: "<<(*it_dedx)->numTotalHits()
1089 << " usedhits: " << (*it_dedx)->numGoodHits() << endreq
1090 << "Electron: expect: " << (*it_dedx)->getDedxExpect(0)
1091 << " sigma: " << (*it_dedx)->getSigmaDedx(0)
1092 << " chi: " << (*it_dedx)->chi(0)
1093 << " prob: " << (*it_dedx)->getPidProb(0) << endreq
1094 << "Muon: expect: " << (*it_dedx)->getDedxExpect(1)
1095 << " sigma: " << (*it_dedx)->getSigmaDedx(1)
1096 << " chi: " << (*it_dedx)->chi(1)
1097 << " prob: " << (*it_dedx)->getPidProb(1) << endreq
1098 << "Pion: expect: " << (*it_dedx)->getDedxExpect(2)
1099 << " sigma: " << (*it_dedx)->getSigmaDedx(2)
1100 << " chi: " << (*it_dedx)->chi(2)
1101 << " prob: " << (*it_dedx)->getPidProb(2) << endreq
1102 << "Kaon: expect: " << (*it_dedx)->getDedxExpect(3)
1103 << " sigma: " << (*it_dedx)->getSigmaDedx(3)
1104 << " chi: " << (*it_dedx)->chi(3)
1105 << " prob: " << (*it_dedx)->getPidProb(3) << endreq
1106 << "Proton: expect: " << (*it_dedx)->getDedxExpect(4)
1107 << " sigma: " << (*it_dedx)->getSigmaDedx(4)
1108 << " chi: " << (*it_dedx)->chi(4)
1109 << " prob: " << (*it_dedx)->getPidProb(4) << endreq;
1110 }
1111 return StatusCode::SUCCESS;
1112}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
Iner_DriftDistCut
#define Iner_DriftDistCut
Definition: DedxCalibParameters.h:58
RPhi_PathMaxCut
#define RPhi_PathMaxCut
Definition: DedxCalibParameters.h:55
Out_DriftDistCut
#define Out_DriftDistCut
Definition: DedxCalibParameters.h:59
MaxHistValue
#define MaxHistValue
Definition: DedxCalibParameters.h:47
pid_dedx
pid_dedx
Definition: DstMdcDedx.h:9
electron
@ electron
Definition: DstMdcDedx.h:9
trackNO3
int trackNO3
Definition: MdcDedxRecon.cxx:51
eventNo
int eventNo
Definition: MdcDedxRecon.cxx:48
RunNO1
int RunNO1
Definition: MdcDedxRecon.cxx:45
trackNO2
int trackNO2
Definition: MdcDedxRecon.cxx:50
RunNO2
int RunNO2
Definition: MdcDedxRecon.cxx:46
trackNO1
int trackNO1
Definition: MdcDedxRecon.cxx:49
RecMdcDedxHitCol
ObjectVector< RecMdcDedxHit > RecMdcDedxHitCol
Definition: RecMdcDedxHit.h:78
DedxHitRefVec
SmartRefVector< RecMdcDedxHit > DedxHitRefVec
Definition: RecMdcDedx.h:27
RecMdcDedxCol
ObjectVector< RecMdcDedx > RecMdcDedxCol
Definition: RecMdcDedx.h:132
HelixSegRefVec
SmartRefVector< RecMdcKalHelixSeg > HelixSegRefVec
Definition: RecMdcKalTrack.h:27
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:22
M_PI
#define M_PI
Definition: TConstant.h:4
NULL
#define NULL
Definition: VertexExtrapolate.h:14
Dedx_Helix
Helix parameter class.
Definition: Dedx_Helix.h:33
DstMdcDedx::setChi
void setChi(double *chi)
Definition: DstMdcDedx.h:77
DstMdcDedx::setTruncAlg
void setTruncAlg(int trunc_alg)
Definition: DstMdcDedx.h:75
DstMdcDedx::setStatus
void setStatus(int status)
Definition: DstMdcDedx.h:74
DstMdcDedx::setNumGoodHits
void setNumGoodHits(int numGoodHits)
Definition: DstMdcDedx.h:81
DstMdcDedx::setProbPH
void setProbPH(double probPH)
Definition: DstMdcDedx.h:83
DstMdcDedx::setNormPH
void setNormPH(double normPH)
Definition: DstMdcDedx.h:84
DstMdcDedx::setParticleId
void setParticleId(int particleId)
Definition: DstMdcDedx.h:73
DstMdcDedx::setTrackId
void setTrackId(int trackId)
Definition: DstMdcDedx.h:72
DstMdcDedx::setNumTotalHits
void setNumTotalHits(int numTotalHits)
Definition: DstMdcDedx.h:82
DstMdcKalTrack::getFError
const HepSymMatrix & getFError(const int pid) const
Definition: DstMdcKalTrack.h:213
DstMdcKalTrack::getFHelix
const HepVector & getFHelix(const int pid) const
Definition: DstMdcKalTrack.h:210
IDedxCorrecSvc::StandardHitCorrec
virtual double StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0
IDedxCorrecSvc::PathL
virtual double PathL(int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const =0
IDedxCorrecSvc::StandardTrackCorrec
virtual double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const =0
IDedxCurSvc::getSigmaSize
virtual const int getSigmaSize()=0
IDedxCurSvc::getCurveSize
virtual const int getCurveSize()=0
IDedxCurSvc::getCurve
virtual const double getCurve(int i)=0
IDedxCurSvc::getSigma
virtual const double getSigma(int i)=0
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
MdcDedxParam::RUN0
static long int RUN0
Definition: MdcDedxParam.h:36
MdcDedxParam::RUN2
static long int RUN2
Definition: MdcDedxParam.h:38
MdcDedxParam::RUN1
static long int RUN1
Definition: MdcDedxParam.h:37
MdcDedxParam::RUN4
static long int RUN4
Definition: MdcDedxParam.h:40
MdcDedxParam::RUN3
static long int RUN3
Definition: MdcDedxParam.h:39
MdcDedxRecon::switchtomdctrack
void switchtomdctrack(int trkid, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1487
MdcDedxRecon::kaltrackrec
void kaltrackrec(RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1291
MdcGeoLayer::Wirst
int Wirst(void) const
Definition: MdcGeoLayer.h:157
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
RecMdcDedxHit::setMdcId
void setMdcId(Identifier mdcid)
Definition: RecMdcDedxHit.h:64
RecMdcDedxHit::setDedx
void setDedx(double dedx)
Definition: RecMdcDedxHit.h:62
RecMdcDedxHit::setFlagLR
void setFlagLR(int lr)
Definition: RecMdcDedxHit.h:61
RecMdcDedxHit::setMdcKalHelixSeg
void setMdcKalHelixSeg(const RecMdcKalHelixSeg *mdcKalHelixSeg)
Definition: RecMdcDedxHit.h:57
RecMdcDedxHit::setTrkId
void setTrkId(int trkid)
Definition: RecMdcDedxHit.h:60
RecMdcDedxHit::setMdcHit
void setMdcHit(const RecMdcHit *mdcHit)
Definition: RecMdcDedxHit.h:58
RecMdcDedxHit::setPathLength
void setPathLength(double pathlength)
Definition: RecMdcDedxHit.h:63
RecMdcDedx::setMdcTrack
void setMdcTrack(RecMdcTrack *trk)
Definition: RecMdcDedx.h:107
RecMdcDedx::setVecDedxHits
void setVecDedxHits(const DedxHitRefVec &vecdedxhit)
Definition: RecMdcDedx.h:76
RecMdcDedx::setDedxNoRun
void setDedxNoRun(double dedx_norun)
Definition: RecMdcDedx.h:80
RecMdcDedx::setMdcKalTrack
void setMdcKalTrack(RecMdcKalTrack *trk)
Definition: RecMdcDedx.h:108
RecMdcDedx::setDedxMoment
void setDedxMoment(double dedx_momentum)
Definition: RecMdcDedx.h:81
RecMdcDedx::setDedxExpect
void setDedxExpect(double *dedx_exp)
Definition: RecMdcDedx.h:90
RecMdcDedx::setSigmaDedx
void setSigmaDedx(double *sigma_dedx)
Definition: RecMdcDedx.h:94
RecMdcDedx::setDedxEsat
void setDedxEsat(double dedx_esat)
Definition: RecMdcDedx.h:79
RecMdcDedx::setDedxHit
void setDedxHit(double dedx_hit)
Definition: RecMdcDedx.h:78
RecMdcDedx::setPidProb
void setPidProb(double *pid_prob)
Definition: RecMdcDedx.h:98
EventModel::Recon::RecMdcDedxCol
_EXTERN_ std::string RecMdcDedxCol
Definition: EventModel.h:88
EventModel::Recon::RecMdcDedxHitCol
_EXTERN_ std::string RecMdcDedxHitCol
Definition: EventModel.h:89
costheta
float costheta
Definition: prepare_data.cxx:28
charge
float charge
Definition: prepare_data.cxx:28
dEdx_meas
float dEdx_meas
Definition: prepare_data.cxx:28
ptrk
float ptrk
Definition: prepare_data.cxx:28

◆ finalize()

StatusCode MdcDedxRecon::finalize ( )

Definition at line 259 of file MdcDedxRecon.cxx.

260{
261 MsgStream log(msgSvc(), name());
262 log << MSG::INFO << "in finalize()" << endreq;
263
264 ex_trks.clear();
265 m_trkalgs.clear();
266
267 if(m_rootfile != "no rootfile")
268 {
269 TFile *f1 = new TFile(m_rootfile.c_str(),"recreate");
270 m_hlist->Write();
271 f1->Close();
272 delete f1;
273 delete m_hitNo_wire;
274 delete m_hitlevel;
275 delete m_hlist;
276 }
277 delete ex_calib;
278
279 cout<<"total event number is : "<<eventNo<<endl;
280 cout<<"total track number is : "<<trackNO1
281 <<" RecMdcTrack number is : "<<trackNO2
282 <<" RecMdcKalTrack number is :"<<trackNO3<<endl;
283 log << MSG::DEBUG <<"MdcDedxRecon terminated!!!"<< endreq;
284 return StatusCode::SUCCESS;
285}
f1
TFile * f1
Definition: DataBase/tau_mode.c:4

◆ initialize()

StatusCode MdcDedxRecon::initialize ( )

Definition at line 77 of file MdcDedxRecon.cxx.

78{
79 MsgStream log(msgSvc(), name());
80 log << MSG::INFO << "in initialize()" << endreq;
81
82 log << MSG::INFO <<"--------------------< MdcDedxRecon V2.1 >---------------------"<< endreq;
83 log << MSG::INFO <<"####### correction for the wire current dependence #######"<< endreq;
84 log << MSG::INFO <<"####### correction for the new z dependence #######"<< endreq;
85 log << MSG::INFO <<"-------------------------------------------------------------"<< endreq;
86 log << MSG::INFO <<"++++++++++++++++++++[ initialization ]+++++++++++++++++++++"<< endreq;
87 log << MSG::INFO << "MdcDedxRecon has been initialized with calib_flag: "<< calib_flag <<endreq;
88 log << MSG::INFO << "MdcDedxRecon has been initialized with landau: "<< landau << endreq;
89 if( landau == 0 ) {truncate = 0.7;}
90 log << MSG::INFO << "MdcDedxRecon has been initialized with ntpFlag: "<< ntpFlag << endreq;
91 log << MSG::INFO << "MdcDedxRecon has been initialized with recalg: "<< recalg << endreq;
92 log << MSG::INFO << "MdcDedxRecon has been initialized with dE/dx cal method m_alg: "<< m_alg << endreq;
93 log << MSG::INFO << "MdcDedxRecon has been initialized with truncate: "<< truncate <<endreq;
94 log << MSG::INFO << "MdcDedxRecon has been initialized with doNewGlobal: "<<doNewGlobal<<endreq;
95 log << MSG::DEBUG << "+++++++++++ MdcDedxRecon initialization end ++++++++++++ "<< endreq;
96 if( truncate <= 0.0 || 1.0 < truncate )
97 {
98 log << MSG::FATAL <<" Initialization ERROR of truncate = "<<truncate<< endreq;
99 log << MSG::FATAL <<" truncate must be within 0.0 to 1.0 ! "<< endreq;
100 log << MSG::FATAL <<" Please stop exec. "<< endreq;
101 }
102 log << MSG::DEBUG <<"-------------------------------------------------------------"<< endreq;
103 log << MSG::DEBUG <<"MdcDedxRecon init 2nd part!!!"<< endreq;
104
105 StatusCode scex = service("DedxCorrecSvc", exsvc);
106 if (scex == StatusCode::SUCCESS)
107 {
108 log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endreq;
109 }
110 else
111 {
112 log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endreq;
113 return StatusCode::SUCCESS;
114 }
115 exsvc->set_flag( calib_flag );
116
117 StatusCode cursvc = service("DedxCurSvc", dedxcursvc);
118 if (cursvc == StatusCode::SUCCESS)
119 {
120 log << MSG::INFO <<"DedxCurSvc is running"<<endreq;
121 }
122 else
123 {
124 log << MSG::ERROR <<"DedxCurSvc is failed"<<endreq;
125 return StatusCode::SUCCESS;
126 }
127
128
129 if( !ex_calib ) ex_calib = new MdcDedxCorrection;
130 //#ifdef DBCALIB
131 // ex_calib->getCalib(); //cleate MdcDedxWire and MdcDedxLayer.
132 //#endif
133
134 //------------------------------produce histogram root files-----------------------------//
135 if(m_rootfile!="no rootfile")
136 {
137 const char* preDir = gDirectory->GetPath();
138 m_hlist = new TObjArray(0);
139 m_hitlevel = new TFolder("dedx_hitlevel","hitlevel");
140 m_hlist -> Add(m_hitlevel);
141 m_hitNo_wire = new TH1F("dedx_HitNo_wire","dedx hitNo vs wire",6797, -0.5, 6796.5);
142 m_hitlevel -> Add(m_hitNo_wire);
143 gDirectory->cd(preDir);
144 }
145
146 //------------------------------produce ntuple files-------------------------------------//
147 if( ntpFlag ==2 )
148 {
149 StatusCode status;
150 NTuplePtr nt1(ntupleSvc(),"FILE103/n103");
151 if ( nt1 )
152 m_nt1 = nt1;
153 else
154 {
155 m_nt1= ntupleSvc()->book("FILE103/n103",CLID_ColumnWiseTuple,"dEdx vs momentum");
156 if ( m_nt1 )
157 {
158 status = m_nt1->addItem("phtm",m_phtm);
159 //status = m_nt1->addItem("median",m_median);
160 //status = m_nt1->addItem("geom",m_geometric);
161 //status = m_nt1->addItem("geom_tr",m_geometric_trunc);
162 //status = m_nt1->addItem("harm",m_harmonic);
163 //status = m_nt1->addItem("harm_tr",m_harmonic_trunc);
164 //status = m_nt1->addItem("transf",m_transform);
165 //status = m_nt1->addItem("logex",m_logtrunc);
166 status = m_nt1->addItem("dEdx_meas", m_dEdx_meas);
167
168 status = m_nt1->addItem("ptrk",m_ptrk);
169 status = m_nt1->addItem("sintheta",m_sintheta);
170 status = m_nt1->addItem("costheta",m_costheta);
171 status = m_nt1->addItem("charge",m_charge);
172 status = m_nt1->addItem("runNO",m_runNO);
173 status = m_nt1->addItem("evtNO",m_evtNO);
174 status = m_nt1->addItem("t0",m_t0);
175 status = m_nt1->addItem("trackId",m_trackId);
176 status = m_nt1->addItem("poca_x",m_poca_x);
177 status = m_nt1->addItem("poca_y",m_poca_y);
178 status = m_nt1->addItem("poca_z",m_poca_z);
179 status = m_nt1->addItem("recalg",m_recalg);
180
181 status = m_nt1->addItem("nhit",m_nhit);
182 status = m_nt1->addItem("usedhit",m_usedhit);
183 status = m_nt1->addItem("ndedxhit",m_nphlisthit,0,100);
184 status = m_nt1->addIndexedItem("dEdx_hit",m_nphlisthit,m_dEdx_hit);
185
186 status = m_nt1->addItem("type",m_parttype);
187 status = m_nt1->addItem("prob",m_prob);
188 status = m_nt1->addItem("expect",m_expect);
189 status = m_nt1->addItem("sigma",m_sigma);
190 status = m_nt1->addItem("chidedx",m_chidedx);
191 status = m_nt1->addItem("chidedx_e",m_chidedxe);
192 status = m_nt1->addItem("chidedx_mu",m_chidedxmu);
193 status = m_nt1->addItem("chidedx_pi",m_chidedxpi);
194 status = m_nt1->addItem("chidedx_k",m_chidedxk);
195 status = m_nt1->addItem("chidedx_p",m_chidedxp);
196 status = m_nt1->addItem("partid",5,m_probpid);
197 status = m_nt1->addItem("expectid",5,m_expectid);
198 status = m_nt1->addItem("sigmaid",5,m_sigmaid);
199 }
200 }
201
202 NTuplePtr nt2(ntupleSvc(),"FILE103/n102");
203 if ( nt2 ) m_nt2 = nt2;
204 else
205 {
206 m_nt2= ntupleSvc()->book("FILE103/n102",CLID_RowWiseTuple,"pulae height raw");
207 if ( m_nt2 )
208 {
209 status = m_nt2->addItem("charge",m_charge1);
210 status = m_nt2->addItem("adc_raw",m_phraw);
211 status = m_nt2->addItem("exraw",m_exraw);
212 status = m_nt2->addItem("runNO1",m_runNO1);
213 status = m_nt2->addItem("evtNO1",m_evtNO1);
214 status = m_nt2->addItem("nhit_hit", m_nhit_hit);
215 status = m_nt2->addItem("wire",m_wire);
216 //status = m_nt2->addItem("doca",m_doca);
217 status = m_nt2->addItem("doca_in",m_doca_in);
218 status = m_nt2->addItem("doca_ex",m_doca_ex);
219 status = m_nt2->addItem("driftdist",m_driftdist);
220 status = m_nt2->addItem("eangle",m_eangle);
221 status = m_nt2->addItem("costheta1",m_costheta1);
222 status = m_nt2->addItem("path_rphi",m_pathL);
223 status = m_nt2->addItem("layer",m_layer);
224 status = m_nt2->addItem("ptrk1",m_ptrk1);
225 status = m_nt2->addItem("ptrk_hit",m_ptrk_hit);
226 status = m_nt2->addItem("t01",m_t01);
227 status = m_nt2->addItem("tdc_raw",m_tdc_raw);
228 status = m_nt2->addItem("driftT",m_driftT);
229 status = m_nt2->addItem("localwid",m_localwid);
230 status = m_nt2->addItem("trackId1",m_trackId1);
231 }
232 }
233 }
234
235 return StatusCode::SUCCESS;
236}
Add
mg Add(gr3)
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
IDedxCorrecSvc::set_flag
virtual void set_flag(int calib_F)=0

◆ kaltrackrec()

void MdcDedxRecon::kaltrackrec ( RecMdcKalTrackCol::iterator  trk,
Identifier  mdcid,
double  tes,
int  RunNO,
int  eventNO,
int  runFlag,
MsgStream  log 
)

Definition at line 1291 of file MdcDedxRecon.cxx.

1292{
1293 vector<double> phlist; //dE/dx only after hit level correction
1294 vector<double> phlist_hit; //dE/dx after hit and track level correction
1295 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
1296 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
1297 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
1298 int Nhits=0;
1299 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
1300
1301 MdcDedxTrk trk_ex( **trk, ParticleFlag);
1302 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
1303 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
1304
1305 HepVector a;
1306 HepSymMatrix tkErrM;
1307 if(ParticleFlag>-1&&ParticleFlag<5)
1308 {
1309 DstMdcKalTrack* dstTrk = *trk;
1310 a = dstTrk->getFHelix(ParticleFlag);
1311 tkErrM = dstTrk->getFError(ParticleFlag);
1312 }
1313 else
1314 {
1315 a = (*trk)->getFHelix();
1316 tkErrM = (*trk)->getFError();
1317 }
1318 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
1319 //cout<<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl; //getFHelix is first layer of MdcKalTrack;
1320
1321 m_d0E = tkErrM[0][0];
1322 m_phi0E = tkErrM[1][1];
1323 m_cpaE = tkErrM[2][2];
1324 m_z0E = tkErrM[3][3];
1325
1326 phi0= a[1];
1327 costheta = cos(M_PI_2-atan(a[4]));
1328 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl; //theta in trk_ex is got by getFHelix();
1329 sintheta = sin(M_PI_2-atan(a[4]));
1330 m_Pt = 1.0/fabs( a[2] );
1331 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
1332 //cout<<"track ptrk: "<<m_P<<" extrack ptrk: "<<trk_ex.P()<<endl;
1333 charge = ( a[2] > 0 )? 1 : -1;
1334 //cout<<"track charge: "<<charge<<" extrack charge: "<<(*trk)->charge()<<endl;
1335 dedxhitrefvec = new DedxHitRefVec;
1336
1337
1338 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
1339 //HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
1340 //if not set ParticleFlag, we will get the last succefully hypothesis;
1341 Nhits = gothits.size();
1342 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
1343 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
1344
1345 RecMdcHit* mdcHit = 0;
1346 HelixSegRefVec::iterator it_gothit = gothits.begin();
1347 for( ; it_gothit != gothits.end(); it_gothit++)
1348 {
1349 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
1350 //HepVector a_hit_ex = (*it_gothit)->getHelixExcl(); //same with getHelixIncl()
1351 HepPoint3D IP(0,0,0);
1352 Dedx_Helix exhel_hit_in(IP, a_hit_in);
1353 //Dedx_Helix exhel_hit_ex(IP, a_hit_ex);
1354 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
1355 //cout<<"getHelixIncl 5 para: "<<a_hit_in[0]<<" "<<a_hit_in[1] <<" "<<a_hit_in[2]<<" "<<a_hit_in[3]<<" "<<a_hit_in[4]<<endl;
1356
1357 mdcid = (*it_gothit)->getMdcId();
1358 layid = MdcID::layer(mdcid);
1359 localwid = MdcID::wire(mdcid);
1360 w0id = geosvc->Layer(layid)->Wirst();
1361 wid = w0id + localwid;
1362 adc_raw = (*it_gothit)->getAdc();
1363 tdc_raw = (*it_gothit)->getTdc();
1364 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
1365 zhit = (*it_gothit)->getZhit();
1366 lr = (*it_gothit)->getFlagLR();
1367 if(lr == 2) cout<<"lr = "<<lr<<endl;
1368 driftD = (*it_gothit)->getDD();
1369 driftd = abs(driftD);
1370 driftT = (*it_gothit)->getDT();
1371 dd_in = (*it_gothit)->getDocaIncl(); //getDocaIncl() include fit unused hit
1372 dd_ex = (*it_gothit)->getDocaExcl(); //getDocaExcl() exclude fit unused hit
1373
1374 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
1375 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
1376 //dd = dd/doca_norm[layid];
1377 //cout<<"lr: "<<lr<<" dd_in: "<<dd_in<<" dd_ex: "<<dd_ex<<endl;
1378
1379 eangle = (*it_gothit)->getEntra();
1380 eangle = eangle/M_PI;
1381 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
1382 rphi_path=pathlength * sintheta;
1383
1384 //cout<<"ph para check: "<<"runFlag: "<<runFlag<<" runNO: "<<runNO<<" pathlength: "<<pathlength<<" wid: "<<wid<<" layid: "<<layid<<" adc_raw: "<<adc_raw<<" dd_in: "<<dd_in<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
1385 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
1386 //cout<<"tes= "<<tes<<endl;
1387 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
1388 //cout<<"adc_raw= "<<adc_raw<<" ph= "<<ph<<endl;
1389 //cout<<"adc_raw = "<<adc_raw<<" ph_hit= "<<ph_hit<<endl;
1390 //if(pathlength == -1)
1391 //cout<<"parameter1: "<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
1392
1393 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
1394 if(runNO<0)
1395 {
1396 if (layid<8)
1397 {
1398 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1399 }
1400 else if(layid >= 8)
1401 {
1402 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1403 }
1404 }
1405 else if(runNO>=0)
1406 {
1407 if(layid <8)
1408 {
1409 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1410 }
1411 else if(layid >= 8)
1412 {
1413 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1414 }
1415 }
1416 //cout<<"recAlg=2 para kal: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd_in<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
1417
1418 if (ph<0.0||ph==0) continue;
1419 else
1420 {
1421 //-----------------------put data into TDS-----------------------------//
1422 dedxhit = new RecMdcDedxHit;
1423 dedxhit->setMdcHit(mdcHit);
1424 dedxhit->setMdcKalHelixSeg(*it_gothit);
1425 dedxhit->setMdcId(mdcid);
1426 dedxhit->setFlagLR(lr);
1427 dedxhit->setTrkId(trk_ex.trk_id());
1428 dedxhit->setDedx(ph_hit);
1429 dedxhit->setPathLength(pathlength);
1430
1431 //---------------------------Fill root file----------------------------//
1432 if(m_rootfile!= "no rootfile")
1433 {
1434 m_hitNo_wire->Fill(wid);
1435 }
1436
1437 //-------------------------Fill ntuple n102---------------------------//
1438 if ( ntpFlag ==2 )
1439 {
1440 m_charge1 = (*trk)->charge();
1441 m_t01 = tes;
1442 m_driftT = driftT;
1443 m_tdc_raw = tdc_raw;
1444 m_phraw = adc_raw;
1445 m_exraw = ph_hit;
1446 m_localwid = localwid;
1447 m_wire = wid;
1448 m_runNO1 = runNO;
1449 m_evtNO1 = eventNO;
1450 m_nhit_hit = Nhits;
1451 m_doca_in = dd_in;
1452 m_doca_ex = dd_ex;
1453 m_driftdist = driftD;
1454 m_eangle = eangle;
1455 m_costheta1 = costheta;
1456 m_pathL = pathlength;
1457 m_layer = layid;
1458 m_ptrk1 = m_P;
1459 m_ptrk_hit = p_hit;
1460 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
1461 m_trackId1 = trk_ex.trk_id();
1462 StatusCode status2 = m_nt2->write();
1463 if ( status2.isFailure() )
1464 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
1465 }
1466 if(layid>3)
1467 {
1468 phlist.push_back(ph);
1469 phlist_hit.push_back(ph_hit);
1470 }
1471 }
1472 dedxhitList->push_back( dedxhit );
1473 dedxhitrefvec->push_back( dedxhit );
1474 }//end of hit loop
1475 trk_ex.set_phlist( phlist );
1476 trk_ex.set_phlist_hit( phlist_hit );
1477 trk_ex.setVecDedxHits( *dedxhitrefvec );
1478 ex_trks.push_back(trk_ex );
1479 m_trkalgs.push_back(m_trkalg);
1480
1481 delete dedxhitrefvec;
1482 phlist.clear();
1483 phlist_hit.clear();
1484 if(ntpFlag>0) trackNO3++;
1485}

Referenced by execute().

◆ mdctrackrec()

void MdcDedxRecon::mdctrackrec ( RecMdcTrackCol::iterator  trk,
Identifier  mdcid,
double  tes,
int  RunNO,
int  eventNO,
int  runFlag,
MsgStream  log 
)

Definition at line 1121 of file MdcDedxRecon.cxx.

1122{
1123 vector<double> phlist; //dE/dx only after hit level correction
1124 vector<double> phlist_hit; //dE/dx after hit and track level correction
1125 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
1126 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
1127 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
1128 int Nhits=0;
1129 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
1130
1131 MdcDedxTrk trk_ex( **trk);
1132 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
1133 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
1134
1135 HepVector a = (*trk)->helix();
1136 HepSymMatrix tkErrM = (*trk)->err();
1137 m_d0E = tkErrM[0][0];
1138 m_phi0E = tkErrM[1][1];
1139 m_cpaE = tkErrM[2][2];
1140 m_z0E = tkErrM[3][3];
1141
1142 HepPoint3D IP(0,0,0);
1143 Dedx_Helix exhel(IP, a); //initialize class Dedx_Helix for DedxCorrecSvc
1144 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
1145 //cout<<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl;
1146 phi0= a[1];
1147 costheta = cos(M_PI_2-atan(a[4]));
1148 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl;
1149 sintheta = sin(M_PI_2-atan(a[4]));
1150 m_Pt = 1.0/fabs( a[2] );
1151 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
1152 charge = ( a[2] > 0 )? 1 : -1;
1153 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
1154 dedxhitrefvec = new DedxHitRefVec;
1155
1156
1157 HitRefVec gothits= (*trk)->getVecHits();
1158 Nhits = gothits.size();
1159 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
1160 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
1161
1162 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
1163 HitRefVec::iterator it_gothit = gothits.begin();
1164 for( ; it_gothit != gothits.end(); it_gothit++)
1165 {
1166 mdcid = (*it_gothit)->getMdcId();
1167 layid = MdcID::layer(mdcid);
1168 localwid = MdcID::wire(mdcid);
1169 w0id = geosvc->Layer(layid)->Wirst();
1170 wid = w0id + localwid;
1171 adc_raw = (*it_gothit)->getAdc();
1172 tdc_raw = (*it_gothit)->getTdc();
1173 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
1174 zhit = (*it_gothit)->getZhit();
1175 lr = (*it_gothit)->getFlagLR();
1176 if(lr == 2) cout<<"lr = "<<lr<<endl;
1177 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
1178 else driftD = (*it_gothit)->getDriftDistRight();
1179 //cout<<"lr: "<<lr<<" driftD: "<<driftD<<endl;
1180 driftd = abs(driftD);
1181 dd = (*it_gothit)->getDoca();
1182 if(lr == 0 || lr == 2 ) dd = -abs(dd);
1183 if(lr == 1 ) dd = abs(dd);
1184 driftT = (*it_gothit)->getDriftT();
1185
1186 eangle = (*it_gothit)->getEntra();
1187 eangle = eangle/M_PI;
1188 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
1189 rphi_path=pathlength * sintheta;
1190
1191 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd: "<<dd<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
1192 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
1193 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
1194 //if(pathlength == -1)
1195 //cout<<"parameter0: "<<"eventNO: "<<eventNO<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
1196
1197 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
1198 if(runNO<0)
1199 {
1200 if (layid<8)
1201 {
1202 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1203 }
1204 else if(layid >= 8)
1205 {
1206 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1207 }
1208 }
1209 else if(runNO>=0)
1210 {
1211 if(layid <8)
1212 {
1213 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1214 }
1215 else if(layid >= 8)
1216 {
1217 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1218 }
1219 }
1220 //cout<<"recAlg=2 para mdc: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
1221
1222 if (ph<0.0||ph==0) continue;
1223 else
1224 {
1225 //-----------------------put data into TDS-----------------------------//
1226 dedxhit = new RecMdcDedxHit;
1227 dedxhit->setMdcHit(*it_gothit);
1228 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg);
1229 dedxhit->setMdcId(mdcid);
1230 dedxhit->setFlagLR(lr);
1231 dedxhit->setTrkId(trk_ex.trk_id());
1232 dedxhit->setDedx(ph_hit);
1233 dedxhit->setPathLength(pathlength);
1234
1235 //---------------------------Fill root file----------------------------//
1236 if(m_rootfile!= "no rootfile")
1237 {
1238 m_hitNo_wire->Fill(wid);
1239 }
1240
1241 //-------------------------Fill ntuple n102---------------------------//
1242 if ( ntpFlag ==2 )
1243 {
1244 m_charge1 = (*trk)->charge();
1245 m_t01 = tes;
1246 m_driftT = driftT;
1247 m_tdc_raw = tdc_raw;
1248 m_phraw = adc_raw;
1249 m_exraw = ph_hit;
1250 m_localwid = localwid;
1251 m_wire = wid;
1252 m_runNO1 = runNO;
1253 m_evtNO1 = eventNO;
1254 m_nhit_hit = Nhits;
1255 m_doca_in = dd;
1256 m_doca_ex = dd;
1257 m_driftdist = driftD;
1258 m_eangle = eangle;
1259 m_costheta1 = costheta;
1260 m_pathL = pathlength;
1261 m_layer = layid;
1262 m_ptrk1 = m_P;
1263 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
1264 m_trackId1 = trk_ex.trk_id();
1265 StatusCode status2 = m_nt2->write();
1266 if ( status2.isFailure() )
1267 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
1268 }
1269 if(layid>3)
1270 {
1271 phlist.push_back(ph);
1272 phlist_hit.push_back(ph_hit);
1273 }
1274 }
1275 dedxhitList->push_back( dedxhit );
1276 dedxhitrefvec->push_back( dedxhit );
1277 }//end of hit loop
1278 trk_ex.set_phlist( phlist );
1279 trk_ex.set_phlist_hit( phlist_hit );
1280 trk_ex.setVecDedxHits( *dedxhitrefvec );
1281 ex_trks.push_back(trk_ex );
1282 m_trkalgs.push_back(m_trkalg);
1283
1284 delete dedxhitrefvec;
1285 phlist.clear();
1286 phlist_hit.clear();
1287 if(ntpFlag>0) trackNO2++;
1288}

Referenced by switchtomdctrack().

◆ switchtomdctrack()

void MdcDedxRecon::switchtomdctrack ( int  trkid,
Identifier  mdcid,
double  tes,
int  RunNO,
int  eventNO,
int  runFlag,
MsgStream  log 
)

Definition at line 1487 of file MdcDedxRecon.cxx.

1488{
1489 //retrieve MdcTrackCol from TDS
1490 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
1491 if (!newtrkCol)
1492 {
1493 log << MSG::WARNING << "Could not find RecMdcTrackCol in switchtomdctrack" << endreq;
1494 return ;
1495 }
1496
1497 RecMdcTrackCol::iterator trk = newtrkCol->begin();
1498 for( ; trk != newtrkCol->end(); trk++)
1499 {
1500 if( (*trk)->trackId()==trkid)
1501 {
1502 HitRefVec gothits= (*trk)->getVecHits();
1503 if(gothits.size()>=2)
1504 mdctrackrec(trk,mdcid,tes,runNO,eventNO,runFlag,log);
1505 }
1506 }
1507}
MdcDedxRecon::mdctrackrec
void mdctrackrec(RecMdcTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1121

Referenced by execute().

◆ tracks()

const std::vector< MdcDedxTrk > & MdcDedxRecon::tracks ( void  ) const

Definition at line 1114 of file MdcDedxRecon.cxx.

1115{
1116 return ex_trks;
1117}
The documentation for this class was generated from the following files: