Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MdcDedxRecon Class Reference

#include <MdcDedxRecon.h>

List of all members.

Public Member Functions

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

Public Attributes

int calib_flag
int cosmic
float cut_sigma
int dbgFlag
int doNewGlobal
int for_calib
int landau
int m_alg
float m_bitrunc1
float m_chidedx1
float m_costheta1
float m_expect1
int m_nhit1
int m_parttype1
float m_pch1
float m_phtm1
float m_probpid1 [5]
float m_sigma1
float m_sintheta1
int m_trkalg
vector< int > m_trkalgs
vector< int > m_trkalgs
float mincut_dedx
int minhit
int ntpFlag
int ParticleFlag
int recalg
int stsw
int timing_check
float truncate
int x_rej

Private Attributes

int begin_run_stat
IDedxCurSvc * dedxcursvc
IDedxCurSvc * dedxcursvc
MdcDedxCorrection * ex_calib
MdcDedxCorrection * ex_calib
std::vector< MdcDedxTrk > ex_trks
std::vector< MdcDedxTrk > ex_trks
IDedxCorrecSvc * exsvc
IDedxCorrecSvc * exsvc
IMdcGeomSvc * geosvc
IMdcGeomSvc * geosvc
IHistogram1D * his_aver
IHistogram1D * his_aver
IHistogram1D * his_ph
IHistogram1D * his_ph
NTuple::Item< float > m_bitrunc
NTuple::Item< float > m_bitrunc
NTuple::Item< float > m_charge
NTuple::Item< float > m_charge
NTuple::Item< double > m_charge1
NTuple::Item< double > m_charge1
NTuple::Item< float > m_chidedx
NTuple::Item< float > m_chidedx
NTuple::Item< float > m_chidedxe
NTuple::Item< float > m_chidedxe
NTuple::Item< float > m_chidedxk
NTuple::Item< float > m_chidedxk
NTuple::Item< float > m_chidedxmu
NTuple::Item< float > m_chidedxmu
NTuple::Item< float > m_chidedxp
NTuple::Item< float > m_chidedxp
NTuple::Item< float > m_chidedxpi
NTuple::Item< float > m_chidedxpi
NTuple::Item< float > m_costhe
NTuple::Item< float > m_costhe
NTuple::Item< float > m_costheta
NTuple::Item< float > m_costheta
NTuple::Item< float > m_dEdx_meas
NTuple::Item< float > m_dEdx_meas
NTuple::Item< float > m_doca
NTuple::Item< float > m_doca
NTuple::Item< float > m_driftdist
NTuple::Item< float > m_driftdist
NTuple::Item< float > m_driftT
NTuple::Item< float > m_driftT
NTuple::Item< float > m_evtNO
NTuple::Item< float > m_evtNO
NTuple::Item< float > m_expect
NTuple::Item< float > m_expect
NTuple::Array< float > m_expectid
NTuple::Array< float > m_expectid
NTuple::Item< double > m_exraw
NTuple::Item< double > m_exraw
NTuple::Item< float > m_layer
NTuple::Item< float > m_layer
NTuple::Item< float > m_localwid
NTuple::Item< float > m_localwid
NTuple::Item< float > m_ndedxhit
NTuple::Item< float > m_ndedxhit
unsigned m_nevt
NTuple::Item< float > m_nhit
NTuple::Item< float > m_nhit
NTuple::Item< float > m_nhit_hit
NTuple::Item< float > m_nhit_hit
NTuple::Tuple * m_nt1
NTuple::Tuple * m_nt1
NTuple::Tuple * m_nt2
NTuple::Tuple * m_nt2
NTuple::Item< float > m_parttype
NTuple::Item< float > m_parttype
NTuple::Item< float > m_pathL
NTuple::Item< float > m_pathL
NTuple::Item< float > m_pch
NTuple::Item< float > m_pch
NTuple::Item< double > m_phraw
NTuple::Item< double > m_phraw
NTuple::Item< float > m_phtm
NTuple::Item< float > m_phtm
NTuple::Array< float > m_probpid
NTuple::Array< float > m_probpid
NTuple::Item< float > m_ptrk1
NTuple::Item< float > m_ptrk1
float m_rec_sw
NTuple::Item< float > m_runNO
NTuple::Item< float > m_runNO
NTuple::Item< double > m_runNO1
NTuple::Item< double > m_runNO1
NTuple::Item< float > m_sigma
NTuple::Item< float > m_sigma
NTuple::Array< float > m_sigmaid
NTuple::Array< float > m_sigmaid
NTuple::Item< float > m_sinenta
NTuple::Item< float > m_sinenta
NTuple::Item< float > m_sintheta
NTuple::Item< float > m_sintheta
NTuple::Item< float > m_t0
NTuple::Item< float > m_t0
NTuple::Item< float > m_t01
NTuple::Item< float > m_t01
NTuple::Item< float > m_tdc_raw
NTuple::Item< float > m_tdc_raw
NTuple::Item< float > m_trackId
NTuple::Item< float > m_trackId
NTuple::Item< float > m_trackId1
NTuple::Item< float > m_trackId1
NTuple::Item< float > m_wire
NTuple::Item< float > m_wire

Constructor & Destructor Documentation

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

>0 : debug on ; 0 : off

truncation rate ( lower side )

mininum number of hits resuired

truncation rate

temporily for q-r correction

xhit rejection 1 : on; 0 : off

0 : gauss distribution

for calibration flag

all calib on

------------------------------------------------------------- 

00063                                                                           :  Algorithm(name, pSvcLocator)
00064 { 
00065    
00066    ex_calib=0;  // MdcDedxCorrection  
00067    dbgFlag = 2;  
00068    mincut_dedx = 0.0; 
00069    minhit = 1;     
00070    truncate = 0.70; 
00071    stsw = 0;  
00072    x_rej = 1;  
00073    landau = 1;  
00074 
00075    calib_flag = 0x7F; 
00076    m_alg = 1;
00077    recalg = 0;   
00078    doNewGlobal = 0;
00079    ntpFlag = 1;
00081    m_rec_sw = 1.8;
00082    m_nevt = 0;
00083    m_nt1 = 0;
00084    m_nt2 = 0;
00085    ParticleFlag = 0;
00086    for_calib = 0;
00087    // Declare the properties  
00088    declareProperty("RecAlg",recalg);  
00089    declareProperty("dEdxMethod",m_alg); 
00090    declareProperty("TruncRate",truncate); 
00091    declareProperty("CalibFlag",calib_flag); 
00092    declareProperty("DebugFlag",dbgFlag); 
00093    declareProperty("NTupleFlag",ntpFlag);
00094    declareProperty("ParticleType",ParticleFlag);
00095    declareProperty("ForCalib",for_calib); 
00096 }

MdcDedxRecon::~MdcDedxRecon (   )  [inline]
  

00031 {};

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

MdcDedxRecon::~MdcDedxRecon (   )  [inline]
  

00031 {};

Member Function Documentation

StatusCode MdcDedxRecon::beginRun (   ) 
 

StatusCode MdcDedxRecon::beginRun (   ) 
  

00260                                  {
00261    MsgStream log(msgSvc(), name());
00262    log << MSG::DEBUG <<"MdcDedxRecon begin_run!!!"<< endreq;
00263    //   begin_run_stat = ex_calib->read_data();
00264    if( dbgFlag > 0 ){
00265            log << MSG::DEBUG <<"begin_run : RunNo="<< endreq;
00266    }
00267   
00268    StatusCode gesc = service("MdcGeomSvc", geosvc);
00269    if (gesc ==  StatusCode::SUCCESS) {
00270       log << MSG::INFO <<"MdcGeomSvc is running"<<endl;
00271    } else {
00272       log << MSG::ERROR <<"MdcGeomSvc is failed"<<endl;
00273       return StatusCode::SUCCESS;
00274    }  
00275 
00276    return StatusCode::SUCCESS;
00277 }

double MdcDedxRecon::cal_dedx (  int  ,
int  ,
float 
) 
 

double MdcDedxRecon::cal_dedx (  int  ,
int  ,
float 
) 
 

void MdcDedxRecon::clearTables (   ) 
 

void MdcDedxRecon::clearTables (   ) 
  

01100                                 {
01101 }

StatusCode MdcDedxRecon::execute (   ) 
 

StatusCode MdcDedxRecon::execute (   ) 
 

remove last event information.

dbg

cal expect 

00292                                  {  
00293 
00294   MsgStream log(msgSvc(), name());
00295   log << MSG::INFO << "in execute()" << endreq;
00296  
00297   // for 652a psip data
00298 /* double DedxCurve_Parameter[5]; 
00299    double DedxSigma_Parameter[14];
00300    for(int i=0;i<5;i++){
00301       DedxCurve_Parameter[i] = dedxcursvc->getCurve(i);  
00302       //std::cout<<"vtxsvc->getCurve(i)"<<dedxcursvc->getCurve(i)<<std::endl;
00303   }
00304   
00305   for(int i=0;i<14;i++){
00306       DedxSigma_Parameter[i] = dedxcursvc->getSigma(i);   
00307       //std::cout<<"vtxsvc->getSigma(i)"<<dedxcursvc->getSigma(i)<<std::endl;
00308       } */
00309  
00310 
00311   vector<double> Curve_Para, Sigma_Para;
00312   //flag definition: vFlag[0] is dedx curve version  
00313   //                 vFlag[1] is dedx sigma version  2:psip 3:jpsi
00314   //                 vFlag[2] is MC ,when = 1
00315   int vFlag[3];
00316   
00317   for(int i=0; i< dedxcursvc->getCurveSize(); i++){
00318     if(i==0) 
00319       vFlag[0] = (int) dedxcursvc->getCurve(i);     // get the dedx curve version
00320     else 
00321       Curve_Para.push_back( dedxcursvc->getCurve(i) ); 
00322        // get the dedx curve parameters
00323   }
00324   for(int i=0; i< dedxcursvc->getSigmaSize(); i++){
00325     if (i==0) 
00326       vFlag[1] = (int) dedxcursvc->getSigma(i);   // get the sigma version 2: psip 3:jpsi
00327     else 
00328       Sigma_Para.push_back( dedxcursvc->getSigma(i) );
00329       //get the dedx sigma parameters
00330   }
00331   
00332  
00333   //check if size of parameters is right
00334   if(vFlag[0] ==1 && Curve_Para.size() != 5){             //version 1: 5 parameters 652a psip data 
00335     cout<<" size of dedx curve paramters for version 1 is not right!"<<endl;
00336   }
00337   if(vFlag[0] ==2 && Curve_Para.size() != 16){            // version 2: 16 parameters  from 652a parameters add to 16
00338     cout<<" size of dedx curve paramters for version 2 is not right!"<<endl;
00339   }
00340 
00341 
00342   if(vFlag[1] ==1 && Sigma_Para.size() != 14){           //vesion 1: 14 parameters 652a psip data (old way)
00343     cout<<" size of dedx sigma paramters for version 1 is not right!"<<endl;
00344   }
00345   if(vFlag[1] ==2 && Sigma_Para.size() != 21){          //version 2: include t0 correction ( for psip09 data) 
00346     cout<<" size of dedx sigma paramters for version 2 is not right!"<<endl;
00347   }
00348   if(vFlag[1] ==3 && Sigma_Para.size() != 18){          //version 3: no t0 correction (for jpsi09 data and beyond)
00349     cout<<" size of dedx sigma paramters for version 3 is not right!"<<endl;
00350   }
00351   if(vFlag[1] ==4 && Sigma_Para.size() != 19){          //version 4: data with mdc track defaulted when kal track not ok , (no t0 correction) 
00352     cout<<" size of dedx sigma paramters for version 4 is not right!"<<endl;
00353   }
00354   if(vFlag[1] ==5 && Sigma_Para.size() != 22){          //version 5: data with mdc track defaulted when kal track not ok , (with t0 correction) 
00355     cout<<" size of dedx sigma paramters for version 5 is not right!"<<endl;
00356   }
00357     
00358   
00359   
00360   Identifier mdcid;
00361   log << MSG::DEBUG <<"MdcDedxRecon event start!!!"<< endreq;
00362   if( begin_run_stat == -999 ) {
00363     log << MSG::ERROR <<" begin_run() did not succeed "<< endreq;
00364      return StatusCode::SUCCESS;
00365   }
00366   //---------- register RecMdcDedxCol and RecMdcDedxHitCol to tds--------- 
00367   //
00368   DataObject *aReconEvent;
00369   eventSvc()->findObject("/Event/Recon",aReconEvent);
00370   if(aReconEvent==NULL) {
00371      aReconEvent = new ReconEvent();
00372      StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00373      if(sc!=StatusCode::SUCCESS) {
00374              log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00375              return( StatusCode::FAILURE);
00376      } 
00377   }
00378 
00379   IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc()); 
00380   DataObject *aDedxcol;
00381   eventSvc()->findObject("/Event/Recon/RecMdcDedxCol",aDedxcol);
00382   if(aDedxcol != NULL) {
00383      //IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00384      dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxCol");
00385      eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxCol");
00386   }
00387   //MdcDedxCol* dedxlist= new MdcDedxCol;
00388   RecMdcDedxCol* dedxList = new RecMdcDedxCol;
00389   StatusCode dedxsc;
00390   dedxsc =  eventSvc()->registerObject(EventModel::Recon::RecMdcDedxCol, dedxList);
00391   if(!dedxsc.isSuccess()) {
00392      log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
00393      return ( StatusCode::FAILURE);
00394   }
00395 
00396   DataObject *aDedxhitcol;
00397   eventSvc()->findObject("/Event/Recon/RecMdcDedxHitCol",aDedxhitcol);
00398   if(aDedxhitcol != NULL) {
00399      //IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00400      dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxHitCol");
00401      eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxHitCol");
00402   }
00403 
00404   // RecMdcDedxHitCol* dedxhitlist= new RecMdcDedxHitCol;
00405   RecMdcDedxHitCol* dedxhitList = new RecMdcDedxHitCol;
00406   StatusCode dedxhitsc;
00407   dedxhitsc =  eventSvc()->registerObject(EventModel::Recon::RecMdcDedxHitCol, dedxhitList);
00408   if(!dedxhitsc.isSuccess()) {
00409      log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
00410      return ( StatusCode::FAILURE);
00411   }
00412 
00413 
00414   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00415   if (!eventHeader) {
00416     log << MSG::INFO << "Could not find Event Header" << endreq;
00417     return( StatusCode::FAILURE);
00418   }
00419   int event = eventHeader->eventNumber();
00420   int runNO = eventHeader->runNumber();
00421   //if(ntpFlag>0) cout<<"event  ----- : "<<event<<endl;
00422   RunNO2= runNO;
00423   if(RunNO1==0) RunNO1=runNO;
00424   if(RunNO2 != RunNO1){
00425     cout<<"RunNO2 = "<<RunNO2 <<"      RunNO1 = "<<RunNO1<<endl;
00426   }
00427   RunNO1 = runNO;
00428   //if(ntpFlag>0) cout<<"event 1 ----- = "<<event<<"        runNO 1 = "<<runNO <<endl;
00429   int runFlag; //data type flag: 0:MC data,   1: 2150,      2:2200
00430   if(runNO<MdcDedxParam::RUN0) runFlag =0;
00431   if(runNO>=MdcDedxParam::RUN1 && runNO<MdcDedxParam::RUN2) runFlag =1;
00432   if(runNO>=MdcDedxParam::RUN2 && runNO<MdcDedxParam::RUN3) runFlag =2;
00433   if(runNO>=MdcDedxParam::RUN3 && runNO<MdcDedxParam::RUN4) runFlag =3;
00434   if(runNO>=MdcDedxParam::RUN4) runFlag =4;
00435   
00436   
00437   if(runNO < 0)     // identify the runno , and for mc data version flag[2] set to 1 !
00438     vFlag[2]=1;
00439   else
00440     vFlag[2]=0;
00441 
00442   
00443   double tes;
00444   int esTimeflag = -1;
00445   SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00446   if( ! estimeCol) {
00447           log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
00448   }else{
00449           RecEsTimeCol::iterator iter_evt = estimeCol->begin();
00450           for(; iter_evt!=estimeCol->end(); iter_evt++){
00451                   tes =  (*iter_evt)->getTest();
00452                   esTimeflag = (*iter_evt)->getStat();
00453           }
00454   }
00455  //cout<<"estime : "<<tes<<endl;
00456 
00457   if(for_calib == 1){
00458      if(runFlag ==2) {if( tes>1000 ) return StatusCode::SUCCESS;;}
00459      else if(runFlag ==3 ){if (tes>700 ) return StatusCode::SUCCESS;}
00460      else if(runFlag ==4 ){if (tes>1400 ) return StatusCode::SUCCESS;}
00461      //if(runNO>0 && runNO < 6202) {if( tes>1000 ) return StatusCode::SUCCESS;}
00462      //else if( tes>1400 ) return StatusCode::SUCCESS;
00463   } 
00464     
00465   if(for_calib == 1 && runFlag ==3 && tes>700) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"      tes : "<<tes<<endl;
00466   if(for_calib == 1 && runFlag ==4 && tes>1400) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"      tes : "<<tes<<endl; 
00467 
00468 
00469   log << MSG::INFO << "now begin the event  " << event << endreq;
00470   //if(ntpFlag>0) cout<<"run number is : "<<runNO<<"       event NO= "<<event<<endl;
00471  
00472   if(dbgFlag > 0) {
00473     ++m_nevt;
00474     log << MSG::DEBUG <<" now calling MdcDedxRecon start : #event "<<m_nevt
00475         << endreq;
00476   }
00477   
00478   int ntrk;
00479   ex_trks.clear();  // clear the vector of MdcDedxTrk,when read a new event  
00480   m_trkalgs.clear();
00481   if( !doNewGlobal ) {
00482     //if(ntpFlag>0) cout<<"recalg = "<<recalg<<endl; 
00483     
00484     if(recalg ==2 ){
00485       //cout<<"using kal algorithm by default, switch to seed track if no kaltrack hits."<<endl;
00486       if(for_calib == 1){
00487         cout<<"warning! calibration should set recalg to 1.0"<<endl;
00488         return( StatusCode::SUCCESS);
00489       }
00490       //retrieve  MdcKalTrackCol from TDS
00491       SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00492       if (!newtrkCol) {
00493         //cout<<"Could not find Kal trk -----------------------"<<endl;
00494         log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
00495         if(ntpFlag>0) return( StatusCode::SUCCESS);
00496       }
00497       if(ntpFlag>0) eventNO++;
00498       ntrk = newtrkCol->size();
00499       log << MSG::DEBUG << "trkCol size:  " << newtrkCol->size() <<endreq;
00500       
00501       
00502       RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
00503       for( ; trk != newtrkCol->end(); trk++) {
00504         
00505         HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
00506         
00507         if(gothits.size()==0){
00508           m_trkalg=0;
00509           int trkid=(*trk)->trackId();
00510           switchtomdctrack(trkid, mdcid,runNO, runFlag, log);
00511         }
00512         else{
00513           m_trkalg=1;
00514           if(gothits.size()<2) continue;
00515           kaltrackrec(trk, mdcid, runNO, runFlag, log );
00516         }
00517       }//end of track loop
00518       
00519     }//end of recalg==2
00520     
00521     else if(recalg ==0 ){
00522       m_trkalg=0;
00523       //retrieve  MdcTrackCol from TDS     
00524       //cout<<"using rec algorithm ................"<<endl;
00525        SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00526        if (!newtrkCol) {
00527           log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
00528           return( StatusCode::SUCCESS);
00529        }
00530        ntrk = newtrkCol->size();    // track numbers of one event
00531        if(ntpFlag>0) eventNO ++;
00532        double poca_x, poca_y, poca_z;
00533        log << MSG::DEBUG << "trkCol size:  " << newtrkCol->size() <<endreq;
00534        //if(ntpFlag>0) cout<< "trkCol size:  " << newtrkCol->size() <<endl;
00535        if(for_calib == 1){if(newtrkCol->size() !=2 ) return( StatusCode::SUCCESS);}  // calib use bhabha event?
00536        RecMdcTrackCol::iterator trk = newtrkCol->begin();              // circulate tracks 
00537        for( ; trk != newtrkCol->end(); trk++) {
00538           if(ntpFlag>0) trackNO1++;
00539           RecMdcDedxHit* dedxhit = new RecMdcDedxHit;
00540           MdcDedxTrk trk_ex( **trk); 
00541           HepVector a = (*trk)->helix();
00542           poca_x = (*trk)->x();
00543           poca_y = (*trk)->y();
00544           poca_z = (*trk)->z();
00545           
00546           HepSymMatrix tkErrM = (*trk)->err();
00547           double m_d0E   = tkErrM[0][0];
00548           double m_phi0E = tkErrM[1][1];
00549           double m_cpaE  = tkErrM[2][2];
00550           double m_z0E   = tkErrM[3][3]; 
00551 
00552           HepPoint3D IP(0,0,0);
00553           Dedx_Helix exhel(IP, a);
00554           for(int i=0;i<5;i++) {
00555                log << MSG::DEBUG <<" helix["<<i<<"] = "<<(*trk)->helix()[i] << endreq;
00556           }
00557 
00558           double theta_temp = M_PI_2-atan(a[4]);
00559           double costhe_temp = cos(M_PI_2-atan(a[4]));
00560           if(for_calib == 1 && abs(costhe_temp)>0.9) continue;
00561  
00562           double m_Pt = 1.0/fabs( a(3) );
00563           double m_P = m_Pt*sqrt(1 + a(5)*a(5));
00564           if(for_calib == 1 && runFlag ==3 &&(m_P>1.95||m_P<1.8)) continue;
00565           if(for_calib == 1 && runFlag ==4 && runNO < MdcDedxParam::RUN5 &&(m_P>1.72||m_P<1.45)) continue;
00566           // for new data psipp data we need to modify the momentum scope
00567           if(for_calib == 1 && runFlag ==4 && runNO>= MdcDedxParam::RUN5 &&(m_P>2.0||m_P<1.7)) continue;
00568 
00569           if(for_calib ==1&&runFlag ==3&&m_P>1.95&&m_P<1.8) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"       m_P : "<<m_P<<endl;
00570           if(for_calib ==1&&runFlag ==4 &&  runNO < MdcDedxParam::RUN5&&m_P>1.72&&m_P<1.45) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"      m_P : "<<m_P<<endl;
00571           if(for_calib ==1&&runFlag ==4 &&  runNO>= MdcDedxParam::RUN5&&m_P>2.0&&m_P<1.7) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"        m_P : "<<m_P<<endl;
00572 
00573           //int charge = ( ((*trk)->helix()[2]) >= 0 )? 1 : -1; 
00574           //int charge = (*trk)->charge();
00575           //cout<<"charge = " <<charge<<endl;     
00576 
00577           log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
00578           log << MSG::DEBUG <<" %%%%%initial id = "<<(*trk)->trackId() <<endreq;
00579           // 2005/5 make phlist here... 
00580           vector<double> phlist;
00581           int layid, localwid,lr, wid, w0id;
00582           double theta, costhe, costhe_pol, sinenta, eangle, dd, dd_pol, ph, adc_raw, ddr, ddl, tdc_raw, driftT;
00583           double zhit, zhit_pol,  m_driftd, sintheta, rphi_path, driftd; 
00584           double pathlength,sindip;
00585           HitRefVec gothits= (*trk)->getVecHits();
00586           log << MSG::DEBUG << "hits size  = " <<gothits.size()<<endreq;
00587           //if(ntpFlag>0) cout<<"hits size  = " <<gothits.size()<<endl;
00588           if(for_calib == 1){if(gothits.size()<20) return StatusCode::SUCCESS;}
00589           if(gothits.size()<2) continue;
00590           //cout<<"nhits = "<<gothits.size()<<endl;
00591           HitRefVec::iterator it_gothit = gothits.begin();
00592           for( ; it_gothit != gothits.end(); it_gothit++) {
00593              mdcid = (*it_gothit)->getMdcId();
00594              layid =  MdcID::layer(mdcid);
00595              localwid = MdcID::wire(mdcid);
00596              w0id = geosvc->Layer(layid)->Wirst();
00597              wid = w0id + localwid;        
00598              adc_raw = (*it_gothit)->getAdc();
00599              //adc_raw = (*it_gothit)->getAdc()*1000000;
00600              tdc_raw = (*it_gothit)->getTdc();
00601              driftT = (*it_gothit)->getDriftT();
00602              //cout<<"adc = "<<adc_raw<<endl;
00603              log << MSG::INFO << "hit layer id  " << layid 
00604                      << "   wire id  " << localwid
00605                      << "   adc_raw  " << adc_raw <<endreq; 
00606              zhit = (*it_gothit)->getZhit();
00607              zhit_pol = fabs(zhit); 
00608              lr = (*it_gothit)->getFlagLR();
00609              if(lr == 2) cout<<"lr = "<<lr<<endl;
00610              if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft();
00611              if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight();
00612              driftd = abs(m_driftd);
00613              dd = (*it_gothit)->getDoca();
00614              if(lr == 0 || lr == 2 ) dd = -abs(dd);
00615              if(lr == 1 ) dd = abs(dd);              
00616              sinenta = (*it_gothit)->getEntra();
00617              eangle = sinenta/M_PI;
00618              theta = M_PI_2-atan(a[4]);
00619              costhe = cos(M_PI_2-atan(a[4]));
00620              sintheta = sin(M_PI_2-atan(a[4]));
00621              //ph = exsvc->StandardCorrec( runFlag, ntpFlag, runNO, exhel, mdcid, adc_raw, dd, eangle, zhit, costhe);
00622              ph = exsvc->StandardHitCorrec(0, runFlag,ntpFlag,runNO,exhel,mdcid,adc_raw,dd,eangle,zhit,costhe);
00623              //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,  tes ); 
00624              //ph = adc_raw;
00625              pathlength = exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit_pol); 
00626              rphi_path = pathlength * sintheta;
00627              //if(ntpFlag>0)
00628              //cout<<event<<setw(10)<<layid<<setw(8)<<localwid<<setw(15)<<driftd<<setw(15)<<rphi_path<<setw(15)<<pathlength<<setw(15)<<ph<<endl; 
00629 
00630              if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
00631              if(runNO<0){  
00632                 if(layid < 4) continue;
00633                 else if(layid>=4&&layid <8){
00634                   if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || dd>Iner_DriftDistCut) continue;
00635                 }
00636                 else if(layid >= 8){
00637                   if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00638                 }
00639              }
00640              else if(runNO>=0){
00641                 if(layid <4) continue;
00642                 else if (layid>=4 && layid<8){
00643                   if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || dd>Iner_DriftDistCut) continue;
00644                 }
00645                 else if(layid >= 8){
00646                   if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00647                 }
00648              }
00649              dedxhit->setMdcId(mdcid);
00650              dedxhit->setFlagLR(lr);
00651              dedxhit->setTrkId(trk_ex.trk_id());
00652              dedxhit->setPathLength(pathlength);
00653              if (ph>0.0) {
00654                 if ( ntpFlag ==2 ) {
00655                    m_charge1 = (*trk)->charge();
00656                    m_t01 = tes;
00657                    //m_driftT = driftT;
00658                    //m_tdc_raw = tdc_raw;
00659                    m_phraw = adc_raw;
00660                    m_exraw = ph;     
00661                    //m_localwid = localwid; 
00662                    //m_wire = wid;
00663                    //m_runNO1 = runNO;
00664                    m_doca = dd;
00665                    m_driftdist = m_driftd;
00666                    m_sinenta = eangle;
00667                    m_costhe = costhe; 
00668                    //m_pathL = pathlength;
00669                    //m_layer = layid;
00670                    m_ptrk1 = m_P;
00671                    StatusCode status2 = m_nt2->write();
00672                    if ( status2.isFailure() )  {
00673                         log << MSG::ERROR << "Cannot fill Ntuple2" << endreq;
00674                    }
00675                 }
00676                 phlist.push_back(ph);
00677              }
00678           }
00679           trk_ex.set_phlist( phlist ); 
00680           ex_trks.push_back(trk_ex );
00681           m_trkalgs.push_back(m_trkalg);
00682           dedxhitList->push_back( dedxhit );
00683           phlist.clear();
00684           if(ntpFlag>0) trackNO3++;
00685           log << MSG::DEBUG << "size in processing:  " << ex_trks.size() <<endreq;
00686        }
00687     }
00688 
00689 
00690     //------------------------ Use information of MDC KAL track rec --------------------------//       
00691     else if(recalg ==1 ){
00692       m_trkalg=1;
00693         //cout<<"using kal algorithm ................"<<endl;
00694         double poca_x, poca_y, poca_z;
00695         //retrieve  MdcKalTrackCol from TDS
00696         SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00697         if (!newtrkCol) {
00698            //cout<<"Could not find Kal trk -----------------------"<<endl;
00699            log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
00700            if(ntpFlag>0) return( StatusCode::SUCCESS);
00701         }
00702         if(ntpFlag>0) eventNO++;
00703         ntrk = newtrkCol->size();
00704         log << MSG::DEBUG << "trkCol size:  " << newtrkCol->size() <<endreq;
00705         if(for_calib == 1){if(newtrkCol->size() !=2 ) return( StatusCode::SUCCESS);}        
00706         RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
00707         for( ; trk != newtrkCol->end(); trk++) {
00708            if(ntpFlag>0) trackNO1++;
00709            vector<double> phlist;
00710            int layid, localwid, lr, wid, w0id;
00711            double theta, costhe, costhe_pol, sinenta, eangle, dd, dd_pol, ph, adc_raw, ddr, ddl, tdc_raw, driftT;
00712            double zhit, zhit_pol,  m_driftd, sintheta, rphi_path, driftd;
00713            double pathlength,sindip;
00714 
00715 
00716            RecMdcDedxHit* dedxhit = new RecMdcDedxHit;
00717            MdcDedxTrk trk_ex( **trk); 
00718            HepVector a_helix_trk = (*trk)->getTHelix();
00719            HepVector a = (*trk)->getFHelix();
00720            HepSymMatrix tkErrM = (*trk)->getFError();
00721            double m_d0E   = tkErrM[0][0];
00722            double m_phi0E = tkErrM[1][1];
00723            double m_cpaE  = tkErrM[2][2];
00724            double m_z0E   = tkErrM[3][3];
00725            HepPoint3D IP(0,0,0);
00726            Dedx_Helix exhel_trk(IP, a_helix_trk);
00727            Dedx_Helix exhel(IP, a);
00728  
00729            double phi0= a_helix_trk[1];
00730            double costhe_temp = cos(M_PI_2-atan(a_helix_trk[4]));     
00731            if(for_calib == 1 && abs(costhe_temp)>0.9 ) continue;
00732            costhe = cos(M_PI_2-atan(a_helix_trk[4]));    
00733            sintheta = sin(M_PI_2-atan(a_helix_trk[4]));
00734            //cout<<a_helix_trk[0]<<"    "<<a_helix_trk[1]<<"    "<<a_helix_trk[2]<<"    "<<a_helix_trk[3]<<"    "<<a_helix_trk[4]<<endl; 
00735            //cout<<a(1)<<"     "<<a(2) <<"      "<<a(3)<<"      "<<a(4)<<"     "<<a(5)<<endl; 
00736            float m_Pt = 1.0/fabs( a(3) );
00737            double m_P = m_Pt*sqrt(1 + a(5)*a(5)); 
00738            int charge = ( a(3) > 0 )? 1 : -1; 
00739       
00740            if(for_calib == 1 && runFlag ==3 &&(m_P>1.88||m_P<1.8)) continue; 
00741            
00742            if(for_calib == 1 && runFlag ==4 && runNO <MdcDedxParam::RUN5 &&(m_P>1.72||m_P<1.45)) continue;
00743            // for new data psipp data we need to modify the momentum scope
00744            if(for_calib == 1 && runFlag==4 && runNO>= MdcDedxParam::RUN5 &&(m_P>2.0||m_P<1.7)) continue;
00745 
00746           if(for_calib ==1&&runFlag ==4 && runNO < MdcDedxParam::RUN5&&m_P>1.72&&m_P<1.45) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"      m_P : "<<m_P<<endl;
00747           if(for_calib ==1&&runFlag ==4 && runNO>= MdcDedxParam::RUN5&&m_P>2.0&&m_P<1.7) cout<<"run : "<<runNO<<"      event : "<<event<<"      runFlag : "<<runFlag<<"        m_P : "<<m_P<<endl;
00748 
00749            log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
00750            log << MSG::DEBUG <<" %%%%%initial id = "<<(*trk)->trackId() <<endreq;
00751            
00752            HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
00753            log << MSG::DEBUG << "hits size  = " <<gothits.size()<<endreq;
00754            //if(ntpFlag>0) cout<<"hits size  = " <<gothits.size()<<endl;
00755            if(for_calib == 1){if(gothits.size()<20) return StatusCode::SUCCESS;}
00756            if(gothits.size()<2) continue;
00757            int Nhits = gothits.size();
00758            
00759            HelixSegRefVec::iterator it_gothit = gothits.begin();
00760            int N_hit=0;
00761            //if(ntpFlag>0&&gothits.size()<15)cout<<"event : "<<event<<"  hits : "<<gothits.size()<<endl;
00762            for( ; it_gothit != gothits.end(); it_gothit++) {
00763                 HepVector a_hit = (*it_gothit)->getHelixIncl();
00764                 //HepVector a_hit = (*it_gothit)->getHelixExcl();    
00765                 HepPoint3D IP(0,0,0);
00766                 Dedx_Helix exhel_hit(IP, a_hit);
00767  
00768                 mdcid = (*it_gothit)->getMdcId();
00769                 layid =  MdcID::layer(mdcid);
00770                 localwid = MdcID::wire(mdcid);
00771                 w0id = geosvc->Layer(layid)->Wirst();
00772                 wid = w0id + localwid; 
00773                 adc_raw = (*it_gothit)->getAdc();
00774                 //adc_raw = (*it_gothit)->getAdc();
00775                 tdc_raw = (*it_gothit)->getTdc();
00776                 log << MSG::INFO << "hit layer id  " << layid
00777                      << "   wire id  " << localwid
00778                      << "   adc_raw  " << adc_raw <<endreq;
00779                 zhit = (*it_gothit)->getZhit();
00780                 lr = (*it_gothit)->getFlagLR();
00781                 m_driftd = (*it_gothit)->getDD();
00782                 driftd = abs(m_driftd);  
00783                 dd = (*it_gothit)->getDocaIncl();//in which: getDocaIncl() include fit unused hit, 
00784                 //getDocaExcl() exclude hit unused in track fit;
00785                 if(lr == 2) cout<<"lr = "<<lr<<endl;
00786                 if(lr==-1 || lr == 2) dd =dd;
00787                 else if(lr ==1) dd =-dd;
00788                 //dd = dd/doca_norm[layid];     
00789                 sinenta = (*it_gothit)->getEntra();
00790                 //costhe = cos(M_PI_2-atan(a[4]));
00791                 //sintheta = sin(M_PI_2-atan(a[4]));  
00792                 eangle = (*it_gothit)->getEntra();
00793                 eangle = eangle/M_PI;
00794                 driftT = (*it_gothit)->getDT();
00795                 
00796                 //int ntpFlag = 1;              
00797                 //ph=exsvc->StandardCorrec(runFlag,ntpFlag,runNO,exhel,mdcid,adc_raw,dd,eangle,zhit,costhe);
00798                 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,exhel_hit,mdcid,adc_raw,dd,eangle,zhit,costhe); 
00799                 //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,tes );
00800                 pathlength=exsvc->PathL( ntpFlag, exhel_hit, layid, localwid, zhit_pol);
00801                 rphi_path=pathlength * sintheta;
00802                 //ph = adc_raw;
00803                 //pathlength=exsvc->PathL( ntpFlag, exhel_hit, layid, localwid, zhit_pol);
00804                 //rphi_path=pathlength * sintheta;        
00805                 //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,tes );
00806                 //if(ntpFlag>0 && pathlength == -1)
00807                 //cout<<"parameter1: "<<layid<<setw(8)<<localwid<<setw(15)<<driftd<<setw(15)<<rphi_path<<setw(15)<<pathlength<<setw(15)<<ph<<endl;
00808                 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
00809                 if(runNO<0){
00810                    if (layid<4) continue;
00811                    else if(layid >= 4 && layid <8){ 
00812                       if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00813                    }  
00814                    else if(layid >= 8){
00815                       if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00816                    }  
00817                 }     
00818                 else if(runNO>=0){
00819                    if(layid <4) continue;
00820                    else if (layid>=4 && layid<8){
00821                       if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00822                    }  
00823                    else if(layid >= 8){
00824                       if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00825                    }
00826                 }
00827                 //if(ntpFlag>0) cout<<layid<<"      "<<localwid<<"       "<<adc_raw<<"      "<<ph<<"        "<<rphi_path<<"      "<<driftd<<endl; 
00828                 dedxhit->setMdcId(mdcid);
00829                 dedxhit->setFlagLR(lr);
00830                 dedxhit->setTrkId(trk_ex.trk_id());
00831                 dedxhit->setPathLength(pathlength);
00832                 if (ph>0.0) {
00833                  if ( ntpFlag ==2 ) { 
00834                     
00835                     //m_charge1 = (*trk)->charge();
00836                     m_t01 = tes;
00837                     //m_driftT = driftT;
00838                     //m_tdc_raw = tdc_raw;
00839                     m_phraw = adc_raw;
00840                     m_exraw = ph;
00841                     //m_localwid = localwid;
00842                     m_wire = wid;
00843                     m_runNO1 = runNO;
00844                     m_nhit_hit = Nhits;
00845                     m_doca = dd;
00846                     m_driftdist = m_driftd;
00847                     m_sinenta = eangle;
00848                     m_costhe = costhe;
00849                     m_pathL = pathlength;
00850                     m_layer = layid;
00851                     m_ptrk1 = m_P;
00852                     //cout<<"adc_raw : "<<adc_raw<<"     "<<ph<<"     "<<dd<<"     "<<sinenta<<endl;
00853                     StatusCode status2 = m_nt2->write();
00854                     if ( status2.isFailure() )  {
00855                          log << MSG::ERROR << "Cannot fill Ntuple2" << endreq;
00856                     }
00857                  }
00858                  phlist.push_back(ph);
00859               }
00860            }
00861            trk_ex.set_phlist( phlist );
00862            ex_trks.push_back(trk_ex );
00863            m_trkalgs.push_back(m_trkalg);
00864            dedxhitList->push_back( dedxhit );
00865            phlist.clear();
00866            if(ntpFlag>0) trackNO3++;
00867            log << MSG::DEBUG << "size in processing:  " << ex_trks.size() <<endreq;
00868         }
00869     }
00870 
00871 
00872      //------------------------------- Rec algorithm finished ---------------------------------//
00873     int trk_cnt = ex_trks.size();
00874 
00875     if( ntrk != ex_trks.size() && dbgFlag > 0 )
00876         log << MSG::DEBUG <<" ntrk = "<<ntrk<<" : trk_cnt = "<<trk_cnt<< endreq;
00877     log << MSG::DEBUG << "dedxtrk size  " << ex_trks.size() <<endreq;
00878     
00879     int idedxid = 0;             // it: MdcDedxTrk* 
00880     for(std::vector<MdcDedxTrk>::iterator it = ex_trks.begin(); it != ex_trks.end(); it++, idedxid++ ) {
00882       if(dbgFlag > 0) {
00883         log << MSG::DEBUG <<"-------------------------------------------------------"<< endreq;
00884         log << MSG::DEBUG <<" trk id ="<< it->trk_id()<<" : P ="<<it->P()   //MdcDedxTrk::trk_id() 
00885             <<"  Pt ="<<it->Pt()<<" : theta ="<<it->theta()
00886             <<" : phi ="<<it->phi()<< " : charge = "<<it->charge()<<endreq;
00887       }
00888       
00889       //if(for_calib == 1){if(it->get_phlist().size()<10) return StatusCode::SUCCESS;}
00890       if(for_calib == 1){if(it->get_phlist().size()<10) continue;}
00891       
00892  
00893       RecMdcDedx*  resdedx = new RecMdcDedx;
00894       int usedhit=0;
00895       double dEdx_meas_hit = it-> cal_dedx_bitrunc(truncate, 0, usedhit);     
00896       if(m_alg==2) {
00897         dEdx_meas_hit = it-> cal_dedx_transform( 1);
00898       }
00899       if(m_alg==3) {
00900         dEdx_meas_hit = it-> cal_dedx_log( 1.0, 1);
00901       }
00902 
00903       
00904       double dEdx_meas_esat=0;
00905       double dEdx_meas_norun=0;
00906       double dEdx_meas=0;
00907       if(usedhit>0){
00908         dEdx_meas_esat = exsvc->StandardTrackCorrec(1,runFlag, ntpFlag, runNO, dEdx_meas_hit,  cos(it->theta()),  tes ); 
00909         dEdx_meas_norun = exsvc->StandardTrackCorrec(2,runFlag, ntpFlag, runNO, dEdx_meas_hit,  cos(it->theta()),  tes ); 
00910         dEdx_meas = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, dEdx_meas_hit,  cos(it->theta()),  tes );    
00911         //double dEdx_meas = dEdx_meas_temp;
00912         m_phtm1 = it->cal_dedx( truncate );
00913         m_pch1 = it->P();
00914         //double m_bitrunc_temp = it->cal_dedx_bitrunc(truncate,0,usedhit );
00915         m_bitrunc1 = exsvc->StandardTrackCorrec(0, runFlag, ntpFlag, runNO, dEdx_meas_hit,  cos(it->theta()),  tes );
00916         //m_bitrunc1 = m_bitrunc_temp;
00917         
00918         //cout<<"m_bitrunc_temp = "<<m_bitrunc_temp<<"      m_bitrunc1 = "<<m_bitrunc1<<"      dEdx_meas_temp= "<<dEdx_meas_temp<<"      dEdx_meas= "<<dEdx_meas<<endl;
00919         //m_median1 = it->cal_dedx_median( truncate );
00920         //m_geometric1 = it->cal_dedx_geometric( truncate );
00921         //m_geometric_trunc1 = it->cal_dedx_geometric_trunc( truncate );
00922         //m_harmonic1 = it->cal_dedx_harmonic( truncate );
00923         //m_harmonic_trunc1 = it->cal_dedx_harmonic_trunc( truncate );
00924         //m_transform1 = it->cal_dedx_transform( 0 );
00925         //m_logtrunc1 = it->cal_dedx_log( 1.0, 0);
00926         //double dEdx_meas = 1000;
00927         //cout<<"nsample = "<<it->nsample()<<"         ph size = "<< it->get_phlist().size()<<endl;     
00928       }
00929 
00930      if( dEdx_meas > 0 ) {
00931        //  ----- fill Mdc dEdx results  ------  
00932        //cout<<"check here --------------------------"<<endl;
00933        int trk_recalg=m_trkalgs[idedxid];
00934        
00935        it->set_dEdx( landau , dEdx_meas, trk_recalg, runFlag, vFlag , tes, Curve_Para, Sigma_Para); 
00936        
00937        //set MdcDedx
00938        //int dedxhit_used = int (it->get_phlist().size()*0.7); 
00939        resdedx->setDedxHit(dEdx_meas_hit);
00940        resdedx->setDedxEsat(dEdx_meas_esat);
00941        resdedx->setDedxNoRun(dEdx_meas_norun);
00942        resdedx->setDedxMoment(it->P());
00943        //if(ntpFlag>0) cout<<dEdx_meas_hit<<"      "<<dEdx_meas_esat<<"      "<<dEdx_meas_norun<<"      "<<it->P()<<endl;
00944        resdedx->setTrackId( it->trk_id() );
00945        resdedx->setProbPH( dEdx_meas );
00946        resdedx->setNormPH( dEdx_meas/550.0 );
00947        resdedx->setDedxExpect( it->pexpect() );
00948        resdedx->setSigmaDedx( it->pexp_sigma() );
00949        resdedx->setPidProb( it->pprob() );
00950        resdedx->setChi( it->pchi_dedx() );
00951        resdedx->setNumTotalHits(it->nsample() );
00952        resdedx->setNumGoodHits(usedhit);
00953        //resdedx->setStatus( it->quality() );
00954        resdedx->setStatus(trk_recalg );
00955        resdedx->setTruncAlg( m_alg );
00956        
00957          //        resdedx->setTrack(trk_ex );
00958        float maxProb=-0.01;
00959        m_parttype1=0.0;
00960        for( int ai = 0; ai < 5; ai++ ) {
00961          m_probpid1[ai]=resdedx->getPidProb(ai+1);
00962          if(m_probpid1[ai]>maxProb) {
00963            maxProb=m_probpid1[ai];
00964            m_parttype1 = ai+1;
00965          }
00966          log<< MSG::INFO<<"test dedxcpl: exp: "<<resdedx->getDedxExpect(ai+1)
00967             <<"  sigma:  "<<resdedx->getSigmaDedx(ai+1)
00968             <<"    prob:  "<< resdedx->getPidProb(ai+1)<<endreq;
00969        }
00970        
00971        m_expect1 = resdedx->getDedxExpect(m_parttype1);
00972        m_sigma1 = resdedx->getSigmaDedx(m_parttype1); 
00973        //m_expect1 = resdedx->getDedxExpect(1);
00974        //m_sigma1 = resdedx->getSigmaDedx(1);
00975        
00976        m_chidedx1  = resdedx->chi(m_parttype1-1); 
00977        m_nhit1 = it->nsample();
00978        m_sintheta1 = sin(it->theta());
00979        //m_sintheta1 = it->theta();
00980        m_costheta1 = cos(it->theta());
00981        
00982        resdedx->setParticleId(m_parttype1);
00983        log << MSG::DEBUG << "Pid before register"
00984            <<resdedx->particleId()<<endreq;
00985        if(ntpFlag ==2 ){    
00986          m_pch=m_pch1;
00987          m_phtm=m_phtm1;
00988          m_bitrunc=m_bitrunc1;
00989          m_dEdx_meas = dEdx_meas;
00990          //cout<<"m_bitrunc = "<<m_bitrunc<<endl;
00991          //m_median=m_median1;
00992          //m_geometric= m_geometric1;
00993          //m_geometric_trunc= m_geometric_trunc1;
00994          //m_harmonic=m_harmonic1;
00995          //m_harmonic_trunc=m_harmonic_trunc1;
00996          //m_transform=m_transform1;
00997          //m_logtrunc=m_logtrunc1;
00998          for(int i=0;i<5;i++)
00999            {m_probpid[i]= m_probpid1[i];
01000            m_expectid[i]= resdedx->getDedxExpect(i+1);
01001            m_sigmaid[i]=resdedx->getSigmaDedx(i+1); 
01002            }
01003          m_parttype=m_parttype1;
01004          m_expect= m_expect1;
01005          m_sigma=m_sigma1;
01006          m_chidedx=m_chidedx1;
01007          
01008          
01009          m_chidedxe=resdedx->chi(0); 
01010          m_chidedxmu=resdedx->chi(1);
01011          m_chidedxpi=resdedx->chi(2);
01012          m_chidedxk=resdedx->chi(3);
01013          m_chidedxp=resdedx->chi(4);        
01014          
01015          //cout<<"m_chidedxe = "<<m_chidedxe<<"     "<<resdedx->chiE()<<endl;    
01016          //m_chidedxe=resdedx->chiE();
01017          //m_chidedxmu=resdedx->chiMu();
01018          //m_chidedxpi=resdedx->chiPi();
01019          //m_chidedxk=resdedx->chiK();
01020          //m_chidedxp=resdedx->chiP();
01021          m_sintheta= m_sintheta1;
01022          m_costheta= m_costheta1;
01023          m_nhit=m_nhit1;
01024          m_ndedxhit=it->get_phlist().size();
01025          m_charge = it->charge();
01026          m_runNO = runNO;
01027          m_evtNO = event;
01028          m_t0 = tes;
01029          //cout<<"event : "<<event<<"       m_charge : "<<m_charge<<"      m_pch : "<<m_pch<<endl;
01030          //m_poca_x = poca_x;
01031          //m_poca_y = poca_y;
01032          //m_poca_z = poca_z;
01033          //cout<<"trk id =   "<<it->trk_id()<<"             m_bitrunc = "<<m_bitrunc<<endl; 
01034          StatusCode status = m_nt1->write();
01035          if ( status.isFailure() )  {
01036            log << MSG::ERROR << "Cannot fill Ntuple1" << endreq;
01037          }
01038        }     
01039      } 
01040      
01041      log << MSG::INFO << "-----------------Summary of this ExTrack---"
01042          << endreq <<"dEdx_mean  "<<dEdx_meas <<"   expct(4)   "
01043          <<it->expect(4) << endreq <<"  sigma(4)   "<<it->exp_sigma(4)
01044          <<"   chi(4)   "<<it->chi_dedx(4)
01045          <<"   prob(4)   "<<it->prob(4)<<endreq;
01046      
01047      dedxList->push_back( resdedx );
01048       
01049     } //ExTrk loop end
01050     
01051   }   //  doNewGlobal==0 END . . . 
01052   
01053   // check MdcDedxCol registered
01054   SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
01055   if (!newexCol) { 
01056     log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
01057     return( StatusCode::SUCCESS);
01058   }
01059   log << MSG::DEBUG << "Begin to check RecMdcDedxCol"<<endreq; 
01060   RecMdcDedxCol::iterator it_dedx = newexCol->begin();
01061   for( ; it_dedx != newexCol->end(); it_dedx++) {
01062            log << MSG::INFO << "retrieved MDC dE/dx:"
01063                << "dEdx Id: " << (*it_dedx)->trackId()
01064                << "part Id: " << (*it_dedx)->particleId()
01065                << "   measured dEdx:  " << (*it_dedx)->probPH()
01066                << "   dEdx std:  " << (*it_dedx)->normPH()
01067                << "   nhits:  " << (*it_dedx)->numTotalHits()
01068                << endreq
01069                << "   Electron: expect: " << (*it_dedx)->getDedxExpect(1)
01070                << "   sigma: " << (*it_dedx)->getSigmaDedx(1)
01071                << "   chi: " << (*it_dedx)->chi(1)
01072                << "   prob: " << (*it_dedx)->getPidProb(1) << endreq
01073                << "   Muon: expect: " << (*it_dedx)->getDedxExpect(2)
01074                << "   sigma: " << (*it_dedx)->getSigmaDedx(2)
01075                << "   chi: " << (*it_dedx)->chi(2)
01076                << "   prob: " << (*it_dedx)->getPidProb(2) << endreq
01077                << "   Pion: expect: " << (*it_dedx)->getDedxExpect(3)
01078                << "   sigma: " << (*it_dedx)->getSigmaDedx(3)
01079                << "   chi: " << (*it_dedx)->chi(3)
01080                << "   prob: " << (*it_dedx)->getPidProb(3) << endreq
01081                << "   Kaon: expect: " << (*it_dedx)->getDedxExpect(4)
01082                << "   sigma: " << (*it_dedx)->getSigmaDedx(4)
01083                << "   chi: " << (*it_dedx)->chi(4)
01084                << "   prob: " << (*it_dedx)->getPidProb(4) << endreq
01085                << "   Proton: expect: " << (*it_dedx)->getDedxExpect(5)
01086                << "   sigma: " << (*it_dedx)->getSigmaDedx(5)
01087                << "   chi: " << (*it_dedx)->chi(5)
01088                << "   prob: " << (*it_dedx)->getPidProb(5) << endreq;
01089   }
01090   
01091   return StatusCode::SUCCESS;
01092 
01093   }

