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