Definition at line 306 of file EsTimeAlg.cxx.
References Helix::a(), abs, alpha, MdcGeoWire::Backward(), IEstTofCaliSvc::BTime1(), IEstTofCaliSvc::BTime2(), Helix::center(), cos(), Bes_Common::DEBUG, check_raw_filter::dist, MdcCalibFunSvc::distToDriftTime(), MdcUtilitySvc::doca(), ELMAS2, Emc_helix::Emc_Get(), TofID::endcap(), calibUtil::ERROR, EST_Trimmer(), TofFz_helix::Etfid, IEstTofCaliSvc::EtfTime(), IEstTofCaliSvc::ETime(), eventNo, Bes_Common::FATAL, MdcGeoWire::Forward(), g_abmom, g_bunchtime, g_difftof_b, g_difftof_e, g_EstimeMdc, g_eventNo, g_mcTestime, g_meant0, g_meantdc, g_meantev, g_meantevdown, g_meantevup, g_ndriftt, g_nmatch, g_nmatch_tot, g_nmatchbarrel, g_nmatchbarrel_1, g_nmatchbarrel_2, g_nmatchend, g_nmatchmdc, g_ntofdown, g_ntofdown1, g_ntofup, g_ntofup1, g_ntrk, g_ntrkMC, g_phi0MC, g_pid, g_ptMC, g_pxMC, g_pyMC, g_pzMC, g_runNo, g_T0, g_t0, g_t0back, g_t0barrelTof, g_t0for, g_t0mean, g_t0offset_b, g_t0offset_e, g_Testime, g_Testime_fzisan, g_theta0MC, g_ttof, g_vel, MdcCalibFunSvc::getT0(), MdcCalibFunSvc::getTimeWalk(), genRecEmupikp::i, Helix::ignoreErrorMatrix(), Bes_Common::INFO, ganga-rec::j, TofFz_helix::Kappa, IMdcGeomSvc::Layer(), MdcID::layer(), m_beforrec, m_bunchtime_MC, m_debug, m_ebeam, m_estFlag, m_estTime, m_evtNo, m_flag, m_mdcCalibFunSvc, m_mdcopt, m_mdcUtilitySvc, m_ntupleflag, m_optCosmic, m_particleTable, m_pass, m_pCalibDataSvc, m_phase, M_PI, m_rawDataProviderSvc, m_TofOpt, m_TofOpt_Cut, m_tuple2, m_tuple3, m_tuple9, m_userawtime_opt, m_useSw, m_useT0cal, m_useTimeFactor, m_useXT, EstParameter::MDC_Debug(), EstParameter::MDC_drCut(), EstParameter::MDC_dzCut(), EstParameter::MDC_Inner(), EstParameter::MDC_Prop(), EstParameter::MDC_Tof(), RawDataUtil::MdcTime(), Helix::momentum(), momentum, msgSvc(), MUMAS2, MXTKHIT, MXTRK, MXWIRE, MdcGeoLayer::NCell(), ndriftt, ntrk, ntrkMC, offset_dt, offset_dt_end, Opt_new(), TofFz_helix::Path_etf, TofFz_helix::Pathl, TofFz_helix::pathlCut(), EstParameter::pathlCut(), Emc_helix::pathlCut(), MdcGeoLayer::PCSiz(), phi0, Emc_helix::phi_emc, pid, PIMAS2, Helix::pivot(), boss::pos, proton, PROTONMAS2, q, TofFz_helix::r_endtof, TofFz_helix::r_etf, Helix::radius(), MdcGeoLayer::Radius(), RCEMC2, RCTOF2, check_raw_filter::runno, runNo, MdcGeoWire::Slant(), storeTDS(), deljobs::string, t(), TofFz_helix::Tanl, Emc_helix::theta_emc, tofCaliSvc, IRawDataProviderSvc::tofDataVectorEstime(), toffset_rawtime, toffset_rawtime_e, TofFz_helix::TofFz_Get(), TofFz_helix::Tofid, RawDataUtil::TofTime(), IEstTofCaliSvc::ValidInfo(), EstParameter::vdrift(), VELPROP, VLIGHT, Bes_Common::WARNING, weight, IMdcGeomSvc::Wire(), MdcID::wire(), Helix::x(), x, Emc_helix::Z_emc, TofFz_helix::Z_etf, TofFz_helix::Z_tof, TofFz_helix::ztofCut(), EstParameter::ztofCutmax(), and EstParameter::ztofCutmin().
00306 {
00307
00308 MsgStream log(msgSvc(), name());
00309 log << MSG::INFO << " tof " << endreq;
00310
00311
00312 EstParameter Estparam;
00313 double offset=0, t_quality=0, tOffset_b=0, tOffset_e=0;
00314 int idtof , tofid_helix[30]={-9},idmatch[3][88]={0},idmatch_emc[3][88]={0} ,idt[15]={0},particleId[30]={0}, tofid_emc[2]={0}, module[20]={0};
00315 int idetf, etfid_helix[30]={-9}, idetfmatch[3][36]={-9}, idmatch_etf_emc[3][36]={0}, etfid_emc[2]={0};
00316 int ntot=0,in=-1,out=-1, emcflag1=0, emcflag2=0, tof_flag=0; double bunchtime=m_bunchtime_MC;
00317 double meant[15]={0.},adc[15]={0.}, momentum[15]={0.},r_endtof[10]={0.};
00318 double ttof[30]={0.},helztof[30]={0.0},mcztof=0.0,forevtime=0.0,backevtime=0.0,meantev[200]={0.},meantevup[200]={0.0},meantevdown[200]={0.0};
00319 double t0forward=0,t0backward=0,t0forward_trk=0,t0backward_trk=0;
00320 double t0forward_add=0,t0backward_add=0,t_Est=-999;
00321 double thetaemc_rec[20]={0.},phiemc_rec[20]={0.},energy_rec[20]={0.},xemc_rec[20]={0.},yemc_rec[20]={0.},zemc_rec[20]={0.};
00322 double r_endetf[10]={0.}, tetf[30]={0.}, helzetf[30]={0.}, helpathetf[36]={0.}, abmom2etf[36]={0.};
00323
00324 int nmatch1=0,nmatch2=0,nmatch_barrel=0,nmatch_end=0,nmatch_mdc=0, nmatch_barrel_1=0, nmatch_barrel_2=0,nmatch=0,ntofup=0,ntofdown=0;
00325 double sum_EstimeMdc=0,sum_EstimeMdcMC=0;
00326 int nbunch=0,tEstFlag=0, runNo=0;
00327 double helpath[88]={0.},helz[88]={0.},abmom2[88]={0.};
00328 double mcTestime=0,trigtiming=0;
00329 double mean_tdc_btof[2][88]={0.}, mean_tdc_etof[3][48]={0.}, mean_tdc_etf[3][36][12]={0.};
00330 int optCosmic=m_optCosmic;
00331 int m_userawtime=m_userawtime_opt;
00332 double Testime_fzisan= -999.;
00333 int TestimeFlag_fzisan= -999;
00334 double TestimeQuality_fzisan= -999.;
00335 double Tof_t0Arr[500]={-999.};
00336
00337 bool useEtofScin = ( m_phase<3 );
00338 bool useEtofMRPC = ( m_phase>2 );
00339
00340
00341 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00342 if (!eventHeader) {
00343 log << MSG::FATAL << "Could not find Event Header" << endreq;
00344 return StatusCode::FAILURE;
00345 }
00346 log << MSG::INFO << "EsTimeAlg: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
00347 int eventNo=eventHeader->eventNumber();
00348 runNo=eventHeader->runNumber();
00349
00350 if(m_ntupleflag && m_tuple2){g_runNo=runNo; g_eventNo=eventNo;}
00351 if(runNo<0) {
00352 offset_dt=0;
00353 offset_dt_end=0;
00354 toffset_rawtime=0;
00355 toffset_rawtime_e=0;
00356 if( useEtofMRPC ) { toffset_rawtime_e = -20.0; }
00357 }
00358 if( runNo>0 && useEtofMRPC ) { toffset_rawtime_e = 1.7; }
00359
00360 if(runNo>0 && runNo<5336) offset_dt=5.8;
00361 if(runNo>5462 && runNo<5466){ offset_dt=1.6; offset_dt_end=0.4;}
00362 if(runNo>5538 && runNo<5920){ offset_dt=-0.2; offset_dt_end=-1.2;}
00363 if(runNo>5924 && runNo<6322){ offset_dt=-0.2; offset_dt_end=-0.1;}
00364 if(runNo>6400 && runNo<6423){ offset_dt=-2.4; offset_dt_end=-1.9;}
00365 if(runNo>6514 && runNo<6668){ offset_dt=0; offset_dt_end=3.4;}
00366 if(runNo>6667 && runNo<6715){ offset_dt=4; offset_dt_end=7.2;}
00367 if(runNo>6715 && runNo<6720){ offset_dt=8; offset_dt_end=7.5;}
00368 if(runNo>6879 && runNo<6909){ offset_dt=0.2; offset_dt_end=-0.1;}
00369 if(m_ntupleflag && m_tuple3){
00370 g_bunchtime=8;
00371 g_t0offset_b=2.0;
00372 g_t0offset_e=2.0;
00373 }
00374
00375 m_pass[0]++;
00376 if(m_evtNo==1 && m_pass[0]%1000 ==0){
00377 cout<<"------------------- Events-----------------: "<<m_pass[0]<<endl;
00378 }
00379 if(m_debug==4) cout<<"m_userawtime: "<<m_userawtime<<endl;
00380 if(m_debug==4) cout<<"EstTofCalibSvc est flag: "<<tofCaliSvc->ValidInfo()<<endl;
00381 if(tofCaliSvc->ValidInfo()==0&&m_userawtime_opt==0)
00382 {
00383 log << MSG::ERROR << "EstTof Calibration Const is Invalid! " << endreq;
00384 return StatusCode::FAILURE;
00385 }
00386 if(tofCaliSvc->ValidInfo()==0||m_userawtime_opt==1) m_userawtime=1;
00387 else m_userawtime=0;
00388
00389 SmartDataPtr<RecEsTimeCol> aRecestimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00390 if (!aRecestimeCol || aRecestimeCol->size()==0) {
00391 if(m_debug==4) log << MSG::INFO << "Could not find RecEsTimeCol from fzsian" << endreq;
00392 }else{
00393 RecEsTimeCol::iterator it_evt = aRecestimeCol->begin();
00394 for(; it_evt!=aRecestimeCol->end(); it_evt++){
00395 Testime_fzisan = (*it_evt)->getTest();
00396 TestimeFlag_fzisan = (*it_evt)->getStat();
00397 TestimeQuality_fzisan = (*it_evt)->getQuality();
00398
00399 log << MSG::INFO << "fzisan : Test = "<<(*it_evt)->getTest()
00400 << ", Status = "<<(*it_evt)->getStat() <<endreq;
00401
00402 if(m_ntupleflag && m_tuple2) g_Testime_fzisan = Testime_fzisan;
00403 }}
00404
00405
00406 std::string fullPath = "/Calib/EsTimeCal";
00407 SmartDataPtr<CalibData::EsTimeCalibData> TEst(m_pCalibDataSvc, fullPath);
00408 if(!TEst){ cout<<"ERROR EsTimeCalibData"<<endl;}
00409 else{
00410 int no = TEst->getTestCalibConstNo();
00411
00412
00413
00414
00415 log<<MSG::INFO<<"offset barrel t0="<< TEst->getToffsetb()
00416 <<", offset endcap t0="<< TEst->getToffsete()
00417 <<", bunch time ="<<TEst->getBunchTime()
00418 <<endreq;
00419 tOffset_b=TEst->getToffsetb();
00420 tOffset_e=TEst->getToffsete();
00421 bunchtime=TEst->getBunchTime();
00422 }
00423 if(m_userawtime) {
00424 tOffset_b=0;
00425 tOffset_e=0;
00426 }
00427 else
00428 {
00429 toffset_rawtime=0;
00430 toffset_rawtime_e=0;
00431 offset_dt=0;
00432 offset_dt_end=0;
00433 }
00434
00435 if(runNo>0&&m_useTimeFactor)
00436 {
00437 bunchtime=RawDataUtil::TofTime(8*1024)*bunchtime/24.0;
00438 }
00439
00440 if(runNo<0) bunchtime=m_bunchtime_MC;
00441
00442
00443 int digiId;
00444 if(runNo<0){
00445 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00446 if (!mcParticleCol) {
00447 log << MSG::INFO<< "Could not find McParticle" << endreq;
00448 }
00449 else{
00450 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00451 digiId = 0;
00452 ntrkMC = 0;
00453 int charge = 0;
00454
00455 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
00456 int statusFlags = (*iter_mc)->statusFlags();
00457 int pid = (*iter_mc)->particleProperty();
00458 log << MSG::INFO
00459 << " MC ParticleId = " << pid
00460 << " statusFlags = " << statusFlags
00461 << " PrimaryParticle = " <<(*iter_mc)->primaryParticle()
00462 << endreq;
00463 if(m_ntupleflag && m_tuple2){
00464 g_theta0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().theta();
00465 g_phi0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().phi();
00466 g_pxMC[ntrkMC] = (*iter_mc)->initialFourMomentum().px()/1000;
00467 g_pyMC[ntrkMC] = (*iter_mc)->initialFourMomentum().py()/1000;
00468 g_pzMC[ntrkMC] = (*iter_mc)->initialFourMomentum().pz()/1000;
00469 g_ptMC[ntrkMC] = sqrt(((*iter_mc)->initialFourMomentum().px())*((*iter_mc)->initialFourMomentum().px())+((*iter_mc)->initialFourMomentum().py())*((*iter_mc)->initialFourMomentum().py()))/1000;
00470 }
00471 if( pid >0 ) {
00472 if(m_particleTable->particle( pid )) charge = m_particleTable->particle( pid )->charge();
00473 } else if ( pid <0 ) {
00474 if(m_particleTable->particle( -pid )) {
00475 charge = m_particleTable->particle( -pid )->charge();
00476 charge *= -1;
00477 }
00478 } else {
00479 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00480 }
00481 log << MSG::DEBUG
00482 << "MC ParticleId = " << pid << " charge = " << charge
00483 << endreq;
00484 if(charge !=0 && abs(cos((*iter_mc)->initialFourMomentum().theta()))<0.93){
00485 ntrkMC++;
00486 }
00487 if(((*iter_mc)->primaryParticle())&& m_ntupleflag && m_tuple2){
00488 g_mcTestime=(*iter_mc)->initialPosition().t();
00489 mcTestime=(*iter_mc)->initialPosition().t();
00490 }
00491 }
00492 if(m_ntupleflag && m_tuple2) g_ntrkMC = ntrkMC;
00493 }
00494 }
00495 if (m_debug) cout<<"bunchtime: "<<bunchtime<<endl;
00496
00497 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00498 if (!newtrkCol || newtrkCol->size()==0) {
00499 log << MSG::INFO<< "Could not find RecMdcTrackCol" << endreq;
00500 } else{
00501 log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq;
00502 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
00503 for( ; iter_trk != newtrkCol->end(); iter_trk++){
00504 log << MSG::DEBUG << "retrieved MDC track:"
00505 << " Track Id: " << (*iter_trk)->trackId()
00506 << " Phi0: " << (*iter_trk)->helix(0)
00507 << " kappa: " << (*iter_trk)->helix(2)
00508 << " Tanl: " << (*iter_trk)->helix(4)
00509 << " Phi terminal: "<< (*iter_trk)->getFiTerm()
00510 << endreq
00511 << "Number of hits: "<< (*iter_trk)->getNhits()
00512 << " Number of stereo hits " << (*iter_trk)->nster()
00513 << endreq;
00514 double kappa = (*iter_trk)->helix(2);
00515 double tanl = (*iter_trk)->helix(4);
00516 momentum[ntot] = 1./fabs(kappa) * sqrt(1.+tanl*tanl);
00517 if((*iter_trk)->helix(3)>50.0) continue;
00518 ntot++;
00519 if (ntot>15) break;
00520 }
00521 }
00522
00523
00524 int emctrk=0;
00525 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00526 if (!aShowerCol || aShowerCol->size()==0) {
00527 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00528 } else{
00529 RecEmcShowerCol::iterator iShowerCol=aShowerCol->begin();
00530 for(;iShowerCol!=aShowerCol->end();iShowerCol++,emctrk++){
00531 if(emctrk>20) break;
00532 phiemc_rec[emctrk]=(*iShowerCol)->position().phi();
00533 thetaemc_rec[emctrk]=(*iShowerCol)->position().theta();
00534 energy_rec[emctrk]=(*iShowerCol)->energy();
00535 xemc_rec[emctrk]=(*iShowerCol)->x();
00536 yemc_rec[emctrk]=(*iShowerCol)->y();
00537 zemc_rec[emctrk]=(*iShowerCol)->z();
00538 module[emctrk]=(*iShowerCol)->module();
00539 }
00540 }
00541 for(int i=0; i<2; i++){
00542 double fi_endtof = atan2(yemc_rec[i],xemc_rec[i] );
00543 if( fi_endtof<0 ) { fi_endtof=2*3.141593+fi_endtof; }
00544 if( module[i]==1 ) {
00545 int Tofid = (int)(fi_endtof/(3.141593/44));
00546 if(Tofid>87) Tofid=Tofid-88;
00547 tofid_emc[i]=Tofid;
00548 idmatch_emc[1][Tofid]=1;
00549 }
00550 else{
00551 if( useEtofScin ) {
00552 int Tofid =(int)(fi_endtof/(3.141593/24));
00553 if( Tofid>47) Tofid=Tofid-48;
00554 tofid_emc[i]= Tofid;
00555 if(module[i]==2) idmatch_emc[2][Tofid]=1;
00556 if(module[i]==0) idmatch_emc[0][Tofid]=1;
00557 }
00558 if( useEtofMRPC ) {
00559 int Tofid = (int)(fi_endtof/(3.141593/18));
00560 if( Tofid>35) Tofid=Tofid-36;
00561 etfid_emc[i]= Tofid;
00562 if(module[i]==2) idmatch_etf_emc[2][Tofid]=1;
00563 if(module[i]==0) idmatch_etf_emc[0][Tofid]=1;
00564 }
00565 }
00566 }
00567
00568 ntrk=0;
00569 if (ntot > 0) {
00570 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
00571 for( int idfztrk=0; iter_trk != newtrkCol->end(); iter_trk++,idfztrk++){
00572 double mdcftrk[5];
00573 mdcftrk[0] = (*iter_trk)->helix(0);
00574 mdcftrk[1] = (*iter_trk)->helix(1);
00575 mdcftrk[2] =-(*iter_trk)->helix(2);
00576 mdcftrk[3] = (*iter_trk)->helix(3);
00577 mdcftrk[4] = (*iter_trk)->helix(4);
00578
00579 if(optCosmic==0){
00580 Emc_helix EmcHit;
00581
00582 EmcHit.pathlCut(Estparam.pathlCut());
00583
00584
00585
00586 if (EmcHit.Emc_Get(sqrt(RCEMC2),idfztrk,mdcftrk) > 0){
00587 double z_emc = EmcHit.Z_emc;
00588 double thetaemc_ext= EmcHit.theta_emc;
00589 double phiemc_ext = EmcHit.phi_emc;
00590
00591 for(int t=0; t<emctrk;t++){
00592 if((thetaemc_ext>=(thetaemc_rec[t]-0.1)) && (thetaemc_ext<=(thetaemc_rec[t]+0.1)) && (phiemc_ext>=(phiemc_rec[t]-0.1)) && (phiemc_ext<=(phiemc_rec[t]+0.1))){
00593 if((energy_rec[t])>=(0.85*momentum[idfztrk]))
00594 particleId[idfztrk]=1;
00595 }
00596 }
00597 }
00598
00599 if(particleId[idfztrk]!=1){
00600
00601 SmartDataPtr<RecMdcDedxCol> newdedxCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
00602 if (!newdedxCol || newdedxCol->size()==0) {
00603 log << MSG::WARNING<< "Could not find RecMdcDedxCol" << endreq;
00604 }
00605 else{
00606 RecMdcDedxCol::iterator it_dedx = newdedxCol->begin();
00607 for( int npid=0; it_dedx != newdedxCol->end(); it_dedx++,npid++) {
00608 log << MSG::INFO << "retrieved MDC dE/dx: "
00609 << "dEdx Id: " << (*it_dedx)->trackId()
00610 << " particle Id: "<< (*it_dedx)->particleType() <<endreq;
00611 if((*it_dedx)->particleType() == proton){
00612 particleId[npid]= 5;
00613 }
00614 if(m_debug==4) cout<<"dedx pid: "<<particleId[npid]<<endl;
00615 }
00616 }
00617 }
00618 }
00619
00620 idtof = -100;
00621 idetf = -100;
00622 TofFz_helix TofHit;
00623
00624
00625 TofHit.pathlCut(Estparam.pathlCut());
00626
00627
00628 TofHit.ztofCut(Estparam.ztofCutmin(),Estparam.ztofCutmax());
00629
00630 int tofpart = TofHit.TofFz_Get(sqrt(RCTOF2),idfztrk,mdcftrk);
00631 if(tofpart < 0) continue;
00632
00633
00634 bool useBarrelScin = ( tofpart==1 );
00635 bool useEndcapScin = ( ( tofpart==0 || tofpart==2 ) && useEtofScin );
00636 bool useEndcapMRPC = ( ( tofpart==0 || tofpart==2 ) && useEtofMRPC );
00637
00638 if( useBarrelScin || useEndcapScin ) {
00639 idtof = TofHit.Tofid;
00640 tofid_helix[idfztrk] = TofHit.Tofid;
00641 }
00642 if( useEndcapMRPC ) {
00643 idetf = TofHit.Etfid;
00644 etfid_helix[idfztrk] = TofHit.Etfid;
00645 }
00646
00647 log << MSG::INFO << "helix to tof hit part: "<<tofpart<<" tof id: "<< idtof << " etf id:" << idetf << endreq;
00648 if( m_debug==4 ) cout << "helix to tof hit part, Id: "<<tofpart<<" , "<< idtof <<endl;
00649 if( ( useBarrelScin && idtof>=0 && idtof<=87 ) || ( useEndcapScin && idtof>=0 && idtof<=47 ) || ( useEndcapMRPC && idetf>=0 && idetf<=35 ) ) {
00650
00651 double abmom = 0.0;
00652 if( useEndcapMRPC ) {
00653 idetfmatch[tofpart][idetf]= 1;
00654 helpathetf[idetf]= TofHit.Path_etf;
00655 helz[idetf] = TofHit.Z_etf;
00656 abmom = 1./fabs(TofHit.Kappa) * sqrt(1.+TofHit.Tanl*TofHit.Tanl);
00657 if(abmom<0.1) continue;
00658 abmom2etf[idetf] = abmom*abmom;
00659 r_endetf[idfztrk]= TofHit.r_etf;
00660 helzetf[idfztrk] = helz[idetf];
00661 }
00662
00663 if( useBarrelScin || useEndcapScin ) {
00664 idmatch[tofpart][idtof] = 1;
00665 helpath[idtof] = TofHit.Pathl;
00666 helz[idtof] = TofHit.Z_tof;
00667 abmom = 1./fabs(TofHit.Kappa) * sqrt(1.+TofHit.Tanl*TofHit.Tanl);
00668 if(abmom<0.1) continue;
00669 abmom2[idtof] = abmom*abmom;
00670 r_endtof[idfztrk]= TofHit.r_endtof;
00671 helztof[idfztrk] = helz[idtof];
00672 }
00673
00674 if( m_debug==4 ) {
00675 cout << "Scintillator info trk number=" << idfztrk << " tofpart=" << tofpart << " idtof=" << idtof << " helpath=" << helpath[idtof] << " helz=" << helz[idtof] << " abmom=" << abmom2[idtof] << " r=" << r_endtof[idfztrk] << " helztof=" << helz[idtof] << endl;
00676 cout << "MRPC info trk number=" << idfztrk << " tofpart=" << tofpart << " idetf=" << idetf << " helpath=" << helpathetf[idetf] << " helz=" << helzetf[idetf] << " abmom=" << abmom2etf[idetf] << " r=" << r_endetf[idfztrk] << " helztof=" << helzetf[idetf] << endl;
00677 }
00678
00679 double vel = 1.0e-6;
00680 if( optCosmic==0 ) {
00681
00682 if( useEndcapMRPC ) {
00683 if( particleId[idfztrk] == 1 ) {
00684 tetf[idfztrk] = sqrt(ELMAS2/abmom2etf[idetf]+1)* helpathetf[idetf]/VLIGHT;
00685 vel = VLIGHT/sqrt(ELMAS2/abmom2etf[idetf]+1);
00686 }
00687 else if( particleId[idfztrk] == 5 ) {
00688 tetf[idfztrk] = sqrt(PROTONMAS2/abmom2etf[idetf]+1)* helpathetf[idetf]/VLIGHT;
00689 vel = VLIGHT/sqrt(PROTONMAS2/abmom2etf[idtof]+1);
00690 }
00691 else {
00692 tetf[idfztrk] = sqrt(PIMAS2/abmom2etf[idetf]+1)* helpathetf[idetf]/VLIGHT;
00693 vel = VLIGHT/sqrt(PIMAS2/abmom2etf[idtof]+1);
00694 }
00695 }
00696
00697 if( useBarrelScin || useEndcapScin ) {
00698 if( particleId[idfztrk] == 1 ) {
00699 ttof[idfztrk] = sqrt(ELMAS2/abmom2[idtof]+1)* helpath[idtof]/VLIGHT;
00700 vel = VLIGHT/sqrt(ELMAS2/abmom2[idtof]+1);
00701 }
00702 else if( particleId[idfztrk] == 5 ) {
00703 ttof[idfztrk] = sqrt(PROTONMAS2/abmom2[idtof]+1)* helpath[idtof]/VLIGHT;
00704 vel = VLIGHT/sqrt(PROTONMAS2/abmom2[idtof]+1);
00705 }
00706 else {
00707 ttof[idfztrk] = sqrt(PIMAS2/abmom2[idtof]+1)* helpath[idtof]/VLIGHT;
00708 vel = VLIGHT/sqrt(PIMAS2/abmom2[idtof]+1);
00709 }
00710 }
00711
00712 }
00713 else{
00714
00715 if( useEndcapMRPC ) {
00716 tetf[idfztrk] = sqrt(MUMAS2/abmom2etf[idetf]+1)* helpathetf[idetf]/VLIGHT;
00717 vel = VLIGHT/sqrt(MUMAS2/abmom2etf[idtof]+1);
00718 }
00719
00720 if( useBarrelScin || useEndcapMRPC ) {
00721 ttof[idfztrk] = sqrt(MUMAS2/abmom2[idtof]+1)* helpath[idtof]/VLIGHT;
00722 vel = VLIGHT/sqrt(MUMAS2/abmom2[idtof]+1);
00723 }
00724 }
00725
00726 if(m_ntupleflag && m_tuple2){
00727 g_vel[idfztrk] = vel;
00728 g_abmom[idfztrk] = abmom;
00729 if( useEndcapMRPC ) {
00730 g_ttof[idfztrk] = tetf[idfztrk];
00731 }
00732 if( useBarrelScin || useEndcapScin ) {
00733 g_ttof[idfztrk] = ttof[idfztrk];
00734 }
00735 g_pid[idfztrk] = particleId[idfztrk];
00736 }
00737 }
00738 ntrk++;
00739 }
00740 if(m_ntupleflag && m_tuple2) g_ntrk= ntrk;
00741 }
00742
00743
00744
00745 if(runNo<0){
00746 SmartDataPtr<TofMcHitCol> tofmcHitCol(eventSvc(),"/Event/MC/TofMcHitCol");
00747 if (!tofmcHitCol) {
00748 log << MSG::ERROR<< "Could not find McParticle" << endreq;
00749
00750 }
00751 else{
00752 TofMcHitCol::iterator iter_mctof = tofmcHitCol->begin();
00753
00754 for (;iter_mctof !=tofmcHitCol->end(); iter_mctof++, digiId++) {
00755 log << MSG::INFO
00756 << " TofMcHit Flight Time = " << (*iter_mctof)->getFlightTime()
00757 << " zPosition = " << ((*iter_mctof)->getPositionZ())/10
00758 << " xPosition = " <<((*iter_mctof)->getPositionX())/10
00759 << " yPosition = " <<((*iter_mctof)->getPositionY())/10
00760 << endreq;
00761 }
00762 }
00763 }
00764
00766 TofDataVector tofDigiVec = m_rawDataProviderSvc->tofDataVectorEstime();
00767
00768 for(TofDataVector::iterator iter2 = tofDigiVec.begin();iter2 != tofDigiVec.end(); iter2++) {
00769 int barrelid;
00770 int layerid;
00771 int tofid;
00772 int strip;
00773
00774 if( !( (*iter2)->is_mrpc() ) ) {
00775 if( (*iter2)->barrel() ) {
00776 barrelid = 1;
00777 tofid = (*iter2)->tofId();
00778 layerid = (*iter2)->layer();
00779 if(layerid==1) tofid=tofid-88;
00780 if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()==1 ) {
00781 double ftdc = (*iter2)->tdc1();
00782 double btdc = (*iter2)->tdc2();
00783 mean_tdc_btof[layerid][tofid]=(ftdc+btdc)/2;
00784 }
00785 else if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()>1 ) {
00786 double ftdc = (*iter2)->tdc1();
00787 double btdc = (*iter2)->tdc2();
00788 mean_tdc_btof[layerid][tofid]=(ftdc+btdc)/2;
00789 }
00790 }
00791 else{
00792 tofid = (*iter2)->tofId();
00793 if(tofid<48) barrelid=0;
00794 if(tofid>47) barrelid=2;
00795 if(barrelid==2) tofid=tofid-48;
00796
00797 if((*iter2)->times()==1){
00798 double ftdc= (*iter2)->tdc();
00799 mean_tdc_etof[barrelid][tofid]=ftdc;
00800 }
00801 else if(((*iter2)->times()>1) && ((*iter2)->times()<5)){
00802 double ftdc= (*iter2)->tdc();
00803 mean_tdc_etof[barrelid][tofid]=ftdc;
00804 }
00805 }
00806 }
00807 else {
00808 tofid = (*iter2)->tofId();
00809 strip = (*iter2)->strip();
00810 if( tofid<36 ) barrelid=0;
00811 if( tofid>35 ) barrelid=2;
00812 if(barrelid==2) tofid=tofid-36;
00813 if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()==1 ) {
00814 double ftdc = (*iter2)->tdc1();
00815 double btdc = (*iter2)->tdc2();
00816 mean_tdc_etf[barrelid][tofid][strip]=(ftdc+btdc)/2;
00817 }
00818 else if( ((*iter2)->quality()&0x5)==0x5 && (*iter2)->times()>1 ) {
00819 double ftdc = (*iter2)->tdc1();
00820 double btdc = (*iter2)->tdc2();
00821 mean_tdc_etf[barrelid][tofid][strip]=(ftdc+btdc)/2;
00822 }
00823 }
00824 }
00825
00826 double difftof_b=100, difftof_e=100;
00827 int tofid1=tofid_emc[0];
00828 int tofid2=tofid_emc[1];
00829 if( module[0]==1 && module[1]==1 ) {
00830 for(int i=0; i<2; i++){
00831 for(int m=0; m<2; m++){
00832 for(int j=-2; j<3; j++){
00833 for(int k=-2; k<3; k++){
00834 int p=tofid1+j;
00835 int q=tofid2+k;
00836 if(p<0) p=p+88;
00837 if(p>87) p=p-88;
00838 if(q<0) q=q+88;
00839 if(q>87) q=q+88;
00840 if(mean_tdc_btof[i][p]==0 || mean_tdc_btof[m][q]==0) continue;
00841 double difftof_b_temp = mean_tdc_btof[i][p]-mean_tdc_btof[m][q];
00842 if(abs(difftof_b_temp)<abs(difftof_b )) difftof_b =difftof_b_temp;
00843 if(m_ntupleflag && m_tuple2) g_difftof_b=difftof_b;
00844 }
00845 }
00846 }
00847 }
00848 }
00849
00850 if( useEtofMRPC ) {
00851 if( module[0]!=1 && module[1]!=1 ) {
00852 tofid1 = etfid_emc[0];
00853 tofid2 = etfid_emc[1];
00854 for(int i=-1; i<2; i++){
00855 for(int j=-1; j<2; j++){
00856 int m=tofid1+i;
00857 int n=tofid2+j;
00858 if(m<0) m=36+m;
00859 if(m>35) m=m-36;
00860 if(n<0) n=36+n;
00861 if(n>35) n=n-36;
00862 if( mean_tdc_etf[0][m] && mean_tdc_etf[2][n]){
00863 double difftof_e_temp= mean_tdc_etf[0][m]-mean_tdc_etf[2][n];
00864 if(abs(difftof_e_temp) < abs(difftof_e)) difftof_e= difftof_e_temp;
00865 if(m_ntupleflag && m_tuple2) g_difftof_e=difftof_e;
00866 }
00867 }
00868 }
00869 }
00870 }
00871
00872 if( useEtofScin ) {
00873 if( module[0]!=1 && module[1]!=1 ) {
00874 for(int i=-1; i<2; i++){
00875 for(int j=-1; j<2; j++){
00876 int m=tofid1+i;
00877 int n=tofid2+j;
00878 if(m<0) m=48+m;
00879 if(m>47) m=m-48;
00880 if(n<0) n=48+n;
00881 if(n>47) n=n-48;
00882 if( mean_tdc_etof[0][m] && mean_tdc_etof[2][n]){
00883 double difftof_e_temp= mean_tdc_etof[0][m]-mean_tdc_etof[2][n];
00884 if(abs(difftof_e_temp) < abs(difftof_e)) difftof_e= difftof_e_temp;
00885 if(m_ntupleflag && m_tuple2) g_difftof_e=difftof_e;
00886 }
00887 }
00888 }
00889 }
00890 }
00891
00892 if(m_optCosmic==1) optCosmic=1;
00893 else optCosmic=0;
00894
00895 digiId = 0;
00896 unsigned int tofid;
00897 unsigned int barrelid;
00898 unsigned int layerid;
00899 t0forward_add = 0;
00900 t0backward_add = 0;
00901 TofDataVector::iterator iter2 = tofDigiVec.begin();
00902 for (;iter2 != tofDigiVec.end(); iter2++, digiId++){
00903 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
00904 barrelid=(*iter2)->barrel();
00905 if((*iter2)->barrel()==0) continue;
00906 if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()==1 ) {
00907 tofid = (*iter2)->tofId();
00908 layerid = (*iter2)->layer();
00909 if(layerid==1) tofid=tofid-88;
00910 log<< MSG::INFO
00911 <<" TofId = "<<tofid
00912 <<" barrelid = "<<barrelid
00913 <<" layerid = "<<layerid
00914 <<" ForwordADC = "<<(*iter2)->adc1()
00915 <<" ForwordTDC = "<<(*iter2)->tdc1()
00916 <<" BackwordADC = "<<(*iter2)->adc2()
00917 <<" BackwordTDC = "<<(*iter2)->tdc2()
00918 << endreq;
00919
00920 double ftdc = (*iter2)->tdc1();
00921 double btdc = (*iter2)->tdc2();
00922 if(m_debug==4) cout<<"barrel 1 ::layer, barrel, tofid, ftdc, btdc: "<<layerid<<" , "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<" , "<<btdc<<endl;
00923 double fadc = (*iter2)->adc1();
00924 double badc = (*iter2)->adc2();
00925 int idptof = ((tofid-1) == -1) ? 87 : tofid-1;
00926 int idntof = ((tofid+1) == 88) ? 0 : tofid+1;
00927 double ztof = fabs((ftdc-btdc)/2)*17.7 , ztof2 = ztof*ztof;
00928 double mean_tdc = 0.5*(btdc+ftdc);
00929 double cor_path = sqrt(RCTOF2+ztof2)/VLIGHT;
00930
00931 if(idmatch[barrelid][tofid]==1||idmatch[barrelid][idptof]==1||idmatch[barrelid][idntof]==1){
00932 for(int i=0;i<=ntot;i++){
00933 if(ttof[i]!=0 && ftdc>0){
00934 if((tofid_helix[i] == tofid) || (tofid_helix[i] == idntof) ||(tofid_helix[i] == idptof)) {
00935 if(barrelid==1 && helztof[i]<117 && helztof[i]>-117 ){
00936 if (optCosmic && tofid<44) {
00937 backevtime = -ttof[i] + (115 + helztof[i])*0.0566 + 12.;
00938 forevtime = -ttof[i] + (115 - helztof[i])*0.0566 + 12.;
00939 meantevup[ntofup]=(backevtime+forevtime)/2;
00940 ntofup++;
00941 }
00942 else{
00943 backevtime = ttof[i] + (115 + helztof[i])*0.0566 + 12.;
00944 forevtime = ttof[i] + (115 - helztof[i])*0.0566 + 12.;
00945 meantevdown[ntofdown]=(backevtime+forevtime)/2;
00946 ntofdown++;
00947 }
00948 if( (*iter2)->adc1()<0.0 || (*iter2)->adc2()<0.0 || m_userawtime){
00949 t0forward_trk = ftdc - forevtime ;
00950 t0backward_trk = btdc - backevtime ;
00951 }
00952 else{
00953 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())-ttof[i];
00954 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())-ttof[i];
00955 if (optCosmic && tofid<44) {
00956 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())+ttof[i];
00957 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())+ttof[i];
00958 }
00959 }
00960
00961 if(t0forward_trk<-3 || t0backward_trk<-3 || fabs(t0forward_trk-t0backward_trk)>15.0) continue;
00962 if(!m_TofOpt&& nmatch_barrel!=0 && fabs((t0forward_trk+t0backward_trk)/2-(t0backward_add+t0forward_add)/2/nmatch_barrel)>11) continue;
00963 if(m_debug ==4 ) cout<<" t0forward_trk, t0backward_trk: "<<t0forward_trk<<" , "<<t0backward_trk<<endl;
00964 if(m_ntupleflag && m_tuple2){
00965 g_t0for[nmatch1] = t0forward_trk ;
00966 g_t0back[nmatch2] = t0backward_trk ;
00967 g_meantdc=(ftdc+btdc)/2;
00968 if(tofid<44) g_ntofup1++;
00969 else g_ntofdown1++;
00970 }
00971 meantev[nmatch]=(backevtime+forevtime)/2;
00972 t0forward_add += t0forward_trk;
00973 t0backward_add += t0backward_trk;
00974 if(nmatch>499) break;
00975 Tof_t0Arr[nmatch]=(t0forward_trk+t0backward_trk)/2.0;
00976 nmatch++;
00977 nmatch1=nmatch1+1;
00978 nmatch2=nmatch2+1;
00979 nmatch_barrel++;
00980 }
00981 }
00982 }
00983 }
00984 }
00985 }
00986 }
00987
00988 if(nmatch_barrel != 0 ){
00989 if(m_ntupleflag && m_tuple2){
00990 g_t0barrelTof=( t0forward_add/nmatch_barrel + t0backward_add/nmatch_barrel)/2;
00991 }
00992 tof_flag = 1;
00993 t_quality=1;
00994 }
00995
00996 if(nmatch_barrel==0){
00997 digiId = 0;
00998 t0forward_add = 0;
00999 t0backward_add = 0;
01000 for (TofDataVector::iterator iter2 = tofDigiVec.begin();iter2 != tofDigiVec.end(); iter2++, digiId++) {
01001 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
01002 barrelid=(*iter2)->barrel();
01003 if((*iter2)->barrel()==0) continue;
01004 if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()>1 ) {
01005 tofid = (*iter2)->tofId();
01006 layerid = (*iter2)->layer();
01007 if(layerid==1) tofid=tofid-88;
01008 log<< MSG::INFO
01009 <<" TofId = "<<tofid
01010 <<" barrelid = "<<barrelid
01011 <<" layerid = "<<layerid
01012 <<" ForwordADC = "<<(*iter2)->adc1()
01013 <<" ForwordTDC = "<<(*iter2)->tdc1()
01014 <<endreq;
01015 double ftdc= (*iter2)->tdc1();
01016 double btdc= (*iter2)->tdc2();
01017 double fadc= (*iter2)->adc1();
01018 double badc= (*iter2)->adc2();
01019 if(m_debug==4) cout<<"barrel 2 ::layer, barrel, tofid, ftdc, btdc: "<<layerid<<" , "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<" , "<<btdc<<endl;
01020 int idptof = ((tofid-1) == -1) ? 87 : tofid-1;
01021 int idntof = ((tofid+1) == 88) ? 0 : tofid+1;
01022 if(idmatch[barrelid][tofid]==1||idmatch[barrelid][idptof]==1||idmatch[barrelid][idntof]==1){
01023 for(int i=0;i<=ntot;i++){
01024 if(ttof[i]!=0 && ftdc>0){
01025 if(tofid_helix[i] == tofid ||(tofid_helix[i] == idptof)||(tofid_helix[i] == idntof)){
01026 if(barrelid==1 && helztof[i]<117 && helztof[i]>-117 ){
01027 if (optCosmic && tofid<44) {
01028 backevtime = -ttof[i] + (115 + helztof[i])*0.0566 + 12.;
01029 forevtime = -ttof[i] + (115 - helztof[i])*0.0566 + 12.;
01030 meantevup[ntofup]=(backevtime+forevtime)/2;
01031 ntofup++;
01032 }
01033 else{
01034 backevtime = ttof[i] + (115 + helztof[i])*0.0566 + 12.;
01035 forevtime = ttof[i] + (115 - helztof[i])*0.0566 + 12.;
01036 meantevdown[ntofdown]=(backevtime+forevtime)/2;
01037 ntofdown++;
01038 }
01039 if( (*iter2)->adc1()<0.0 || (*iter2)->adc2()<0.0 || m_userawtime){
01040 t0forward_trk = ftdc - forevtime ;
01041 t0backward_trk = btdc - backevtime ;
01042 }
01043 else{
01044 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())-ttof[i];
01045 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())-ttof[i];
01046 if (optCosmic && tofid<44) {
01047 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())+ttof[i];
01048 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())+ttof[i];
01049 }
01050 }
01051
01052 if(t0forward_trk<-3 || t0backward_trk<-3 || fabs(t0forward_trk-t0backward_trk)>15.0) continue;
01053 if(!m_TofOpt&&nmatch_barrel!=0 && fabs((t0forward_trk+t0backward_trk)/2-(t0backward_add+t0forward_add)/2/nmatch_barrel)>11) continue;
01054 if(m_debug == 4) cout<<"t0forward_trk, t0backward_trk: "<<t0forward_trk<<" , "<<t0backward_trk<<endl;
01055 if(m_ntupleflag && m_tuple2){
01056 g_t0for[nmatch1] = t0forward_trk ;
01057 g_t0back[nmatch2] = t0backward_trk ;
01058 g_meantdc=(ftdc+btdc)/2;
01059 if(tofid<44) g_ntofup1++;
01060 else g_ntofdown1++;
01061 }
01062 meantev[nmatch]=(backevtime+forevtime)/2;
01063 t0forward_add += t0forward_trk;
01064 t0backward_add += t0backward_trk;
01065 if(nmatch>499) break;
01066 Tof_t0Arr[nmatch]=(t0forward_trk+t0backward_trk)/2.0;
01067 nmatch++;
01068 nmatch1=nmatch1+1;
01069 nmatch2=nmatch2+1;
01070 nmatch_barrel++;
01071 }
01072 }
01073 }
01074 }
01075 }
01076 }
01077 }
01078 if(nmatch_barrel) tof_flag=2;
01079 }
01080
01081 if(ntot==0 || nmatch_barrel==0) {
01082 for (TofDataVector::iterator iter2 = tofDigiVec.begin();iter2 != tofDigiVec.end(); iter2++, digiId++) {
01083 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
01084 barrelid=(*iter2)->barrel();
01085 if((*iter2)->barrel()==0) continue;
01086 if( ((*iter2)->quality() & 0x5)==0x5 && (*iter2)->times()==1 ) {
01087 tofid = (*iter2)->tofId();
01088 layerid = (*iter2)->layer();
01089 if(layerid==1) tofid=tofid-88;
01090 log<< MSG::INFO
01091 <<" TofId = "<<tofid
01092 <<" barrelid = "<<barrelid
01093 <<" layerid = "<<layerid
01094 <<" ForwordADC = "<<(*iter2)->adc1()
01095 <<" ForwordTDC = "<<(*iter2)->tdc1()
01096 <<endreq;
01097 double ftdc= (*iter2)->tdc1();
01098 double btdc= (*iter2)->tdc2();
01099 double fadc= (*iter2)->adc1();
01100 double badc= (*iter2)->adc2();
01101 if(m_debug==4) cout<<"barrel 3 ::layer, barrel, tofid, ftdc, btdc: "<<layerid<<" , "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<" , "<<btdc<<endl;
01102 int idptof = ((tofid-1) == -1) ? 87 : tofid-1;
01103 int idntof = ((tofid+1) == 88) ? 0 : tofid+1;
01104 for(int i=0; i<2; i++){
01105 if(tofid_emc[i] == tofid || tofid_emc[i] == idptof || tofid_emc[i] == idntof){
01106 if(zemc_rec[0]||zemc_rec[1]){
01107 if(tofid ==tofid_emc[i] || tofid_emc[i]==idntof || tofid_emc[i]==idptof){
01108 if(ftdc>2000.|| module[i]!=1) continue;
01109 ttof[i]= sqrt(ELMAS2/(m_ebeam*m_ebeam)+1)* sqrt(xemc_rec[i]*xemc_rec[i]+yemc_rec[i]*yemc_rec[i]+zemc_rec[i]*zemc_rec[i])/VLIGHT;
01110 if(optCosmic==1 && tofid<44){
01111 backevtime = -ttof[i] + (115 + zemc_rec[i])*0.0566 + 12.;
01112 forevtime = -ttof[i] + (115 - zemc_rec[i])*0.0566 + 12.;
01113 meantevup[ntofup]=(backevtime+forevtime)/2;
01114 ntofup++;
01115 }
01116 else{
01117 backevtime = ttof[i] + (115 + zemc_rec[i])*0.0566 + 12.;
01118 forevtime = ttof[i] + (115 - zemc_rec[i])*0.0566 + 12.;
01119 meantevdown[ntofdown]=(backevtime+forevtime)/2;
01120 ntofdown++;
01121 }
01122 if( (*iter2)->adc1()<0.0 || (*iter2)->adc2()<0.0 || m_userawtime){
01123 t0forward_trk=ftdc-forevtime;
01124 t0backward_trk=btdc-backevtime;
01125 }
01126 else{
01127 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())-ttof[i];
01128 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())-ttof[i];
01129 if (optCosmic && tofid<44) {
01130 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())+ttof[i];
01131 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())+ttof[i];
01132 }
01133 }
01134
01135 if(t0forward_trk<-1 || t0backward_trk<-1 || fabs(t0forward_trk-t0backward_trk)>15.0) continue;
01136 if(!m_TofOpt&&nmatch_barrel!=0 && fabs((t0forward_trk+t0backward_trk)/2-(t0backward_add+t0forward_add)/2/nmatch_barrel)>11) continue;
01137 if(m_debug == 4) cout<<"t0forward_trk, t0backward_trk: "<<t0forward_trk<<" , "<<t0backward_trk<<endl;
01138 meantev[nmatch]=(backevtime+forevtime)/2;
01139 t0forward_add += t0forward_trk;
01140 t0backward_add += t0backward_trk;
01141 if(nmatch>499) break;
01142 Tof_t0Arr[nmatch]=(t0forward_trk+t0backward_trk)/2.0;
01143 nmatch++;
01144 nmatch_barrel++;
01145 emcflag1=1;
01146 }
01147 }
01148 }
01149 }
01150 }
01151 }
01152 if(nmatch_barrel) tof_flag=3;
01153 }
01154
01155 if( nmatch_barrel != 0 ) {
01156 t0forward = t0forward_add/nmatch_barrel;
01157 t0backward = t0backward_add/nmatch_barrel;
01158 if(optCosmic==0){
01159 if(m_TofOpt)
01160 {
01161 t_Est=EST_Trimmer(Opt_new(Tof_t0Arr,nmatch,m_TofOpt_Cut),tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01162 }
01163 else t_Est=EST_Trimmer((t0forward+t0backward)/2,tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01164 if(t_Est<0) t_Est=0;
01165 if(tof_flag==1) tEstFlag=111;
01166 else if(tof_flag==2) tEstFlag=121;
01167 else if(tof_flag==3) tEstFlag=131;
01168 }
01169 if(optCosmic){
01170 t_Est=(t0forward+t0backward)/2;
01171 if(tof_flag==1) tEstFlag=211;
01172 else if(tof_flag==2) tEstFlag=221;
01173 else if(tof_flag==3) tEstFlag=231;
01174 }
01175 if(m_ntupleflag && m_tuple2) g_meant0=(t0forward+t0backward)/2;
01176 }
01177
01178
01179 digiId = 0;
01180 t0forward_add = 0;
01181 t0backward_add = 0;
01182
01183 if(nmatch_barrel==0){
01184 for (TofDataVector::iterator iter2 = tofDigiVec.begin();iter2 != tofDigiVec.end(); iter2++, digiId++) {
01185 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
01186 barrelid=(*iter2)->barrel();
01187 if((*iter2)->barrel()==0) continue;
01188 if(((*iter2)->quality() & 0x5) ==0x4){
01189 tofid = (*iter2)->tofId();
01190 layerid = (*iter2)->layer();
01191 if(layerid==1) tofid=tofid-88;
01192 log<< MSG::INFO
01193 <<" TofId = "<<tofid
01194 <<" barrelid = "<<barrelid
01195 <<" layerid = "<<layerid
01196 <<" ForwordADC = "<<(*iter2)->adc1()
01197 <<" ForwordTDC = "<<(*iter2)->tdc1()
01198 <<endreq;
01199
01200 double ftdc= (*iter2)->tdc1();
01201 double fadc= (*iter2)->adc1();
01202 if(m_debug==4) cout<<"barrel 4 ::layer, barrel, tofid, ftdc: "<<layerid<<" , "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<endl;
01203 int idptof = ((tofid-1) == -1) ? 87 : tofid-1;
01204 int idntof = ((tofid+1) == 88) ? 0 : tofid+1;
01205 if(idmatch[barrelid][tofid]==1||idmatch[barrelid][idptof]==1||idmatch[barrelid][idntof]==1){
01206 for(int i=0;i<=ntot;i++){
01207 if(ttof[i]!=0 && ftdc>0){
01208 if(tofid_helix[i] == tofid ||(tofid_helix[i] == idptof) || (tofid_helix[i] == idntof)){
01209 if(barrelid==1 && helztof[i]<117 && helztof[i]>-117 ){
01210 if (optCosmic && tofid<44) {
01211 forevtime = -ttof[i] + (115 - helztof[i])*0.0566 + 12.;
01212 meantevup[ntofup]=forevtime;
01213 ntofup++;
01214 }
01215 else{
01216 forevtime = ttof[i] + (115 - helztof[i])*0.0566 + 12.;
01217 meantevdown[ntofdown]=forevtime;
01218 ntofdown++;
01219 }
01220 if( (*iter2)->adc1()<0.0 || m_userawtime){
01221 t0forward_trk = ftdc - forevtime ;
01222 }
01223 else{
01224 t0forward_trk = tofCaliSvc->BTime1((*iter2)->adc1(), (*iter2)->tdc1(),helztof[i], (*iter2)->tofId())-ttof[i];
01225 }
01226
01227 if(t0forward_trk<-1) continue;
01228 if(!m_TofOpt&&nmatch_barrel_1!=0 && fabs((t0forward_trk)-(t0forward_add)/nmatch_barrel_1)>11) continue;
01229 if(m_debug == 4) cout<<"t0forward_trk: "<<t0forward_trk<<endl;
01230 if(m_ntupleflag && m_tuple2){
01231 g_t0for[nmatch1] = t0forward_trk ;
01232 g_meantdc=ftdc;
01233 if(tofid<44) g_ntofup1++;
01234 else g_ntofdown1++;
01235 }
01236 meantev[nmatch]=forevtime;
01237 t0forward_add += t0forward_trk;
01238
01239 if(nmatch>499) break;
01240 Tof_t0Arr[nmatch]=t0forward_trk;
01241 nmatch++;
01242 nmatch1++;
01243 nmatch_barrel_1++;
01244 }
01245 }
01246 }
01247 }
01248 }
01249 }
01250 else if(((*iter2)->quality() & 0x5) ==0x1){
01251 tofid = (*iter2)->tofId();
01252 layerid = (*iter2)->layer();
01253 if(layerid==1) tofid=tofid-88;
01254 log<< MSG::INFO
01255 <<" TofId = "<<tofid
01256 <<" barrelid = "<<barrelid
01257 <<" layerid = "<<layerid
01258 <<" BackwordADC = "<<(*iter2)->adc2()
01259 <<" BackwordTDC = "<<(*iter2)->tdc2()
01260 << endreq;
01261
01262 double btdc= (*iter2)->tdc2();
01263 if(m_debug==4) cout<<"barrel 5 ::layer, barrel, tofid, btdc: "<<layerid<<" , "<<barrelid<<" , "<<tofid<<" , "<<btdc<<endl;
01264 double badc= (*iter2)->adc2();
01265 int idptof = ((tofid-1) == -1) ? 87 : tofid-1;
01266 int idntof = ((tofid+1) == 88) ? 0 : tofid+1;
01267 if(idmatch[barrelid][tofid]==1||idmatch[barrelid][idptof]==1||idmatch[barrelid][idntof]==1){
01268 for(int i=0;i<=ntot;i++){
01269 if(ttof[i]!=0 && btdc>0){
01270 if((tofid_helix[i] == tofid) || (tofid_helix[i] == idntof) ||(tofid_helix[i] == idptof)){
01271 if(barrelid==1 && helztof[i]<117 && helztof[i]>-117 ){
01272 if (optCosmic && tofid<44) {
01273 backevtime = -ttof[i] + (115 + helztof[i])*0.0566 + 12.;
01274 meantevup[ntofup]=backevtime;
01275 ntofup++;
01276 }
01277 else{
01278 backevtime = ttof[i] + (115 + helztof[i])*0.0566 + 12.;
01279 meantevdown[ntofdown]=backevtime;
01280 ntofdown++;
01281 }
01282
01283 if( (*iter2)->adc2()<0.0 || m_userawtime){
01284 t0backward_trk = btdc - backevtime ;
01285 }
01286 else{
01287 t0backward_trk = tofCaliSvc->BTime2((*iter2)->adc2(), (*iter2)->tdc2(),helztof[i], (*iter2)->tofId())-ttof[i];
01288 }
01289
01290 if(t0backward_trk<-1) continue;
01291 if(!m_TofOpt&&nmatch_barrel_2!=0 && fabs((t0backward_trk)-(t0backward_add)/nmatch_barrel_2)>11) continue;
01292 if(m_debug == 4) cout<<"t0backward_trk: "<<t0backward_trk<<endl;
01293 if(m_ntupleflag && m_tuple2){
01294 g_t0back[nmatch2] = t0backward_trk ;
01295 g_meantdc=btdc;
01296 if(tofid<44) g_ntofup1++;
01297 else g_ntofdown1++;
01298 }
01299 meantev[nmatch]=backevtime;
01300 t0backward_add += t0backward_trk;
01301 if(nmatch>499) break;
01302 Tof_t0Arr[nmatch]=t0backward_trk;
01303 nmatch++;
01304 nmatch2++;
01305 nmatch_barrel_2++;
01306 }
01307 }
01308 }
01309 }
01310 }
01311 }
01312 }
01313 }
01314 if(nmatch_barrel_1 != 0 ){
01315 t0forward = t0forward_add/nmatch_barrel_1;
01316 if(optCosmic==0){
01317 if(m_TofOpt)
01318 {
01319 t_Est=EST_Trimmer(Opt_new(Tof_t0Arr,nmatch,m_TofOpt_Cut),tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01320 }
01321 else t_Est=EST_Trimmer(t0forward,tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01322 if(t_Est<0) t_Est=0;
01323 tEstFlag=141;
01324 }
01325 else{
01326 t_Est=t0forward;
01327 tEstFlag=241;
01328 }
01329 if(m_ntupleflag && m_tuple2) g_meant0=t0forward;
01330 }
01331 if(nmatch_barrel_2 != 0 ){
01332 t0backward = t0backward_add/nmatch_barrel_2;
01333 if(optCosmic==0){
01334 if(m_TofOpt)
01335 {
01336 t_Est=EST_Trimmer(Opt_new(Tof_t0Arr,nmatch,m_TofOpt_Cut),tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01337 }
01338 else t_Est=EST_Trimmer(t0backward,tOffset_b,toffset_rawtime,offset_dt,bunchtime);
01339 if(t_Est<0) t_Est=0;
01340 tEstFlag=141;
01341 }
01342 else{
01343 t_Est=t0backward;
01344 tEstFlag=241;
01345 }
01346 if(m_ntupleflag && m_tuple2) g_meant0=t0backward;
01347 }
01348
01349 digiId = 0;
01350 t0forward_add = 0;
01351 t0backward_add = 0;
01352 if(nmatch_barrel==0){
01353 for (TofDataVector::iterator iter2 = tofDigiVec.begin();iter2 != tofDigiVec.end(); iter2++, digiId++) {
01354 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
01355 barrelid=(*iter2)->barrel();
01356 if((*iter2)->barrel()!=0) continue;
01357 if((*iter2)->times()!=1) continue;
01358 tofid = (*iter2)->tofId();
01359
01360 if( !( (*iter2)->is_mrpc() ) ) {
01361 if( tofid<48 ) { barrelid=0; }
01362 if( tofid>47 ) { barrelid=2; }
01363 if( barrelid==2 ) { tofid=tofid-48; }
01364 }
01365
01366 else if( (*iter2)->is_mrpc() ) {
01367 if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 0 ) { barrelid=0; }
01368 else if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 1 ) { barrelid=2; }
01369 if( barrelid==2 ) { tofid=tofid-36; }
01370 }
01371
01372 log<< MSG::INFO
01373 <<" is_mrpc = " << (*iter2)->is_mrpc()
01374 <<" TofId = "<< tofid
01375 <<" barrelid = "<< barrelid
01376 <<endreq
01377 <<" ForwordADC = "<< (*iter2)->adc()
01378 <<" ForwordTDC = "<< (*iter2)->tdc()
01379 << endreq;
01380 double ftdc= (*iter2)->tdc();
01381 double fadc= (*iter2)->adc();
01382 if(m_debug ==4) cout<<"endcap::single hit,barrelid,tofid,tdc: "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<endl;
01383
01384
01385 if( !( (*iter2)->is_mrpc() ) && useEtofScin ) {
01386 int idptof = ((tofid-1) == -1) ? 47 : tofid-1;
01387 int idntof = ((tofid+1) == 48) ? 0 : tofid+1;
01388
01389 if(idmatch[barrelid][tofid]==1||idmatch[barrelid][idptof]==1||idmatch[barrelid][idntof]==1){
01390 for(int i=0;i<=ntot;i++){
01391 if(ttof[i]!=0 && ftdc>0 && ftdc<2000.){
01392 if((tofid_helix[i] == tofid) ||(tofid_helix[i] == idntof) ||(tofid_helix[i] == idptof)){
01393 if( barrelid==0 || barrelid==2 ){
01394 if( r_endtof[i]>=41 && r_endtof[i]<=90 ) {
01395 if( optCosmic && ( tofid<24 || ( tofid>48 && tofid<71 ) ) ) {
01396 forevtime = -ttof[i] + r_endtof[i]*0.09 + 12.2;
01397 meantevup[ntofup] = forevtime;
01398 ntofup++;
01399 }
01400 else {
01401 forevtime = ttof[i] + r_endtof[i]*0.09 + 12.2;
01402 meantevdown[ntofdown] = forevtime;
01403 ntofdown++;
01404 }
01405 if( (*iter2)->adc()<0.0 || m_userawtime){
01406 t0forward_trk = ftdc - forevtime;
01407 }
01408 else{
01409 t0forward_trk = tofCaliSvc->ETime((*iter2)->adc(), (*iter2)->tdc(),r_endtof[i], (*iter2)->tofId())-ttof[i];
01410 }
01411
01412 if(t0forward_trk<-1.) continue;
01413 if( !m_TofOpt && nmatch_end!=0 && fabs( t0forward_trk - t0forward_add/nmatch_end)>9 ) continue;
01414 t0forward_add += t0forward_trk;
01415 if(nmatch>499) break;
01416 Tof_t0Arr[nmatch]=t0forward_trk;
01417 meantev[nmatch]=forevtime/2;
01418 nmatch++;
01419 nmatch_end++;
01420 }
01421 }
01422 }
01423 if( m_debug==4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
01424 }
01425 }
01426 }
01427 }
01428
01429 if( (*iter2)->is_mrpc() && useEtofMRPC ) {
01430 if( ((*iter2)->quality() & 0x5)!=0x5 ) continue;
01431 double btdc= (*iter2)->tdc2();
01432 double badc= (*iter2)->adc2();
01433 int idptof = ((tofid-1) == -1) ? 35 : tofid-1;
01434 int idntof = ((tofid+1) == 36) ? 0 : tofid+1;
01435
01436 if( idetfmatch[barrelid][tofid]==1 || idetfmatch[barrelid][idptof]==1 || idetfmatch[barrelid][idntof]==1 ) {
01437 for( int i=0; i<=ntot; i++ ) {
01438 if( tetf[i]!=0 && ftdc>0 && ftdc<2000.) {
01439 if( etfid_helix[i]==tofid || etfid_helix[i]==idntof || etfid_helix[i] == idptof ) {
01440 if( barrelid==0 || barrelid==2 ) {
01441 if( r_endetf[i]>=41 && r_endetf[i]<=90 ) {
01442 if( optCosmic && ( tofid<18 || ( tofid>35 && tofid<54 ) ) ) {
01443 forevtime = -tetf[i] + 17.7;
01444 meantevup[ntofup] = forevtime;
01445 ntofup++;
01446 }
01447 else {
01448 forevtime = tetf[i] + 17.7;
01449 meantevdown[ntofdown] = forevtime;
01450 ntofdown++;
01451 }
01452 if( m_userawtime ) {
01453 double fbtdc = ( ftdc + btdc )/2.0;
01454 if( ftdc>0 && btdc<0 ) { fbtdc = ftdc; }
01455 else if( ftdc<0 && btdc>0 ) { fbtdc = btdc; }
01456 else if( ftdc<0 && btdc<0 ) continue;
01457 t0forward_trk = fbtdc - forevtime;
01458 }
01459 else {
01460 t0forward_trk = tofCaliSvc->EtfTime( (*iter2)->tdc1(), (*iter2)->tdc2(), (*iter2)->tofId(), (*iter2)->strip() )-tetf[i];
01461 }
01462
01463 if( t0forward_trk<-1 ) continue;
01464 if( m_TofOpt && nmatch_end!=0 && fabs(t0forward_trk-t0forward_add/nmatch_end)>9 ) continue;
01465 if( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
01466 t0forward_add += t0forward_trk;
01467 if(nmatch>499) break;
01468 Tof_t0Arr[nmatch] = t0forward_trk;
01469 meantev[nmatch] = forevtime;
01470 nmatch++;
01471 nmatch_end++;
01472 }
01473 }
01474 }
01475 }
01476 }
01477 }
01478 }
01479 }
01480 if( nmatch_end ) { tof_flag=5; }
01481 }
01482
01483 if( nmatch_barrel==0 && nmatch_end==0 ) {
01484 for( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end(); iter2++, digiId++ ) {
01485 barrelid=(*iter2)->barrel();
01486 if( (*iter2)->barrel()!=0 ) continue;
01487 if( (*iter2)->times()!=1 ) continue;
01488 tofid = (*iter2)->tofId();
01489
01490 if( !( (*iter2)->is_mrpc() ) ) {
01491 if( tofid<48 ) { barrelid=0; }
01492 if( tofid>47 ) { barrelid=2; }
01493 if( barrelid==2 ) { tofid=tofid-48; }
01494 }
01495
01496 else if( (*iter2)->is_mrpc() ) {
01497 if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 0 ) { barrelid=0; }
01498 else if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 1 ) { barrelid=2; }
01499 if( barrelid==2 ) { tofid=tofid-36; }
01500 }
01501
01502 log<< MSG::INFO
01503 <<" is_mrpc = " << (*iter2)->is_mrpc()
01504 <<" TofId = "<< tofid
01505 <<" barrelid = "<< barrelid
01506 <<endreq
01507 <<" ForwordADC = "<< (*iter2)->adc()
01508 <<" ForwordTDC = "<< (*iter2)->tdc()
01509 << endreq;
01510 double ftdc= (*iter2)->tdc();
01511 double fadc= (*iter2)->adc();
01512 if(m_debug ==4) cout<<"endcap::single hit,barrelid,tofid,tdc: "<<barrelid<<" , "<<tofid<<" , "<<ftdc<<endl;
01513
01514
01515 if( !( (*iter2)->is_mrpc() ) && useEtofScin ) {
01516 int idptof = ((tofid-1) == -1) ? 47 : tofid-1;
01517 int idntof = ((tofid+1) == 48) ? 0 : tofid+1;
01518 for( int i=0; i<2; i++ ) {
01519 if( zemc_rec[0] || zemc_rec[1] ) {
01520 if( tofid==tofid_emc[i] || tofid_emc[i]==idntof || tofid_emc[i]==idptof ) {
01521 if( ftdc>2000. || module[i]==1 ) continue;
01522 ttof[i]= sqrt(ELMAS2/(m_ebeam*m_ebeam)+1)*sqrt(xemc_rec[i]*xemc_rec[i]+yemc_rec[i]*yemc_rec[i]+133*133)/VLIGHT;
01523 r_endtof[i]=sqrt(xemc_rec[i]*xemc_rec[i]+yemc_rec[i]*yemc_rec[i]);
01524 if( optCosmic && ( tofid<24 || ( tofid>48 && tofid<71 ) ) ) {
01525 forevtime = -ttof[i] + r_endtof[i]*0.09 + 12.2;
01526 meantevup[ntofup] = forevtime;
01527 ntofup++;
01528 }
01529 else {
01530 forevtime = ttof[i] + r_endtof[i]*0.09 + 12.2;
01531 meantevdown[ntofdown] = forevtime;
01532 ntofdown++;
01533 }
01534 if( (*iter2)->adc()<0.0 || m_userawtime){
01535 t0forward_trk = ftdc - forevtime;
01536 }
01537 else{
01538 t0forward_trk = tofCaliSvc->ETime((*iter2)->adc(), (*iter2)->tdc(),r_endtof[i], (*iter2)->tofId())-ttof[i];
01539 if(m_debug ==4) cout<<"emc flag t0forward_trk: "<<t0forward_trk<<endl;
01540 }
01541
01542 if( t0forward_trk<-1.) continue;
01543 if( !m_TofOpt && nmatch_end!=0 && fabs( t0forward_trk - t0forward_add/nmatch_end)>9 ) continue;
01544 meantev[nmatch] = forevtime;
01545 t0forward_add += t0forward_trk;
01546 if(nmatch>499) break;
01547 Tof_t0Arr[nmatch] = t0forward_trk;
01548 nmatch++;
01549 nmatch_end++;
01550 emcflag2=1;
01551 }
01552 }
01553 }
01554 }
01555
01556 if( (*iter2)->is_mrpc() && useEtofMRPC ) {
01557 double btdc= (*iter2)->tdc2();
01558 double badc= (*iter2)->adc2();
01559 int idptof = ((tofid-1) == -1) ? 35 : tofid-1;
01560 int idntof = ((tofid+1) == 36) ? 0 : tofid+1;
01561 for( int i=0; i<2; i++ ) {
01562 if( zemc_rec[0] || zemc_rec[1] ) {
01563 if( tofid==etfid_emc[i] || etfid_emc[i]==idntof || etfid_emc[i]==idptof ) {
01564
01565 if( ftdc>2000.|| module[i]==1 ) continue;
01566 tetf[i]= sqrt(ELMAS2/(m_ebeam*m_ebeam)+1)* sqrt(xemc_rec[i]*xemc_rec[i]+yemc_rec[i]*yemc_rec[i]+133*133)/VLIGHT;
01567 r_endetf[i] = sqrt(xemc_rec[i]*xemc_rec[i]+yemc_rec[i]*yemc_rec[i]);
01568 if( optCosmic && ( tofid<18 || ( tofid>35 && tofid<54 ) ) ) {
01569 forevtime = -tetf[i] + 17.7;
01570 meantevup[ntofup] = forevtime;
01571 ntofup++;
01572 }
01573 else {
01574 forevtime = tetf[i] + 17.7;
01575 meantevdown[ntofdown] = forevtime;
01576 ntofdown++;
01577 }
01578
01579 if( m_userawtime ) {
01580 double fbtdc = ( ftdc + btdc )/2.0;
01581 if( ftdc>0 && btdc<0 ) { fbtdc = ftdc; }
01582 else if( ftdc<0 && btdc>0 ) { fbtdc = btdc; }
01583 else if( ftdc<0 && btdc<0 ) continue;
01584 t0forward_trk = fbtdc - forevtime;
01585 }
01586 else {
01587 t0forward_trk = tofCaliSvc->EtfTime( (*iter2)->tdc1(), (*iter2)->tdc2(), (*iter2)->tofId(), (*iter2)->strip() )-tetf[i];
01588 }
01589
01590 if( t0forward_trk<-1 ) continue;
01591 if( !m_TofOpt && nmatch_end!=0 && fabs( t0forward_trk - t0forward_add/nmatch_end)>9 ) continue;
01592 if( m_debug==4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
01593 t0forward_add += t0forward_trk;
01594 if(nmatch>499) break;
01595 Tof_t0Arr[nmatch]=t0forward_trk;
01596 nmatch++;
01597 nmatch_end++;
01598 emcflag2=1;
01599 }
01600 }
01601 }
01602 }
01603 }
01604 if( nmatch_end ) { tof_flag=5; }
01605 }
01606
01607 if( nmatch_barrel==0 && nmatch_end==0 ) {
01608 for( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end(); iter2++, digiId++ ) {
01609 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
01610 barrelid = (*iter2)->barrel();
01611 if( (*iter2)->barrel()!=0 ) continue;
01612 if( (*iter2)->times()>1 && (*iter2)->times()<5 ) {
01613 tofid = (*iter2)->tofId();
01614
01615 if( !( (*iter2)->is_mrpc() ) ) {
01616 if( tofid<48 ) { barrelid=0; }
01617 if( tofid>47 ) { barrelid=2; }
01618 if( barrelid==2 ) { tofid=tofid-48; }
01619 }
01620
01621 else if( (*iter2)->is_mrpc() ) {
01622 if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 0 ) { barrelid=0; }
01623 else if( ( TofID::endcap( Identifier( (*iter2)->identify() ) ) ) == 1 ) { barrelid=2; }
01624 if( barrelid==2 ) { tofid=tofid-36; }
01625 }
01626 log<< MSG::INFO
01627 <<" TofId = "<<tofid
01628 <<" barrelid = "<<barrelid
01629 <<endreq
01630 <<" ForwordADC = "<< (*iter2)->adc()
01631 <<" ForwordTDC = "<< (*iter2)->tdc()
01632 << endreq;
01633 double ftdc = (*iter2)->tdc();
01634 double fadc = (*iter2)->adc();
01635 if( m_debug==4 ) { cout << "endcap::multi hit,barrelid,tofid,tdc: " << barrelid << " , " << tofid << " , " << ftdc << endl; }
01636
01637
01638 if( !( (*iter2)->is_mrpc() ) && useEtofScin ) {
01639 int idptof = ((tofid-1) == -1) ? 47 : tofid-1;
01640 int idntof = ((tofid+1) == 48) ? 0 : tofid+1;
01641
01642 if( idmatch[barrelid][tofid]==1 || idmatch[barrelid][idptof]==1 || idmatch[barrelid][idntof]==1 ) {
01643 for( int i=0; i<=ntot; i++ ) {
01644 if( ttof[i]!=0 && ftdc>0 ) {
01645 if( (tofid_helix[i]==tofid) || (tofid_helix[i]==idntof) || (tofid_helix[i]==idptof) ) {
01646 if( barrelid==0 || barrelid==2 ) {
01647 if( r_endtof[i]>=41 && r_endtof[i]<=90 ) {
01648 if( optCosmic && ( tofid<24 || ( tofid>48 && tofid<71 ) ) ) {
01649 forevtime = -ttof[i] + r_endtof[i]*0.09 + 12.2;
01650 meantevup[ntofup]=forevtime;
01651 ntofup++;
01652 }
01653 else{
01654 forevtime = ttof[i] + r_endtof[i]*0.09 + 12.2;
01655 meantevdown[ntofdown]=forevtime;
01656 ntofdown++;
01657 }
01658 if( (*iter2)->adc()<0.0 || m_userawtime){
01659 t0forward_trk=ftdc-forevtime;
01660 }
01661 else{
01662 t0forward_trk = tofCaliSvc->ETime((*iter2)->adc(), (*iter2)->tdc(),r_endtof[i], (*iter2)->tofId())-ttof[i];
01663 }
01664
01665 if( t0forward_trk<-1.) continue;
01666 if( !m_TofOpt && nmatch_end!=0 && fabs( t0forward_trk - t0forward_add/nmatch_end)>9 ) continue;
01667 meantev[nmatch] = forevtime;
01668 t0forward_add += t0forward_trk;
01669 if( nmatch>499 ) break;
01670 Tof_t0Arr[nmatch] = t0forward_trk;
01671 nmatch++;
01672 nmatch_end++;
01673 }
01674 }
01675 if( m_debug==4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
01676 }
01677 }
01678 }
01679 }
01680 }
01681
01682 if( (*iter2)->is_mrpc() && useEtofMRPC ) {
01683 double btdc= (*iter2)->tdc2();
01684 double badc= (*iter2)->adc2();
01685 int idptof = ((tofid-1) == -1) ? 35 : tofid-1;
01686 int idntof = ((tofid+1) == 36) ? 0 : tofid+1;
01687
01688 if( idetfmatch[barrelid][tofid]==1 || idetfmatch[barrelid][idptof]==1 || idetfmatch[barrelid][idntof]==1 ) {
01689 for( int i=0; i<=ntot; i++ ) {
01690 if( tetf[i]!=0 && ftdc>0 && ftdc<2000.) {
01691 if( etfid_helix[i]==tofid || etfid_helix[i]==idntof || etfid_helix[i] == idptof ) {
01692 if( barrelid==0 || barrelid==2 ) {
01693 if( r_endetf[i]>=41 && r_endetf[i]<=90 ) {
01694 if( optCosmic && ( tofid<18 || ( tofid>35 && tofid<54 ) ) ) {
01695 forevtime = -tetf[i] + 17.7;
01696 meantevup[ntofup] = forevtime;
01697 ntofup++;
01698 }
01699 else {
01700 forevtime = tetf[i] + 17.7;
01701 meantevdown[ntofdown] = forevtime;
01702 ntofdown++;
01703 }
01704 if( m_userawtime ) {
01705 double fbtdc = ( ftdc + btdc )/2.0;
01706 if( ftdc>0 && btdc<0 ) { fbtdc = ftdc; }
01707 else if( ftdc<0 && btdc>0 ) { fbtdc = btdc; }
01708 else if( ftdc<0 && btdc<0 ) continue;
01709 t0forward_trk = fbtdc - forevtime;
01710 }
01711 else {
01712 t0forward_trk = tofCaliSvc->EtfTime( (*iter2)->tdc1(), (*iter2)->tdc2(), (*iter2)->tofId(), (*iter2)->strip() )-tetf[i];
01713 }
01714
01715 if( t0forward_trk<-1 ) continue;
01716 if( !m_TofOpt && nmatch_end!=0 && fabs( t0forward_trk - t0forward_add/nmatch_end )>9 ) continue;
01717 if( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
01718 t0forward_add += t0forward_trk;
01719 if(nmatch>499) break;
01720 Tof_t0Arr[nmatch] = t0forward_trk;
01721 meantev[nmatch] = forevtime;
01722 nmatch++;
01723 nmatch_end++;
01724 }
01725 }
01726 }
01727 }
01728 }
01729 }
01730 }
01731 }
01732 }
01733 if( nmatch_end ) { tof_flag=7; }
01734 }
01735
01736 if(m_ntupleflag && m_tuple2){
01737 g_nmatchbarrel = nmatch_barrel;
01738 g_nmatchbarrel_1 = nmatch_barrel_1;
01739 g_nmatchbarrel_2 = nmatch_barrel_2;
01740 g_nmatchend = nmatch_end;
01741 }
01742
01743 if( nmatch_end !=0 ) {
01744 t0forward = t0forward_add/nmatch_end;
01745 if( optCosmic==0 ) {
01746 if( m_TofOpt ) {
01747 t_Est=EST_Trimmer(Opt_new(Tof_t0Arr,nmatch,m_TofOpt_Cut),tOffset_e,toffset_rawtime_e,offset_dt_end,bunchtime);
01748 }
01749 else { t_Est=EST_Trimmer(t0forward,tOffset_e,toffset_rawtime_e,offset_dt_end,bunchtime); }
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767 if(t_Est<0) t_Est=0;
01768 if(tof_flag==5) tEstFlag=151;
01769 else if(tof_flag==7) tEstFlag=171;
01770 if(emcflag2==1) tEstFlag=161;
01771
01772
01773
01774
01775
01776
01777
01778
01779 }
01780 if( optCosmic ) {
01781 t_Est=t0forward;
01782 if(tof_flag==5) tEstFlag=251;
01783 else if(tof_flag==7) tEstFlag=271;
01784 if(emcflag2==1) tEstFlag=261;
01785 }
01786 if(m_ntupleflag && m_tuple2) g_meant0=t0forward;
01787 }
01788
01789 double t0_comp=-999;
01790 double T0=-999;
01791
01792 if(nmatch_barrel==0 && nmatch_end==0 && m_flag==1){
01793 double mhit[43][300]={0.};
01794 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
01795 if (!mdcDigiCol) {
01796 log << MSG::INFO<< "Could not find MDC digi" << endreq;
01797 return StatusCode::FAILURE;
01798 }
01799
01800 IMdcGeomSvc* mdcGeomSvc;
01801 StatusCode sc = service("MdcGeomSvc", mdcGeomSvc);
01802 if (sc != StatusCode::SUCCESS) {
01803 return StatusCode::FAILURE;
01804 }
01805
01806 MdcDigiCol::iterator iter1 = mdcDigiCol->begin();
01807 digiId = 0;
01808 Identifier mdcId;
01809 int layerId;
01810 int wireId;
01811
01812 for (;iter1 != mdcDigiCol->end(); iter1++, digiId++) {
01813 mdcId = (*iter1)->identify();
01814 layerId = MdcID::layer(mdcId);
01815 wireId = MdcID::wire(mdcId);
01816
01817 mhit[layerId][wireId]=RawDataUtil::MdcTime((*iter1)->getTimeChannel());
01818 mhit[layerId][wireId]-=1.28*(mdcGeomSvc->Layer(layerId)->Radius())/299.8;
01819
01820 mdcGeomSvc->Wire(layerId,wireId);
01821
01822 double tof;
01823
01824 }
01825
01826 int Iused[43][300]={0},tmp=0;
01827 bool Lpat,Lpat11,Lpat12,Lpat2,Lpat31,Lpat32;
01828 double t0_all=0,t0_mean=0;
01829 double r[4]={0.};
01830 double chi2=999.0;
01831 double phi[4]={0.},corr[4]={0.},driftt[4]={0.};
01832 double ambig=1;
01833 double mchisq=50000;
01834
01835
01836 for(int i=5;i<10;i++){
01837
01838 double T1=0.5*0.1*(mdcGeomSvc->Layer(4*i+0)->PCSiz())/0.004;
01839 double T2=0.5*0.1*(mdcGeomSvc->Layer(4*i+1)->PCSiz())/0.004;
01840 double T3=0.5*0.1*(mdcGeomSvc->Layer(4*i+2)->PCSiz())/0.004;
01841 double T4=0.5*0.1*(mdcGeomSvc->Layer(4*i+3)->PCSiz())/0.004;
01842 r[0]=(mdcGeomSvc->Layer(4*i+0)->Radius())*0.1;
01843 r[1]=(mdcGeomSvc->Layer(4*i+1)->Radius())*0.1;
01844 r[2]=(mdcGeomSvc->Layer(4*i+2)->Radius())*0.1;
01845 r[3]=(mdcGeomSvc->Layer(4*i+3)->Radius())*0.1;
01846 double r0=r[0]-r[1]-(r[2]-r[1])/2;
01847 double r1=-(r[2]-r[1])/2;
01848 double r2=(r[2]-r[1])/2;
01849 double r3=r[3]-r[2]+(r[2]-r[1])/2;
01850
01851 for(int j=0;j<mdcGeomSvc->Layer(i*4+3)->NCell();j++){
01852
01853 int Icp=0;
01854 Icp=j-1;
01855 if(Icp<0) Icp=mdcGeomSvc->Layer(i*4+3)->NCell();
01856
01857 Lpat=(mhit[4*i][j]!=0) && (mhit[4*i][Icp]==0) &&( mhit[4*i][j+1]==0) && (Iused[4*i][j]==0);
01858 if(Lpat==1){
01859
01860 }
01861 if(Lpat){
01862 Lpat11=(mhit[4*i+1][Icp]==0)&&(Iused[4*i+1][j]==0)&&(mhit[4*i+1][j]!=0)&&(mhit[4*i+1][j+1]==0);
01863 Lpat12=(mhit[4*i+1][j]==0)&&(Iused[4*i+1][j+1]==0)&&(mhit[4*i+1][j+1]!=0)&&(mhit[4*i+1][j+2]==0);
01864 Lpat2=(mhit[4*i+2][Icp]==0)&&(Iused[4*i+2][j]==0)&&(mhit[4*i+2][j]!=0)&&(mhit[4*i+2][j+1]==0);
01865 Lpat31=(mhit[4*i+3][Icp]==0)&&(Iused[4*i+3][j]==0)&&(mhit[4*i+3][j]!=0)&&(mhit[4*i+3][j+1]==0);
01866 Lpat32=(mhit[4*i+3][j]==0)&&(Iused[4*i+3][j+1]==0)&&(mhit[4*i+3][j+1]!=0)&&(mhit[4*i+3][j+2]==0);
01867
01868 if(Lpat11 && Lpat2 && Lpat31 ){
01869
01870 Iused[4*i+0][j]=1;
01871 Iused[4*i+1][j]=1;
01872 Iused[4*i+2][j]=1;
01873 Iused[4*i+3][j]=1;
01874 double t_i = mhit[4*i+0][j]+mhit[4*i+2][j];
01875 double t_j = mhit[4*i+1][j]+mhit[4*i+3][j];
01876 double l_j = T2+T4;
01877 double r_i = r0+r2;
01878 double r_j = r1+r3;
01879 double r_2k= r0*r0+r1*r1+r2*r2+r3*r3;
01880 double rt_i = r0*mhit[4*i+0][j]+r2*mhit[4*i+2][j];
01881 double rt_j = r1*mhit[4*i+1][j]+r3*mhit[4*i+3][j];
01882 double rl_j = r1*T2+r3*T4;
01883
01884 double deno= 4*r_2k-2*(r_i*r_i+r_j*r_j);
01885
01886 if (deno!=0){
01887 double Pa=(4*(rt_i-rt_j+rl_j)-(t_i+t_j-l_j)*(r_i-r_j)-(t_i-t_j+l_j)*(r_i+r_j))/deno;
01888 double Pb= 0.25*(t_i-t_j+l_j-(r_i+r_j) * Pa);
01889 double Ang=fabs(90.-fabs(atan(Pa)*180./3.141593));
01890
01891 t0_all+= (-0.25*((r_i-r_j)*Pa-t_i-t_j+l_j));
01892
01893 double chi2_tmp;
01894 for(int t0c=0;t0c<17;t0c+=8){
01895 chi2_tmp=(mhit[4*i+0][j]-t0c-r0 * Pa-Pb)*(mhit[4*i+0][j]-t0c-r0 * Pa-Pb)+(T2-mhit[4*i+1][j]+t0c-r1 * Pa-Pb)*(T2-mhit[4*i+1][j]+t0c-r1 * Pa-Pb)+(mhit[4*i+2][j]-t0c-r2 * Pa-Pb)*(mhit[4*i+2][j]-t0c-r2 * Pa-Pb) + (T4-mhit[4*i+3][j]+t0c-r3 * Pa-Pb)*(T4-mhit[4*i+3][j]+t0c-r3 * Pa-Pb);
01896 if(chi2_tmp<chi2){
01897 chi2=chi2_tmp;
01898 t0_comp=t0c;
01899 }
01900 }
01901 tmp+=1;
01902 }
01903
01904
01905 for(int tmpT0=0;tmpT0<17;tmpT0+=8){
01906 driftt[0]=mhit[4*i+0][j]-tmpT0;
01907 driftt[1]=mhit[4*i+1][j]-tmpT0;
01908 driftt[2]=mhit[4*i+2][j]-tmpT0;
01909 driftt[3]=mhit[4*i+3][j]-tmpT0;
01910
01911 phi[0]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4)->NCell())+2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell())/2;
01912 phi[1]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell());
01913 phi[2]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+2)->NCell())+2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell())/2;
01914 phi[3]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+3)->NCell());
01915 phi[0]-=ambig*driftt[0]*0.004/r[0];
01916 phi[1]+=ambig*driftt[1]*0.004/r[1];
01917 phi[2]-=ambig*driftt[2]*0.004/r[2];
01918 phi[3]+=ambig*driftt[3]*0.004/r[3];
01919 double s1, sx, sy, sxx, sxy;
01920 double delinv, denom;
01921 double weight;
01922 double sigma;
01923 double x[4];
01924 x[0]=r0;
01925 x[1]=r1;
01926 x[2]=r2;
01927 x[3]=r3;
01928 sigma=0.12;
01929 s1 = sx = sy = sxx = sxy = 0.0;
01930 double chisq = 0.0;
01931 for (int ihit = 0; ihit < 4; ihit++) {
01932 weight = 1. / (sigma * sigma);
01933 s1 += weight;
01934 sx += x[ihit] * weight;
01935 sy += phi[ihit] * weight;
01936 sxx += x[ihit] * (x[ihit] * weight);
01937 sxy += phi[ihit] * (x[ihit] * weight);
01938 }
01939 double resid[4]={0.};
01940
01941 denom = s1 * sxx - sx * sx;
01942 delinv = (denom == 0.0) ? 1.e20 : 1. / denom;
01943 double intercept = (sy * sxx - sx * sxy) * delinv;
01944 double slope = (s1 * sxy - sx * sy) * delinv;
01945
01946
01947 for (int ihit = 0; ihit < 4; ihit++) {
01948 resid[ihit] = ( phi[ihit] - intercept - slope * x[ihit] );
01949 chisq += resid[ihit] * resid[ihit]/(sigma*sigma) ;
01950 }
01951 if(chisq<mchisq){
01952 mchisq=chisq;
01953 T0=tmpT0;
01954 }
01955 }
01956 }
01957 if(Lpat12 && Lpat2 && Lpat32){
01958 Iused[4*i+0][j]=1;
01959 Iused[4*i+1][j+1]=1;
01960 Iused[4*i+2][j]=1;
01961 Iused[4*i+3][j+1]=1;
01962
01963 double t_i = mhit[4*i+0][j]+mhit[4*i+2][j];
01964 double t_j = mhit[4*i+1][j+1]+mhit[4*i+3][j+1];
01965 double l_j = T2+T4;
01966 double r_i = r0+r2;
01967 double r_j = r1+r3;
01968 double r_2k= r0*r0+r1*r1+r2*r2+r3*r3;
01969 double rt_i = r0*mhit[4*i+0][j]+r2*mhit[4*i+2][j];
01970 double rt_j = r1*mhit[4*i+1][j+1]+r3*mhit[4*i+3][j+1];
01971 double rl_j = r1*T2+r3*T4;
01972 double deno= 4*r_2k-2*(r_i*r_i+r_j*r_j);
01973
01974 if (deno!=0){
01975 double Pa=(4*(rt_i-rt_j+rl_j)-(t_i+t_j-l_j)*(r_i-r_j)-(t_i-t_j+l_j)*(r_i+r_j))/deno;
01976 double Pb= 0.25*(t_i-t_j+l_j-(r_i+r_j) * Pa);
01977 double Ang=fabs(90.-fabs(atan(Pa)*180./3.141593));
01978 t0_all+= (-0.25*((r_i-r_j)*Pa-t_i-t_j+l_j));
01979 tmp+=1;
01980 double chi2_tmp;
01981
01982 for(int t0c=0;t0c<17;t0c+=8){
01983 chi2_tmp=(mhit[4*i+0][j]-t0c-r0 * Pa-Pb)*(mhit[4*i+0][j]-t0c-r0 * Pa-Pb)+(T2-mhit[4*i+1][j+1]+t0c-r1 * Pa-Pb)*(T2-mhit[4*i+1][j+1]+t0c-r1 * Pa-Pb)+(mhit[4*i+2][j]-t0c-r2 * Pa-Pb)*(mhit[4*i+2][j]-t0c-r2 * Pa-Pb) + (T4-mhit[4*i+3][j+1]+t0c-r3 * Pa-Pb)*(T4-mhit[4*i+3][j+1]+t0c-r3 * Pa-Pb);
01984
01985 if(chi2_tmp<chi2){
01986 chi2=chi2_tmp;
01987 t0_comp=t0c;
01988 }
01989 }
01990 }
01991
01992
01993
01994 for(int tmpT0=0;tmpT0<17;tmpT0+=8){
01995 driftt[0]=mhit[4*i+0][j]-tmpT0;
01996 driftt[1]=mhit[4*i+1][j+1]-tmpT0;
01997 driftt[2]=mhit[4*i+2][j]-tmpT0;
01998 driftt[3]=mhit[4*i+3][j+1]-tmpT0;
01999
02000 phi[0]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4)->NCell())+2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell())/2;
02001 phi[1]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell());
02002 phi[2]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+2)->NCell())+2*3.14159265/(mdcGeomSvc->Layer(i*4+1)->NCell())/2;
02003 phi[3]=j*2*3.14159265/(mdcGeomSvc->Layer(i*4+3)->NCell());
02004 phi[0]-=ambig*driftt[0]*0.004/r[0];
02005 phi[1]+=ambig*driftt[1]*0.004/r[1];
02006 phi[2]-=ambig*driftt[2]*0.004/r[2];
02007 phi[3]+=ambig*driftt[3]*0.004/r[3];
02008 double s1, sx, sy, sxx, sxy;
02009 double delinv, denom;
02010 double weight;
02011 double sigma;
02012 double x[4];
02013 x[0]=r0;
02014 x[1]=r1;
02015 x[2]=r2;
02016 x[3]=r3;
02017 sigma=0.12;
02018 s1 = sx = sy = sxx = sxy = 0.0;
02019 double chisq = 0.0;
02020 for (int ihit = 0; ihit < 4; ihit++) {
02021 weight = 1. / (sigma * sigma);
02022 s1 += weight;
02023 sx += x[ihit] * weight;
02024 sy += phi[ihit] * weight;
02025 sxx += x[ihit] * (x[ihit] * weight);
02026 sxy += phi[ihit] * (x[ihit] * weight);
02027 }
02028 double resid[4]={0.};
02029
02030 denom = s1 * sxx - sx * sx;
02031 delinv = (denom == 0.0) ? 1.e20 : 1. / denom;
02032 double intercept = (sy * sxx - sx * sxy) * delinv;
02033 double slope = (s1 * sxy - sx * sy) * delinv;
02034
02035
02036 for (int ihit = 0; ihit < 4; ihit++) {
02037 resid[ihit] = ( phi[ihit] - intercept - slope * x[ihit] );
02038 chisq += resid[ihit] * resid[ihit]/(sigma*sigma) ;
02039 }
02040
02041 if(chisq<mchisq){
02042 mchisq=chisq;
02043 T0=tmpT0;
02044 }
02045 }
02046 }
02047 }
02048 }
02049 }
02050
02051 if(tmp!=0){
02052 t0_mean=t0_all/tmp;
02053 }
02054 if(m_ntupleflag && m_tuple2) g_t0mean=t0_mean;
02055
02056 t_Est=T0 + tOffset_b;
02057 tEstFlag=2;
02058 }
02059 if(m_ntupleflag && m_tuple2){
02060 g_t0=t0_comp;
02061 g_T0=T0;
02062 }
02063 if(nmatch_barrel==0 && nmatch_end==0 && nmatch_barrel_1==0&&nmatch_barrel_2==0&&m_mdcCalibFunSvc&&m_flag==2){
02064
02065 log << MSG::INFO << " mdc " << endreq;
02066
02067 #define MXWIRE 6860
02068 #define MXTKHIT 6860
02069 #define MXTRK 15
02070
02071
02072 double VELPROP=33.33;
02073
02074
02075 int nhits_ax = 0;
02076 int nhits_ax2 = 0;
02077 int nhits_st = 0;
02078 int nhits_st2 = 0;
02079
02080 double tev_ax[MXTKHIT];
02081 double tev_st[MXTKHIT];
02082 double tevv[MXTKHIT];
02083 double toft;
02084 double prop;
02085 double t0_minus_TDC[MXWIRE];
02086
02087 double T0_cal=-230;
02088 double Mdc_t0Arr[500];
02089
02090 int expmc=1;
02091 double scale=1.;
02092 int expno, runno;
02093 ndriftt=0;
02094
02095
02096
02097
02098 if(ntot > MXTRK) {
02099 ntot=MXTRK;
02100 }
02101
02102 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
02103 if (!newtrkCol || newtrkCol->size()==0) {
02104 log << MSG::INFO<< "Could not find RecMdcTrackCol" << endreq;
02105 return StatusCode::SUCCESS;
02106 }
02107 log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq;
02108
02109 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
02110
02111 for( int tempntrack=0; iter_trk != newtrkCol->end(); iter_trk++,tempntrack++) {
02112 log << MSG::DEBUG << "retrieved MDC track:"
02113 << " Track Id: " << (*iter_trk)->trackId()
02114 << " Dr: " <<(*iter_trk)->helix(0)
02115 << " Phi0: " << (*iter_trk)->helix(1)
02116 << " kappa: " << (*iter_trk)->helix(2)
02117 << " Dz: " << (*iter_trk)->helix(3)
02118 << " Tanl: " << (*iter_trk)->helix(4)
02119 << " Phi terminal: "<< (*iter_trk)->getFiTerm()
02120 << endreq
02121 << "Number of hits: "<< (*iter_trk)->getNhits()
02122 << " Number of stereo hits " << (*iter_trk)->nster()
02123 << endreq;
02124
02125
02126 const HepPoint3D pivot0(0.,0.,0.);
02127 HepVector a(5,0);
02128
02129 a[0] = (*iter_trk)->helix(0);
02130 a[1] = (*iter_trk)->helix(1);
02131 a[2] = (*iter_trk)->helix(2);
02132 a[3] = (*iter_trk)->helix(3);
02133 a[4] = (*iter_trk)->helix(4);
02134
02135
02136 if (abs(a[0])>Estparam.MDC_drCut() || abs(a[3])>Estparam.MDC_dzCut() || abs(a[4])>500.) continue;
02137
02138 double phi0 = a[1];
02139 double kappa = abs(a[2]);
02140 double dirmag = sqrt(1. + a[4]*a[4]);
02141
02142 double mom = abs(dirmag/kappa);
02143 double beta=mom/sqrt(mom*mom+PIMAS2);
02144 if (particleId[tempntrack]== 1) { beta=mom/sqrt(mom*mom+ELMAS2);}
02145 if(particleId[tempntrack]== 5) { beta=mom/sqrt(mom*mom+PROTONMAS2);}
02146
02147
02148 Helix helix0(pivot0,a);
02149 double rho = helix0.radius();
02150 double unit_s = abs(rho * dirmag);
02151
02152 helix0.ignoreErrorMatrix();
02153 HepPoint3D hcen = helix0.center();
02154 double xc= hcen(0);
02155 double yc= hcen(1);
02156
02157 if( xc==0.0 && yc==0.0 ) continue;
02158
02159 double direction =1.;
02160 if(optCosmic!=0) {
02161 double phi = atan2(helix0.momentum(0.).y(),helix0.momentum(0.).x());
02162 if(phi> 0. && phi <= M_PI) direction=-1.;
02163 }
02164
02165 IMdcGeomSvc* mdcGeomSvc;
02166 StatusCode sc = service("MdcGeomSvc", mdcGeomSvc);
02167 double zeast;
02168 double zwest;
02169 double m_vp[43]={0.}, m_zst[43]={0.};
02170 for(int lay=0; lay<43; lay++){
02171 zwest = mdcGeomSvc->Wire(lay, 0)->Forward().z();
02172 zeast = mdcGeomSvc->Wire(lay, 0)->Backward().z();
02173
02174
02175 if(lay < 8) m_vp[lay] = 220.0;
02176 else m_vp[lay] = 240.0;
02177
02178 if( 0 == (lay % 2) ){
02179 m_zst[lay] = zwest;
02180 } else{
02181 m_zst[lay] = zeast;
02182 }
02183 }
02184
02185
02186 log << MSG::DEBUG << "hitList of this track:" << endreq;
02187 HitRefVec gothits = (*iter_trk)->getVecHits();
02188 HitRefVec::iterator it_gothit = gothits.begin();
02189 for( ; it_gothit != gothits.end(); it_gothit++){
02190
02191 log << MSG::DEBUG << "hits Id: "<<(*it_gothit)->getId()
02192 << " hits MDC layerId wireId " <<MdcID::layer((*it_gothit)->getMdcId())
02193 << " "<<MdcID::wire((*it_gothit)->getMdcId())
02194 << endreq
02195 << " hits TDC " <<(*it_gothit)->getTdc()
02196 << endreq;
02197
02198 int layer=MdcID::layer((*it_gothit)->getMdcId());
02199 int wid=MdcID::wire((*it_gothit)->getMdcId());
02200 double tdc=(*it_gothit)->getTdc() ;
02201
02202 double trkchi2=(*iter_trk)->chi2();
02203 if(trkchi2>100)continue;
02204 double hitChi2=(*it_gothit)->getChisqAdd();
02205 HepVector helix_par = (*iter_trk)->helix();
02206 HepSymMatrix helixErr=(*iter_trk)->err();
02207
02208 if((layer>=8&&layer<=19) ||(layer>=36&&layer<=42)){
02209
02211
02212
02213 const MdcGeoWire* GeoRef = mdcGeomSvc->Wire(layer,wid);
02214
02217
02218
02219
02220
02222
02223
02224 if(Estparam.MDC_Inner()==0 && layer <=3) continue;
02225
02226 double xw = GeoRef->Forward().x()/10;
02227 double yw = GeoRef->Forward().y()/10;
02228
02229 HepPoint3D pivot1(xw,yw,0.);
02230 helix0.pivot(pivot1);
02231 double zw=helix0.a()[3];
02232
02233
02234 double dphi = (-xc*(xw-xc)-yc*(yw-yc)) /
02235 sqrt((xc*xc+yc*yc)*((xw-xc)*(xw-xc)+(yw-yc)*(yw-yc)));
02236 dphi = acos(dphi);
02237 double pathtof = abs(unit_s * dphi);
02238 if (kappa!=0) {
02239 toft = pathtof/VLIGHT/beta;
02240 } else {
02241 toft = pathtof/VLIGHT;
02242 }
02243
02244
02247
02248 if (zw >(GeoRef->Backward().z())/10) zw =(GeoRef->Backward().z())/10;
02249 if (zw <(GeoRef->Forward().z())/10) zw =(GeoRef->Forward().z())/10;
02250
02251 double slant = GeoRef ->Slant();
02252 prop = zw / cos(slant) / VELPROP;
02253
02254 double Tw = 0;
02255
02256
02257
02258
02259
02260 double driftt;
02261 double dummy;
02262 int lr=2;
02263
02264 double p[3];
02265 p[0]=helix0.momentum(0.).x();
02266 p[1]=helix0.momentum(0.).y();
02267 double pos[2];
02268 pos[0]=xw; pos[1]=yw;
02269 double alpha;
02271
02272
02273 double dist = 0.;
02274
02275 dist=(m_mdcUtilitySvc->doca(layer, wid, helix_par, helixErr))*10.0;
02276
02277 if(dist<0.) lr=1;
02278 else lr=0;
02279 dist=fabs(dist);
02280 if(dist> 0.4*(mdcGeomSvc->Layer(layer))->PCSiz()) continue;
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298 int idummy;
02299
02300 if(!m_useXT) {driftt = dist/Estparam.vdrift();}
02301 else {
02302 double entrance=(*it_gothit)->getEntra();
02303 driftt = m_mdcCalibFunSvc->distToDriftTime(dist, layer, wid,lr,entrance);
02304 }
02305 if(m_useT0cal)
02306 {
02307 T0_cal=m_mdcCalibFunSvc->getT0(layer, wid)+m_mdcCalibFunSvc->getTimeWalk(layer,tdc);
02308 }
02309
02310 double zprop = fabs(zw - m_zst[layer]);
02311 double tp = zprop / m_vp[layer];
02312
02313 if(driftt>tdc) continue;
02314 double difft=tdc-driftt-toft-tp-T0_cal;
02315 if(ndriftt>=500) break;
02316 if(difft<-10) continue;
02317 Mdc_t0Arr[ndriftt]=difft;
02318
02319 sum_EstimeMdc=sum_EstimeMdc+difft;
02320 ndriftt++;
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330 double tev= -t0_minus_TDC[wid]+ driftt;
02331 if(Estparam.MDC_Tof() !=0) tev += direction*toft;
02332 if(Estparam.MDC_Prop()!=0) tev += prop;
02334
02335
02336
02337 nhits_ax++;
02338 tev_ax[nhits_ax-1]=tev;
02339
02340 if(Estparam.MDC_Debug()!=0) log << MSG::INFO << "*** tev ***" <<tev <<endreq;
02341 double driftt_mea = t0_minus_TDC[wid];
02342
02343 if(abs(driftt - driftt_mea)<75.) {
02344
02345 nhits_ax2++;
02346 if(Estparam.MDC_Debug()!=0) log << MSG::INFO << "*** tev2 ***" <<tev <<endreq;
02347 }
02348 }
02349
02350
02351 else if(((layer>=4&&layer<=7)||(layer>=20&&layer<=35))&&m_useSw){
02352
02353 IMdcGeomSvc* mdcGeomSvc;
02354 StatusCode sc = service("MdcGeomSvc", mdcGeomSvc);
02355 const MdcGeoWire* GeoRef = mdcGeomSvc->Wire(layer,wid);
02356
02357
02359
02360 double bx= GeoRef->Backward().x()/10;
02361 double by= GeoRef->Backward().y()/10;
02362 double bz= GeoRef->Backward().z()/10;
02363 double fx= GeoRef->Forward().x()/10;
02364 double fy= GeoRef->Forward().y()/10;
02365 double fz= GeoRef->Forward().z()/10;
02366
02369 HepPoint3D fwd(fx,fy,fz);
02370 HepPoint3D bck(bx,by,bz);
02371
02372 Hep3Vector wire = (CLHEP::Hep3Vector)bck - (CLHEP::Hep3Vector)fwd;
02373 HepPoint3D try1=(fwd + bck) * .5;
02374 helix0.pivot(try1);
02375 HepPoint3D try2 = (helix0.x(0).z() - bck.z())/ wire.z() * wire + bck;
02376 helix0.pivot(try2);
02377 HepPoint3D try3 = (helix0.x(0).z() - bck.z())/ wire.z() * wire + bck;
02378 helix0.pivot(try3);
02379
02380 double xh = helix0.x(0.).x();
02381 double yh = helix0.x(0.).y();
02382 double z = helix0.x(0.).z();
02383
02384
02385 double dphi = (-xc*(xh-xc)-yc*(yh-yc)) /
02386 sqrt((xc*xc+yc*yc)*((xh-xc)*(xh-xc)+(yh-yc)*(yh-yc)));
02387 dphi = acos(dphi);
02388 double pathtof = abs(unit_s * dphi);
02389 if (kappa!=0) {
02390 toft = pathtof/VLIGHT/beta;
02391 } else {
02392 toft = pathtof/VLIGHT;
02393 }
02394
02395
02396
02397 if (z != 9999.) {
02398
02399 if(z < fz ) z = fz;
02400
02401 if(z > bz ) z = bz;
02402 double slant = GeoRef->Slant();
02403 prop = z / cos(slant) / VELPROP;
02404 } else {
02405 prop = 0.;
02406 }
02407
02408
02409 double Tw = 0;
02410
02411
02412
02413
02414
02415 double driftt=0;
02416 double dummy;
02417
02418 double xw = fx + (bx-fx)/(bz-fz)*(z-fz);
02419 double yw = fy + (by-fy)/(bz-fz)*(z-fz);
02420
02421 HepPoint3D pivot1(xw,yw,z);
02422 helix0.pivot(pivot1);
02423
02424 double zw=helix0.a()[3];
02425
02426 int lr=2;
02427
02428 double p[3];
02429 p[0]=helix0.momentum(0.).x();
02430 p[1]=helix0.momentum(0.).y();
02431 double pos[2];
02432 pos[0]=xw; pos[1]=yw;
02433 double alpha;
02435
02436
02437 double dist=0.;
02438
02439 dist=(m_mdcUtilitySvc->doca(layer, wid, helix_par, helixErr))*10.0;
02440
02441 if(dist<0.) lr=1;
02442 else lr=0;
02443 dist=fabs(dist);
02444 if(dist> 0.4*(mdcGeomSvc->Layer(layer))->PCSiz()) continue;
02445
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462 if(!m_useXT){driftt = dist/(Estparam.vdrift());}
02463 else {
02464 double entrance=(*it_gothit)->getEntra();
02465 driftt = m_mdcCalibFunSvc->distToDriftTime(dist, layer, wid,lr,entrance);
02466 }
02467 if(m_useT0cal)
02468 {
02469 T0_cal=m_mdcCalibFunSvc->getT0(layer, wid)+m_mdcCalibFunSvc->getTimeWalk(layer,tdc);
02470 }
02471
02472 double zprop = fabs(zw - m_zst[layer]);
02473 double tp = zprop / m_vp[layer];
02474
02475 if(driftt>tdc) continue;
02476 double difft=tdc-driftt-toft-tp-T0_cal;
02477 if(difft<-10) continue;
02478 if(ndriftt>=500) break;
02479 Mdc_t0Arr[ndriftt]=difft;
02480
02481
02482
02483 sum_EstimeMdc=sum_EstimeMdc+difft;
02484 ndriftt++;
02485
02486
02487
02488
02489 double tev= -t0_minus_TDC[wid]+ driftt;
02490 if(Estparam.MDC_Tof() !=0) tev += direction*toft;
02491 if(Estparam.MDC_Prop()!=0) tev += prop;
02493
02494
02495
02496
02497
02498 nhits_st++;
02499 tev_st[nhits_st-1]= tev;
02500
02501 if(Estparam.MDC_Debug()!=0) log << MSG::INFO << "*** tev_st ***" <<tev <<endreq;
02502 double driftt_mea = t0_minus_TDC[wid];
02503
02504 if(abs(driftt - driftt_mea) <75.) {
02505
02506 nhits_st2++;
02507 if(Estparam.MDC_Debug()!=0) log << MSG::INFO << "*** tev_st2 ***" <<tev <<endreq;
02508 }
02509 }
02510
02511 }
02512 nmatch_mdc++;
02513 }
02514
02515
02516 if(m_ntupleflag && m_tuple2) g_nmatchmdc = nmatch_mdc;
02517 if(ndriftt!=0){
02518 if(m_mdcopt) {
02519 sum_EstimeMdc=Opt_new(Mdc_t0Arr,ndriftt,400.0);
02520 }
02521 else { sum_EstimeMdc= sum_EstimeMdc/ndriftt;}
02522 if(m_ntupleflag && m_tuple2) g_EstimeMdc=sum_EstimeMdc;
02523 t_Est=sum_EstimeMdc + tOffset_b;
02524 if(t_Est<0) t_Est=0;
02525 if(optCosmic==0){
02526 tEstFlag=102;
02527 nbunch=((int)(t_Est-offset))/bunchtime;
02528
02529 if((t_Est-offset-nbunch*bunchtime)>(bunchtime/2)) nbunch=nbunch+1;
02530 t_Est=nbunch*bunchtime+offset + tOffset_b;
02531
02532
02533
02534
02535
02536
02537
02538
02539 }
02540 if(optCosmic){
02541 t_Est=sum_EstimeMdc;
02542 tEstFlag=202;
02543 }
02544 }
02545 if(m_ntupleflag && m_tuple2) g_ndriftt=ndriftt;
02546 }
02547
02548 if(t_Est!=-999){
02549
02550 if((!m_beforrec) && (Testime_fzisan != t_Est) ){
02551 if(tEstFlag==211) tEstFlag=213;
02552 if(tEstFlag==212) tEstFlag=216;
02553 if(tEstFlag==111) tEstFlag=113;
02554 if(tEstFlag==112) tEstFlag=116;
02555 }
02556
02557 if( optCosmic ){
02558 StatusCode scStoreTds = storeTDS(t_Est,tEstFlag,t_quality);
02559 if (scStoreTds!=StatusCode::SUCCESS){ return scStoreTds; }
02560 }else if(!optCosmic){
02561 StatusCode scStoreTds = storeTDS(t_Est,tEstFlag,t_quality);
02562 if (scStoreTds!=StatusCode::SUCCESS){ return scStoreTds; }
02563 }
02564 }else{
02565
02566 if(m_beforrec){
02567
02568
02569 double segTest = Testime_fzisan + tOffset_b;
02570 int segFlag = TestimeFlag_fzisan;
02571 double segQuality = TestimeQuality_fzisan;
02572 StatusCode scStoreTds = storeTDS(segTest,segFlag,segQuality);
02573 if (scStoreTds!=StatusCode::SUCCESS){ return scStoreTds; }
02574 }
02575 }
02576
02577
02578
02579 SmartDataPtr<RecEsTimeCol> arecestimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
02580 if (!arecestimeCol) {
02581 if(m_debug==4) log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq;
02582 return( StatusCode::SUCCESS);
02583 }
02584 RecEsTimeCol::iterator iter_evt = arecestimeCol->begin();
02585 for(; iter_evt!=arecestimeCol->end(); iter_evt++){
02586 log << MSG::INFO
02587 << "Test = "<<(*iter_evt)->getTest()
02588 << ", Status = "<<(*iter_evt)->getStat()
02589 <<endreq;
02590 if(m_ntupleflag && m_tuple2){
02591 g_Testime=(*iter_evt)->getTest();
02592 }
02593
02594 }
02595
02596 if(m_ntupleflag){
02597 if(m_tuple2){
02598 for(g_ntofdown=0;g_ntofdown<ntofdown;g_ntofdown++){ g_meantevdown[g_ntofdown]=meantevdown[g_ntofdown];}
02599 for(g_ntofup=0;g_ntofup<ntofup;g_ntofup++){ g_meantevup[g_ntofup]=meantevup[g_ntofup];}
02600 g_nmatch_tot=nmatch;
02601 m_estFlag=tEstFlag;
02602 m_estTime=t_Est;
02603 StatusCode status = m_tuple2->write();
02604 if (!status.isSuccess()) {
02605 log << MSG::ERROR << "Can't fill ntuple!" << endreq;
02606 }
02607 }
02608 if(m_tuple9){
02609 for(g_nmatch=0;g_nmatch<nmatch;g_nmatch++)
02610 {
02611 g_meantev[g_nmatch]=meantev[g_nmatch];
02612 }
02613 StatusCode status = m_tuple9->write();
02614 if (!status.isSuccess()) {
02615 log << MSG::ERROR << "Can't fill ntuple!" << endreq;
02616 }
02617 }
02618 }
02619 return StatusCode::SUCCESS;
02620 }