StatusCode MdcDedxRecon::finalize (   ) 
 

StatusCode MdcDedxRecon::finalize (   ) 
  

00279                                  {
00280    MsgStream log(msgSvc(), name());
00281    ex_trks.clear();
00282    m_trkalgs.clear();
00283    log << MSG::DEBUG <<"MdcDedxRecon terminated!!!"<< endreq;
00284    //cout<<"total event number is : "<<eventNO<<endl;
00285    //cout<<"total track number is : "<<trackNO1<<"        "<<trackNO2<<"        "<<trackNO3<<endl;
00286    //   SmartDataPtr<IHistogram1D> hist1D( histoSvc(), "/stat/ph/10");
00287    // if(hist1D != 0) histoSvc()->print( hist1D );
00288    return StatusCode::SUCCESS;
00289 }

StatusCode MdcDedxRecon::initialize (   ) 
 

StatusCode MdcDedxRecon::initialize (   ) 
 

cleate MdcDedxWire and MdcDedxLayer. 

00100 { 
00101  
00102    MsgStream log(msgSvc(), name());
00103    eventNO =0;
00104    trackNO1 =0;  
00105    trackNO2 =0;
00106    trackNO3 =0;
00107 
00108    //IDedxCurSvc*  vtxsvc;
00109    StatusCode cursvc = service("DedxCurSvc", dedxcursvc); 
00110    //Gaudi::svcLocator()->service("DedxCurSvc", vtxsvc);
00111    if (cursvc ==  StatusCode::SUCCESS) {
00112       log << MSG::INFO <<"DedxCurSvc is running"<<endl;
00113    } else {
00114       log << MSG::ERROR <<"DedxCurSvc is failed"<<endl;
00115       return StatusCode::SUCCESS;
00116    }
00117  
00118    StatusCode scex = service("DedxCorrecSvc", exsvc);
00119    if (scex ==  StatusCode::SUCCESS) {              
00120       log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endl;
00121    } else {
00122       log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endl;
00123       return StatusCode::SUCCESS;
00124    }
00125    exsvc->set_flag( calib_flag );
00126  
00127    log << MSG::INFO <<"--------------------<  MdcDedxRecon  V2.1   >---------------------"<< endreq;
00128    log << MSG::INFO <<"#######  correction for the wire current dependence   #######"<< endreq;
00129    log << MSG::INFO <<"#######      correction for the new z dependence      #######"<< endreq;
00130    log << MSG::INFO <<"-------------------------------------------------------------"<< endreq;
00131    log << MSG::INFO <<"++++++++++++++++++++[  initialization  ]+++++++++++++++++++++"<< endreq;
00132    log << MSG::INFO << "MdcDedxRecon has been initialized with dbgFlag: "
00133            << dbgFlag <<  endreq;
00134 
00135    log << MSG::INFO << "MdcDedxRecon has been initialized with ntpFlag: "
00136            << ntpFlag <<  endreq;   
00137 
00138    log << MSG::INFO << "MdcDedxRecon has been initialized with mincut_dedx: "
00139            << mincut_dedx <<  endreq;
00140    log << MSG::INFO << "MdcDedxRecon has been initialized with minhit: "
00141            << minhit <<  endreq;
00142    log << MSG::INFO << "MdcDedxRecon has been initialized with landau: "
00143            << landau <<  endreq;
00144    if( landau == 0 ){
00145            truncate = 0.7;
00146            x_rej = 0;
00147    }
00148    log << MSG::INFO << "MdcDedxRecon has been initialized with truncate: "
00149            << truncate <<  endreq;
00150    log << MSG::INFO << "MdcDedxRecon has been initialized with stsw: "
00151            << stsw <<  endreq;
00152    log << MSG::INFO << "MdcDedxRecon has been initialized with x_rej: "
00153            << x_rej <<  endreq;
00154    m_rec_sw = x_rej + truncate;
00155    log << MSG::INFO << "MdcDedxRecon has been initialized with calib_flag: "
00156            << calib_flag <<  endreq;
00157    log << MSG::INFO << "MdcDedxRecon has been initialized with doNewGlobal: "
00158            << doNewGlobal <<  endreq;
00159 
00161    if( !ex_calib ) ex_calib = MdcDedxCorrection::getMdcDedxCorrection();
00163 
00164    log << MSG::DEBUG  << "+++++++++++ MdcDedxRecon initialization end ++++++++++++ "<< endreq;
00165    if( truncate <= 0.0 || 1.0 < truncate ){
00166            log << MSG::FATAL <<" Initialization ERROR of truncate  = "<<truncate<< endreq;
00167            log << MSG::FATAL <<" truncate must be within 0.0 to 1.0 ! "<< endreq;
00168            log << MSG::FATAL <<" Please stop exec. "<< endreq;
00169    }
00170    log << MSG::DEBUG <<"-------------------------------------------------------------"<< endreq;
00171    log << MSG::DEBUG <<"MdcDedxRecon init 2nd part!!!"<< endreq;
00172    //   begin_run_stat = ex_calib->read_data();
00173    if( dbgFlag > 0 ){
00174            log << MSG::DEBUG <<"init 2nd part : RunNo="<< endreq;
00175    }
00176 #ifdef DBCALIB
00177    ex_calib->getCalib();
00178 #endif
00179    //   his_curv= histoSvc()->book("/stat/ph",30,"dEdx vs momentum",100,0.01,10.6,50,0.0,0.05); 
00180 
00181    if( ntpFlag ==2 ){
00182        NTuplePtr  nt1(ntupleSvc(),"FILE103/n103");
00183        StatusCode status;
00184        if ( nt1 ) 
00185             m_nt1 = nt1;
00186        else {
00187              m_nt1= ntupleSvc()->book("FILE103/n103",CLID_ColumnWiseTuple,"dEdx vs momentum");
00188              if ( m_nt1 )  {
00189                   status = m_nt1->addItem("ptrk",m_pch);
00190                   status = m_nt1->addItem("sintheta",m_sintheta);          
00191                   status = m_nt1->addItem("costheta",m_costheta);
00192                   status = m_nt1->addItem("dEdx_meas", m_dEdx_meas);
00193                   status = m_nt1->addItem("nhit",m_nhit);
00194                   status = m_nt1->addItem("ndedxhit",m_ndedxhit);
00195                   status = m_nt1->addItem("charge",m_charge);        
00196                   status = m_nt1->addItem("runNO",m_runNO);  
00197                   status = m_nt1->addItem("evtNO",m_evtNO); 
00198                   status = m_nt1->addItem("t0",m_t0); 
00199           
00200                   status = m_nt1->addItem("dedx_tm",m_phtm);
00201                   status = m_nt1->addItem("bitrunc",m_bitrunc);
00202                   //status = m_nt1->addItem("median",m_median);
00203                   //status = m_nt1->addItem("geom",m_geometric);
00204                   //status = m_nt1->addItem("geom_tr",m_geometric_trunc);
00205                   //status = m_nt1->addItem("harm",m_harmonic);
00206                   //status = m_nt1->addItem("harm_tr",m_harmonic_trunc);
00207                   //status = m_nt1->addItem("transf",m_transform);
00208                   //status = m_nt1->addItem("logex",m_logtrunc);          
00209 
00210                   status = m_nt1->addItem("trackId",m_trackId);                  
00211                   status = m_nt1->addItem("partid",5,m_probpid);
00212                   status = m_nt1->addItem("expectid",5,m_expectid);
00213                   status = m_nt1->addItem("sigmaid",5,m_sigmaid); 
00214 
00215                   status = m_nt1->addItem("type",m_parttype);
00216                   status = m_nt1->addItem("expect",m_expect);
00217                   status = m_nt1->addItem("sigma",m_sigma);          
00218                   status = m_nt1->addItem("chidedx",m_chidedx);          
00219                   status = m_nt1->addItem("chidedx_e",m_chidedxe);
00220                   status = m_nt1->addItem("chidedx_mu",m_chidedxmu);
00221                   status = m_nt1->addItem("chidedx_pi",m_chidedxpi);
00222                   status = m_nt1->addItem("chidedx_k",m_chidedxk);
00223                   status = m_nt1->addItem("chidedx_p",m_chidedxp);
00224                   
00225                   //status = m_nt1->addItem("poca_x",m_poca_x);
00226                   //status = m_nt1->addItem("poca_y",m_poca_y);
00227                   //status = m_nt1->addItem("poca_z",m_poca_z); 
00228              }
00229        }
00230 
00231        NTuplePtr  nt2(ntupleSvc(),"FILE103/n102");
00232        if ( nt2 )  m_nt2 = nt2;
00233        else {
00234              m_nt2= ntupleSvc()->book("FILE103/n102",CLID_RowWiseTuple,"pulae height raw");
00235              if ( m_nt2 )  {
00236                   status = m_nt2->addItem("charge",m_charge1);
00237                   status = m_nt2->addItem("adc_raw",m_phraw);
00238                   status = m_nt2->addItem("exraw",m_exraw);       
00239                   status = m_nt2->addItem("runNO1",m_runNO1);
00240                   status = m_nt2->addItem("nhit_hit", m_nhit_hit);
00241                   status = m_nt2->addItem("wire",m_wire);
00242                   status = m_nt2->addItem("doca",m_doca);
00243                   status = m_nt2->addItem("driftdist",m_driftdist);
00244                   status = m_nt2->addItem("sinenta",m_sinenta);
00245                   status = m_nt2->addItem("costhe",m_costhe);
00246                   status = m_nt2->addItem("path_rphi",m_pathL);       
00247                   status = m_nt2->addItem("layer",m_layer);
00248                   status = m_nt2->addItem("ptrk1",m_ptrk1);
00249                   status = m_nt2->addItem("t01",m_t01); 
00250                   //status = m_nt2->addItem("tdc_raw",m_tdc_raw);  
00251                   //status = m_nt2->addItem("driftT",m_driftT);
00252                   //status = m_nt1->addItem("trackId1",m_trackId1);
00253                   //status = m_nt2->addItem("localwid",m_localwid);
00254              }
00255        } 
00256    }
00257    return StatusCode::SUCCESS;
00258 }

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

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

