00001
00002
00003
00004
00005 #include "stdio.h"
00006 #include "math.h"
00007 #include <iostream>
00008 #include <fstream>
00009 #include <string>
00010 #include <algorithm>
00011 #include "GaudiKernel/MsgStream.h"
00012 #include "GaudiKernel/AlgFactory.h"
00013 #include "GaudiKernel/ISvcLocator.h"
00014 #include "GaudiKernel/SmartDataPtr.h"
00015 #include "GaudiKernel/IDataProviderSvc.h"
00016 #include "GaudiKernel/IDataManagerSvc.h"
00017 #include "GaudiKernel/PropertyMgr.h"
00018 #include "GaudiKernel/IHistogramSvc.h"
00019 #include "AIDA/IHistogramFactory.h"
00020 #include "GaudiKernel/INTupleSvc.h"
00021 #include "EventModel/EventHeader.h"
00022 #include "ReconEvent/ReconEvent.h"
00023 #include "EvTimeEvent/RecEsTime.h"
00024 #include "MdcRecEvent/RecMdcHit.h"
00025 #include "MdcRecEvent/RecMdcTrack.h"
00026 #include "MdcRecEvent/RecMdcKalTrack.h"
00027 #include "MdcRecEvent/RecMdcKalHelixSeg.h"
00028 #include "MdcRecEvent/RecMdcDedx.h"
00029 #include "MdcRecEvent/RecMdcDedxHit.h"
00030 #include "TTree.h"
00031 #include "Identifier/Identifier.h"
00032 #include "Identifier/MdcID.h"
00033
00034 #include "MdcDedxAlg/MdcDedxRecon.h"
00035 #include "MdcDedxAlg/MdcDedxCorrection.h"
00036 #include "MdcDedxAlg/MdcDedxTrk.h"
00037 #include "MdcDedxAlg/MdcDedxParam.h"
00038
00039
00040 using namespace std;
00041 using CLHEP::HepVector;
00042
00043 extern "C" {float prob_ (float *, int*);}
00044
00045 int RunNO1 = 0;
00046 int RunNO2;
00047
00048 int eventNo = 0;
00049 int trackNO1 = 0;
00050 int trackNO2 = 0;
00051 int trackNO3 = 0;
00052 MdcDedxRecon::MdcDedxRecon(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
00053 {
00054 ex_calib=0;
00055 calib_flag = 0x7F;
00056 landau = 1;
00057 ntpFlag = 1;
00058 doNewGlobal = 0;
00059 recalg = 0;
00060
00061 ParticleFlag = -1;
00062 m_alg = 1;
00063 truncate = 0.70;
00064
00065
00066 declareProperty("CalibFlag",calib_flag);
00067 declareProperty("NTupleFlag",ntpFlag);
00068 declareProperty("RecAlg",recalg);
00069 declareProperty("ParticleType",ParticleFlag);
00070 declareProperty("dEdxMethod",m_alg);
00071 declareProperty("TruncRate",truncate);
00072 declareProperty("RootFile",m_rootfile = std::string("no rootfile"));
00073 }
00074
00075
00076 StatusCode MdcDedxRecon::initialize()
00077 {
00078 MsgStream log(msgSvc(), name());
00079 log << MSG::INFO << "in initialize()" << endreq;
00080
00081 log << MSG::INFO <<"--------------------< MdcDedxRecon V2.1 >---------------------"<< endreq;
00082 log << MSG::INFO <<"####### correction for the wire current dependence #######"<< endreq;
00083 log << MSG::INFO <<"####### correction for the new z dependence #######"<< endreq;
00084 log << MSG::INFO <<"-------------------------------------------------------------"<< endreq;
00085 log << MSG::INFO <<"++++++++++++++++++++[ initialization ]+++++++++++++++++++++"<< endreq;
00086 log << MSG::INFO << "MdcDedxRecon has been initialized with calib_flag: "<< calib_flag <<endreq;
00087 log << MSG::INFO << "MdcDedxRecon has been initialized with landau: "<< landau << endreq;
00088 if( landau == 0 ) {truncate = 0.7;}
00089 log << MSG::INFO << "MdcDedxRecon has been initialized with ntpFlag: "<< ntpFlag << endreq;
00090 log << MSG::INFO << "MdcDedxRecon has been initialized with recalg: "<< recalg << endreq;
00091 log << MSG::INFO << "MdcDedxRecon has been initialized with dE/dx cal method m_alg: "<< m_alg << endreq;
00092 log << MSG::INFO << "MdcDedxRecon has been initialized with truncate: "<< truncate <<endreq;
00093 log << MSG::INFO << "MdcDedxRecon has been initialized with doNewGlobal: "<<doNewGlobal<<endreq;
00094 log << MSG::DEBUG << "+++++++++++ MdcDedxRecon initialization end ++++++++++++ "<< endreq;
00095 if( truncate <= 0.0 || 1.0 < truncate )
00096 {
00097 log << MSG::FATAL <<" Initialization ERROR of truncate = "<<truncate<< endreq;
00098 log << MSG::FATAL <<" truncate must be within 0.0 to 1.0 ! "<< endreq;
00099 log << MSG::FATAL <<" Please stop exec. "<< endreq;
00100 }
00101 log << MSG::DEBUG <<"-------------------------------------------------------------"<< endreq;
00102 log << MSG::DEBUG <<"MdcDedxRecon init 2nd part!!!"<< endreq;
00103
00104 StatusCode scex = service("DedxCorrecSvc", exsvc);
00105 if (scex == StatusCode::SUCCESS)
00106 {
00107 log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endreq;
00108 }
00109 else
00110 {
00111 log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endreq;
00112 return StatusCode::SUCCESS;
00113 }
00114 exsvc->set_flag( calib_flag );
00115
00116 StatusCode cursvc = service("DedxCurSvc", dedxcursvc);
00117 if (cursvc == StatusCode::SUCCESS)
00118 {
00119 log << MSG::INFO <<"DedxCurSvc is running"<<endreq;
00120 }
00121 else
00122 {
00123 log << MSG::ERROR <<"DedxCurSvc is failed"<<endreq;
00124 return StatusCode::SUCCESS;
00125 }
00126
00127
00128 if( !ex_calib ) ex_calib = new MdcDedxCorrection;
00129
00130
00131
00132
00133
00134 if(m_rootfile!="no rootfile")
00135 {
00136 const char* preDir = gDirectory->GetPath();
00137 m_hlist = new TObjArray(0);
00138 m_hitlevel = new TFolder("dedx_hitlevel","hitlevel");
00139 m_hlist -> Add(m_hitlevel);
00140 m_hitNo_wire = new TH1F("dedx_HitNo_wire","dedx hitNo vs wire",6797, -0.5, 6796.5);
00141 m_hitlevel -> Add(m_hitNo_wire);
00142 gDirectory->cd(preDir);
00143 }
00144
00145
00146 if( ntpFlag ==2 )
00147 {
00148 StatusCode status;
00149 NTuplePtr nt1(ntupleSvc(),"FILE103/n103");
00150 if ( nt1 )
00151 m_nt1 = nt1;
00152 else
00153 {
00154 m_nt1= ntupleSvc()->book("FILE103/n103",CLID_ColumnWiseTuple,"dEdx vs momentum");
00155 if ( m_nt1 )
00156 {
00157 status = m_nt1->addItem("phtm",m_phtm);
00158
00159
00160
00161
00162
00163
00164
00165 status = m_nt1->addItem("dEdx_meas", m_dEdx_meas);
00166
00167 status = m_nt1->addItem("ptrk",m_ptrk);
00168 status = m_nt1->addItem("sintheta",m_sintheta);
00169 status = m_nt1->addItem("costheta",m_costheta);
00170 status = m_nt1->addItem("charge",m_charge);
00171 status = m_nt1->addItem("runNO",m_runNO);
00172 status = m_nt1->addItem("evtNO",m_evtNO);
00173 status = m_nt1->addItem("t0",m_t0);
00174 status = m_nt1->addItem("trackId",m_trackId);
00175 status = m_nt1->addItem("poca_x",m_poca_x);
00176 status = m_nt1->addItem("poca_y",m_poca_y);
00177 status = m_nt1->addItem("poca_z",m_poca_z);
00178 status = m_nt1->addItem("recalg",m_recalg);
00179
00180 status = m_nt1->addItem("nhit",m_nhit);
00181 status = m_nt1->addItem("usedhit",m_usedhit);
00182 status = m_nt1->addItem("ndedxhit",m_nphlisthit,0,100);
00183 status = m_nt1->addIndexedItem("dEdx_hit",m_nphlisthit,m_dEdx_hit);
00184
00185 status = m_nt1->addItem("type",m_parttype);
00186 status = m_nt1->addItem("prob",m_prob);
00187 status = m_nt1->addItem("expect",m_expect);
00188 status = m_nt1->addItem("sigma",m_sigma);
00189 status = m_nt1->addItem("chidedx",m_chidedx);
00190 status = m_nt1->addItem("chidedx_e",m_chidedxe);
00191 status = m_nt1->addItem("chidedx_mu",m_chidedxmu);
00192 status = m_nt1->addItem("chidedx_pi",m_chidedxpi);
00193 status = m_nt1->addItem("chidedx_k",m_chidedxk);
00194 status = m_nt1->addItem("chidedx_p",m_chidedxp);
00195 status = m_nt1->addItem("partid",5,m_probpid);
00196 status = m_nt1->addItem("expectid",5,m_expectid);
00197 status = m_nt1->addItem("sigmaid",5,m_sigmaid);
00198 }
00199 }
00200
00201 NTuplePtr nt2(ntupleSvc(),"FILE103/n102");
00202 if ( nt2 ) m_nt2 = nt2;
00203 else
00204 {
00205 m_nt2= ntupleSvc()->book("FILE103/n102",CLID_RowWiseTuple,"pulae height raw");
00206 if ( m_nt2 )
00207 {
00208 status = m_nt2->addItem("charge",m_charge1);
00209 status = m_nt2->addItem("adc_raw",m_phraw);
00210 status = m_nt2->addItem("exraw",m_exraw);
00211 status = m_nt2->addItem("runNO1",m_runNO1);
00212 status = m_nt2->addItem("nhit_hit", m_nhit_hit);
00213 status = m_nt2->addItem("wire",m_wire);
00214
00215 status = m_nt2->addItem("doca_in",m_doca_in);
00216 status = m_nt2->addItem("doca_ex",m_doca_ex);
00217 status = m_nt2->addItem("driftdist",m_driftdist);
00218 status = m_nt2->addItem("eangle",m_eangle);
00219 status = m_nt2->addItem("costheta1",m_costheta1);
00220 status = m_nt2->addItem("path_rphi",m_pathL);
00221 status = m_nt2->addItem("layer",m_layer);
00222 status = m_nt2->addItem("ptrk1",m_ptrk1);
00223 status = m_nt2->addItem("ptrk_hit",m_ptrk_hit);
00224 status = m_nt2->addItem("t01",m_t01);
00225 status = m_nt2->addItem("tdc_raw",m_tdc_raw);
00226 status = m_nt2->addItem("driftT",m_driftT);
00227 status = m_nt2->addItem("localwid",m_localwid);
00228 status = m_nt2->addItem("trackId1",m_trackId1);
00229 }
00230 }
00231 }
00232
00233 return StatusCode::SUCCESS;
00234 }
00235
00236
00237 StatusCode MdcDedxRecon::beginRun()
00238 {
00239 MsgStream log(msgSvc(), name());
00240 log << MSG::DEBUG <<"in MdcDedxRecon::beginrun()!!!"<< endreq;
00241
00242 StatusCode gesc = service("MdcGeomSvc", geosvc);
00243 if (gesc == StatusCode::SUCCESS)
00244 {
00245 log << MSG::INFO <<"MdcGeomSvc is running"<<endreq;
00246 }
00247 else
00248 {
00249 log << MSG::ERROR <<"MdcGeomSvc is failed"<<endreq;
00250 return StatusCode::SUCCESS;
00251 }
00252
00253 return StatusCode::SUCCESS;
00254 }
00255
00256
00257 StatusCode MdcDedxRecon::finalize()
00258 {
00259 MsgStream log(msgSvc(), name());
00260 log << MSG::INFO << "in finalize()" << endreq;
00261
00262 ex_trks.clear();
00263 m_trkalgs.clear();
00264
00265 if(m_rootfile != "no rootfile")
00266 {
00267 TFile *f1 = new TFile(m_rootfile.c_str(),"recreate");
00268 m_hlist->Write();
00269 f1->Close();
00270 delete f1;
00271 delete m_hitNo_wire;
00272 delete m_hitlevel;
00273 delete m_hlist;
00274 }
00275 delete ex_calib;
00276
00277 cout<<"total event number is : "<<eventNo<<endl;
00278 cout<<"total track number is : "<<trackNO1
00279 <<" RecMdcTrack number is : "<<trackNO2
00280 <<" RecMdcKalTrack number is :"<<trackNO3<<endl;
00281 log << MSG::DEBUG <<"MdcDedxRecon terminated!!!"<< endreq;
00282 return StatusCode::SUCCESS;
00283 }
00284
00285
00286 StatusCode MdcDedxRecon::execute()
00287 {
00288 MsgStream log(msgSvc(), name());
00289 log << MSG::INFO << "in execute()" << endreq;
00290
00291 vector<double> Curve_Para, Sigma_Para;
00292 int vFlag[3];
00293
00294 for(int i=0; i< dedxcursvc->getCurveSize(); i++)
00295 {
00296 if(i==0) vFlag[0] = (int) dedxcursvc->getCurve(i);
00297 else Curve_Para.push_back( dedxcursvc->getCurve(i) );
00298 }
00299 for(int i=0; i< dedxcursvc->getSigmaSize(); i++)
00300 {
00301 if(i==0) vFlag[1] = (int) dedxcursvc->getSigma(i);
00302 else Sigma_Para.push_back( dedxcursvc->getSigma(i) );
00303 }
00304
00305
00306 if(vFlag[0] ==1 && Curve_Para.size() != 5)
00307 cout<<" size of dedx curve paramters for version 1 is not right!"<<endl;
00308 if(vFlag[0] ==2 && Curve_Para.size() != 16)
00309 cout<<" size of dedx curve paramters for version 2 is not right!"<<endl;
00310
00311 if(vFlag[1] ==1 && Sigma_Para.size() != 14)
00312 cout<<" size of dedx sigma paramters for version 1 is not right!"<<endl;
00313 if(vFlag[1] ==2 && Sigma_Para.size() != 21)
00314 cout<<" size of dedx sigma paramters for version 2 is not right!"<<endl;
00315 if(vFlag[1] ==3 && Sigma_Para.size() != 18)
00316 cout<<" size of dedx sigma paramters for version 3 is not right!"<<endl;
00317 if(vFlag[1] ==4 && Sigma_Para.size() != 19)
00318 cout<<" size of dedx sigma paramters for version 4 is not right!"<<endl;
00319 if(vFlag[1] ==5 && Sigma_Para.size() != 22)
00320 cout<<" size of dedx sigma paramters for version 5 is not right!"<<endl;
00321
00322
00323
00324 DataObject *aReconEvent;
00325 eventSvc()->findObject("/Event/Recon",aReconEvent);
00326 if(aReconEvent==NULL)
00327 {
00328 aReconEvent = new ReconEvent();
00329 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00330 if(sc!=StatusCode::SUCCESS)
00331 {
00332 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00333 return( StatusCode::FAILURE);
00334 }
00335 }
00336
00337 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00338
00339 DataObject *aDedxcol;
00340 eventSvc()->findObject("/Event/Recon/RecMdcDedxCol",aDedxcol);
00341 if(aDedxcol != NULL)
00342 {
00343 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxCol");
00344 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxCol");
00345 }
00346 dedxList = new RecMdcDedxCol;
00347 StatusCode dedxsc;
00348 dedxsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxCol, dedxList);
00349 if(!dedxsc.isSuccess())
00350 {
00351 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
00352 return ( StatusCode::FAILURE);
00353 }
00354
00355 DataObject *aDedxhitcol;
00356 eventSvc()->findObject("/Event/Recon/RecMdcDedxHitCol",aDedxhitcol);
00357 if(aDedxhitcol != NULL)
00358 {
00359 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxHitCol");
00360 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxHitCol");
00361 }
00362 dedxhitList = new RecMdcDedxHitCol;
00363 StatusCode dedxhitsc;
00364 dedxhitsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxHitCol, dedxhitList);
00365 if(!dedxhitsc.isSuccess())
00366 {
00367 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
00368 return ( StatusCode::FAILURE);
00369 }
00370
00371 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00372 if (!eventHeader)
00373 {
00374 log << MSG::INFO << "Could not find Event Header" << endreq;
00375 return( StatusCode::FAILURE);
00376 }
00377 int eventNO = eventHeader->eventNumber();
00378 int runNO = eventHeader->runNumber();
00379 log << MSG::INFO << "now begin the event: runNO= "<<runNO<<" eventNO= "<< eventNO<< endreq;
00380 RunNO2= runNO;
00381 if(RunNO1==0) RunNO1=runNO;
00382 if(RunNO2 != RunNO1)
00383 {
00384 cout<<"RunNO2 = "<<RunNO2 <<" RunNO1 = "<<RunNO1<<endl;
00385 }
00386 RunNO1 = runNO;
00387 int runFlag;
00388 if(runNO<MdcDedxParam::RUN0) runFlag =0;
00389 else if(runNO>=MdcDedxParam::RUN1 && runNO<MdcDedxParam::RUN2) runFlag =1;
00390 else if(runNO>=MdcDedxParam::RUN2 && runNO<MdcDedxParam::RUN3) runFlag =2;
00391 else if(runNO>=MdcDedxParam::RUN3 && runNO<MdcDedxParam::RUN4) runFlag =3;
00392 else runFlag =4;
00393
00394
00395 if(runNO < 0) vFlag[2]=1;
00396 else vFlag[2]=0;
00397
00398 double tes = -999.0;
00399 int esTimeflag = -1;
00400 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00401 if( ! estimeCol)
00402 {
00403 log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
00404 }
00405 else
00406 {
00407 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
00408 for(; iter_evt!=estimeCol->end(); iter_evt++)
00409 {
00410 tes = (*iter_evt)->getTest();
00411 esTimeflag = (*iter_evt)->getStat();
00412 }
00413 }
00414
00415
00416
00417 Identifier mdcid;
00418 int ntrk;
00419 ex_trks.clear();
00420 m_trkalgs.clear();
00421 if( !doNewGlobal )
00422 {
00423 log << MSG::DEBUG << "recalg: "<<recalg<<endreq;
00424
00425
00426 if(recalg ==2 )
00427 {
00428
00429 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00430 if (!newtrkCol)
00431 {
00432 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
00433 return StatusCode::SUCCESS;
00434 }
00435 if(ntpFlag>0) eventNo++;
00436 ntrk = newtrkCol->size();
00437 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
00438
00439 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
00440
00441 for( ; trk != newtrkCol->end(); trk++)
00442 {
00443 if(ntpFlag>0) trackNO1++;
00444
00445 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
00446
00447 if(gothits.size()==0)
00448 {
00449 m_trkalg=0;
00450 int trkid=(*trk)->trackId();
00451 switchtomdctrack(trkid, mdcid, tes, runNO, runFlag, log);
00452 }
00453 else
00454 {
00455 m_trkalg=1;
00456 if(gothits.size()<2) continue;
00457 kaltrackrec(trk, mdcid, tes, runNO, runFlag, log );
00458 }
00459 }
00460 }
00461
00462
00463 else if(recalg ==0 )
00464 {
00465 m_trkalg=0;
00466
00467
00468 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00469 if (!newtrkCol)
00470 {
00471 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
00472 return StatusCode::SUCCESS;
00473 }
00474 if(ntpFlag>0) eventNo++;
00475 ntrk = newtrkCol->size();
00476 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
00477
00478 vector<double> phlist;
00479 vector<double> phlist_hit;
00480 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
00481 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
00482 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
00483 int Nhits=0;
00484 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
00485
00486 RecMdcTrackCol::iterator trk = newtrkCol->begin();
00487 for( ; trk != newtrkCol->end(); trk++)
00488 {
00489 if(ntpFlag>0) trackNO1++;
00490
00491 MdcDedxTrk trk_ex( **trk);
00492 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
00493 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
00494
00495 HepVector a = (*trk)->helix();
00496 HepSymMatrix tkErrM = (*trk)->err();
00497 m_d0E = tkErrM[0][0];
00498 m_phi0E = tkErrM[1][1];
00499 m_cpaE = tkErrM[2][2];
00500 m_z0E = tkErrM[3][3];
00501
00502 HepPoint3D IP(0,0,0);
00503 Dedx_Helix exhel(IP, a);
00504 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
00505
00506 phi0= a[1];
00507 costheta = cos(M_PI_2-atan(a[4]));
00508
00509 sintheta = sin(M_PI_2-atan(a[4]));
00510 m_Pt = 1.0/fabs( a[2] );
00511 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
00512 charge = ( a[2] > 0 )? 1 : -1;
00513
00514 dedxhitrefvec = new DedxHitRefVec;
00515
00516
00517 HitRefVec gothits= (*trk)->getVecHits();
00518 Nhits = gothits.size();
00519 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
00520 if(gothits.size()<2)
00521 {
00522 delete dedxhitrefvec;
00523 continue;
00524 }
00525
00526
00527 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
00528 HitRefVec::iterator it_gothit = gothits.begin();
00529 for( ; it_gothit != gothits.end(); it_gothit++)
00530 {
00531 mdcid = (*it_gothit)->getMdcId();
00532 layid = MdcID::layer(mdcid);
00533 localwid = MdcID::wire(mdcid);
00534 w0id = geosvc->Layer(layid)->Wirst();
00535 wid = w0id + localwid;
00536 adc_raw = (*it_gothit)->getAdc();
00537 tdc_raw = (*it_gothit)->getTdc();
00538 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
00539 zhit = (*it_gothit)->getZhit();
00540 lr = (*it_gothit)->getFlagLR();
00541 if(lr == 2) cout<<"lr = "<<lr<<endl;
00542 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
00543 else driftD = (*it_gothit)->getDriftDistRight();
00544
00545 driftd = abs(driftD);
00546 dd = (*it_gothit)->getDoca();
00547 if(lr == 0 || lr == 2 ) dd = -abs(dd);
00548 if(lr == 1 ) dd = abs(dd);
00549 driftT = (*it_gothit)->getDriftT();
00550
00551 eangle = (*it_gothit)->getEntra();
00552 eangle = eangle/M_PI;
00553 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
00554 rphi_path=pathlength * sintheta;
00555
00556
00557 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
00558 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, ph, costheta,tes );
00559
00560
00561
00562 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
00563 if(runNO<0)
00564 {
00565 if (layid<8)
00566 {
00567 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00568 }
00569 else if(layid >= 8)
00570 {
00571 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00572 }
00573 }
00574 else if(runNO>=0)
00575 {
00576 if(layid <8)
00577 {
00578 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00579 }
00580 else if(layid >= 8)
00581 {
00582 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00583 }
00584 }
00585
00586
00587 if (ph<0.0||ph==0) continue;
00588 else
00589 {
00590
00591 dedxhit = new RecMdcDedxHit;
00592 dedxhit->setMdcHit(*it_gothit);
00593 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg );
00594 dedxhit->setMdcId(mdcid);
00595 dedxhit->setFlagLR(lr);
00596 dedxhit->setTrkId(trk_ex.trk_id());
00597 dedxhit->setDedx(ph_hit);
00598 dedxhit->setPathLength(pathlength);
00599
00600
00601 if(m_rootfile!= "no rootfile")
00602 {
00603 m_hitNo_wire->Fill(wid);
00604 }
00605
00606
00607 if ( ntpFlag ==2 )
00608 {
00609 m_charge1 = (*trk)->charge();
00610 m_t01 = tes;
00611 m_driftT = driftT;
00612 m_tdc_raw = tdc_raw;
00613 m_phraw = adc_raw;
00614 m_exraw = ph_hit;
00615 m_localwid = localwid;
00616 m_wire = wid;
00617 m_runNO1 = runNO;
00618 m_nhit_hit = Nhits;
00619 m_doca_in = dd;
00620 m_doca_ex = dd;
00621 m_driftdist = driftD;
00622 m_eangle = eangle;
00623 m_costheta1 = costheta;
00624 m_pathL = pathlength;
00625 m_layer = layid;
00626 m_ptrk1 = m_P;
00627
00628 m_trackId1 = trk_ex.trk_id();
00629 StatusCode status2 = m_nt2->write();
00630 if ( status2.isFailure() )
00631 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
00632 }
00633 if(layid>3)
00634 {
00635 phlist.push_back(ph);
00636 phlist_hit.push_back(ph_hit);
00637 }
00638 }
00639 dedxhitList->push_back( dedxhit );
00640 dedxhitrefvec->push_back( dedxhit );
00641 }
00642 trk_ex.set_phlist( phlist );
00643 trk_ex.set_phlist_hit( phlist_hit );
00644 trk_ex.setVecDedxHits( *dedxhitrefvec );
00645 ex_trks.push_back(trk_ex );
00646 m_trkalgs.push_back(m_trkalg);
00647
00648 delete dedxhitrefvec;
00649 phlist.clear();
00650 phlist_hit.clear();
00651 if(ntpFlag>0) trackNO2++;
00652 }
00653 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
00654 }
00655
00656
00657 else if(recalg ==1 )
00658 {
00659 m_trkalg=1;
00660
00661
00662 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00663 if (!newtrkCol)
00664 {
00665 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
00666 return StatusCode::SUCCESS;
00667 }
00668 if(ntpFlag>0) eventNo++;
00669 ntrk = newtrkCol->size();
00670 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
00671
00672 vector<double> phlist;
00673 vector<double> phlist_hit;
00674 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
00675 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
00676 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
00677 int Nhits=0;
00678 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
00679
00680
00681 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
00682 for( ; trk != newtrkCol->end(); trk++)
00683 {
00684 if(ntpFlag>0) trackNO1++;
00685
00686 MdcDedxTrk trk_ex( **trk, ParticleFlag);
00687 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
00688 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
00689
00690 HepVector a;
00691 HepSymMatrix tkErrM;
00692 if(ParticleFlag>-1&&ParticleFlag<5)
00693 {
00694 DstMdcKalTrack* dstTrk = *trk;
00695 a = dstTrk->getFHelix(ParticleFlag);
00696 tkErrM = dstTrk->getFError(ParticleFlag);
00697 }
00698 else
00699 {
00700 a = (*trk)->getFHelix();
00701 tkErrM = (*trk)->getFError();
00702 }
00703
00704 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
00705
00706 m_d0E = tkErrM[0][0];
00707 m_phi0E = tkErrM[1][1];
00708 m_cpaE = tkErrM[2][2];
00709 m_z0E = tkErrM[3][3];
00710
00711 phi0= a[1];
00712 costheta = cos(M_PI_2-atan(a[4]));
00713
00714 sintheta = sin(M_PI_2-atan(a[4]));
00715 m_Pt = 1.0/fabs( a[2] );
00716 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
00717
00718 charge = ( a[2] > 0 )? 1 : -1;
00719
00720 dedxhitrefvec = new DedxHitRefVec;
00721
00722
00723 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
00724
00725
00726 Nhits = gothits.size();
00727 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
00728 if(gothits.size()<2)
00729 {
00730 delete dedxhitrefvec;
00731 continue;
00732 }
00733
00734
00735 RecMdcHit* mdcHit = 0;
00736 HelixSegRefVec::iterator it_gothit = gothits.begin();
00737 for( ; it_gothit != gothits.end(); it_gothit++)
00738 {
00739 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
00740
00741 HepPoint3D IP(0,0,0);
00742 Dedx_Helix exhel_hit_in(IP, a_hit_in);
00743
00744 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
00745
00746
00747 mdcid = (*it_gothit)->getMdcId();
00748 layid = MdcID::layer(mdcid);
00749 localwid = MdcID::wire(mdcid);
00750 w0id = geosvc->Layer(layid)->Wirst();
00751 wid = w0id + localwid;
00752 adc_raw = (*it_gothit)->getAdc();
00753 tdc_raw = (*it_gothit)->getTdc();
00754 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
00755 zhit = (*it_gothit)->getZhit();
00756 lr = (*it_gothit)->getFlagLR();
00757 if(lr == 2) cout<<"lr = "<<lr<<endl;
00758 driftD = (*it_gothit)->getDD();
00759 driftd = abs(driftD);
00760 driftT = (*it_gothit)->getDT();
00761 dd_in = (*it_gothit)->getDocaIncl();
00762 dd_ex = (*it_gothit)->getDocaExcl();
00763
00764 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
00765 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
00766
00767
00768
00769 eangle = (*it_gothit)->getEntra();
00770 eangle = eangle/M_PI;
00771 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
00772 rphi_path=pathlength * sintheta;
00773
00774 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
00775 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, ph, costheta,tes );
00776
00777
00778
00779 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
00780 if(runNO<0)
00781 {
00782 if (layid<8)
00783 {
00784 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00785 }
00786 else if(layid >= 8)
00787 {
00788 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00789 }
00790 }
00791 else if(runNO>=0)
00792 {
00793 if(layid <8)
00794 {
00795 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
00796 }
00797 else if(layid >= 8)
00798 {
00799 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
00800 }
00801 }
00802
00803
00804 if (ph<0.0||ph==0) continue;
00805 else
00806 {
00807
00808 dedxhit = new RecMdcDedxHit;
00809 dedxhit->setMdcHit(mdcHit);
00810 dedxhit->setMdcKalHelixSeg(*it_gothit);
00811 dedxhit->setMdcId(mdcid);
00812 dedxhit->setFlagLR(lr);
00813 dedxhit->setTrkId(trk_ex.trk_id());
00814 dedxhit->setDedx(ph_hit);
00815 dedxhit->setPathLength(pathlength);
00816
00817
00818 if(m_rootfile!= "no rootfile")
00819 {
00820 m_hitNo_wire->Fill(wid);
00821 }
00822
00823
00824 if ( ntpFlag ==2 )
00825 {
00826 m_charge1 = (*trk)->charge();
00827 m_t01 = tes;
00828 m_driftT = driftT;
00829 m_tdc_raw = tdc_raw;
00830 m_phraw = adc_raw;
00831 m_exraw = ph_hit;
00832 m_localwid = localwid;
00833 m_wire = wid;
00834 m_runNO1 = runNO;
00835 m_nhit_hit = Nhits;
00836 m_doca_in = dd_in;
00837 m_doca_ex = dd_ex;
00838 m_driftdist = driftD;
00839 m_eangle = eangle;
00840 m_costheta1 = costheta;
00841 m_pathL = pathlength;
00842 m_layer = layid;
00843 m_ptrk1 = m_P;
00844 m_ptrk_hit = p_hit;
00845
00846 m_trackId1 = trk_ex.trk_id();
00847 StatusCode status2 = m_nt2->write();
00848 if ( status2.isFailure() )
00849 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
00850 }
00851 if(layid>3)
00852 {
00853 phlist.push_back(ph);
00854 phlist_hit.push_back(ph_hit);
00855 }
00856 }
00857 dedxhitList->push_back( dedxhit );
00858 dedxhitrefvec->push_back( dedxhit );
00859 }
00860 trk_ex.set_phlist( phlist );
00861 trk_ex.set_phlist_hit( phlist_hit );
00862 trk_ex.setVecDedxHits( *dedxhitrefvec );
00863 ex_trks.push_back(trk_ex );
00864 m_trkalgs.push_back(m_trkalg);
00865
00866 delete dedxhitrefvec;
00867 phlist.clear();
00868 phlist_hit.clear();
00869 if(ntpFlag>0) trackNO3++;
00870 }
00871 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
00872 }
00873
00874
00875
00876 if( ntrk != ex_trks.size())
00877 log << MSG::DEBUG <<"ntrkCol size: "<<ntrk<<" dedxtrk size: "<<ex_trks.size()<< endreq;
00878
00879 double poca_x=0,poca_y=0,poca_z=0;
00880 float sintheta=0,costheta=0,ptrk=0,charge=0;
00881 int Nhit=0,Nphlisthit=0;
00882 int usedhit = 0;
00883 double phtm=0,median=0,geometric=0,geometric_trunc=0,harmonic=0,harmonic_trunc=0,transform=0,logtrunc=0;
00884
00885 enum pid_dedx parType_temp(electron);
00886 float probpid_temp=-0.01,expect_temp=-0.01,sigma_temp=-0.01,chidedx_temp=-0.01;
00887
00888 double dEdx_meas_hit=0;
00889 double dEdx_meas=0,dEdx_meas_esat=0,dEdx_meas_norun=0;
00890 int trk_recalg=-1;
00891
00892 int idedxid = 0;
00893 for(std::vector<MdcDedxTrk>::iterator it = ex_trks.begin(); it != ex_trks.end(); it++, idedxid++ )
00894 {
00895 log << MSG::DEBUG <<"-------------------------------------------------------"<< endreq;
00896 log << MSG::DEBUG <<" trk id ="<< it->trk_id()<<" : P ="<<it->P() <<" Pt ="<<it->Pt()<<" : theta ="<<it->theta() <<" : phi ="<<it->phi()<< " : charge = "<<it->charge()<<endreq;
00897 log << MSG::DEBUG <<"all hits on track: "<<it->nsample()<<" phlist size: "<<it->get_phlist().size()<<endreq;
00898
00899 if(it->trk_ptr_kal()!=0)
00900 {
00901 poca_x = it->trk_ptr_kal()->x();
00902 poca_y = it->trk_ptr_kal()->y();
00903 poca_z = it->trk_ptr_kal()->z();
00904 }
00905 else if(it->trk_ptr()!=0)
00906 {
00907 poca_x = it->trk_ptr()->x();
00908 poca_y = it->trk_ptr()->y();
00909 poca_z = it->trk_ptr()->z();
00910 }
00911
00912
00913 sintheta = sin(it->theta());
00914 costheta = cos(it->theta());
00915 ptrk = it->P();
00916 charge = it->charge();
00917 Nhit = it->nsample();
00918 Nphlisthit = it->get_phlist_hit().size();
00919
00920
00921
00922 phtm = it->cal_dedx( truncate );
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932 if(m_alg==1) dEdx_meas_hit = it->cal_dedx_bitrunc(truncate, 0, usedhit);
00933 else if(m_alg==2) dEdx_meas_hit = it-> cal_dedx_transform(1);
00934 else if(m_alg==3) dEdx_meas_hit = it->cal_dedx_log(1.0, 1);
00935 else cout<<"-------------Truncate Algorithm Error!!!------------------------"<<endl;
00936 if(m_alg==1 && usedhit==0)
00937 {
00938 cout<<"***************bad extrk with no hits!!!!!******************"<<endl;
00939 continue;
00940 }
00941
00942 usedhit = (int)(Nphlisthit*truncate);
00943
00944
00945 dEdx_meas = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, dEdx_meas_hit, cos(it->theta()), tes );
00946 dEdx_meas_esat = exsvc->StandardTrackCorrec(1,runFlag, ntpFlag, runNO, dEdx_meas_hit, cos(it->theta()), tes );
00947 dEdx_meas_norun = exsvc->StandardTrackCorrec(2,runFlag, ntpFlag, runNO, dEdx_meas_hit, cos(it->theta()), tes );
00948 log << MSG::INFO << "===================== " << endreq << endreq;
00949 log << MSG::DEBUG <<"dEdx_meas_hit: "<<dEdx_meas_hit<<" dEdx_meas: "<<dEdx_meas<<" dEdx_meas_esat: "<<dEdx_meas_esat<<" dEdx_meas_norun: "<<dEdx_meas_norun<<" ptrk: "<<it->P()<<endreq;
00950 log << MSG::INFO << "ptrk " << ptrk << " costhe " << costheta << " runno " << runNO << " evtno " << eventNO << endreq;
00951
00952
00953 trk_recalg = m_trkalgs[idedxid];
00954
00955 if(dEdx_meas<0 || dEdx_meas==0) continue;
00956 it->set_dEdx( landau , dEdx_meas, trk_recalg, runFlag, vFlag , tes, Curve_Para, Sigma_Para, ex_calib);
00957 parType_temp = electron;
00958 probpid_temp=-0.01;
00959 expect_temp=-0.01;
00960 sigma_temp=-0.01;
00961 chidedx_temp=-0.01;
00962 for(int i=0;i<5;i++)
00963 {
00964 if(it->pprob()[i]>probpid_temp)
00965 {
00966 parType_temp = (enum pid_dedx) i;
00967 probpid_temp = it->pprob()[i];
00968 expect_temp = it->pexpect()[i];
00969 sigma_temp = it->pexp_sigma()[i];
00970 chidedx_temp = it->pchi_dedx()[i];
00971 }
00972 }
00973 log<< MSG::INFO<<"expect dE/dx: type: "<<parType_temp<<" probpid: "<<probpid_temp<<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq;
00974
00975
00976
00977 resdedx = new RecMdcDedx;
00978 resdedx->setDedxHit(dEdx_meas_hit);
00979 resdedx->setDedxEsat(dEdx_meas_esat);
00980 resdedx->setDedxNoRun(dEdx_meas_norun);
00981 resdedx->setDedxMoment(it->P());
00982 resdedx->setTrackId( it->trk_id() );
00983 resdedx->setProbPH( dEdx_meas );
00984 resdedx->setNormPH( dEdx_meas/550.0 );
00985 resdedx->setDedxExpect( it->pexpect() );
00986 resdedx->setSigmaDedx( it->pexp_sigma() );
00987 resdedx->setPidProb( it->pprob() );
00988 resdedx->setChi( it->pchi_dedx() );
00989
00990 resdedx->setNumTotalHits(it->nsample() );
00991 resdedx->setNumGoodHits(usedhit);
00992
00993
00994 resdedx->setStatus(trk_recalg );
00995
00996 resdedx->setTruncAlg( m_alg );
00997 resdedx->setParticleId(parType_temp);
00998
00999 resdedx->setVecDedxHits(it->getVecDedxHits());
01000 resdedx->setMdcTrack(it->trk_ptr());
01001 resdedx->setMdcKalTrack(it->trk_ptr_kal());
01002
01003
01004 if(ntpFlag ==2)
01005 {
01006 m_phtm=phtm;
01007
01008
01009
01010
01011
01012
01013
01014 m_dEdx_meas = dEdx_meas;
01015
01016 m_poca_x = poca_x;
01017 m_poca_y = poca_y;
01018 m_poca_z = poca_z;
01019 m_ptrk=ptrk;
01020 m_sintheta=sintheta;
01021 m_costheta=costheta;
01022 m_charge=charge;
01023 m_runNO = runNO;
01024 m_evtNO = eventNO;
01025
01026 m_t0 = tes;
01027 m_trackId = it->trk_id();
01028 m_recalg = trk_recalg;
01029
01030 m_nhit=Nhit;
01031 m_nphlisthit=Nphlisthit;
01032 m_usedhit=usedhit;
01033 for(int i=0; i<Nphlisthit; i++) m_dEdx_hit[i] = it->get_phlist_hit()[i];
01034
01035 m_parttype = parType_temp;
01036 m_prob = probpid_temp;
01037 m_expect = expect_temp;
01038 m_sigma = sigma_temp;
01039 m_chidedx = chidedx_temp;
01040 m_chidedxe=it->pchi_dedx()[0];
01041 m_chidedxmu=it->pchi_dedx()[1];
01042 m_chidedxpi=it->pchi_dedx()[2];
01043 m_chidedxk=it->pchi_dedx()[3];
01044 m_chidedxp=it->pchi_dedx()[4];
01045 for(int i=0;i<5;i++)
01046 {
01047 m_probpid[i]= it->pprob()[i];
01048 m_expectid[i]= it->pexpect()[i];
01049 m_sigmaid[i]= it->pexp_sigma()[i];
01050 }
01051
01052 StatusCode status = m_nt1->write();
01053 if ( status.isFailure() )
01054 {
01055 log << MSG::ERROR << "Cannot fill Ntuple n103" << endreq;
01056 }
01057 }
01058
01059 log<< MSG::INFO<<"-----------------Summary of this ExTrack----------------"<<endreq
01060 <<"dEdx_mean: "<<dEdx_meas<<" type: "<<parType_temp<<" probpid: "<<probpid_temp
01061 <<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq<<endreq;
01062
01063 dedxList->push_back( resdedx );
01064 }
01065 }
01066
01067
01068
01069 SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
01070 if (!newexCol)
01071 {
01072 log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
01073 return( StatusCode::SUCCESS);
01074 }
01075 log << MSG::DEBUG << "----------------Begin to check RecMdcDedxCol-----------------"<<endreq;
01076 RecMdcDedxCol::iterator it_dedx = newexCol->begin();
01077 for( ; it_dedx != newexCol->end(); it_dedx++)
01078 {
01079 log << MSG::INFO << "retrieved MDC dE/dx:" << endreq
01080 << "dEdx Id: " << (*it_dedx)->trackId()
01081 << " part Id: " << (*it_dedx)->particleType()
01082 << " measured dEdx: " << (*it_dedx)->probPH()
01083 << " dEdx std: " << (*it_dedx)->normPH() << endreq
01084 << "hits on track: "<<(*it_dedx)->numTotalHits()
01085 << " usedhits: " << (*it_dedx)->numGoodHits() << endreq
01086 << "Electron: expect: " << (*it_dedx)->getDedxExpect(0)
01087 << " sigma: " << (*it_dedx)->getSigmaDedx(0)
01088 << " chi: " << (*it_dedx)->chi(0)
01089 << " prob: " << (*it_dedx)->getPidProb(0) << endreq
01090 << "Muon: expect: " << (*it_dedx)->getDedxExpect(1)
01091 << " sigma: " << (*it_dedx)->getSigmaDedx(1)
01092 << " chi: " << (*it_dedx)->chi(1)
01093 << " prob: " << (*it_dedx)->getPidProb(1) << endreq
01094 << "Pion: expect: " << (*it_dedx)->getDedxExpect(2)
01095 << " sigma: " << (*it_dedx)->getSigmaDedx(2)
01096 << " chi: " << (*it_dedx)->chi(2)
01097 << " prob: " << (*it_dedx)->getPidProb(2) << endreq
01098 << "Kaon: expect: " << (*it_dedx)->getDedxExpect(3)
01099 << " sigma: " << (*it_dedx)->getSigmaDedx(3)
01100 << " chi: " << (*it_dedx)->chi(3)
01101 << " prob: " << (*it_dedx)->getPidProb(3) << endreq
01102 << "Proton: expect: " << (*it_dedx)->getDedxExpect(4)
01103 << " sigma: " << (*it_dedx)->getSigmaDedx(4)
01104 << " chi: " << (*it_dedx)->chi(4)
01105 << " prob: " << (*it_dedx)->getPidProb(4) << endreq;
01106 }
01107 return StatusCode::SUCCESS;
01108 }
01109
01110 const std::vector<MdcDedxTrk>&MdcDedxRecon::tracks(void) const
01111 {
01112 return ex_trks;
01113 }
01114
01115 void MdcDedxRecon::clearTables() {}
01116
01117 void MdcDedxRecon::mdctrackrec(RecMdcTrackCol::iterator trk, Identifier mdcid, double tes, int runNO,int runFlag, MsgStream log )
01118 {
01119 vector<double> phlist;
01120 vector<double> phlist_hit;
01121 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
01122 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
01123 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
01124 int Nhits=0;
01125 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
01126
01127 MdcDedxTrk trk_ex( **trk);
01128 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
01129 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
01130
01131 HepVector a = (*trk)->helix();
01132 HepSymMatrix tkErrM = (*trk)->err();
01133 m_d0E = tkErrM[0][0];
01134 m_phi0E = tkErrM[1][1];
01135 m_cpaE = tkErrM[2][2];
01136 m_z0E = tkErrM[3][3];
01137
01138 HepPoint3D IP(0,0,0);
01139 Dedx_Helix exhel(IP, a);
01140 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
01141
01142 phi0= a[1];
01143 costheta = cos(M_PI_2-atan(a[4]));
01144
01145 sintheta = sin(M_PI_2-atan(a[4]));
01146 m_Pt = 1.0/fabs( a[2] );
01147 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
01148 charge = ( a[2] > 0 )? 1 : -1;
01149
01150 dedxhitrefvec = new DedxHitRefVec;
01151
01152
01153 HitRefVec gothits= (*trk)->getVecHits();
01154 Nhits = gothits.size();
01155 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
01156
01157
01158 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
01159 HitRefVec::iterator it_gothit = gothits.begin();
01160 for( ; it_gothit != gothits.end(); it_gothit++)
01161 {
01162 mdcid = (*it_gothit)->getMdcId();
01163 layid = MdcID::layer(mdcid);
01164 localwid = MdcID::wire(mdcid);
01165 w0id = geosvc->Layer(layid)->Wirst();
01166 wid = w0id + localwid;
01167 adc_raw = (*it_gothit)->getAdc();
01168 tdc_raw = (*it_gothit)->getTdc();
01169 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
01170 zhit = (*it_gothit)->getZhit();
01171 lr = (*it_gothit)->getFlagLR();
01172 if(lr == 2) cout<<"lr = "<<lr<<endl;
01173 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
01174 else driftD = (*it_gothit)->getDriftDistRight();
01175
01176 driftd = abs(driftD);
01177 dd = (*it_gothit)->getDoca();
01178 if(lr == 0 || lr == 2 ) dd = -abs(dd);
01179 if(lr == 1 ) dd = abs(dd);
01180 driftT = (*it_gothit)->getDriftT();
01181
01182 eangle = (*it_gothit)->getEntra();
01183 eangle = eangle/M_PI;
01184 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
01185 rphi_path=pathlength * sintheta;
01186
01187
01188 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
01189 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, ph, costheta,tes );
01190
01191
01192
01193 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
01194 if(runNO<0)
01195 {
01196 if (layid<8)
01197 {
01198 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01199 }
01200 else if(layid >= 8)
01201 {
01202 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01203 }
01204 }
01205 else if(runNO>=0)
01206 {
01207 if(layid <8)
01208 {
01209 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01210 }
01211 else if(layid >= 8)
01212 {
01213 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01214 }
01215 }
01216
01217
01218 if (ph<0.0||ph==0) continue;
01219 else
01220 {
01221
01222 dedxhit = new RecMdcDedxHit;
01223 dedxhit->setMdcHit(*it_gothit);
01224 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg);
01225 dedxhit->setMdcId(mdcid);
01226 dedxhit->setFlagLR(lr);
01227 dedxhit->setTrkId(trk_ex.trk_id());
01228 dedxhit->setDedx(ph_hit);
01229 dedxhit->setPathLength(pathlength);
01230
01231
01232 if(m_rootfile!= "no rootfile")
01233 {
01234 m_hitNo_wire->Fill(wid);
01235 }
01236
01237
01238 if ( ntpFlag ==2 )
01239 {
01240 m_charge1 = (*trk)->charge();
01241 m_t01 = tes;
01242 m_driftT = driftT;
01243 m_tdc_raw = tdc_raw;
01244 m_phraw = adc_raw;
01245 m_exraw = ph_hit;
01246 m_localwid = localwid;
01247 m_wire = wid;
01248 m_runNO1 = runNO;
01249 m_nhit_hit = Nhits;
01250 m_doca_in = dd;
01251 m_doca_ex = dd;
01252 m_driftdist = driftD;
01253 m_eangle = eangle;
01254 m_costheta1 = costheta;
01255 m_pathL = pathlength;
01256 m_layer = layid;
01257 m_ptrk1 = m_P;
01258
01259 m_trackId1 = trk_ex.trk_id();
01260 StatusCode status2 = m_nt2->write();
01261 if ( status2.isFailure() )
01262 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
01263 }
01264 if(layid>3)
01265 {
01266 phlist.push_back(ph);
01267 phlist_hit.push_back(ph_hit);
01268 }
01269 }
01270 dedxhitList->push_back( dedxhit );
01271 dedxhitrefvec->push_back( dedxhit );
01272 }
01273 trk_ex.set_phlist( phlist );
01274 trk_ex.set_phlist_hit( phlist_hit );
01275 trk_ex.setVecDedxHits( *dedxhitrefvec );
01276 ex_trks.push_back(trk_ex );
01277 m_trkalgs.push_back(m_trkalg);
01278
01279 delete dedxhitrefvec;
01280 phlist.clear();
01281 phlist_hit.clear();
01282 if(ntpFlag>0) trackNO2++;
01283 }
01284
01285
01286 void MdcDedxRecon::kaltrackrec(RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int runNO,int runFlag, MsgStream log )
01287 {
01288 vector<double> phlist;
01289 vector<double> phlist_hit;
01290 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
01291 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
01292 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
01293 int Nhits=0;
01294 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
01295
01296 MdcDedxTrk trk_ex( **trk, ParticleFlag);
01297 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
01298 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
01299
01300 HepVector a;
01301 HepSymMatrix tkErrM;
01302 if(ParticleFlag>-1&&ParticleFlag<5)
01303 {
01304 DstMdcKalTrack* dstTrk = *trk;
01305 a = dstTrk->getFHelix(ParticleFlag);
01306 tkErrM = dstTrk->getFError(ParticleFlag);
01307 }
01308 else
01309 {
01310 a = (*trk)->getFHelix();
01311 tkErrM = (*trk)->getFError();
01312 }
01313 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
01314
01315
01316 m_d0E = tkErrM[0][0];
01317 m_phi0E = tkErrM[1][1];
01318 m_cpaE = tkErrM[2][2];
01319 m_z0E = tkErrM[3][3];
01320
01321 phi0= a[1];
01322 costheta = cos(M_PI_2-atan(a[4]));
01323
01324 sintheta = sin(M_PI_2-atan(a[4]));
01325 m_Pt = 1.0/fabs( a[2] );
01326 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
01327
01328 charge = ( a[2] > 0 )? 1 : -1;
01329
01330 dedxhitrefvec = new DedxHitRefVec;
01331
01332
01333 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
01334
01335
01336 Nhits = gothits.size();
01337 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
01338
01339
01340 RecMdcHit* mdcHit = 0;
01341 HelixSegRefVec::iterator it_gothit = gothits.begin();
01342 for( ; it_gothit != gothits.end(); it_gothit++)
01343 {
01344 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
01345
01346 HepPoint3D IP(0,0,0);
01347 Dedx_Helix exhel_hit_in(IP, a_hit_in);
01348
01349 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
01350
01351
01352 mdcid = (*it_gothit)->getMdcId();
01353 layid = MdcID::layer(mdcid);
01354 localwid = MdcID::wire(mdcid);
01355 w0id = geosvc->Layer(layid)->Wirst();
01356 wid = w0id + localwid;
01357 adc_raw = (*it_gothit)->getAdc();
01358 tdc_raw = (*it_gothit)->getTdc();
01359 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
01360 zhit = (*it_gothit)->getZhit();
01361 lr = (*it_gothit)->getFlagLR();
01362 if(lr == 2) cout<<"lr = "<<lr<<endl;
01363 driftD = (*it_gothit)->getDD();
01364 driftd = abs(driftD);
01365 driftT = (*it_gothit)->getDT();
01366 dd_in = (*it_gothit)->getDocaIncl();
01367 dd_ex = (*it_gothit)->getDocaExcl();
01368
01369 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
01370 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
01371
01372
01373
01374 eangle = (*it_gothit)->getEntra();
01375 eangle = eangle/M_PI;
01376 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
01377 rphi_path=pathlength * sintheta;
01378
01379
01380 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
01381
01382 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, ph, costheta,tes );
01383
01384
01385
01386
01387
01388 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
01389 if(runNO<0)
01390 {
01391 if (layid<8)
01392 {
01393 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01394 }
01395 else if(layid >= 8)
01396 {
01397 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01398 }
01399 }
01400 else if(runNO>=0)
01401 {
01402 if(layid <8)
01403 {
01404 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
01405 }
01406 else if(layid >= 8)
01407 {
01408 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
01409 }
01410 }
01411
01412
01413 if (ph<0.0||ph==0) continue;
01414 else
01415 {
01416
01417 dedxhit = new RecMdcDedxHit;
01418 dedxhit->setMdcHit(mdcHit);
01419 dedxhit->setMdcKalHelixSeg(*it_gothit);
01420 dedxhit->setMdcId(mdcid);
01421 dedxhit->setFlagLR(lr);
01422 dedxhit->setTrkId(trk_ex.trk_id());
01423 dedxhit->setDedx(ph_hit);
01424 dedxhit->setPathLength(pathlength);
01425
01426
01427 if(m_rootfile!= "no rootfile")
01428 {
01429 m_hitNo_wire->Fill(wid);
01430 }
01431
01432
01433 if ( ntpFlag ==2 )
01434 {
01435 m_charge1 = (*trk)->charge();
01436 m_t01 = tes;
01437 m_driftT = driftT;
01438 m_tdc_raw = tdc_raw;
01439 m_phraw = adc_raw;
01440 m_exraw = ph_hit;
01441 m_localwid = localwid;
01442 m_wire = wid;
01443 m_runNO1 = runNO;
01444 m_nhit_hit = Nhits;
01445 m_doca_in = dd_in;
01446 m_doca_ex = dd_ex;
01447 m_driftdist = driftD;
01448 m_eangle = eangle;
01449 m_costheta1 = costheta;
01450 m_pathL = pathlength;
01451 m_layer = layid;
01452 m_ptrk1 = m_P;
01453 m_ptrk_hit = p_hit;
01454
01455 m_trackId1 = trk_ex.trk_id();
01456 StatusCode status2 = m_nt2->write();
01457 if ( status2.isFailure() )
01458 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
01459 }
01460 if(layid>3)
01461 {
01462 phlist.push_back(ph);
01463 phlist_hit.push_back(ph_hit);
01464 }
01465 }
01466 dedxhitList->push_back( dedxhit );
01467 dedxhitrefvec->push_back( dedxhit );
01468 }
01469 trk_ex.set_phlist( phlist );
01470 trk_ex.set_phlist_hit( phlist_hit );
01471 trk_ex.setVecDedxHits( *dedxhitrefvec );
01472 ex_trks.push_back(trk_ex );
01473 m_trkalgs.push_back(m_trkalg);
01474
01475 delete dedxhitrefvec;
01476 phlist.clear();
01477 phlist_hit.clear();
01478 if(ntpFlag>0) trackNO3++;
01479 }
01480
01481 void MdcDedxRecon::switchtomdctrack(int trkid,Identifier mdcid, double tes, int runNO,int runFlag, MsgStream log)
01482 {
01483
01484 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
01485 if (!newtrkCol)
01486 {
01487 log << MSG::WARNING << "Could not find RecMdcTrackCol in switchtomdctrack" << endreq;
01488 return ;
01489 }
01490
01491 RecMdcTrackCol::iterator trk = newtrkCol->begin();
01492 for( ; trk != newtrkCol->end(); trk++)
01493 {
01494 if( (*trk)->trackId()==trkid)
01495 {
01496 HitRefVec gothits= (*trk)->getVecHits();
01497 if(gothits.size()>=2)
01498 mdctrackrec(trk,mdcid,tes,runNO,runFlag,log);
01499 }
01500 }
01501 }