BOSS 6.6.4.p01
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 runFlag, MsgStream log)
 
void kaltrackrec (RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int runFlag, MsgStream log)
 
void switchtomdctrack (int trkid, Identifier mdcid, double tes, int RunNO, 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 52 of file MdcDedxRecon.cxx.

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

◆ ~MdcDedxRecon()

MdcDedxRecon::~MdcDedxRecon ( )
inline

Definition at line 33 of file MdcDedxRecon.h.

33{};

Member Function Documentation

◆ beginRun()

StatusCode MdcDedxRecon::beginRun ( )

Definition at line 237 of file MdcDedxRecon.cxx.

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

◆ clearTables()

void MdcDedxRecon::clearTables ( )

Definition at line 1115 of file MdcDedxRecon.cxx.

1115{}

◆ execute()

StatusCode MdcDedxRecon::execute ( )

Definition at line 286 of file MdcDedxRecon.cxx.

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

◆ finalize()

StatusCode MdcDedxRecon::finalize ( )

Definition at line 257 of file MdcDedxRecon.cxx.

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

◆ initialize()

StatusCode MdcDedxRecon::initialize ( )

Definition at line 76 of file MdcDedxRecon.cxx.

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

Definition at line 1286 of file MdcDedxRecon.cxx.

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

Referenced by execute().

◆ mdctrackrec()

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

Definition at line 1117 of file MdcDedxRecon.cxx.

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

Referenced by switchtomdctrack().

◆ switchtomdctrack()

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

Definition at line 1481 of file MdcDedxRecon.cxx.

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

Referenced by execute().

◆ tracks()

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

Definition at line 1110 of file MdcDedxRecon.cxx.

1111{
1112 return ex_trks;
1113}
The documentation for this class was generated from the following files: