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