01216                                                                                                                      {
01217 
01218   vector<double> phlist;
01219   int layid, localwid, lr, wid, w0id;
01220   double theta, costhe, costhe_pol, sinenta, eangle, dd, dd_pol, ph, adc_raw, ddr, ddl, tdc_raw, driftT;
01221   double zhit, zhit_pol,  m_driftd, sintheta, rphi_path, driftd;
01222   double pathlength,sindip;
01223   
01224   
01225   RecMdcDedxHit* dedxhit = new RecMdcDedxHit;
01226   MdcDedxTrk trk_ex( **trk); 
01227   HepVector a_helix_trk = (*trk)->getTHelix();
01228   HepVector a = (*trk)->getFHelix();
01229   HepSymMatrix tkErrM = (*trk)->getFError();
01230   double m_d0E   = tkErrM[0][0];
01231   double m_phi0E = tkErrM[1][1];
01232   double m_cpaE  = tkErrM[2][2];
01233   double m_z0E   = tkErrM[3][3];
01234   HepPoint3D IP(0,0,0);
01235   Dedx_Helix exhel_trk(IP, a_helix_trk);
01236   Dedx_Helix exhel(IP, a);
01237   
01238   double phi0= a_helix_trk[1];
01239   double costhe_temp = cos(M_PI_2-atan(a_helix_trk[4]));     
01240   
01241   costhe = cos(M_PI_2-atan(a_helix_trk[4]));    
01242   sintheta = sin(M_PI_2-atan(a_helix_trk[4]));
01243   //cout<<a_helix_trk[0]<<"    "<<a_helix_trk[1]<<"    "<<a_helix_trk[2]<<"    "<<a_helix_trk[3]<<"    "<<a_helix_trk[4]<<endl; 
01244   //cout<<a(1)<<"     "<<a(2) <<"      "<<a(3)<<"      "<<a(4)<<"     "<<a(5)<<endl; 
01245   float m_Pt = 1.0/fabs( a(3) );
01246   double m_P = m_Pt*sqrt(1 + a(5)*a(5)); 
01247   int charge = ( a(3) > 0 )? 1 : -1; 
01248   
01249   log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
01250   log << MSG::DEBUG <<" %%%%%initial id = "<<(*trk)->trackId() <<endreq;
01251   
01252   HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
01253   log << MSG::DEBUG << "hits size  = " <<gothits.size()<<endreq;
01254   //if(ntpFlag>0) cout<<"hits size  = " <<gothits.size()<<endl;
01255   
01256 
01257   HelixSegRefVec::iterator it_gothit = gothits.begin();
01258   int N_hit=0;
01259   
01260   for( ; it_gothit != gothits.end(); it_gothit++) {
01261     HepVector a_hit = (*it_gothit)->getHelixIncl();
01262     //HepVector a_hit = (*it_gothit)->getHelixExcl();    
01263     HepPoint3D IP(0,0,0);
01264     Dedx_Helix exhel_hit(IP, a_hit);
01265     
01266     mdcid = (*it_gothit)->getMdcId();
01267     layid =  MdcID::layer(mdcid);
01268     localwid = MdcID::wire(mdcid);
01269     w0id = geosvc->Layer(layid)->Wirst();
01270     wid = w0id + localwid; 
01271     adc_raw = (*it_gothit)->getAdc();
01272     //adc_raw = (*it_gothit)->getAdc();
01273     tdc_raw = (*it_gothit)->getTdc();
01274     log << MSG::INFO << "hit layer id  " << layid
01275         << "   wire id  " << localwid
01276         << "   adc_raw  " << adc_raw <<endreq;
01277     zhit = (*it_gothit)->getZhit();
01278     lr = (*it_gothit)->getFlagLR();
01279     m_driftd = (*it_gothit)->getDD();
01280     driftd = abs(m_driftd);  
01281     dd = (*it_gothit)->getDocaIncl();//in which: getDocaIncl() include fit unused hit, 
01282     //getDocaExcl() exclude hit unused in track fit;
01283     if(lr == 2) cout<<"lr = "<<lr<<endl;
01284     if(lr==-1 || lr == 2) dd =dd;
01285     else if(lr ==1) dd =-dd;
01286     //dd = dd/doca_norm[layid];     
01287     sinenta = (*it_gothit)->getEntra();
01288     //costhe = cos(M_PI_2-atan(a[4]));
01289     //sintheta = sin(M_PI_2-atan(a[4]));  
01290     eangle = (*it_gothit)->getEntra();
01291     eangle = eangle/M_PI;
01292     driftT = (*it_gothit)->getDT();
01293     
01294     //int ntpFlag = 1;              
01295     //ph=exsvc->StandardCorrec(runFlag,ntpFlag,runNO,exhel,mdcid,adc_raw,dd,eangle,zhit,costhe);
01296     ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,exhel_hit,mdcid,adc_raw,dd,eangle,zhit,costhe); 
01297     //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,tes );
01298     pathlength=exsvc->PathL( ntpFlag, exhel_hit, layid, localwid, zhit_pol);
01299     rphi_path=pathlength * sintheta;
01300     //ph = adc_raw;
01301     //pathlength=exsvc->PathL( ntpFlag, exhel_hit, layid, localwid, zhit_pol);
01302     //rphi_path=pathlength * sintheta;        
01303     //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,tes );
01304     //if(ntpFlag>0 && pathlength == -1)
01305     //cout<<"parameter1: "<<layid<<setw(8)<<localwid<<setw(15)<<driftd<<setw(15)<<rphi_path<<setw(15)<<pathlength<<setw(15)<<ph<<endl;
01306     if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
01307     if(runNO<0){
01308       if (layid<4) continue;
01309       else if(layid >= 4 && layid <8){ 
01310         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01311       }  
01312       else if(layid >= 8){
01313         if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01314       }  
01315     }     
01316     else if(runNO>=0){
01317       if(layid <4) continue;
01318       else if (layid>=4 && layid<8){
01319         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01320       }  
01321       else if(layid >= 8){
01322         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01323       }
01324     }
01325     //if(ntpFlag>0) cout<<layid<<"      "<<localwid<<"       "<<adc_raw<<"      "<<ph<<"        "<<rphi_path<<"      "<<driftd<<endl; 
01326     dedxhit->setMdcId(mdcid);
01327     dedxhit->setFlagLR(lr);
01328     dedxhit->setTrkId(trk_ex.trk_id());
01329     dedxhit->setPathLength(pathlength);
01330     if (ph>0.0) {
01331       phlist.push_back(ph);
01332     }
01333   }
01334   trk_ex.set_phlist( phlist );
01335   ex_trks.push_back(trk_ex );
01336   m_trkalgs.push_back(m_trkalg);
01337   //dedxhitList->push_back( dedxhit );
01338   phlist.clear();
01339           
01340   log << MSG::DEBUG << "size in processing:  " << ex_trks.size() <<endreq;  
01341 }

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

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

