BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcDedxRecon.cxx
Go to the documentation of this file.
1// BesIII MDC dE/dx Reconstruction Module
2// Class: MdcDedxRecon
3// Created by Dayong WANG 2003/11
4
5#include "stdio.h"
6#include "math.h"
7#include <iostream>
8#include <fstream>
9#include <string>
10#include <algorithm>
11#include "GaudiKernel/MsgStream.h"
12#include "GaudiKernel/AlgFactory.h"
13#include "GaudiKernel/ISvcLocator.h"
14#include "GaudiKernel/SmartDataPtr.h"
15#include "GaudiKernel/IDataProviderSvc.h"
16#include "GaudiKernel/IDataManagerSvc.h"
17#include "GaudiKernel/PropertyMgr.h"
18#include "GaudiKernel/IHistogramSvc.h"
19#include "AIDA/IHistogramFactory.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "EventModel/EventHeader.h"
22#include "ReconEvent/ReconEvent.h"
23#include "EvTimeEvent/RecEsTime.h"
24#include "MdcRecEvent/RecMdcHit.h"
25#include "MdcRecEvent/RecMdcTrack.h"
26#include "MdcRecEvent/RecMdcKalTrack.h"
27#include "MdcRecEvent/RecMdcKalHelixSeg.h"
28#include "MdcRecEvent/RecMdcDedx.h"
29#include "MdcRecEvent/RecMdcDedxHit.h"
30#include "TTree.h"
31#include "Identifier/Identifier.h"
32#include "Identifier/MdcID.h"
33
34#include "MdcDedxAlg/MdcDedxRecon.h"
35#include "MdcDedxAlg/MdcDedxCorrection.h"
36#include "MdcDedxAlg/MdcDedxTrk.h"
37#include "MdcDedxAlg/MdcDedxParam.h"
38
39
40using namespace std;
41using CLHEP::HepVector;
42
43extern "C" {float prob_ (float *, int*);}
44
45int RunNO1 = 0;
46int RunNO2;
47
48int eventNo = 0;
49int trackNO1 = 0;
50int trackNO2 = 0;
51int trackNO3 = 0;
52DECLARE_COMPONENT(MdcDedxRecon)
53MdcDedxRecon::MdcDedxRecon(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
54{
55 ex_calib=0; // get MdcDedxCorrection
56 calib_flag = 0x7F; //all calib on
57 landau = 1; //0: gauss distribution; 1:landau distribution;
58 ntpFlag = 1;
59 doNewGlobal = 0;
60 recalg = 0; //reconstruction method: 0:RecMdcTrack; 1:RecMdcKalTrack;
61 //2:use RecMdcTrack when no RecMdcKalTrack
62 ParticleFlag = -1; //e:0, mu:1, pi:2, K:3, p:4, cosmic:5
63 m_alg = 1; //calculation method of dE/dx of one track; 1: truncation and get average;
64 truncate = 0.70; // truncation rate
65
66 // Declare the properties
67 declareProperty("CalibFlag",calib_flag);
68 declareProperty("NTupleFlag",ntpFlag);
69 declareProperty("RecAlg",recalg);
70 declareProperty("ParticleType",ParticleFlag);
71 declareProperty("dEdxMethod",m_alg);
72 declareProperty("TruncRate",truncate);
73 declareProperty("RootFile",m_rootfile = std::string("no rootfile"));
74}
75
76
77StatusCode MdcDedxRecon::initialize()
78{
79 MsgStream log(msgSvc(), name());
80 log << MSG::INFO << "in initialize()" << endreq;
81
82 log << MSG::INFO <<"--------------------< MdcDedxRecon V2.1 >---------------------"<< endreq;
83 log << MSG::INFO <<"####### correction for the wire current dependence #######"<< endreq;
84 log << MSG::INFO <<"####### correction for the new z dependence #######"<< endreq;
85 log << MSG::INFO <<"-------------------------------------------------------------"<< endreq;
86 log << MSG::INFO <<"++++++++++++++++++++[ initialization ]+++++++++++++++++++++"<< endreq;
87 log << MSG::INFO << "MdcDedxRecon has been initialized with calib_flag: "<< calib_flag <<endreq;
88 log << MSG::INFO << "MdcDedxRecon has been initialized with landau: "<< landau << endreq;
89 if( landau == 0 ) {truncate = 0.7;}
90 log << MSG::INFO << "MdcDedxRecon has been initialized with ntpFlag: "<< ntpFlag << endreq;
91 log << MSG::INFO << "MdcDedxRecon has been initialized with recalg: "<< recalg << endreq;
92 log << MSG::INFO << "MdcDedxRecon has been initialized with dE/dx cal method m_alg: "<< m_alg << endreq;
93 log << MSG::INFO << "MdcDedxRecon has been initialized with truncate: "<< truncate <<endreq;
94 log << MSG::INFO << "MdcDedxRecon has been initialized with doNewGlobal: "<<doNewGlobal<<endreq;
95 log << MSG::DEBUG << "+++++++++++ MdcDedxRecon initialization end ++++++++++++ "<< endreq;
96 if( truncate <= 0.0 || 1.0 < truncate )
97 {
98 log << MSG::FATAL <<" Initialization ERROR of truncate = "<<truncate<< endreq;
99 log << MSG::FATAL <<" truncate must be within 0.0 to 1.0 ! "<< endreq;
100 log << MSG::FATAL <<" Please stop exec. "<< endreq;
101 }
102 log << MSG::DEBUG <<"-------------------------------------------------------------"<< endreq;
103 log << MSG::DEBUG <<"MdcDedxRecon init 2nd part!!!"<< endreq;
104
105 StatusCode scex = service("DedxCorrecSvc", exsvc);
106 if (scex == StatusCode::SUCCESS)
107 {
108 log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endreq;
109 }
110 else
111 {
112 log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endreq;
113 return StatusCode::SUCCESS;
114 }
115 exsvc->set_flag( calib_flag );
116
117 StatusCode cursvc = service("DedxCurSvc", dedxcursvc);
118 if (cursvc == StatusCode::SUCCESS)
119 {
120 log << MSG::INFO <<"DedxCurSvc is running"<<endreq;
121 }
122 else
123 {
124 log << MSG::ERROR <<"DedxCurSvc is failed"<<endreq;
125 return StatusCode::SUCCESS;
126 }
127
128
129 if( !ex_calib ) ex_calib = new MdcDedxCorrection;
130 //#ifdef DBCALIB
131 // ex_calib->getCalib(); //cleate MdcDedxWire and MdcDedxLayer.
132 //#endif
133
134 //------------------------------produce histogram root files-----------------------------//
135 if(m_rootfile!="no rootfile")
136 {
137 const char* preDir = gDirectory->GetPath();
138 m_hlist = new TObjArray(0);
139 m_hitlevel = new TFolder("dedx_hitlevel","hitlevel");
140 m_hlist -> Add(m_hitlevel);
141 m_hitNo_wire = new TH1F("dedx_HitNo_wire","dedx hitNo vs wire",6797, -0.5, 6796.5);
142 m_hitlevel -> Add(m_hitNo_wire);
143 gDirectory->cd(preDir);
144 }
145
146 //------------------------------produce ntuple files-------------------------------------//
147 if( ntpFlag ==2 )
148 {
149 StatusCode status;
150 NTuplePtr nt1(ntupleSvc(),"FILE103/n103");
151 if ( nt1 )
152 m_nt1 = nt1;
153 else
154 {
155 m_nt1= ntupleSvc()->book("FILE103/n103",CLID_ColumnWiseTuple,"dEdx vs momentum");
156 if ( m_nt1 )
157 {
158 status = m_nt1->addItem("phtm",m_phtm);
159 //status = m_nt1->addItem("median",m_median);
160 //status = m_nt1->addItem("geom",m_geometric);
161 //status = m_nt1->addItem("geom_tr",m_geometric_trunc);
162 //status = m_nt1->addItem("harm",m_harmonic);
163 //status = m_nt1->addItem("harm_tr",m_harmonic_trunc);
164 //status = m_nt1->addItem("transf",m_transform);
165 //status = m_nt1->addItem("logex",m_logtrunc);
166 status = m_nt1->addItem("dEdx_meas", m_dEdx_meas);
167
168 status = m_nt1->addItem("ptrk",m_ptrk);
169 status = m_nt1->addItem("sintheta",m_sintheta);
170 status = m_nt1->addItem("costheta",m_costheta);
171 status = m_nt1->addItem("charge",m_charge);
172 status = m_nt1->addItem("runNO",m_runNO);
173 status = m_nt1->addItem("evtNO",m_evtNO);
174 status = m_nt1->addItem("t0",m_t0);
175 status = m_nt1->addItem("trackId",m_trackId);
176 status = m_nt1->addItem("poca_x",m_poca_x);
177 status = m_nt1->addItem("poca_y",m_poca_y);
178 status = m_nt1->addItem("poca_z",m_poca_z);
179 status = m_nt1->addItem("recalg",m_recalg);
180
181 status = m_nt1->addItem("nhit",m_nhit);
182 status = m_nt1->addItem("usedhit",m_usedhit);
183 status = m_nt1->addItem("ndedxhit",m_nphlisthit,0,100);
184 status = m_nt1->addIndexedItem("dEdx_hit",m_nphlisthit,m_dEdx_hit);
185
186 status = m_nt1->addItem("type",m_parttype);
187 status = m_nt1->addItem("prob",m_prob);
188 status = m_nt1->addItem("expect",m_expect);
189 status = m_nt1->addItem("sigma",m_sigma);
190 status = m_nt1->addItem("chidedx",m_chidedx);
191 status = m_nt1->addItem("chidedx_e",m_chidedxe);
192 status = m_nt1->addItem("chidedx_mu",m_chidedxmu);
193 status = m_nt1->addItem("chidedx_pi",m_chidedxpi);
194 status = m_nt1->addItem("chidedx_k",m_chidedxk);
195 status = m_nt1->addItem("chidedx_p",m_chidedxp);
196 status = m_nt1->addItem("partid",5,m_probpid);
197 status = m_nt1->addItem("expectid",5,m_expectid);
198 status = m_nt1->addItem("sigmaid",5,m_sigmaid);
199 }
200 }
201
202 NTuplePtr nt2(ntupleSvc(),"FILE103/n102");
203 if ( nt2 ) m_nt2 = nt2;
204 else
205 {
206 m_nt2= ntupleSvc()->book("FILE103/n102",CLID_RowWiseTuple,"pulae height raw");
207 if ( m_nt2 )
208 {
209 status = m_nt2->addItem("charge",m_charge1);
210 status = m_nt2->addItem("adc_raw",m_phraw);
211 status = m_nt2->addItem("exraw",m_exraw);
212 status = m_nt2->addItem("runNO1",m_runNO1);
213 status = m_nt2->addItem("evtNO1",m_evtNO1);
214 status = m_nt2->addItem("nhit_hit", m_nhit_hit);
215 status = m_nt2->addItem("wire",m_wire);
216 //status = m_nt2->addItem("doca",m_doca);
217 status = m_nt2->addItem("doca_in",m_doca_in);
218 status = m_nt2->addItem("doca_ex",m_doca_ex);
219 status = m_nt2->addItem("driftdist",m_driftdist);
220 status = m_nt2->addItem("eangle",m_eangle);
221 status = m_nt2->addItem("costheta1",m_costheta1);
222 status = m_nt2->addItem("path_rphi",m_pathL);
223 status = m_nt2->addItem("layer",m_layer);
224 status = m_nt2->addItem("ptrk1",m_ptrk1);
225 status = m_nt2->addItem("ptrk_hit",m_ptrk_hit);
226 status = m_nt2->addItem("t01",m_t01);
227 status = m_nt2->addItem("tdc_raw",m_tdc_raw);
228 status = m_nt2->addItem("driftT",m_driftT);
229 status = m_nt2->addItem("localwid",m_localwid);
230 status = m_nt2->addItem("trackId1",m_trackId1);
231 }
232 }
233 }
234
235 return StatusCode::SUCCESS;
236}
237
238
239StatusCode MdcDedxRecon::beginRun()
240{
241 MsgStream log(msgSvc(), name());
242 log << MSG::DEBUG <<"in MdcDedxRecon::beginrun()!!!"<< endreq;
243
244 StatusCode gesc = service("MdcGeomSvc", geosvc);
245 if (gesc == StatusCode::SUCCESS)
246 {
247 log << MSG::INFO <<"MdcGeomSvc is running"<<endreq;
248 }
249 else
250 {
251 log << MSG::ERROR <<"MdcGeomSvc is failed"<<endreq;
252 return StatusCode::SUCCESS;
253 }
254
255 return StatusCode::SUCCESS;
256}
257
258
259StatusCode MdcDedxRecon::finalize()
260{
261 MsgStream log(msgSvc(), name());
262 log << MSG::INFO << "in finalize()" << endreq;
263
264 ex_trks.clear();
265 m_trkalgs.clear();
266
267 if(m_rootfile != "no rootfile")
268 {
269 TFile *f1 = new TFile(m_rootfile.c_str(),"recreate");
270 m_hlist->Write();
271 f1->Close();
272 delete f1;
273 delete m_hitNo_wire;
274 delete m_hitlevel;
275 delete m_hlist;
276 }
277 delete ex_calib;
278
279 cout<<"total event number is : "<<eventNo<<endl;
280 cout<<"total track number is : "<<trackNO1
281 <<" RecMdcTrack number is : "<<trackNO2
282 <<" RecMdcKalTrack number is :"<<trackNO3<<endl;
283 log << MSG::DEBUG <<"MdcDedxRecon terminated!!!"<< endreq;
284 return StatusCode::SUCCESS;
285}
286
287
288StatusCode MdcDedxRecon::execute()
289{
290 MsgStream log(msgSvc(), name());
291 log << MSG::INFO << "in execute()" << endreq;
292
293 vector<double> Curve_Para, Sigma_Para;
294 int vFlag[3];//vFlag[0]:dedx curve version; vFlag[1]:dedx sigma version; vFlag[2]: 0:data; 1:MC
295
296 for(int i=0; i< dedxcursvc->getCurveSize(); i++) // get the dedx curve parameters
297 {
298 if(i==0) vFlag[0] = (int) dedxcursvc->getCurve(i); //first element is dedx curve version
299 else Curve_Para.push_back( dedxcursvc->getCurve(i) ); //dedx curve parameters
300 }
301 for(int i=0; i< dedxcursvc->getSigmaSize(); i++)
302 {
303 if(i==0) vFlag[1] = (int) dedxcursvc->getSigma(i); //dedx sigma version: 2: psip; 3:jpsi
304 else Sigma_Para.push_back( dedxcursvc->getSigma(i) ); //dedx sigma parameters
305 }
306
307 //check if size of parameters is right
308 if(vFlag[0] ==1 && Curve_Para.size() != 5) //version 1: 5 parameters 652a psip data
309 cout<<" size of dedx curve paramters for version 1 is not right!"<<endl;
310 if(vFlag[0] ==2 && Curve_Para.size() != 16) //version 2: increase from 5 parameters of 652 to 16
311 cout<<" size of dedx curve paramters for version 2 is not right!"<<endl;
312
313 if(vFlag[1] ==1 && Sigma_Para.size() != 14) //vesion 1: 14 parameters 652a psip data (old way)
314 cout<<" size of dedx sigma paramters for version 1 is not right!"<<endl;
315 if(vFlag[1] ==2 && Sigma_Para.size() != 21) //version 2: include t0 correction (for psip09 data)
316 cout<<" size of dedx sigma paramters for version 2 is not right!"<<endl;
317 if(vFlag[1] ==3 && Sigma_Para.size() != 18) //version 3: no t0 correction (for jpsi09 data and beyond)
318 cout<<" size of dedx sigma paramters for version 3 is not right!"<<endl;
319 if(vFlag[1] ==4 && Sigma_Para.size() != 19) //version 4: data with mdc track defaulted when kal track not ok(no t0 correction)
320 cout<<" size of dedx sigma paramters for version 4 is not right!"<<endl;
321 if(vFlag[1] ==5 && Sigma_Para.size() != 22) //version 5: data with mdc track defaulted when kal track not ok(with t0 correction)
322 cout<<" size of dedx sigma paramters for version 5 is not right!"<<endl;
323
324
325 //---------- register RecMdcDedxCol and RecMdcDedxHitCol to tds---------//
326 DataObject *aReconEvent;
327 eventSvc()->findObject("/Event/Recon",aReconEvent);
328 if(aReconEvent==NULL)
329 {
330 aReconEvent = new ReconEvent();
331 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
332 if(sc!=StatusCode::SUCCESS)
333 {
334 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
335 return( StatusCode::FAILURE);
336 }
337 }
338
339 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
340
341 DataObject *aDedxcol;
342 eventSvc()->findObject("/Event/Recon/RecMdcDedxCol",aDedxcol);
343 if(aDedxcol != NULL)
344 {
345 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxCol");
346 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxCol");
347 }
348 dedxList = new RecMdcDedxCol;
349 StatusCode dedxsc;
350 dedxsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxCol, dedxList);
351 if(!dedxsc.isSuccess())
352 {
353 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
354 return ( StatusCode::FAILURE);
355 }
356
357 DataObject *aDedxhitcol;
358 eventSvc()->findObject("/Event/Recon/RecMdcDedxHitCol",aDedxhitcol);
359 if(aDedxhitcol != NULL)
360 {
361 dataManSvc->clearSubTree("/Event/Recon/RecMdcDedxHitCol");
362 eventSvc()->unregisterObject("/Event/Recon/RecMdcDedxHitCol");
363 }
364 dedxhitList = new RecMdcDedxHitCol;
365 StatusCode dedxhitsc;
366 dedxhitsc = eventSvc()->registerObject(EventModel::Recon::RecMdcDedxHitCol, dedxhitList);
367 if(!dedxhitsc.isSuccess())
368 {
369 log << MSG::FATAL << " Could not register Mdc dedx collection" <<endreq;
370 return ( StatusCode::FAILURE);
371 }
372
373 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
374 if (!eventHeader)
375 {
376 log << MSG::INFO << "Could not find Event Header" << endreq;
377 return( StatusCode::FAILURE);
378 }
379 int eventNO = eventHeader->eventNumber();
380 int runNO = eventHeader->runNumber();
381 log << MSG::INFO << "now begin the event: runNO= "<<runNO<<" eventNO= "<< eventNO<< endreq;
382 RunNO2= runNO;
383 if(RunNO1==0) RunNO1=runNO;
384 if(RunNO2 != RunNO1)
385 {
386 cout<<"RunNO2 = "<<RunNO2 <<" RunNO1 = "<<RunNO1<<endl;
387 }
388 RunNO1 = runNO;
389 int runFlag; //data type flag
390 if(runNO<MdcDedxParam::RUN0) runFlag =0; //MC: less than 0
391 else if(runNO>=MdcDedxParam::RUN1 && runNO<MdcDedxParam::RUN2) runFlag =1;
392 else if(runNO>=MdcDedxParam::RUN2 && runNO<MdcDedxParam::RUN3) runFlag =2;
393 else if(runNO>=MdcDedxParam::RUN3 && runNO<MdcDedxParam::RUN4) runFlag =3;
394 else runFlag =4;
395
396 //vFlag[2] identify MC or data: 1:Mc; 0:data
397 if(runNO < 0) vFlag[2]=1;
398 else vFlag[2]=0;
399
400 double tes = -999.0;
401 int esTimeflag = -1;
402 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
403 if( ! estimeCol)
404 {
405 log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
406 }
407 else
408 {
409 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
410 for(; iter_evt!=estimeCol->end(); iter_evt++)
411 {
412 tes = (*iter_evt)->getTest();
413 esTimeflag = (*iter_evt)->getStat();
414 }
415 }
416 //cout<<"estime : "<<tes<<endl;
417
418
419 Identifier mdcid;
420 int ntrk;
421 ex_trks.clear(); // clear the vector of MdcDedxTrk,when read a new event
422 m_trkalgs.clear();
423 if( !doNewGlobal )
424 {
425 log << MSG::DEBUG << "recalg: "<<recalg<<endreq;
426
427 //---------using kal algorithm by default, switch to seek mdc track if no kaltrack hits //
428 if(recalg ==2 )
429 {
430 //retrieve MdcKalTrackCol from TDS
431 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
432 if (!newtrkCol)
433 {
434 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
435 return StatusCode::SUCCESS;
436 }
437 if(ntpFlag>0) eventNo++;
438 ntrk = newtrkCol->size();
439 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
440
441 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
442 //cout<<"MdcDedxRecon newtrkCol.size() "<<newtrkCol->size()<<endl;
443 for( ; trk != newtrkCol->end(); trk++)
444 {
445 if(ntpFlag>0) trackNO1++;
446
447 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
448 //if not set ParticleFlag, we will get the last succefully hypothesis;
449 if(gothits.size()==0)
450 {
451 m_trkalg=0;
452 int trkid=(*trk)->trackId();
453 switchtomdctrack(trkid, mdcid, tes, runNO, eventNO, runFlag, log);
454 }
455 else
456 {
457 m_trkalg=1;
458 if(gothits.size()<2) continue;
459 kaltrackrec(trk, mdcid, tes, runNO, eventNO, runFlag, log );
460 }
461 }//end of track loop
462 }//end of recalg==2
463
464 //------------------------ Use information of MDC track rec --------------------------//
465 else if(recalg ==0 )
466 {
467 m_trkalg=0;
468
469 //retrieve MdcTrackCol from TDS
470 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
471 if (!newtrkCol)
472 {
473 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
474 return StatusCode::SUCCESS;
475 }
476 if(ntpFlag>0) eventNo++;
477 ntrk = newtrkCol->size();
478 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
479
480 vector<double> phlist; //dE/dx only after hit level correction
481 vector<double> phlist_hit; //dE/dx after hit and track level correction
482 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
483 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
484 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
485 int Nhits=0;
486 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
487
488 RecMdcTrackCol::iterator trk = newtrkCol->begin();
489 for( ; trk != newtrkCol->end(); trk++)
490 {
491 if(ntpFlag>0) trackNO1++;
492
493 MdcDedxTrk trk_ex( **trk);
494 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
495 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
496
497 HepVector a = (*trk)->helix();
498 HepSymMatrix tkErrM = (*trk)->err();
499 m_d0E = tkErrM[0][0];
500 m_phi0E = tkErrM[1][1];
501 m_cpaE = tkErrM[2][2];
502 m_z0E = tkErrM[3][3];
503
504 HepPoint3D IP(0,0,0);
505 Dedx_Helix exhel(IP, a); //initialize class Dedx_Helix for DedxCorrecSvc
506 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
507 //cout<<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl;
508 phi0= a[1];
509 costheta = cos(M_PI_2-atan(a[4]));
510 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack costheta: "<<trk_ex.theta()<<endl;
511 sintheta = sin(M_PI_2-atan(a[4]));
512 m_Pt = 1.0/fabs( a[2] );
513 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
514 charge = ( a[2] > 0 )? 1 : -1;
515 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
516 dedxhitrefvec = new DedxHitRefVec;
517
518
519 HitRefVec gothits= (*trk)->getVecHits();
520 Nhits = gothits.size();
521 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
522 if(gothits.size()<2)
523 {
524 delete dedxhitrefvec;
525 continue;
526 }
527 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
528
529 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
530 HitRefVec::iterator it_gothit = gothits.begin();
531 for( ; it_gothit != gothits.end(); it_gothit++)
532 {
533 mdcid = (*it_gothit)->getMdcId();
534 layid = MdcID::layer(mdcid);
535 localwid = MdcID::wire(mdcid);
536 w0id = geosvc->Layer(layid)->Wirst();
537 wid = w0id + localwid;
538 adc_raw = (*it_gothit)->getAdc();
539 tdc_raw = (*it_gothit)->getTdc();
540 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
541 zhit = (*it_gothit)->getZhit();
542 lr = (*it_gothit)->getFlagLR();
543 if(lr == 2) cout<<"lr = "<<lr<<endl;
544 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
545 else driftD = (*it_gothit)->getDriftDistRight();
546 //cout<<"lr: "<<lr<<" driftD: "<<driftD<<endl;
547 driftd = abs(driftD);
548 dd = (*it_gothit)->getDoca();
549 if(lr == 0 || lr == 2 ) dd = -abs(dd);
550 if(lr == 1 ) dd = abs(dd);
551 driftT = (*it_gothit)->getDriftT();
552
553 eangle = (*it_gothit)->getEntra();
554 eangle = eangle/M_PI;
555 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
556 rphi_path=pathlength * sintheta;
557
558 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd: "<<dd<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
559 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
560 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
561 //if(pathlength == -1)
562 //cout<<"parameter0: "<<"eventNO: "<<eventNO<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
563
564 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
565 if(runNO<0)
566 {
567 if (layid<8)
568 {
569 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
570 }
571 else if(layid >= 8)
572 {
573 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
574 }
575 }
576 else if(runNO>=0)
577 {
578 if(layid <8)
579 {
580 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
581 }
582 else if(layid >= 8)
583 {
584 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
585 }
586 }
587 //cout<<"recAlg=0 para mdc: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
588
589 if (ph<0.0||ph==0) continue;
590 else
591 {
592 //-----------------------put data into TDS-----------------------------//
593 dedxhit = new RecMdcDedxHit;
594 dedxhit->setMdcHit(*it_gothit);
595 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg );
596 dedxhit->setMdcId(mdcid);
597 dedxhit->setFlagLR(lr);
598 dedxhit->setTrkId(trk_ex.trk_id());
599 dedxhit->setDedx(ph_hit);
600 dedxhit->setPathLength(pathlength);
601
602 //---------------------------Fill root file----------------------------//
603 if(m_rootfile!= "no rootfile")
604 {
605 m_hitNo_wire->Fill(wid);
606 }
607
608 //-------------------------Fill ntuple n102---------------------------//
609 if ( ntpFlag ==2 )
610 {
611 m_charge1 = (*trk)->charge();
612 m_t01 = tes;
613 m_driftT = driftT;
614 m_tdc_raw = tdc_raw;
615 m_phraw = adc_raw;
616 m_exraw = ph_hit;
617 m_localwid = localwid;
618 m_wire = wid;
619 m_runNO1 = runNO;
620 m_evtNO1 = eventNO;
621 m_nhit_hit = Nhits;
622 m_doca_in = dd;
623 m_doca_ex = dd;
624 m_driftdist = driftD;
625 m_eangle = eangle;
626 m_costheta1 = costheta;
627 m_pathL = pathlength;
628 m_layer = layid;
629 m_ptrk1 = m_P;
630 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
631 m_trackId1 = trk_ex.trk_id();
632 StatusCode status2 = m_nt2->write();
633 if ( status2.isFailure() )
634 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
635 }
636 if(layid>3)
637 {
638 phlist.push_back(ph);
639 phlist_hit.push_back(ph_hit);
640 }
641 }
642 dedxhitList->push_back( dedxhit );
643 dedxhitrefvec->push_back( dedxhit );
644 }//end of hit loop
645 trk_ex.set_phlist( phlist );
646 trk_ex.set_phlist_hit( phlist_hit );
647 trk_ex.setVecDedxHits( *dedxhitrefvec );
648 ex_trks.push_back(trk_ex );
649 m_trkalgs.push_back(m_trkalg);
650
651 delete dedxhitrefvec;
652 phlist.clear();
653 phlist_hit.clear();
654 if(ntpFlag>0) trackNO2++;
655 }//end of track loop
656 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
657 }//end of recalg==0
658
659 //------------------------ Use information of MDC KAL track rec -----------------------//
660 else if(recalg ==1 )
661 {
662 m_trkalg=1;
663
664 //retrieve MdcKalTrackCol from TDS
665 SmartDataPtr<RecMdcKalTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
666 if (!newtrkCol)
667 {
668 log << MSG::WARNING << "Could not find RecMdcKalTrackCol" << endreq;
669 return StatusCode::SUCCESS;
670 }
671 if(ntpFlag>0) eventNo++;
672 ntrk = newtrkCol->size();
673 log << MSG::DEBUG << "track collection size: " << newtrkCol->size() <<endreq;
674
675 vector<double> phlist; //dE/dx only after hit level correction
676 vector<double> phlist_hit; //dE/dx after hit and track level correction
677 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
678 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
679 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
680 int Nhits=0;
681 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
682
683
684 RecMdcKalTrackCol::iterator trk = newtrkCol->begin();
685 for( ; trk != newtrkCol->end(); trk++)
686 {
687 if(ntpFlag>0) trackNO1++;
688
689 MdcDedxTrk trk_ex( **trk, ParticleFlag);
690 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
691 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
692
693 HepVector a;
694 HepSymMatrix tkErrM;
695 if(ParticleFlag>-1&&ParticleFlag<5)
696 {
697 DstMdcKalTrack* dstTrk = *trk;
698 a = dstTrk->getFHelix(ParticleFlag);
699 tkErrM = dstTrk->getFError(ParticleFlag);
700 }
701 else
702 {
703 a = (*trk)->getFHelix();
704 tkErrM = (*trk)->getFError();
705 }
706 //cout<<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl; //getFHelix is first layer of MdcKalTrack;
707 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
708
709 m_d0E = tkErrM[0][0];
710 m_phi0E = tkErrM[1][1];
711 m_cpaE = tkErrM[2][2];
712 m_z0E = tkErrM[3][3];
713
714 phi0= a[1];
715 costheta = cos(M_PI_2-atan(a[4]));
716 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl; //theta in trk_ex is got by getFHelix();
717 sintheta = sin(M_PI_2-atan(a[4]));
718 m_Pt = 1.0/fabs( a[2] );
719 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
720 //cout<<"track ptrk: "<<m_P<<" extrack ptrk: "<<trk_ex.P()<<endl;
721 charge = ( a[2] > 0 )? 1 : -1;
722 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
723 dedxhitrefvec = new DedxHitRefVec;
724
725
726 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
727 //HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
728 //if not set ParticleFlag, we will get the last succefully hypothesis;
729 Nhits = gothits.size();
730 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
731 if(gothits.size()<2)
732 {
733 delete dedxhitrefvec;
734 continue;
735 }
736 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
737
738 RecMdcHit* mdcHit = 0;
739 HelixSegRefVec::iterator it_gothit = gothits.begin();
740 for( ; it_gothit != gothits.end(); it_gothit++)
741 {
742 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
743 //HepVector a_hit_ex = (*it_gothit)->getHelixExcl(); //same with getHelixIncl()
744 HepPoint3D IP(0,0,0);
745 Dedx_Helix exhel_hit_in(IP, a_hit_in);
746 //Dedx_Helix exhel_hit_ex(IP, a_hit_ex);
747 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
748 //cout<<"getHelixIncl 5 para: "<<a_hit_in[0]<<" "<<a_hit_in[1] <<" "<<a_hit_in[2]<<" "<<a_hit_in[3]<<" "<<a_hit_in[4]<<endl;
749
750 mdcid = (*it_gothit)->getMdcId();
751 layid = MdcID::layer(mdcid);
752 localwid = MdcID::wire(mdcid);
753 w0id = geosvc->Layer(layid)->Wirst();
754 wid = w0id + localwid;
755 adc_raw = (*it_gothit)->getAdc();
756 tdc_raw = (*it_gothit)->getTdc();
757 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
758 zhit = (*it_gothit)->getZhit();
759 lr = (*it_gothit)->getFlagLR();
760 if(lr == 2) cout<<"lr = "<<lr<<endl;
761 driftD = (*it_gothit)->getDD();
762 driftd = abs(driftD);
763 driftT = (*it_gothit)->getDT();
764 dd_in = (*it_gothit)->getDocaIncl(); //getDocaIncl() include fit unused hit
765 dd_ex = (*it_gothit)->getDocaExcl(); //getDocaExcl() exclude fit unused hit
766
767 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
768 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
769 //dd = dd/doca_norm[layid];
770 //cout<<"lr "<<lr<<" dd_in: "<<dd_in<<" dd_ex: "<<dd_ex<<endl;
771
772 eangle = (*it_gothit)->getEntra();
773 eangle = eangle/M_PI;
774 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
775 rphi_path=pathlength * sintheta;
776 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd_in: "<<dd_in<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
777 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
778 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
779 //if(pathlength == -1)
780 //cout<<"parameter1: "<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
781
782 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
783 if(runNO<0)
784 {
785 if (layid<8)
786 {
787 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
788 }
789 else if(layid >= 8)
790 {
791 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
792 }
793 }
794 else if(runNO>=0)
795 {
796 if(layid <8)
797 {
798 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
799 }
800 else if(layid >= 8)
801 {
802 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
803 }
804 }
805 //cout<<"recAlg=1 para kal: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd_in<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
806
807 if (ph<0.0||ph==0) continue;
808 else
809 {
810 //-----------------------put data into TDS-----------------------------//
811 dedxhit = new RecMdcDedxHit;
812 dedxhit->setMdcHit(mdcHit);
813 dedxhit->setMdcKalHelixSeg(*it_gothit);
814 dedxhit->setMdcId(mdcid);
815 dedxhit->setFlagLR(lr);
816 dedxhit->setTrkId(trk_ex.trk_id());
817 dedxhit->setDedx(ph_hit);
818 dedxhit->setPathLength(pathlength);
819
820 //---------------------------Fill root file----------------------------//
821 if(m_rootfile!= "no rootfile")
822 {
823 m_hitNo_wire->Fill(wid);
824 }
825
826 //-------------------------Fill ntuple n102---------------------------//
827 if ( ntpFlag ==2 )
828 {
829 m_charge1 = (*trk)->charge();
830 m_t01 = tes;
831 m_driftT = driftT;
832 m_tdc_raw = tdc_raw;
833 m_phraw = adc_raw;
834 m_exraw = ph_hit;
835 m_localwid = localwid;
836 m_wire = wid;
837 m_runNO1 = runNO;
838 m_evtNO1 = eventNO;
839 m_nhit_hit = Nhits;
840 m_doca_in = dd_in;
841 m_doca_ex = dd_ex;
842 m_driftdist = driftD;
843 m_eangle = eangle;
844 m_costheta1 = costheta;
845 m_pathL = pathlength;
846 m_layer = layid;
847 m_ptrk1 = m_P;
848 m_ptrk_hit = p_hit;
849 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
850 m_trackId1 = trk_ex.trk_id();
851 StatusCode status2 = m_nt2->write();
852 if ( status2.isFailure() )
853 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
854 }
855 if(layid>3)
856 {
857 phlist.push_back(ph);
858 phlist_hit.push_back(ph_hit);
859 }
860 }
861 dedxhitList->push_back( dedxhit );
862 dedxhitrefvec->push_back( dedxhit );
863 }//end of hit loop
864 trk_ex.set_phlist( phlist );
865 trk_ex.set_phlist_hit( phlist_hit );
866 trk_ex.setVecDedxHits( *dedxhitrefvec );
867 ex_trks.push_back(trk_ex );
868 m_trkalgs.push_back(m_trkalg);
869
870 delete dedxhitrefvec;
871 phlist.clear();
872 phlist_hit.clear();
873 if(ntpFlag>0) trackNO3++;
874 }//end of track loop
875 log << MSG::DEBUG << "size in processing: " << ex_trks.size() <<endreq;
876 }//end of recalg==1
877 //--------------------------------- Hit level finished ---------------------------------//
878
879 //-------------------------------- Track level begin --------------------------------//
880 if( ntrk != ex_trks.size())
881 log << MSG::DEBUG <<"ntrkCol size: "<<ntrk<<" dedxtrk size: "<<ex_trks.size()<< endreq;
882
883 double poca_x=0,poca_y=0,poca_z=0;
884 float sintheta=0,costheta=0,ptrk=0,charge=0;
885 int Nhit=0,Nphlisthit=0;
886 int usedhit = 0;
887 double phtm=0,median=0,geometric=0,geometric_trunc=0,harmonic=0,harmonic_trunc=0,transform=0,logtrunc=0;
888
889 enum pid_dedx parType_temp(electron);
890 float probpid_temp=-0.01,expect_temp=-0.01,sigma_temp=-0.01,chidedx_temp=-0.01;
891
892 double dEdx_meas_hit=0;
893 double dEdx_meas=0,dEdx_meas_esat=0,dEdx_meas_norun=0;
894 int trk_recalg=-1;
895
896 int idedxid = 0;
897 for(std::vector<MdcDedxTrk>::iterator it = ex_trks.begin(); it != ex_trks.end(); it++, idedxid++ )
898 {
899 log << MSG::DEBUG <<"-------------------------------------------------------"<< endreq;
900 log << MSG::DEBUG <<" trk id ="<< it->trk_id()<<" : P ="<<it->P() <<" Pt ="<<it->Pt()<<" : theta ="<<it->theta() <<" : phi ="<<it->phi()<< " : charge = "<<it->charge()<<endreq;
901 log << MSG::DEBUG <<"all hits on track: "<<it->nsample()<<" phlist size: "<<it->get_phlist().size()<<endreq;
902
903 if(it->trk_ptr_kal()!=0)
904 {
905 poca_x = it->trk_ptr_kal()->x(); //get poca, default pid is pion; change pid using setPidType();
906 poca_y = it->trk_ptr_kal()->y();
907 poca_z = it->trk_ptr_kal()->z();
908 }
909 else if(it->trk_ptr()!=0)
910 {
911 poca_x = it->trk_ptr()->x();
912 poca_y = it->trk_ptr()->y();
913 poca_z = it->trk_ptr()->z();
914 }
915 //cout<<"poca_x: "<<poca_x<<" poca_y: "<<poca_y<<" poca_z: "<<poca_z<<endl;
916
917 sintheta = sin(it->theta());
918 costheta = cos(it->theta());
919 ptrk = it->P();
920 charge = it->charge();
921 Nhit = it->nsample(); //total hits on track used as sample;
922 Nphlisthit = it->get_phlist_hit().size(); //hits in phlist_hit, exclude first 4 layers;
923 //usedhit: hits after truncting phlist and used in cal dE/dx value;
924
925 //--------------------------extrk truncation--------------------------------//
926 phtm = it->cal_dedx( truncate );
927 //cout<<"phtm: "<<phtm<<endl;
928 //median = it->cal_dedx_median( truncate );
929 //geometric = it->cal_dedx_geometric( truncate );
930 //geometric_trunc = it->cal_dedx_geometric_trunc( truncate );
931 //harmonic = it->cal_dedx_harmonic( truncate );
932 //harmonic_trunc = it->cal_dedx_harmonic_trunc( truncate );
933 //transform = it->cal_dedx_transform( 0 );
934 //logtrunc = it->cal_dedx_log( 1.0, 0);
935
936 if(m_alg==1) dEdx_meas_hit = it->cal_dedx_bitrunc(truncate, 0, usedhit);
937 else if(m_alg==2) dEdx_meas_hit = it-> cal_dedx_transform(1);
938 else if(m_alg==3) dEdx_meas_hit = it->cal_dedx_log(1.0, 1);
939 else cout<<"-------------Truncate Algorithm Error!!!------------------------"<<endl;
940 if(m_alg==1 && usedhit==0)
941 {
942 cout<<"***************bad extrk with no hits!!!!!******************"<<endl;
943 continue;
944 }
945 // force to use the same definition of usedhit in TDS and what used in setDedx() function
946 usedhit = (int)(Nphlisthit*truncate);
947
948 //--------------------track level correction of extrk dE/dx Value---------------//
949 dEdx_meas = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
950 dEdx_meas_esat = exsvc->StandardTrackCorrec(1,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
951 dEdx_meas_norun = exsvc->StandardTrackCorrec(2,runFlag, ntpFlag, runNO, eventNO, dEdx_meas_hit, cos(it->theta()), tes );
952 log << MSG::INFO << "===================== " << endreq << endreq;
953 log << MSG::DEBUG <<"dEdx_meas_hit: "<<dEdx_meas_hit<<" dEdx_meas: "<<dEdx_meas<<" dEdx_meas_esat: "<<dEdx_meas_esat<<" dEdx_meas_norun: "<<dEdx_meas_norun<<" ptrk: "<<it->P()<<endreq;
954 log << MSG::INFO << "ptrk " << ptrk << " costhe " << costheta << " runno " << runNO << " evtno " << eventNO << endreq;
955 //cout<<"dEdx_meas_hit: "<<dEdx_meas_hit<<" dEdx_meas: "<<dEdx_meas<<" dEdx_meas_esat: "<<dEdx_meas_esat<<" dEdx_meas_norun: "<<dEdx_meas_norun<<" ptrk: "<<it->P()<<endl;
956
957 trk_recalg = m_trkalgs[idedxid];
958
959 if(dEdx_meas<0 || dEdx_meas==0) continue;
960 it->set_dEdx( landau , dEdx_meas, trk_recalg, runFlag, vFlag , tes, Curve_Para, Sigma_Para, ex_calib); // calculate expect
961 parType_temp = electron;
962 probpid_temp=-0.01;
963 expect_temp=-0.01;
964 sigma_temp=-0.01;
965 chidedx_temp=-0.01;
966 for(int i=0;i<5;i++)
967 {
968 if(it->pprob()[i]>probpid_temp)
969 {
970 parType_temp = (enum pid_dedx) i; //e:0, mu:1, pi:2, K:3, p:4
971 probpid_temp = it->pprob()[i];
972 expect_temp = it->pexpect()[i];
973 sigma_temp = it->pexp_sigma()[i];
974 chidedx_temp = it->pchi_dedx()[i];
975 }
976 }
977 log<< MSG::INFO<<"expect dE/dx: type: "<<parType_temp<<" probpid: "<<probpid_temp<<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq;
978 //cout<<"dEdx_meas: "<<dEdx_meas<<endl;
979
980 //-----------------------put data into TDS-----------------------------//
981 resdedx = new RecMdcDedx;
982 resdedx->setDedxHit(dEdx_meas_hit);
983 resdedx->setDedxEsat(dEdx_meas_esat);
984 resdedx->setDedxNoRun(dEdx_meas_norun);
985 resdedx->setDedxMoment(it->P());
986 resdedx->setTrackId( it->trk_id() );
987 resdedx->setProbPH( dEdx_meas );
988 resdedx->setNormPH( dEdx_meas/550.0 );
989 resdedx->setDedxExpect( it->pexpect() );
990 resdedx->setSigmaDedx( it->pexp_sigma() );
991 resdedx->setPidProb( it->pprob() );
992 resdedx->setChi( it->pchi_dedx() );
993 //cout<<"recdedx: "<<" "<<resdedx->getPidProb(parType_temp)<<" "<<resdedx->getDedxExpect(parType_temp)<<" "<<resdedx->getSigmaDedx(parType_temp)<<" "<<resdedx->chi(parType_temp)<<endl;
994 resdedx->setNumTotalHits(it->nsample() ); //all hits on track;
995 resdedx->setNumGoodHits(usedhit); //hits after truncing the phlist
996 //phlist are all hits on track excluding first 4 layers;
997 //resdedx->setStatus( it->quality() );
998 resdedx->setStatus(trk_recalg );
999 //cout<<"trk_recalg: "<<trk_recalg<<" trk stat: "<<it->quality()<<endl;
1000 resdedx->setTruncAlg( m_alg );
1001 resdedx->setParticleId(parType_temp);
1002 //cout<<"ParticleType: "<<parType_temp<<" "<<resdedx->particleType()<<endl;
1003 resdedx->setVecDedxHits(it->getVecDedxHits());
1004 resdedx->setMdcTrack(it->trk_ptr());
1005 resdedx->setMdcKalTrack(it->trk_ptr_kal());
1006
1007 //-------------------------Fill ntuple n103---------------------------//
1008 if(ntpFlag ==2)
1009 {
1010 m_phtm=phtm;
1011 //m_median=median;
1012 //m_geometric=geometric;
1013 //m_geometric_trunc=geometric_trunc;
1014 //m_harmonic=harmonic;
1015 //m_harmonic_trunc=harmonic_trunc;
1016 //m_transform=transform;
1017 //m_logtrunc=logtrunc;
1018 m_dEdx_meas = dEdx_meas;
1019
1020 m_poca_x = poca_x;
1021 m_poca_y = poca_y;
1022 m_poca_z = poca_z;
1023 m_ptrk=ptrk;
1024 m_sintheta=sintheta;
1025 m_costheta=costheta;
1026 m_charge=charge;
1027 m_runNO = runNO;
1028 m_evtNO = eventNO;
1029
1030 m_t0 = tes;
1031 m_trackId = it->trk_id();
1032 m_recalg = trk_recalg;
1033
1034 m_nhit=Nhit;
1035 m_nphlisthit=Nphlisthit;
1036 m_usedhit=usedhit;
1037 for(int i=0; i<Nphlisthit; i++) m_dEdx_hit[i] = it->get_phlist_hit()[i];
1038
1039 m_parttype = parType_temp;
1040 m_prob = probpid_temp;
1041 m_expect = expect_temp;
1042 m_sigma = sigma_temp;
1043 m_chidedx = chidedx_temp;
1044 m_chidedxe=it->pchi_dedx()[0];
1045 m_chidedxmu=it->pchi_dedx()[1];
1046 m_chidedxpi=it->pchi_dedx()[2];
1047 m_chidedxk=it->pchi_dedx()[3];
1048 m_chidedxp=it->pchi_dedx()[4];
1049 for(int i=0;i<5;i++)
1050 {
1051 m_probpid[i]= it->pprob()[i];
1052 m_expectid[i]= it->pexpect()[i];
1053 m_sigmaid[i]= it->pexp_sigma()[i];
1054 }
1055
1056 StatusCode status = m_nt1->write();
1057 if ( status.isFailure() )
1058 {
1059 log << MSG::ERROR << "Cannot fill Ntuple n103" << endreq;
1060 }
1061 }
1062
1063 log<< MSG::INFO<<"-----------------Summary of this ExTrack----------------"<<endreq
1064 <<"dEdx_mean: "<<dEdx_meas<<" type: "<<parType_temp<<" probpid: "<<probpid_temp
1065 <<" expect: "<<expect_temp<<" sigma: "<<sigma_temp<<" chidedx: "<<chidedx_temp<<endreq<<endreq;
1066
1067 dedxList->push_back( resdedx );
1068 }//ExTrk loop end
1069 } //doNewGlobal==0 END . . .
1070
1071
1072 // check MdcDedxCol registered
1073 SmartDataPtr<RecMdcDedxCol> newexCol(eventSvc(),"/Event/Recon/RecMdcDedxCol");
1074 if (!newexCol)
1075 {
1076 log << MSG::FATAL << "Could not find RecMdcDedxCol" << endreq;
1077 return( StatusCode::SUCCESS);
1078 }
1079 log << MSG::DEBUG << "----------------Begin to check RecMdcDedxCol-----------------"<<endreq;
1080 RecMdcDedxCol::iterator it_dedx = newexCol->begin();
1081 for( ; it_dedx != newexCol->end(); it_dedx++)
1082 {
1083 log << MSG::INFO << "retrieved MDC dE/dx:" << endreq
1084 << "dEdx Id: " << (*it_dedx)->trackId()
1085 << " part Id: " << (*it_dedx)->particleType()
1086 << " measured dEdx: " << (*it_dedx)->probPH()
1087 << " dEdx std: " << (*it_dedx)->normPH() << endreq
1088 << "hits on track: "<<(*it_dedx)->numTotalHits()
1089 << " usedhits: " << (*it_dedx)->numGoodHits() << endreq
1090 << "Electron: expect: " << (*it_dedx)->getDedxExpect(0)
1091 << " sigma: " << (*it_dedx)->getSigmaDedx(0)
1092 << " chi: " << (*it_dedx)->chi(0)
1093 << " prob: " << (*it_dedx)->getPidProb(0) << endreq
1094 << "Muon: expect: " << (*it_dedx)->getDedxExpect(1)
1095 << " sigma: " << (*it_dedx)->getSigmaDedx(1)
1096 << " chi: " << (*it_dedx)->chi(1)
1097 << " prob: " << (*it_dedx)->getPidProb(1) << endreq
1098 << "Pion: expect: " << (*it_dedx)->getDedxExpect(2)
1099 << " sigma: " << (*it_dedx)->getSigmaDedx(2)
1100 << " chi: " << (*it_dedx)->chi(2)
1101 << " prob: " << (*it_dedx)->getPidProb(2) << endreq
1102 << "Kaon: expect: " << (*it_dedx)->getDedxExpect(3)
1103 << " sigma: " << (*it_dedx)->getSigmaDedx(3)
1104 << " chi: " << (*it_dedx)->chi(3)
1105 << " prob: " << (*it_dedx)->getPidProb(3) << endreq
1106 << "Proton: expect: " << (*it_dedx)->getDedxExpect(4)
1107 << " sigma: " << (*it_dedx)->getSigmaDedx(4)
1108 << " chi: " << (*it_dedx)->chi(4)
1109 << " prob: " << (*it_dedx)->getPidProb(4) << endreq;
1110 }
1111 return StatusCode::SUCCESS;
1112}
1113
1114const std::vector<MdcDedxTrk>&MdcDedxRecon::tracks(void) const
1115{
1116 return ex_trks;
1117}
1118
1119void MdcDedxRecon::clearTables() {}
1120
1121void MdcDedxRecon::mdctrackrec(RecMdcTrackCol::iterator trk, Identifier mdcid, double tes, int runNO, int eventNO, int runFlag, MsgStream log )
1122{
1123 vector<double> phlist; //dE/dx only after hit level correction
1124 vector<double> phlist_hit; //dE/dx after hit and track level correction
1125 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
1126 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
1127 double adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd=0,eangle=0;
1128 int Nhits=0;
1129 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
1130
1131 MdcDedxTrk trk_ex( **trk);
1132 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
1133 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
1134
1135 HepVector a = (*trk)->helix();
1136 HepSymMatrix tkErrM = (*trk)->err();
1137 m_d0E = tkErrM[0][0];
1138 m_phi0E = tkErrM[1][1];
1139 m_cpaE = tkErrM[2][2];
1140 m_z0E = tkErrM[3][3];
1141
1142 HepPoint3D IP(0,0,0);
1143 Dedx_Helix exhel(IP, a); //initialize class Dedx_Helix for DedxCorrecSvc
1144 log << MSG::DEBUG <<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
1145 //cout<<"MDC Helix 5 parameters: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl;
1146 phi0= a[1];
1147 costheta = cos(M_PI_2-atan(a[4]));
1148 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl;
1149 sintheta = sin(M_PI_2-atan(a[4]));
1150 m_Pt = 1.0/fabs( a[2] );
1151 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
1152 charge = ( a[2] > 0 )? 1 : -1;
1153 //cout<<"track charge: "<<charge<<" extrack charge: "<<trk_ex.charge()<<endl;
1154 dedxhitrefvec = new DedxHitRefVec;
1155
1156
1157 HitRefVec gothits= (*trk)->getVecHits();
1158 Nhits = gothits.size();
1159 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
1160 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
1161
1162 RecMdcKalHelixSeg* mdcKalHelixSeg = 0;
1163 HitRefVec::iterator it_gothit = gothits.begin();
1164 for( ; it_gothit != gothits.end(); it_gothit++)
1165 {
1166 mdcid = (*it_gothit)->getMdcId();
1167 layid = MdcID::layer(mdcid);
1168 localwid = MdcID::wire(mdcid);
1169 w0id = geosvc->Layer(layid)->Wirst();
1170 wid = w0id + localwid;
1171 adc_raw = (*it_gothit)->getAdc();
1172 tdc_raw = (*it_gothit)->getTdc();
1173 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
1174 zhit = (*it_gothit)->getZhit();
1175 lr = (*it_gothit)->getFlagLR();
1176 if(lr == 2) cout<<"lr = "<<lr<<endl;
1177 if(lr == 0 || lr == 2) driftD = (*it_gothit)->getDriftDistLeft();
1178 else driftD = (*it_gothit)->getDriftDistRight();
1179 //cout<<"lr: "<<lr<<" driftD: "<<driftD<<endl;
1180 driftd = abs(driftD);
1181 dd = (*it_gothit)->getDoca();
1182 if(lr == 0 || lr == 2 ) dd = -abs(dd);
1183 if(lr == 1 ) dd = abs(dd);
1184 driftT = (*it_gothit)->getDriftT();
1185
1186 eangle = (*it_gothit)->getEntra();
1187 eangle = eangle/M_PI;
1188 pathlength=exsvc->PathL( ntpFlag, exhel, layid, localwid, zhit);
1189 rphi_path=pathlength * sintheta;
1190
1191 //cout<<"ph para check: "<<"adc_raw: "<<adc_raw<<" dd: "<<dd<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
1192 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd,eangle,costheta);
1193 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
1194 //if(pathlength == -1)
1195 //cout<<"parameter0: "<<"eventNO: "<<eventNO<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
1196
1197 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
1198 if(runNO<0)
1199 {
1200 if (layid<8)
1201 {
1202 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1203 }
1204 else if(layid >= 8)
1205 {
1206 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1207 }
1208 }
1209 else if(runNO>=0)
1210 {
1211 if(layid <8)
1212 {
1213 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1214 }
1215 else if(layid >= 8)
1216 {
1217 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1218 }
1219 }
1220 //cout<<"recAlg=2 para mdc: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
1221
1222 if (ph<0.0||ph==0) continue;
1223 else
1224 {
1225 //-----------------------put data into TDS-----------------------------//
1226 dedxhit = new RecMdcDedxHit;
1227 dedxhit->setMdcHit(*it_gothit);
1228 dedxhit->setMdcKalHelixSeg(mdcKalHelixSeg);
1229 dedxhit->setMdcId(mdcid);
1230 dedxhit->setFlagLR(lr);
1231 dedxhit->setTrkId(trk_ex.trk_id());
1232 dedxhit->setDedx(ph_hit);
1233 dedxhit->setPathLength(pathlength);
1234
1235 //---------------------------Fill root file----------------------------//
1236 if(m_rootfile!= "no rootfile")
1237 {
1238 m_hitNo_wire->Fill(wid);
1239 }
1240
1241 //-------------------------Fill ntuple n102---------------------------//
1242 if ( ntpFlag ==2 )
1243 {
1244 m_charge1 = (*trk)->charge();
1245 m_t01 = tes;
1246 m_driftT = driftT;
1247 m_tdc_raw = tdc_raw;
1248 m_phraw = adc_raw;
1249 m_exraw = ph_hit;
1250 m_localwid = localwid;
1251 m_wire = wid;
1252 m_runNO1 = runNO;
1253 m_evtNO1 = eventNO;
1254 m_nhit_hit = Nhits;
1255 m_doca_in = dd;
1256 m_doca_ex = dd;
1257 m_driftdist = driftD;
1258 m_eangle = eangle;
1259 m_costheta1 = costheta;
1260 m_pathL = pathlength;
1261 m_layer = layid;
1262 m_ptrk1 = m_P;
1263 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
1264 m_trackId1 = trk_ex.trk_id();
1265 StatusCode status2 = m_nt2->write();
1266 if ( status2.isFailure() )
1267 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
1268 }
1269 if(layid>3)
1270 {
1271 phlist.push_back(ph);
1272 phlist_hit.push_back(ph_hit);
1273 }
1274 }
1275 dedxhitList->push_back( dedxhit );
1276 dedxhitrefvec->push_back( dedxhit );
1277 }//end of hit loop
1278 trk_ex.set_phlist( phlist );
1279 trk_ex.set_phlist_hit( phlist_hit );
1280 trk_ex.setVecDedxHits( *dedxhitrefvec );
1281 ex_trks.push_back(trk_ex );
1282 m_trkalgs.push_back(m_trkalg);
1283
1284 delete dedxhitrefvec;
1285 phlist.clear();
1286 phlist_hit.clear();
1287 if(ntpFlag>0) trackNO2++;
1288}
1289
1290
1291void MdcDedxRecon::kaltrackrec(RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int runNO, int eventNO, int runFlag, MsgStream log )
1292{
1293 vector<double> phlist; //dE/dx only after hit level correction
1294 vector<double> phlist_hit; //dE/dx after hit and track level correction
1295 double m_d0E=0,m_phi0E=0,m_cpaE=0,m_z0E=0,phi0=0,costheta=0,sintheta=0,m_Pt=0,m_P=0;
1296 int charge=0,layid=0,localwid=0,w0id=0,wid=0,lr=-2;
1297 double p_hit=0,adc_raw=0,tdc_raw=0,zhit=0,driftd=0,driftD=0,driftT=0,dd_in=0,dd_ex=0,eangle=0;
1298 int Nhits=0;
1299 double ph=0,ph_hit=0,pathlength=0,rphi_path=0;
1300
1301 MdcDedxTrk trk_ex( **trk, ParticleFlag);
1302 log << MSG::DEBUG <<" %%%%% trk id = "<<trk_ex.trk_id() <<endreq;
1303 log << MSG::DEBUG <<" %%%%% initial id = "<<(*trk)->trackId() <<endreq;
1304
1305 HepVector a;
1306 HepSymMatrix tkErrM;
1307 if(ParticleFlag>-1&&ParticleFlag<5)
1308 {
1309 DstMdcKalTrack* dstTrk = *trk;
1310 a = dstTrk->getFHelix(ParticleFlag);
1311 tkErrM = dstTrk->getFError(ParticleFlag);
1312 }
1313 else
1314 {
1315 a = (*trk)->getFHelix();
1316 tkErrM = (*trk)->getFError();
1317 }
1318 log << MSG::DEBUG <<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endreq;
1319 //cout<<"FHelix 5 parameters: "<<a[0]<<" "<<a[1] <<" "<<a[2]<<" "<<a[3]<<" "<<a[4]<<endl; //getFHelix is first layer of MdcKalTrack;
1320
1321 m_d0E = tkErrM[0][0];
1322 m_phi0E = tkErrM[1][1];
1323 m_cpaE = tkErrM[2][2];
1324 m_z0E = tkErrM[3][3];
1325
1326 phi0= a[1];
1327 costheta = cos(M_PI_2-atan(a[4]));
1328 //cout<<"track theta: "<<M_PI_2-atan(a[4])<<" extrack theta: "<<trk_ex.theta()<<endl; //theta in trk_ex is got by getFHelix();
1329 sintheta = sin(M_PI_2-atan(a[4]));
1330 m_Pt = 1.0/fabs( a[2] );
1331 m_P = m_Pt*sqrt(1 + a[4]*a[4]);
1332 //cout<<"track ptrk: "<<m_P<<" extrack ptrk: "<<trk_ex.P()<<endl;
1333 charge = ( a[2] > 0 )? 1 : -1;
1334 //cout<<"track charge: "<<charge<<" extrack charge: "<<(*trk)->charge()<<endl;
1335 dedxhitrefvec = new DedxHitRefVec;
1336
1337
1338 HelixSegRefVec gothits= (*trk)->getVecHelixSegs(ParticleFlag);
1339 //HelixSegRefVec gothits= (*trk)->getVecHelixSegs();
1340 //if not set ParticleFlag, we will get the last succefully hypothesis;
1341 Nhits = gothits.size();
1342 log << MSG::DEBUG << "hits size = " <<gothits.size()<<endreq;
1343 //if(gothits.size()<15) cout<<"eventNO : "<<eventNO<<" hits : "<<gothits.size()<<endl;
1344
1345 RecMdcHit* mdcHit = 0;
1346 HelixSegRefVec::iterator it_gothit = gothits.begin();
1347 for( ; it_gothit != gothits.end(); it_gothit++)
1348 {
1349 HepVector a_hit_in = (*it_gothit)->getHelixIncl();
1350 //HepVector a_hit_ex = (*it_gothit)->getHelixExcl(); //same with getHelixIncl()
1351 HepPoint3D IP(0,0,0);
1352 Dedx_Helix exhel_hit_in(IP, a_hit_in);
1353 //Dedx_Helix exhel_hit_ex(IP, a_hit_ex);
1354 p_hit = 1.0/fabs(a_hit_in(3))*sqrt(1+a_hit_in(5)*a_hit_in(5));
1355 //cout<<"getHelixIncl 5 para: "<<a_hit_in[0]<<" "<<a_hit_in[1] <<" "<<a_hit_in[2]<<" "<<a_hit_in[3]<<" "<<a_hit_in[4]<<endl;
1356
1357 mdcid = (*it_gothit)->getMdcId();
1358 layid = MdcID::layer(mdcid);
1359 localwid = MdcID::wire(mdcid);
1360 w0id = geosvc->Layer(layid)->Wirst();
1361 wid = w0id + localwid;
1362 adc_raw = (*it_gothit)->getAdc();
1363 tdc_raw = (*it_gothit)->getTdc();
1364 log << MSG::INFO <<"hit layer id " <<layid<< " wire id: " <<localwid<< " adc_raw: " <<adc_raw<<" tdc_raw: "<<tdc_raw<<endreq;
1365 zhit = (*it_gothit)->getZhit();
1366 lr = (*it_gothit)->getFlagLR();
1367 if(lr == 2) cout<<"lr = "<<lr<<endl;
1368 driftD = (*it_gothit)->getDD();
1369 driftd = abs(driftD);
1370 driftT = (*it_gothit)->getDT();
1371 dd_in = (*it_gothit)->getDocaIncl(); //getDocaIncl() include fit unused hit
1372 dd_ex = (*it_gothit)->getDocaExcl(); //getDocaExcl() exclude fit unused hit
1373
1374 if(lr==-1 || lr == 2) {dd_in = dd_in; dd_ex = dd_ex;}
1375 else if(lr ==1) {dd_in = -dd_in; dd_ex = -dd_ex;}
1376 //dd = dd/doca_norm[layid];
1377 //cout<<"lr: "<<lr<<" dd_in: "<<dd_in<<" dd_ex: "<<dd_ex<<endl;
1378
1379 eangle = (*it_gothit)->getEntra();
1380 eangle = eangle/M_PI;
1381 pathlength=exsvc->PathL( ntpFlag, exhel_hit_in, layid, localwid, zhit);
1382 rphi_path=pathlength * sintheta;
1383
1384 //cout<<"ph para check: "<<"runFlag: "<<runFlag<<" runNO: "<<runNO<<" pathlength: "<<pathlength<<" wid: "<<wid<<" layid: "<<layid<<" adc_raw: "<<adc_raw<<" dd_in: "<<dd_in<<" eangle: "<<eangle<<" zhit: "<<zhit<<" costheta: "<<costheta<<endl;
1385 ph=exsvc->StandardHitCorrec(0,runFlag,ntpFlag,runNO,eventNO,pathlength,wid,layid,adc_raw,dd_in,eangle,costheta);
1386 //cout<<"tes= "<<tes<<endl;
1387 ph_hit = exsvc->StandardTrackCorrec(0,runFlag, ntpFlag, runNO, eventNO, ph, costheta,tes );
1388 //cout<<"adc_raw= "<<adc_raw<<" ph= "<<ph<<endl;
1389 //cout<<"adc_raw = "<<adc_raw<<" ph_hit= "<<ph_hit<<endl;
1390 //if(pathlength == -1)
1391 //cout<<"parameter1: "<<" layid: "<<layid<<" localwid: "<<localwid<<" driftd: "<<driftd<<" rphi_path: "<<rphi_path<<" pathlength: "<<pathlength<<" ph: "<<ph<<" ph_hit: "<<ph_hit<<endl;
1392
1393 if(runNO<=5911 && runNO>=5854 && layid ==14) continue;
1394 if(runNO<0)
1395 {
1396 if (layid<8)
1397 {
1398 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1399 }
1400 else if(layid >= 8)
1401 {
1402 if(adc_raw<MinHistValue || adc_raw>MaxHistValue ||rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1403 }
1404 }
1405 else if(runNO>=0)
1406 {
1407 if(layid <8)
1408 {
1409 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;
1410 }
1411 else if(layid >= 8)
1412 {
1413 if(adc_raw<MinHistValue || adc_raw>MaxHistValue || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut ) continue;
1414 }
1415 }
1416 //cout<<"recAlg=2 para kal: m_phraw: "<<adc_raw<<" m_exraw: "<<ph_hit<<" m_doca: "<<dd_in<<" m_eangle: "<<eangle<<" m_costheta: "<<costheta<<endl;
1417
1418 if (ph<0.0||ph==0) continue;
1419 else
1420 {
1421 //-----------------------put data into TDS-----------------------------//
1422 dedxhit = new RecMdcDedxHit;
1423 dedxhit->setMdcHit(mdcHit);
1424 dedxhit->setMdcKalHelixSeg(*it_gothit);
1425 dedxhit->setMdcId(mdcid);
1426 dedxhit->setFlagLR(lr);
1427 dedxhit->setTrkId(trk_ex.trk_id());
1428 dedxhit->setDedx(ph_hit);
1429 dedxhit->setPathLength(pathlength);
1430
1431 //---------------------------Fill root file----------------------------//
1432 if(m_rootfile!= "no rootfile")
1433 {
1434 m_hitNo_wire->Fill(wid);
1435 }
1436
1437 //-------------------------Fill ntuple n102---------------------------//
1438 if ( ntpFlag ==2 )
1439 {
1440 m_charge1 = (*trk)->charge();
1441 m_t01 = tes;
1442 m_driftT = driftT;
1443 m_tdc_raw = tdc_raw;
1444 m_phraw = adc_raw;
1445 m_exraw = ph_hit;
1446 m_localwid = localwid;
1447 m_wire = wid;
1448 m_runNO1 = runNO;
1449 m_evtNO1 = eventNO;
1450 m_nhit_hit = Nhits;
1451 m_doca_in = dd_in;
1452 m_doca_ex = dd_ex;
1453 m_driftdist = driftD;
1454 m_eangle = eangle;
1455 m_costheta1 = costheta;
1456 m_pathL = pathlength;
1457 m_layer = layid;
1458 m_ptrk1 = m_P;
1459 m_ptrk_hit = p_hit;
1460 //cout<<"adc_raw : "<<adc_raw<<" "<<ph_hit<<" "<<dd_in<<" "<<dd_ex<<" "<<eangle<<endl;
1461 m_trackId1 = trk_ex.trk_id();
1462 StatusCode status2 = m_nt2->write();
1463 if ( status2.isFailure() )
1464 log << MSG::ERROR << "Cannot fill Ntuple n102" << endreq;
1465 }
1466 if(layid>3)
1467 {
1468 phlist.push_back(ph);
1469 phlist_hit.push_back(ph_hit);
1470 }
1471 }
1472 dedxhitList->push_back( dedxhit );
1473 dedxhitrefvec->push_back( dedxhit );
1474 }//end of hit loop
1475 trk_ex.set_phlist( phlist );
1476 trk_ex.set_phlist_hit( phlist_hit );
1477 trk_ex.setVecDedxHits( *dedxhitrefvec );
1478 ex_trks.push_back(trk_ex );
1479 m_trkalgs.push_back(m_trkalg);
1480
1481 delete dedxhitrefvec;
1482 phlist.clear();
1483 phlist_hit.clear();
1484 if(ntpFlag>0) trackNO3++;
1485}
1486
1487void MdcDedxRecon::switchtomdctrack(int trkid,Identifier mdcid, double tes, int runNO, int eventNO, int runFlag, MsgStream log)
1488{
1489 //retrieve MdcTrackCol from TDS
1490 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
1491 if (!newtrkCol)
1492 {
1493 log << MSG::WARNING << "Could not find RecMdcTrackCol in switchtomdctrack" << endreq;
1494 return ;
1495 }
1496
1497 RecMdcTrackCol::iterator trk = newtrkCol->begin();
1498 for( ; trk != newtrkCol->end(); trk++)
1499 {
1500 if( (*trk)->trackId()==trkid)
1501 {
1502 HitRefVec gothits= (*trk)->getVecHits();
1503 if(gothits.size()>=2)
1504 mdctrackrec(trk,mdcid,tes,runNO,eventNO,runFlag,log);
1505 }
1506 }
1507}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
Add
mg Add(gr3)
f1
TFile * f1
Definition: DataBase/tau_mode.c:4
Iner_DriftDistCut
#define Iner_DriftDistCut
Definition: DedxCalibParameters.h:58
RPhi_PathMaxCut
#define RPhi_PathMaxCut
Definition: DedxCalibParameters.h:55
Out_DriftDistCut
#define Out_DriftDistCut
Definition: DedxCalibParameters.h:59
MaxHistValue
#define MaxHistValue
Definition: DedxCalibParameters.h:47
pid_dedx
pid_dedx
Definition: DstMdcDedx.h:9
electron
@ electron
Definition: DstMdcDedx.h:9
trackNO3
int trackNO3
Definition: MdcDedxRecon.cxx:51
eventNo
int eventNo
Definition: MdcDedxRecon.cxx:48
RunNO1
int RunNO1
Definition: MdcDedxRecon.cxx:45
trackNO2
int trackNO2
Definition: MdcDedxRecon.cxx:50
RunNO2
int RunNO2
Definition: MdcDedxRecon.cxx:46
prob_
float prob_(float *, int *)
trackNO1
int trackNO1
Definition: MdcDedxRecon.cxx:49
RecMdcDedxHitCol
ObjectVector< RecMdcDedxHit > RecMdcDedxHitCol
Definition: RecMdcDedxHit.h:78
DedxHitRefVec
SmartRefVector< RecMdcDedxHit > DedxHitRefVec
Definition: RecMdcDedx.h:27
RecMdcDedxCol
ObjectVector< RecMdcDedx > RecMdcDedxCol
Definition: RecMdcDedx.h:132
HelixSegRefVec
SmartRefVector< RecMdcKalHelixSeg > HelixSegRefVec
Definition: RecMdcKalTrack.h:27
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:22
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
M_PI
#define M_PI
Definition: TConstant.h:4
NULL
#define NULL
Definition: VertexExtrapolate.h:14
Dedx_Helix
Helix parameter class.
Definition: Dedx_Helix.h:33
DstMdcDedx::setChi
void setChi(double *chi)
Definition: DstMdcDedx.h:77
DstMdcDedx::setTruncAlg
void setTruncAlg(int trunc_alg)
Definition: DstMdcDedx.h:75
DstMdcDedx::setStatus
void setStatus(int status)
Definition: DstMdcDedx.h:74
DstMdcDedx::setNumGoodHits
void setNumGoodHits(int numGoodHits)
Definition: DstMdcDedx.h:81
DstMdcDedx::setProbPH
void setProbPH(double probPH)
Definition: DstMdcDedx.h:83
DstMdcDedx::setNormPH
void setNormPH(double normPH)
Definition: DstMdcDedx.h:84
DstMdcDedx::setParticleId
void setParticleId(int particleId)
Definition: DstMdcDedx.h:73
DstMdcDedx::setTrackId
void setTrackId(int trackId)
Definition: DstMdcDedx.h:72
DstMdcDedx::setNumTotalHits
void setNumTotalHits(int numTotalHits)
Definition: DstMdcDedx.h:82
DstMdcKalTrack::getFError
const HepSymMatrix & getFError(const int pid) const
Definition: DstMdcKalTrack.h:213
DstMdcKalTrack::getFHelix
const HepVector & getFHelix(const int pid) const
Definition: DstMdcKalTrack.h:210
IDedxCorrecSvc::set_flag
virtual void set_flag(int calib_F)=0
IDedxCorrecSvc::StandardHitCorrec
virtual double StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0
IDedxCorrecSvc::PathL
virtual double PathL(int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const =0
IDedxCorrecSvc::StandardTrackCorrec
virtual double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const =0
IDedxCurSvc::getSigmaSize
virtual const int getSigmaSize()=0
IDedxCurSvc::getCurveSize
virtual const int getCurveSize()=0
IDedxCurSvc::getCurve
virtual const double getCurve(int i)=0
IDedxCurSvc::getSigma
virtual const double getSigma(int i)=0
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
MdcDedxParam::RUN0
static long int RUN0
Definition: MdcDedxParam.h:36
MdcDedxParam::RUN2
static long int RUN2
Definition: MdcDedxParam.h:38
MdcDedxParam::RUN1
static long int RUN1
Definition: MdcDedxParam.h:37
MdcDedxParam::RUN4
static long int RUN4
Definition: MdcDedxParam.h:40
MdcDedxParam::RUN3
static long int RUN3
Definition: MdcDedxParam.h:39
MdcDedxRecon::mdctrackrec
void mdctrackrec(RecMdcTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1121
MdcDedxRecon::switchtomdctrack
void switchtomdctrack(int trkid, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1487
MdcDedxRecon::initialize
StatusCode initialize()
Definition: MdcDedxRecon.cxx:77
MdcDedxRecon::execute
StatusCode execute()
Definition: MdcDedxRecon.cxx:288
MdcDedxRecon::tracks
const std::vector< MdcDedxTrk > & tracks(void) const
Definition: MdcDedxRecon.cxx:1114
MdcDedxRecon::finalize
StatusCode finalize()
Definition: MdcDedxRecon.cxx:259
MdcDedxRecon::kaltrackrec
void kaltrackrec(RecMdcKalTrackCol::iterator trk, Identifier mdcid, double tes, int RunNO, int eventNO, int runFlag, MsgStream log)
Definition: MdcDedxRecon.cxx:1291
MdcDedxRecon::beginRun
StatusCode beginRun()
Definition: MdcDedxRecon.cxx:239
MdcDedxTrk::trk_id
int trk_id() const
Definition: MdcDedxTrk.h:58
MdcDedxTrk::setVecDedxHits
void setVecDedxHits(const DedxHitRefVec &vecdedxhit)
Definition: MdcDedxTrk.h:47
MdcDedxTrk::set_phlist_hit
void set_phlist_hit(const vector< double > &phlist)
Definition: MdcDedxTrk.h:46
MdcDedxTrk::set_phlist
void set_phlist(const vector< double > &phlist)
Definition: MdcDedxTrk.h:45
MdcGeoLayer::Wirst
int Wirst(void) const
Definition: MdcGeoLayer.h:157
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
RecMdcDedxHit::setMdcId
void setMdcId(Identifier mdcid)
Definition: RecMdcDedxHit.h:64
RecMdcDedxHit::setDedx
void setDedx(double dedx)
Definition: RecMdcDedxHit.h:62
RecMdcDedxHit::setFlagLR
void setFlagLR(int lr)
Definition: RecMdcDedxHit.h:61
RecMdcDedxHit::setMdcKalHelixSeg
void setMdcKalHelixSeg(const RecMdcKalHelixSeg *mdcKalHelixSeg)
Definition: RecMdcDedxHit.h:57
RecMdcDedxHit::setTrkId
void setTrkId(int trkid)
Definition: RecMdcDedxHit.h:60
RecMdcDedxHit::setMdcHit
void setMdcHit(const RecMdcHit *mdcHit)
Definition: RecMdcDedxHit.h:58
RecMdcDedxHit::setPathLength
void setPathLength(double pathlength)
Definition: RecMdcDedxHit.h:63
RecMdcDedx::setMdcTrack
void setMdcTrack(RecMdcTrack *trk)
Definition: RecMdcDedx.h:107
RecMdcDedx::setVecDedxHits
void setVecDedxHits(const DedxHitRefVec &vecdedxhit)
Definition: RecMdcDedx.h:76
RecMdcDedx::setDedxNoRun
void setDedxNoRun(double dedx_norun)
Definition: RecMdcDedx.h:80
RecMdcDedx::setMdcKalTrack
void setMdcKalTrack(RecMdcKalTrack *trk)
Definition: RecMdcDedx.h:108
RecMdcDedx::setDedxMoment
void setDedxMoment(double dedx_momentum)
Definition: RecMdcDedx.h:81
RecMdcDedx::setDedxExpect
void setDedxExpect(double *dedx_exp)
Definition: RecMdcDedx.h:90
RecMdcDedx::setSigmaDedx
void setSigmaDedx(double *sigma_dedx)
Definition: RecMdcDedx.h:94
RecMdcDedx::setDedxEsat
void setDedxEsat(double dedx_esat)
Definition: RecMdcDedx.h:79
RecMdcDedx::setDedxHit
void setDedxHit(double dedx_hit)
Definition: RecMdcDedx.h:78
RecMdcDedx::setPidProb
void setPidProb(double *pid_prob)
Definition: RecMdcDedx.h:98
EventModel::Recon::RecMdcDedxCol
_EXTERN_ std::string RecMdcDedxCol
Definition: EventModel.h:88
EventModel::Recon::RecMdcDedxHitCol
_EXTERN_ std::string RecMdcDedxHitCol
Definition: EventModel.h:89
costheta
float costheta
Definition: prepare_data.cxx:28
charge
float charge
Definition: prepare_data.cxx:28
dEdx_meas
float dEdx_meas
Definition: prepare_data.cxx:28
ptrk
float ptrk
Definition: prepare_data.cxx:28