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