01103                                                                                                                   {
01104   
01105   RecMdcDedxHit* dedxhit = new RecMdcDedxHit;
01106   MdcDedxTrk trk_ex( **trk); 
01107   HepVector a = (*trk)->helix();
01108     
01109   HepSymMatrix tkErrM = (*trk)->err();
01110   double m_d0E   = tkErrM[0][0];
01111   double m_phi0E = tkErrM[1][1];
01112   double m_cpaE  = tkErrM[2][2];
01113   double m_z0E   = tkErrM[3][3]; 
01114   
01115   HepPoint3D IP(0,0,0);
01116   Dedx_Helix exhel(IP, a);
01117   for(int i=0;i<5;i++) {
01118     log << MSG::DEBUG <<" helix["<<i<<"] = "<<(*trk)->helix()[i] << endreq;
01119   }
01120 
01121   double theta_temp = M_PI_2-atan(a[4]);
01122   double costhe_temp = cos(M_PI_2-atan(a[4]));
01123     
01124   double m_Pt = 1.0/fabs( a(3) );
01125   double m_P = m_Pt*sqrt(1 + a(5)*a(5));
01126     
01127   log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
01128   log << MSG::DEBUG <<" %%%%%initial id = "<<(*trk)->trackId() <<endreq;
01129   // 2005/5 make phlist here... 
01130   vector<double> phlist;
01131   int layid, localwid,lr, wid, w0id;
01132   double theta, costhe, costhe_pol, sinenta, eangle, dd, dd_pol, ph, adc_raw, ddr, ddl, tdc_raw, driftT;
01133   double zhit, zhit_pol,  m_driftd, sintheta, rphi_path, driftd; 
01134   double pathlength,sindip;
01135   HitRefVec gothits= (*trk)->getVecHits();
01136   log << MSG::DEBUG << "hits size  = " <<gothits.size()<<endreq;
01137 
01138     
01139   HitRefVec::iterator it_gothit = gothits.begin();
01140   for( ; it_gothit != gothits.end(); it_gothit++) {
01141     mdcid = (*it_gothit)->getMdcId();
01142     layid =  MdcID::layer(mdcid);
01143     localwid = MdcID::wire(mdcid);
01144     w0id = geosvc->Layer(layid)->Wirst();
01145     wid = w0id + localwid;        
01146     adc_raw = (*it_gothit)->getAdc();
01147     //adc_raw = (*it_gothit)->getAdc()*1000000;
01148     tdc_raw = (*it_gothit)->getTdc();
01149     driftT = (*it_gothit)->getDriftT();
01150     //cout<<"adc = "<<adc_raw<<endl;
01151     log << MSG::INFO << "hit layer id  " << layid 
01152         << "   wire id  " << localwid
01153         << "   adc_raw  " << adc_raw <<endreq; 
01154     zhit = (*it_gothit)->getZhit();
01155     zhit_pol = fabs(zhit); 
01156     lr = (*it_gothit)->getFlagLR();
01157     if(lr == 2) cout<<"lr = "<<lr<<endl;
01158     if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft();
01159     if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight();
01160     driftd = abs(m_driftd);
01161     dd = (*it_gothit)->getDoca();
01162     if(lr == 0 || lr == 2 ) dd = -abs(dd);
01163     if(lr == 1 ) dd = abs(dd);              
01164     sinenta = (*it_gothit)->getEntra();
01165     eangle = sinenta/M_PI;
01166     theta = M_PI_2-atan(a[4]);
01167     costhe = cos(M_PI_2-atan(a[4]));
01168     sintheta = sin(M_PI_2-atan(a[4]));
01169     //ph = exsvc->StandardCorrec( runFlag, ntpFlag, runNO, exhel, mdcid, adc_raw, dd, eangle, zhit, costhe);
01170     ph = exsvc->StandardHitCorrec(0, runFlag,ntpFlag,runNO,exhel,mdcid,adc_raw,dd,eangle,zhit,costhe);
01171     //ph = exsvc->StandardTrackCorrec(runFlag, ntpFlag, runNO, ph, costhe,  tes ); 
01172     //ph = adc_raw;
01173     pathlength = exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit_pol); 
01174     rphi_path = pathlength * sintheta;
01175    
01176     if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
01177     if(runNO<0){  
01178       if(layid < 4) continue;
01179       else if(layid>=4&&layid <8){
01180         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || dd>Iner_DriftDistCut) continue;
01181       }
01182       else if(layid >= 8){
01183         if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01184       }
01185     }
01186     else if(runNO>=0){
01187       if(layid <4) continue;
01188       else if (layid>=4 && layid<8){
01189         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || dd>Iner_DriftDistCut) continue;
01190       }
01191       else if(layid >= 8){
01192         if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01193       }
01194     }
01195     dedxhit->setMdcId(mdcid);
01196     dedxhit->setFlagLR(lr);
01197     dedxhit->setTrkId(trk_ex.trk_id());
01198     dedxhit->setPathLength(pathlength);
01199     if (ph>0.0) {
01200       phlist.push_back(ph);
01201     }
01202   
01203   }//end of hits loop
01204   
01205   trk_ex.set_phlist( phlist ); 
01206   ex_trks.push_back(trk_ex );
01207   m_trkalgs.push_back(m_trkalg);
01208   //dedxhitList->push_back( dedxhit );
01209   phlist.clear();
01210   
01211   log << MSG::DEBUG << "size in processing:  " << ex_trks.size() <<endreq;
01212 
01213 }

