BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcNavigation.cxx
Go to the documentation of this file.
1#include "MdcNavigation/MdcNavigation.h"
2#include "CLHEP/Geometry/Point3D.h"
3#include "EventModel/Event.h"
4#include "EventModel/EventHeader.h"
5#include "EventModel/EventModel.h"
6#include "EventNavigator/EventNavigator.h"
7#include "GaudiKernel/IPartPropSvc.h"
8#include "GaudiKernel/MsgStream.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "HepPDT/ParticleData.hh"
11#include "HepPDT/ParticleDataTable.hh"
12#include "MdcGeom/Constants.h"
13#include "RawEvent/RawDataUtil.h"
14#include <cmath>
15
16#include "CLHEP/Matrix/SymMatrix.h"
17#include "MdcData/MdcHit.h"
18#include "MdcGeom/BesAngle.h"
19#include "MdcGeom/MdcDetector.h"
20#include "MdcGeom/MdcLayer.h"
21#include "TrkBase/HelixTraj.h"
22#include "TrkBase/TrkPoca.h"
23#include "TrkBase/TrkPocaBase.h"
24
25#include "EvTimeEvent/RecEsTime.h"
26#include "Identifier/MdcID.h"
27#include "MdcRawEvent/MdcDigi.h"
28
29#include "GaudiKernel/INTupleSvc.h"
30#include "GaudiKernel/NTuple.h"
31using namespace Event;
32#ifndef ENABLE_BACKWARDS_COMPATIBILITY
33// backwards compatibility will be enabled ONLY in CLHEP 1.9
34typedef HepGeom::Point3D<double> HepPoint3D;
35#endif
36
37using namespace std;
38using namespace Event;
39const double epsilon = 0.000000001;
40DECLARE_COMPONENT( MdcNavigation )
41MdcNavigation::MdcNavigation( const std::string& name, ISvcLocator* pSvcLocator )
42 : Algorithm( name, pSvcLocator ) {
43 declareProperty( "hist", m_hist = 0 );
44 declareProperty( "nMcHit", m_nMcHit = 5 );
45 declareProperty( "mc", m_mc = 1 );
46
47 declareProperty( "maxMdcDigi", m_maxMdcDigi = 0 );
48 declareProperty( "keepBadTdc", m_keepBadTdc = 0 );
49 declareProperty( "dropHot", m_dropHot = 0 );
50 declareProperty( "keepUnmatch", m_keepUnmatch = 0 );
51
52 declareProperty( "poca", m_poca = false );
53 declareProperty( "doSag", m_doSag = false );
54
55 declareProperty( "d0Cut", m_d0Cut = 1. );
56 declareProperty( "z0Cut", m_z0Cut = 10. );
57 declareProperty( "debug", m_debug = 0 );
58}
59
60// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
61
62StatusCode MdcNavigation::initialize() {
63 MsgStream log( msgSvc(), name() );
64 StatusCode sc = StatusCode::SUCCESS;
65
66 t_nTk = 0;
67 // Initailize magnetic filed
68 sc = service( "MagneticFieldSvc", m_pIMF );
69 if ( sc != StatusCode::SUCCESS )
70 {
71 log << MSG::ERROR << "Unable to open Magnetic field service" << endmsg;
72 return StatusCode::FAILURE;
73 }
74
75 // Get the Particle Properties Service
76 IPartPropSvc* p_PartPropSvc;
77 static const bool CREATEIFNOTTHERE( true );
78 sc = service( "PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE );
79 if ( !sc.isSuccess() || 0 == p_PartPropSvc )
80 {
81 log << MSG::ERROR << " Could not initialize PartPropSvc" << endmsg;
82 return sc;
83 }
84
85 m_particleTable = p_PartPropSvc->PDT();
86
87 IRawDataProviderSvc* irawDataProviderSvc;
88 sc = service( "RawDataProviderSvc", irawDataProviderSvc );
89 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>( irawDataProviderSvc );
90 if ( sc.isFailure() )
91 {
92 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endmsg;
93 return StatusCode::FAILURE;
94 }
95
96 if ( m_hist )
97 {
98 sc = bookNTuple();
99 if ( !sc.isSuccess() )
100 {
101 log << MSG::WARNING << " Could not book NTuple" << endmsg;
102 m_hist = 0;
103 }
104 }
105
106 keepedParticles = new int( 211 );
107
108 return StatusCode::SUCCESS;
109}
110
111StatusCode MdcNavigation::beginRun() {
112 MsgStream log( msgSvc(), name() );
113 log << MSG::INFO << "in beginRun()" << endmsg;
114
115 m_gm = MdcDetector::instance( m_doSag );
116 if ( NULL == m_gm ) return StatusCode::FAILURE;
117
118 return StatusCode::SUCCESS;
119}
120
121// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
122StatusCode MdcNavigation::execute() {
123 setFilterPassed( false );
124 MsgStream log( msgSvc(), name() );
125 StatusCode sc = StatusCode::SUCCESS;
126
127 // Get EventNavigator from the TDS
128 SmartDataPtr<EventNavigator> navigator( eventSvc(), "/Event/Navigator" );
129 if ( !navigator )
130 {
131 log << MSG::WARNING << " Unable to retrieve EventNavigator" << endmsg;
132 m_rawData = true;
133 }
134 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
135 SmartDataPtr<RecMdcHitCol> recMdcHitCol( eventSvc(), "/Event/Recon/RecMdcHitCol" );
136
137 // get eventNo, t0 and MdcDigi
138 if ( m_hist )
139 {
140 sc = fillInit();
141 if ( sc != StatusCode::SUCCESS ) { return StatusCode::FAILURE; }
142 }
143 if ( m_mc )
144 {
145 // Get McParticleCol
146 SmartDataPtr<Event::McParticleCol> mcParticles( eventSvc(), "/Event/MC/McParticleCol" );
147 SmartDataPtr<Event::MdcMcHitCol> mcHit( eventSvc(), "/Event/MC/MdcMcHitCol" );
148 if ( !mcParticles )
149 { log << MSG::WARNING << " Unable to retrieve McParticleCol" << endmsg; }
150 else
151 {
152 // For each McParticle ...
153 t_mcTkNum = 0;
154 McParticleCol::iterator it = mcParticles->begin();
155 log << MSG::INFO << "mcParticles size = " << mcParticles->size()
156 << endmsg; // yzhang debug
157 for ( ; it != mcParticles->end(); it++ )
158 {
159 // int tkId = (*it)->trackIndex();
160 t_mcTkNum++;
161 }
162 }
163 }
164 t_mcTkNum = 0;
165 t_recTkNum = 0;
166 // for each rec track
167
168 if ( !recMdcTrackCol )
169 {
170 log << MSG::WARNING << " Unable to retrieve recMdcTrackCol" << endmsg;
171 return StatusCode::SUCCESS;
172 }
173 t_recTkNum = recMdcTrackCol->size();
174
175 //=============loop over tracks==============
176 RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
177 for ( ; it != recMdcTrackCol->end(); it++ )
178 {
179 if ( m_mc && navigator )
180 {
181 McParticleVector particles = navigator->getMcParticles( *it );
182 t_mcTkNum = particles.size();
183 RecMdcTrackVector tracks = navigator->getMdcTracks( particles[0] );
184 // for FIRST parent particle
185 if ( sc != StatusCode::SUCCESS ) return StatusCode::FAILURE;
186 }
187 sc = fillRecTrack( *it, t_mcTkNum, t_recTkNum );
188 t_nTk++;
189 if ( sc != StatusCode::SUCCESS ) return StatusCode::FAILURE;
190 } // end of loop over tracks
191
192 //=============loop over hits==============
193 fillRecHits( *recMdcHitCol );
194
195 if ( m_hist ) { fillEvent(); }
196
197 return StatusCode::SUCCESS;
198}
199
200// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
201
202StatusCode MdcNavigation::finalize() {
203 MsgStream log( msgSvc(), name() );
204 log << MSG::INFO << "in finalize()" << endmsg;
205 std::cout << "nTk == " << t_nTk << std::endl; // yzhang debug
206 delete keepedParticles; // FIXME
207 return StatusCode::SUCCESS;
208}
209
210// routine to calculate momentum using helix parameters of fitted track
211Hep3Vector MdcNavigation::momentum( const RecMdcTrack* trk ) {
212 // double dr = trk->helix(0);
213 double fi0 = trk->helix( 1 );
214 double cpa = trk->helix( 2 );
215 double tanl = trk->helix( 4 );
216
217 double pxy = 0.;
218 if ( cpa != 0 ) pxy = 1 / fabs( cpa );
219
220 double px = pxy * ( -sin( fi0 ) );
221 double py = pxy * cos( fi0 );
222 double pz = pxy * tanl;
223
224 Hep3Vector p;
225 p.set( px, py, pz ); // MeV/c
226 return p;
227}
228
229StatusCode MdcNavigation::fillRecTrack( const RecMdcTrack* tk, int mcTkNum, int recTkNum ) {
230
231 // int mcTkId = m_na_mcTkId;
232 // m_na_nDigi = nDigiTk[mcTkId];//fill # of hit in this mc track
233 m_na_nEvtNoise = nNoise; // fill # of noise hit in this event
234 // m_na_mcTkNum = mcTkNum;// how many mc track are this track include
235 m_na_recTkNum = recTkNum; // how many rec track are this track's mc particle include
236 CLHEP::Hep3Vector rec_mom = momentum( tk );
237 // fill rec momentum
238 m_na_p = rec_mom.mag();
239 m_na_pt = rec_mom.perp();
240 m_na_pz = rec_mom.z();
241 // cout << "MSG::INFO Retrieved McParticles for for RecMdcTrack # " << it->getId()
242 // << " of reconstructed momentum " << rec_mom << " GeV/c" << endl;
243 // five param and it's error matrix
244 m_na_d0 = tk->helix( 0 );
245 m_na_phi0 = tk->helix( 1 );
246 m_na_cpa = tk->helix( 2 );
247 m_na_z0 = tk->helix( 3 );
248 m_na_tanl = tk->helix( 4 );
249
250 if ( m_na_d0 > m_d0Cut )
251 {
252 std::cout << __FILE__ << " " << __LINE__ << " evtNo: " << t_eventNo
253 << " d0:" << std::setw( 5 ) << m_na_d0 << ">" << m_d0Cut << std::endl;
254 setFilterPassed( true );
255 }
256 if ( m_na_z0 > m_z0Cut )
257 {
258 std::cout << __FILE__ << " " << __LINE__ << " evtNo: " << t_eventNo
259 << " z0:" << std::setw( 5 ) << m_na_z0 << ">" << m_z0Cut << std::endl;
260 setFilterPassed( true );
261 }
262 // const CLHEP::HepSymMatrix tkErrM = tk->err();
263 m_na_d0E = tk->err( 0 );
264 m_na_phi0E = tk->err( 2 );
265 m_na_cpaE = tk->err( 5 );
266 m_na_z0E = tk->err( 9 );
267 m_na_tanlE = tk->err( 14 );
268 m_na_q = tk->charge();
269 // fill track about
270 m_na_chi2 = tk->chi2();
271 // m_na_nHit = tk->getNhits();
272 m_na_nDof = tk->ndof();
273 if ( m_na_nDof > 0 ) { m_na_chi2Dof = m_na_chi2 / (float)m_na_nDof; }
274 else { m_na_chi2Dof = 200.; }
275 m_na_chi2Prob = probab( m_na_nDof, m_na_chi2 );
276 m_na_nSt = tk->nster();
277 m_na_fiTerm = tk->getFiTerm();
278
279 if ( tk->stat() == 0 ) { t_trkRecoTk++; }
280 else if ( tk->stat() == 1 ) { t_curlTk++; }
281 else if ( tk->stat() == 2 ) { t_patRecTk++; }
282 else if ( tk->stat() == 3 ) { t_xRecTk++; }
283 m_na_tkStat = tk->stat();
284
285 ////----fill rec Hit
286 // int ihit = 0;
287 // int layerEff[43];
288 // for (int ii=0;ii<43;ii++){
289 // layerEff[ii]=0;
290 // }
291 int noiseHit = 0;
292 int matchHit = 0;
293 int nAct = 0;
294 HitRefVec hl = tk->getVecHits();
295 HitRefVec::iterator hitIt = hl.begin();
296 for ( ; hitIt != hl.end(); ++hitIt )
297 {
298 RecMdcHit* h = *hitIt;
299 if ( !h ) continue;
300 if ( h->getStat() != 0 ) nAct++;
301 // //fill residual
302 // double m_na_lr = h->getFlagLR();
303 // double ddrift = -999;double ddoca = -999;
304 // ddoca = h->getDoca();
305 // if( 0 == m_na_lr ){ ddrift = h->getDriftDistLeft();
306 // }else{ ddrift = h->getDriftDistRight();}
307 // m_na_resid[ihit] = fabs(ddrift) - fabs(ddoca);
308 // if( 0 == m_na_lr ){ m_na_resid[ihit] *= -1.0;}
309 // m_na_driftD[ihit] = ddrift;
310 // m_na_driftT[ihit] = h->getDriftT();
311 // m_na_doca[ihit] = ddoca;
312 // m_na_entra[ihit] = h->getEntra();
313 // m_na_zhit[ihit] = h->getZhit();
314 // m_na_chi2add[ihit] = h->getChisqAdd();
315 // m_na_flaglr[ihit] = h->getFlagLR();
316 // m_na_hitStat[ihit] = h->getStat();
317 // m_na_Tdc[ihit] = h->getTdc();
318 // m_na_Adc[ihit] = h->getAdc();
319 // m_na_act[ihit] = h->getStat();
320
321 int tlayer = MdcID::layer( h->getMdcId() );
322 int twire = MdcID::wire( h->getMdcId() );
323 // m_na_layer[ihit] = tlayer;
324 // m_na_wire[ihit] = twire;
325 // m_na_gwire[ihit] = Constants::nWireBeforeLayer[tlayer] + twire;
326 // if (0==layerEff[tlayer]){
327 // layerEff[tlayer]=1;
328 // g_layerEff->fill(tlayer);
329 // }
330
331 if ( havedigi[tlayer][twire] < 0 ) { ++noiseHit; }
332 // else{
333 // if(havedigi[tlayer][twire] == mcTkId) ++matchHit;
334 // m_na_hitTkId[ihit] = havedigi[tlayer][twire];
335 // }
336 // ++ihit;
337 }
338 m_na_nAct = nAct;
339 m_na_nNoise = noiseHit;
340 m_na_nMatch = matchHit;
341 g_tupleMc->write();
342 //------------------------------------------
343 // fill Rec Track
344 //------------------------------------------
345
346 if ( m_poca )
347 {
348 uint32_t getDigiFlag = 0;
349 getDigiFlag += m_maxMdcDigi;
350 if ( m_dropHot ) getDigiFlag |= MdcRawDataProvider::b_dropHot;
351 if ( m_keepBadTdc ) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
352 if ( m_keepUnmatch ) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
353 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec( getDigiFlag );
354 MdcDigiCol::iterator iter = mdcDigiVec.begin();
355 for ( ; iter != mdcDigiVec.end(); iter++ ) { poca( ( *iter ), tk->helix(), tk->err() ); }
356 }
357
358 return StatusCode::SUCCESS;
359}
360
361StatusCode MdcNavigation::bookNTuple() {
362 MsgStream log( msgSvc(), name() );
363 StatusCode sc = StatusCode::SUCCESS;
364 g_layerEff = histoSvc()->book( "layerEff", "layerEff", 43, -0.5, 42.5 );
365
366 NTuplePtr nt1( ntupleSvc(), "MdcNavigation/rec" );
367 if ( nt1 ) { g_tupleMc = nt1; }
368 else
369 {
370 g_tupleMc = ntupleSvc()->book( "MdcNavigation/rec", CLID_ColumnWiseTuple,
371 "rec and delta with mc truth" );
372 if ( g_tupleMc )
373 {
374 sc = g_tupleMc->addItem( "tkStat", m_na_tkStat );
375 sc = g_tupleMc->addItem( "p", m_na_p );
376 sc = g_tupleMc->addItem( "pt", m_na_pt );
377 sc = g_tupleMc->addItem( "pz", m_na_pz );
378 sc = g_tupleMc->addItem( "d0", m_na_d0 );
379 sc = g_tupleMc->addItem( "phi0", m_na_phi0 );
380 sc = g_tupleMc->addItem( "cpa", m_na_cpa );
381 sc = g_tupleMc->addItem( "z0", m_na_z0 );
382 sc = g_tupleMc->addItem( "tanl", m_na_tanl );
383 sc = g_tupleMc->addItem( "d0E", m_na_d0E );
384 sc = g_tupleMc->addItem( "phi0E", m_na_phi0E );
385 sc = g_tupleMc->addItem( "cpaE", m_na_cpaE );
386 sc = g_tupleMc->addItem( "z0E", m_na_z0E );
387 sc = g_tupleMc->addItem( "tanlE", m_na_tanlE );
388 sc = g_tupleMc->addItem( "q", m_na_q );
389
390 sc = g_tupleMc->addItem( "dP", m_na_dP );
391 sc = g_tupleMc->addItem( "dPt", m_na_dPt );
392 sc = g_tupleMc->addItem( "dPz", m_na_dPz );
393 sc = g_tupleMc->addItem( "dD0", m_na_dD0 );
394 sc = g_tupleMc->addItem( "dPhi0", m_na_dPhi0 );
395 sc = g_tupleMc->addItem( "dCpa", m_na_dCpa );
396 sc = g_tupleMc->addItem( "dZ0", m_na_dZ0 );
397 sc = g_tupleMc->addItem( "dTanl", m_na_dTanl );
398
399 sc = g_tupleMc->addItem( "d0Res", m_na_d0Res );
400 sc = g_tupleMc->addItem( "phiRes", m_na_phi0Res );
401 sc = g_tupleMc->addItem( "z0Res", m_na_z0Res );
402 sc = g_tupleMc->addItem( "tanlRes", m_na_tanlRes );
403 sc = g_tupleMc->addItem( "cpaRes", m_na_cpaRes );
404
405 sc = g_tupleMc->addItem( "recTkNum", m_na_recTkNum );
406 sc = g_tupleMc->addItem( "mcTkId", m_na_mcTkId );
407 sc = g_tupleMc->addItem( "mcTkNum", m_na_mcTkNum );
408 sc = g_tupleMc->addItem( "evtNo", m_na_evtNo );
409 sc = g_tupleMc->addItem( "nSt", m_na_nSt );
410 sc = g_tupleMc->addItem( "nDof", m_na_nDof );
411 sc = g_tupleMc->addItem( "chi2", m_na_chi2 );
412 sc = g_tupleMc->addItem( "chi2Dof", m_na_chi2Dof );
413 sc = g_tupleMc->addItem( "chi2Porb", m_na_chi2Prob );
414 sc = g_tupleMc->addItem( "fiTerm", m_na_fiTerm );
415 sc = g_tupleMc->addItem( "nMatch", m_na_nMatch );
416 sc = g_tupleMc->addItem( "nAct", m_na_nAct );
417 sc = g_tupleMc->addItem( "nTkNoise", m_na_nNoise );
418 sc = g_tupleMc->addItem( "nEvtNoise", m_na_nEvtNoise );
419 sc = g_tupleMc->addItem( "nHit", m_na_nHit, 0, 10000 );
420 sc = g_tupleMc->addItem( "nDigi", m_na_nDigi, 0, 10000 );
421 sc = g_tupleMc->addIndexedItem( "resid", m_na_nHit, m_na_resid );
422 sc = g_tupleMc->addIndexedItem( "driftD", m_na_nHit, m_na_driftD );
423 sc = g_tupleMc->addIndexedItem( "driftT", m_na_nHit, m_na_driftT );
424 sc = g_tupleMc->addIndexedItem( "doca", m_na_nHit, m_na_doca );
425 sc = g_tupleMc->addIndexedItem( "entra", m_na_nHit, m_na_entra );
426 sc = g_tupleMc->addIndexedItem( "zhit", m_na_nHit, m_na_zhit );
427 sc = g_tupleMc->addIndexedItem( "chi2add", m_na_nHit, m_na_chi2add );
428 sc = g_tupleMc->addIndexedItem( "flaglr", m_na_nHit, m_na_flaglr );
429 sc = g_tupleMc->addIndexedItem( "hitStat", m_na_nHit, m_na_hitStat );
430 sc = g_tupleMc->addIndexedItem( "tdc", m_na_nHit, m_na_Tdc );
431 sc = g_tupleMc->addIndexedItem( "adc", m_na_nHit, m_na_Adc );
432 sc = g_tupleMc->addIndexedItem( "act", m_na_nHit, m_na_act );
433 sc = g_tupleMc->addIndexedItem( "layer", m_na_nHit, m_na_layer );
434 sc = g_tupleMc->addIndexedItem( "wire", m_na_nHit, m_na_wire );
435 sc = g_tupleMc->addIndexedItem( "gwire", m_na_nHit, m_na_gwire );
436 sc = g_tupleMc->addIndexedItem( "hitTkId", m_na_nHit, m_na_hitTkId );
437 sc = g_tupleMc->addIndexedItem( "digiTkId", m_na_nDigi, m_na_digiTkId );
438 sc = g_tupleMc->addIndexedItem( "mclayer", m_na_nDigi, m_na_digiLayer );
439 }
440 else
441 {
442 log << MSG::ERROR << " Cannot book N-tuple:" << long( g_tupleMc ) << endmsg;
443 return StatusCode::FAILURE;
444 }
445 } // end book of g_tupleMc
446 NTuplePtr nt3( ntupleSvc(), "MdcNavigation/evt" );
447 if ( nt3 ) { g_tupleEvt = nt3; }
448 else
449 {
450 g_tupleEvt = ntupleSvc()->book( "MdcNavigation/evt", CLID_ColumnWiseTuple, "event" );
451 if ( g_tupleEvt )
452 {
453 sc = g_tupleEvt->addItem( "nTdsTk", m_na_t3recTk );
454 sc = g_tupleEvt->addItem( "trkReco", m_na_t3TrkReco );
455 sc = g_tupleEvt->addItem( "curlFinder", m_na_t3Curl );
456 sc = g_tupleEvt->addItem( "patRec", m_na_t3PatRec );
457 sc = g_tupleEvt->addItem( "xRec", m_na_t3XRec );
458 sc = g_tupleEvt->addItem( "mcTkNum", m_na_t3mcTk );
459 sc = g_tupleEvt->addItem( "evtNo", m_na_t3evtNo );
460 sc = g_tupleEvt->addItem( "t0", m_na_t3t0 );
461 sc = g_tupleEvt->addItem( "t0Truth", m_na_t3t0Truth );
462 sc = g_tupleEvt->addItem( "t0Stat", m_na_t3t0Stat );
463 sc = g_tupleEvt->addItem( "runNo", m_na_t3runNo );
464 sc = g_tupleEvt->addItem( "nDigi", m_na_t3nDigi, 0, 10000 );
465 sc = g_tupleEvt->addIndexedItem( "layer", m_na_t3nDigi, m_na_t3layer );
466 sc = g_tupleEvt->addIndexedItem( "wire", m_na_t3nDigi, m_na_t3wire );
467 sc = g_tupleEvt->addIndexedItem( "gwire", m_na_t3nDigi, m_na_t3gwire );
468 sc = g_tupleEvt->addIndexedItem( "rt", m_na_t3nDigi, m_na_t3rt );
469 sc = g_tupleEvt->addIndexedItem( "rtNot0", m_na_t3nDigi, m_na_t3rtNot0 );
470 sc = g_tupleEvt->addIndexedItem( "rc", m_na_t3nDigi, m_na_t3rc );
471 sc = g_tupleEvt->addIndexedItem( "phi", m_na_t3nDigi, m_na_t3phi );
472 sc = g_tupleEvt->addIndexedItem( "ovfl", m_na_t3nDigi, m_na_t3ovfl );
473 sc = g_tupleEvt->addIndexedItem( "tNum", m_na_t3nDigi, m_na_t3tNum );
474 }
475 }
476}
477
478StatusCode MdcNavigation::fillInit() {
479 MsgStream log( msgSvc(), name() );
480 StatusCode sc = StatusCode::SUCCESS;
481
482 // int mcTkId = m_na_mcTkId;
483 // m_na_nDigi = nDigiTk[mcTkId];//fill # of hit in this mc track
484 t_trkRecoTk = 0;
485 t_curlTk = 0;
486 t_patRecTk = 0;
487 t_xRecTk = 0;
488 //-------------Get event header
489 SmartDataPtr<Event::EventHeader> evtHead( eventSvc(), "/Event/EventHeader" );
490 if ( evtHead )
491 {
492 t_runNo = evtHead->runNumber();
493 t_eventNo = evtHead->eventNumber();
494 }
495 else
496 {
497 log << MSG::WARNING << "Could not retrieve event header" << endmsg;
498 return StatusCode::FAILURE;
499 }
500 if ( m_debug ) std::cout << "evtNo:" << t_eventNo << std::endl;
501
502 //------------Get event time
503 t_t0 = -1;
504 t_t0Stat = -1;
505 SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
506
507 if ( !aevtimeCol || aevtimeCol->size() == 0 )
508 { log << MSG::WARNING << "Could not find RecEsTimeCol" << endmsg; }
509 else
510 {
511 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
512 t_t0 = ( *iter_evt )->getTest();
513 t_t0Stat = ( *iter_evt )->getStat();
514 }
515
516 //------------- Get McDigi collection
517 uint32_t getDigiFlag = 0;
518 getDigiFlag += m_maxMdcDigi;
519 if ( m_dropHot ) getDigiFlag |= MdcRawDataProvider::b_dropHot;
520 if ( m_keepBadTdc ) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
521 if ( m_keepUnmatch ) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
522 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec( getDigiFlag );
523 if ( ( mdcDigiVec.size() == 0 ) )
524 {
525 log << MSG::WARNING << t_eventNo << "No digi or could not find event in MdcDigiVec"
526 << endmsg;
527 }
528
529 for ( int ii = 0; ii < 43; ii++ )
530 {
531 for ( int jj = 0; jj < 288; jj++ )
532 {
533 havedigi[ii][jj] = -99; // no hit or noise
534 }
535 }
536
537 for ( int imc = 0; imc < 100; imc++ )
538 {
539 nDigiTk[imc] = 0;
540 digiLayer[imc] = 0;
541 }
542
543 for ( int ii = 0; ii < 43; ii++ )
544 for ( int jj = 0; jj < 288; jj++ ) multiTdcCount[ii][jj] = 0;
545 MdcDigiCol::iterator iter = mdcDigiVec.begin();
546 for ( ; iter != mdcDigiVec.end(); iter++ )
547 {
548 double t = RawDataUtil::MdcTime( ( *iter )->getTimeChannel() );
549 double c = ( *iter )->getChargeChannel();
550 int l = MdcID::layer( ( *iter )->identify() );
551 int w = MdcID::wire( ( *iter )->identify() );
552 multiTdcCount[l][w]++;
553 }
554
555 int t_i = 0;
556 iter = mdcDigiVec.begin();
557 for ( ; iter != mdcDigiVec.end(); iter++, ++t_i )
558 {
559 double t = RawDataUtil::MdcTime( ( *iter )->getTimeChannel() );
560 double c = ( *iter )->getChargeChannel();
561 int l = MdcID::layer( ( *iter )->identify() );
562 int w = MdcID::wire( ( *iter )->identify() );
563 if ( !m_rawData )
564 {
565 int tkid = ( *iter )->getTrackIndex();
566 havedigi[l][w] = tkid;
567 if ( g_tupleMc )
568 {
569 // m_na_digiTkId[t_i] = tkid;
570 }
571 if ( tkid > -1 ) { ++nDigiTk[tkid]; }
572 else { nNoise++; }
573 }
574 // if (g_tupleMc) m_na_digiLayer[t_i] = l;
575 }
576 return sc;
577}
578
579StatusCode MdcNavigation::skipMcParticle( const McParticle* it, int nKindKeeped, int* pid ) {
580
581 MsgStream log( msgSvc(), name() );
582
583 int pdg_code = ( *it ).particleProperty();
584 if ( fabs( pdg_code ) >= 8 )
585 {
586 const HepPDT::ParticleData* particles = m_particleTable->particle( abs( pdg_code ) );
587 if ( !particles )
588 { log << MSG::WARNING << "Exotic particle found with PDG code " << pdg_code << endmsg; }
589 else
590 {
591 // skip neutrals
592 if ( particles->charge() == 0 )
593 {
594 log << MSG::INFO << "Skip charge == 0 mc particle " << pdg_code << endmsg;
595 return StatusCode::FAILURE;
596 }
597 }
598 }
599
600 int mcTkId = ( *it ).trackIndex();
601 int nMcHit = 0;
602 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol( eventSvc(), "/Event/MC/MdcMcHitCol" );
603 if ( !mcMdcMcHitCol ) { log << MSG::INFO << "Could not find MdcMcHit" << endmsg; }
604 else
605 {
606 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
607 for ( ; iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ )
608 {
609 int iMcTk = ( *iter_mchit )->getTrackIndex();
610 if ( mcTkId == iMcTk ) nMcHit++;
611 }
612 }
613 if ( nMcHit < m_nMcHit ) return StatusCode::FAILURE; // 5 for default
614
615 bool keeped = false;
616 for ( int i = 0; i < nKindKeeped; i++ )
617 {
618 if ( abs( pdg_code ) == pid[i] ) keeped = true;
619 }
620
621 if ( !keeped ) return StatusCode::FAILURE;
622
623 return StatusCode::SUCCESS;
624} // end of skipMcParticle()
625
626StatusCode MdcNavigation::fillEvent() {
627 if ( !g_tupleEvt ) return StatusCode::FAILURE;
628 uint32_t getDigiFlag = 0;
629 getDigiFlag += m_maxMdcDigi;
630 if ( m_dropHot ) getDigiFlag |= MdcRawDataProvider::b_dropHot;
631 if ( m_keepBadTdc ) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
632 if ( m_keepUnmatch ) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
633 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec( getDigiFlag );
634
635 MdcDigiCol::iterator iter = mdcDigiVec.begin();
636 int t_i = 0;
637 for ( ; iter != mdcDigiVec.end(); iter++, ++t_i )
638 {
639 double t = RawDataUtil::MdcTime( ( *iter )->getTimeChannel() );
640 double c = ( *iter )->getChargeChannel();
641 int l = MdcID::layer( ( *iter )->identify() );
642 int w = MdcID::wire( ( *iter )->identify() );
643 if ( g_tupleEvt )
644 {
645 m_na_t3layer[t_i] = l;
646 m_na_t3wire[t_i] = w;
647 m_na_t3gwire[t_i] = Constants::nWireBeforeLayer[l] + w;
648 m_na_t3rt[t_i] = t;
649 m_na_t3rtNot0[t_i] = t - t_t0;
650 m_na_t3rc[t_i] = c;
651 const MdcDigi* digi = const_cast<const MdcDigi*>( *iter );
652 m_na_t3ovfl[t_i] = ( *iter )->getOverflow();
653 m_na_t3tNum[t_i] = ( ( *iter )->getOverflow() & 0xF0 ) >> 4;
654 }
655 }
656 if ( g_tupleEvt ) m_na_t3nDigi = t_i;
657
658 m_na_t3TrkReco = t_trkRecoTk;
659 m_na_t3Curl = t_curlTk;
660 m_na_t3PatRec = t_patRecTk;
661 m_na_t3XRec = t_xRecTk;
662
663 m_na_t3t0 = t_t0;
664 m_na_t3t0Stat = t_t0Stat;
665
666 m_na_t3recTk = t_recTkNum;
667 m_na_t3mcTk = t_mcTkNum;
668
669 m_na_t3evtNo = t_eventNo;
670 m_na_t3runNo = t_runNo;
671 g_tupleEvt->write();
672
673 return StatusCode::SUCCESS;
674}
675
676double MdcNavigation::poca( const MdcDigi* aDigi, const HepVector helixPar,
677 const HepSymMatrix errMat ) {
678 int lay = MdcID::layer( aDigi->identify() );
679 int wire = MdcID::wire( aDigi->identify() );
680 //======from five track parameter to calculate the exact position=====//
681 double ALPHA_loc, rho, pt, kappa, phiIn;
682 int charge;
683 double Bz = m_pIMF->getReferField() * 1000.;
684 ALPHA_loc = 333.567 * Bz; // magnetic field const
685 charge = ( kappa >= 0 ) ? 1 : -1;
686 rho = ALPHA_loc / kappa;
687 pt = fabs( 1.0 / kappa );
688 HepVector helix( helix );
689 helix[0] = -helixPar[0]; // cm
690 helix[1] = helixPar[1] + pi / 2;
691 helix[2] = -1. / helixPar[2];
692 helix[3] = helixPar[3]; // cm
693 helix[4] = helixPar[4];
694 HelixTraj* m_helixTraj;
695 MdcSagTraj* m_wireTraj_I;
696 MdcSWire* m_mdcSWire_I;
697 TrkPoca* m_trkPoca;
698 TrkPoca* m_trkPoca_I;
699 TrkPoca* m_trkPoca_O;
700 MdcSagTraj* m_wireTraj_O;
701 MdcSWire* m_mdcSWire_O;
702 m_helixTraj = new HelixTraj( helix, errMat );
703 //---------------wire in and wire out of this layer---------------------
704 const MdcLayer* layPtr = m_gm->Layer( lay );
705 double fltLenIn = layPtr->rMid();
706 phiIn = helix[1] + helix[2] * fltLenIn; // phi0 + omega * L
707 BesAngle tmp( phiIn - layPtr->phiEPOffset() );
708 int wlow = (int)floor( layPtr->nWires() * tmp.rad() / twopi );
709 int wbig = (int)ceil( layPtr->nWires() * tmp.rad() / twopi );
710 if ( tmp == 0 )
711 {
712 wlow = -1;
713 wbig = 1;
714 }
715 int wireIn, wireOut;
716 wireIn = wbig;
717 wireOut = wlow;
718 std::cout << " in MdcNavigation lay/4 " << lay / 4 << " phi " << tmp << " wire " << wireIn
719 << " " << wireOut << std::endl;
720}
721
722StatusCode MdcNavigation::fillRecHits( RecMdcHitCol& hitCol ) {
723 int ihit = 0;
724 RecMdcHitCol::iterator itHit = hitCol.begin();
725 for ( ; itHit != hitCol.end(); itHit++ )
726 {
727 RecMdcHit* h = *itHit;
728 if ( !h ) continue;
729 double m_na_lr = h->getFlagLR();
730 double ddrift = -999;
731 double ddoca = -999;
732 ddoca = h->getDoca();
733 if ( 0 == m_na_lr ) { ddrift = h->getDriftDistLeft(); }
734 else { ddrift = h->getDriftDistRight(); }
735 m_na_resid[ihit] = fabs( ddrift ) - fabs( ddoca );
736 if ( 0 == m_na_lr ) { m_na_resid[ihit] *= -1.0; }
737 m_na_driftD[ihit] = ddrift;
738 m_na_driftT[ihit] = h->getDriftT();
739 m_na_doca[ihit] = ddoca;
740 m_na_entra[ihit] = h->getEntra();
741 m_na_zhit[ihit] = h->getZhit();
742 m_na_chi2add[ihit] = h->getChisqAdd();
743 m_na_flaglr[ihit] = h->getFlagLR();
744 m_na_hitStat[ihit] = h->getStat();
745 m_na_Tdc[ihit] = h->getTdc();
746 m_na_Adc[ihit] = h->getAdc();
747 m_na_act[ihit] = h->getStat();
748 int tlayer = MdcID::layer( h->getMdcId() );
749 int twire = MdcID::wire( h->getMdcId() );
750 m_na_layer[ihit] = tlayer;
751 m_na_wire[ihit] = twire;
752 m_na_gwire[ihit] = Constants::nWireBeforeLayer[tlayer] + twire;
753 ++ihit;
754 } // end of loop over hits
755 m_na_nHit = ihit;
756 return StatusCode::SUCCESS;
757}
758
759double MdcNavigation::probab( const int& ndof, const double& chisq ) {
760
761 // constants
762 static double srtopi = 2.0 / sqrt( 2.0 * M_PI );
763 static double upl = 100.0;
764
765 double prob = 0.0;
766 if ( ndof <= 0 ) { return prob; }
767 if ( chisq < 0.0 ) { return prob; }
768 if ( ndof <= 60 )
769 {
770 // full treatment
771 if ( chisq > upl ) { return prob; }
772 double sum = exp( -0.5 * chisq );
773 double term = sum;
774
775 int m = ndof / 2;
776 if ( 2 * m == ndof )
777 {
778 if ( m == 1 ) { return sum; }
779 for ( int i = 2; i <= m; i++ )
780 {
781 term = 0.5 * term * chisq / ( i - 1 );
782 sum += term;
783 } //(int i=2; i<=m)
784 return sum;
785 // even
786 }
787 else
788 {
789 // odd
790 double srty = sqrt( chisq );
791 double temp = srty / M_SQRT2;
792 prob = erfc( temp );
793 if ( ndof == 1 ) { return prob; }
794 if ( ndof == 3 ) { return ( srtopi * srty * sum + prob ); }
795 m = m - 1;
796 for ( int i = 1; i <= m; i++ )
797 {
798 term = term * chisq / ( 2 * i + 1 );
799 sum += term;
800 } //(int i=1; i<=m; i++)
801 return ( srtopi * srty * sum + prob );
802
803 } //(2*m==ndof)
804 }
805 else
806 {
807 // asymtotic Gaussian approx
808 double srty = sqrt( chisq ) - sqrt( ndof - 0.5 );
809 if ( srty < 12.0 ) { prob = 0.5 * erfc( srty ); };
810
811 } // ndof<30
812
813 return prob;
814} // endof probab
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition BesBdkRc/src/fortran/HepEvt.h:23
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
RecMdcTrackVector
std::vector< const RecMdcTrack * > RecMdcTrackVector
Definition Event/EventNavigator/include/EventNavigator/EventNavigator.h:50
McParticleVector
std::vector< const Event::McParticle * > McParticleVector
Definition Event/EventNavigator/include/EventNavigator/EventNavigator.h:49
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
pi
double pi
Definition EvtPhokharaDef.hh:17
w
double w
Definition EvtPhokharaDef.hh:28
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
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
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcRecEvent/RecMdcHit.h:92
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcRecEvent/RecMdcTrack.h:22
epsilon
const double epsilon
Definition MdcNavigation.cxx:39
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
histoSvc
IHistogramSvc * histoSvc()
Definition ServiceAccessor.h:23
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
M_PI
#define M_PI
Definition TConstant.h:4
DstMdcTrack::err
const HepSymMatrix err() const
Definition DstMdcTrack.cxx:149
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
MdcDetector::instance
static MdcDetector * instance()
Definition MdcDetector.cxx:18
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
MdcNavigation::initialize
StatusCode initialize()
Definition MdcNavigation.cxx:62
MdcNavigation::execute
StatusCode execute()
Definition MdcNavigation.cxx:122
MdcNavigation::finalize
StatusCode finalize()
Definition MdcNavigation.cxx:202
MdcNavigation::beginRun
StatusCode beginRun()
Definition MdcNavigation.cxx:111
MdcNavigation::MdcNavigation
MdcNavigation(const std::string &name, ISvcLocator *pSvcLocator)
Definition MdcNavigation.cxx:41
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition Event/RawEvent/include/RawEvent/RawDataUtil.h:7
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
t
int t()
Definition t.c:1