int MdcDedxRecon::ReadRaw (  void   ) 
 

int MdcDedxRecon::ReadRaw (  void   ) 
 

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

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

01343                                                                                                      {
01344 
01345   SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
01346   if (!newtrkCol) {
01347     log << MSG::WARNING << "Could not find RecMdcTrackCol in swithcto mdctrack" << endreq;
01348   }
01349   
01350   RecMdcTrackCol::iterator trk = newtrkCol->begin();     
01351   for( ; trk != newtrkCol->end(); trk++) {
01352   
01353     if( (*trk)->trackId()==trkid){
01354       HitRefVec gothits= (*trk)->getVecHits();
01355       if(gothits.size()>=2)
01356         mdctrackrec(trk,mdcid,runNO,runFlag,log);
01357     }
01358 
01359   }
01360   
01361 }

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

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

01096                               {
01097   return ex_trks;
01098 }

Member Data Documentation

int MdcDedxRecon::begin_run_stat [private]
 

int MdcDedxRecon::calib_flag
 

int MdcDedxRecon::cosmic
 

float MdcDedxRecon::cut_sigma
 

int MdcDedxRecon::dbgFlag
 

IDedxCurSvc* MdcDedxRecon::dedxcursvc [private]
 

IDedxCurSvc* MdcDedxRecon::dedxcursvc [private]
 

int MdcDedxRecon::doNewGlobal
 

MdcDedxCorrection* MdcDedxRecon::ex_calib [private]
 

MdcDedxCorrection* MdcDedxRecon::ex_calib [private]
 

std::vector<MdcDedxTrk> MdcDedxRecon::ex_trks [private]
 

std::vector<MdcDedxTrk> MdcDedxRecon::ex_trks [private]
 

IDedxCorrecSvc* MdcDedxRecon::exsvc [private]
 

IDedxCorrecSvc* MdcDedxRecon::exsvc [private]
 

int MdcDedxRecon::for_calib
 

IMdcGeomSvc* MdcDedxRecon::geosvc [private]
 

IMdcGeomSvc* MdcDedxRecon::geosvc [private]
 

IHistogram1D* MdcDedxRecon::his_aver [private]
 

IHistogram1D* MdcDedxRecon::his_aver [private]
 

IHistogram1D* MdcDedxRecon::his_ph [private]
 

IHistogram1D* MdcDedxRecon::his_ph [private]
 

int MdcDedxRecon::landau
 

int MdcDedxRecon::m_alg
 

NTuple::Item<float> MdcDedxRecon::m_bitrunc [private]
 

NTuple::Item<float> MdcDedxRecon::m_bitrunc [private]
 

float MdcDedxRecon::m_bitrunc1
 

NTuple::Item<float> MdcDedxRecon::m_charge [private]
 

NTuple::Item<float> MdcDedxRecon::m_charge [private]
 

NTuple::Item<double> MdcDedxRecon::m_charge1 [private]
 

NTuple::Item<double> MdcDedxRecon::m_charge1 [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedx [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedx [private]
 

float MdcDedxRecon::m_chidedx1
 

NTuple::Item<float> MdcDedxRecon::m_chidedxe [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxe [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxk [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxk [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxmu [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxmu [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxp [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxp [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxpi [private]
 

NTuple::Item<float> MdcDedxRecon::m_chidedxpi [private]
 

NTuple::Item<float> MdcDedxRecon::m_costhe [private]
 

NTuple::Item<float> MdcDedxRecon::m_costhe [private]
 

NTuple::Item<float> MdcDedxRecon::m_costheta [private]
 

NTuple::Item<float> MdcDedxRecon::m_costheta [private]
 

float MdcDedxRecon::m_costheta1
 

NTuple::Item<float> MdcDedxRecon::m_dEdx_meas [private]
 

NTuple::Item<float> MdcDedxRecon::m_dEdx_meas [private]
 

NTuple::Item<float> MdcDedxRecon::m_doca [private]
 

NTuple::Item<float> MdcDedxRecon::m_doca [private]
 

NTuple::Item<float> MdcDedxRecon::m_driftdist [private]
 

NTuple::Item<float> MdcDedxRecon::m_driftdist [private]
 

NTuple::Item<float> MdcDedxRecon::m_driftT [private]
 

NTuple::Item<float> MdcDedxRecon::m_driftT [private]
 

NTuple::Item<float> MdcDedxRecon::m_evtNO [private]
 

NTuple::Item<float> MdcDedxRecon::m_evtNO [private]
 

NTuple::Item<float> MdcDedxRecon::m_expect [private]
 

NTuple::Item<float> MdcDedxRecon::m_expect [private]
 

float MdcDedxRecon::m_expect1
 

NTuple::Array<float> MdcDedxRecon::m_expectid [private]
 

NTuple::Array<float> MdcDedxRecon::m_expectid [private]
 

NTuple::Item<double> MdcDedxRecon::m_exraw [private]
 

NTuple::Item<double> MdcDedxRecon::m_exraw [private]
 

NTuple::Item<float> MdcDedxRecon::m_layer [private]
 

NTuple::Item<float> MdcDedxRecon::m_layer [private]
 

NTuple::Item<float> MdcDedxRecon::m_localwid [private]
 

NTuple::Item<float> MdcDedxRecon::m_localwid [private]
 

NTuple::Item<float> MdcDedxRecon::m_ndedxhit [private]
 

NTuple::Item<float> MdcDedxRecon::m_ndedxhit [private]
 

unsigned MdcDedxRecon::m_nevt [private]
 

NTuple::Item<float> MdcDedxRecon::m_nhit [private]
 

NTuple::Item<float> MdcDedxRecon::m_nhit [private]
 

int MdcDedxRecon::m_nhit1
 

NTuple::Item<float> MdcDedxRecon::m_nhit_hit [private]
 

NTuple::Item<float> MdcDedxRecon::m_nhit_hit [private]
 

NTuple::Tuple* MdcDedxRecon::m_nt1 [private]
 

NTuple::Tuple* MdcDedxRecon::m_nt1 [private]
 

NTuple::Tuple* MdcDedxRecon::m_nt2 [private]
 

NTuple::Tuple* MdcDedxRecon::m_nt2 [private]
 

NTuple::Item<float> MdcDedxRecon::m_parttype [private]
 

NTuple::Item<float> MdcDedxRecon::m_parttype [private]
 

int MdcDedxRecon::m_parttype1
 

NTuple::Item<float> MdcDedxRecon::m_pathL [private]
 

NTuple::Item<float> MdcDedxRecon::m_pathL [private]
 

NTuple::Item<float> MdcDedxRecon::m_pch [private]
 

NTuple::Item<float> MdcDedxRecon::m_pch [private]
 

float MdcDedxRecon::m_pch1
 

NTuple::Item<double> MdcDedxRecon::m_phraw [private]
 

NTuple::Item<double> MdcDedxRecon::m_phraw [private]
 

NTuple::Item<float> MdcDedxRecon::m_phtm [private]
 

NTuple::Item<float> MdcDedxRecon::m_phtm [private]
 

float MdcDedxRecon::m_phtm1
 

NTuple::Array<float> MdcDedxRecon::m_probpid [private]
 

NTuple::Array<float> MdcDedxRecon::m_probpid [private]
 

float MdcDedxRecon::m_probpid1
 

NTuple::Item<float> MdcDedxRecon::m_ptrk1 [private]
 

NTuple::Item<float> MdcDedxRecon::m_ptrk1 [private]
 

float MdcDedxRecon::m_rec_sw [private]
 

NTuple::Item<float> MdcDedxRecon::m_runNO [private]
 

NTuple::Item<float> MdcDedxRecon::m_runNO [private]
 

NTuple::Item<double> MdcDedxRecon::m_runNO1 [private]
 

NTuple::Item<double> MdcDedxRecon::m_runNO1 [private]
 

NTuple::Item<float> MdcDedxRecon::m_sigma [private]
 

NTuple::Item<float> MdcDedxRecon::m_sigma [private]
 

float MdcDedxRecon::m_sigma1
 

NTuple::Array<float> MdcDedxRecon::m_sigmaid [private]
 

NTuple::Array<float> MdcDedxRecon::m_sigmaid [private]
 

NTuple::Item<float> MdcDedxRecon::m_sinenta [private]
 

NTuple::Item<float> MdcDedxRecon::m_sinenta [private]
 

NTuple::Item<float> MdcDedxRecon::m_sintheta [private]
 

NTuple::Item<float> MdcDedxRecon::m_sintheta [private]
 

float MdcDedxRecon::m_sintheta1
 

NTuple::Item<float> MdcDedxRecon::m_t0 [private]
 

NTuple::Item<float> MdcDedxRecon::m_t0 [private]
 

NTuple::Item<float> MdcDedxRecon::m_t01 [private]
 

NTuple::Item<float> MdcDedxRecon::m_t01 [private]
 

NTuple::Item<float> MdcDedxRecon::m_tdc_raw [private]
 

NTuple::Item<float> MdcDedxRecon::m_tdc_raw [private]
 

NTuple::Item<float> MdcDedxRecon::m_trackId [private]
 

NTuple::Item<float> MdcDedxRecon::m_trackId [private]
 

NTuple::Item<float> MdcDedxRecon::m_trackId1 [private]
 

NTuple::Item<float> MdcDedxRecon::m_trackId1 [private]
 

int MdcDedxRecon::m_trkalg
 

vector<int> MdcDedxRecon::m_trkalgs
 

vector<int> MdcDedxRecon::m_trkalgs
 

NTuple::Item<float> MdcDedxRecon::m_wire [private]
 

NTuple::Item<float> MdcDedxRecon::m_wire [private]
 

float MdcDedxRecon::mincut_dedx
 

int MdcDedxRecon::minhit
 

int MdcDedxRecon::ntpFlag
 

int MdcDedxRecon::ParticleFlag
 

int MdcDedxRecon::recalg
 

int MdcDedxRecon::stsw
 

int MdcDedxRecon::timing_check
 

float MdcDedxRecon::truncate
 

int MdcDedxRecon::x_rej
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:24:55 2011 for BOSS6.5.5 by  doxygen 1.3.9.1