130 {
131 MsgStream log(
msgSvc(), name() );
132 log << MSG::INFO << "DedxCalibEvent::genNtuple()" << endmsg;
133
134 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
135 if ( !eventHeader )
136 {
137 log << MSG::INFO << "Could not find Event Header" << endmsg;
138 return;
139 }
140 int eventNO = eventHeader->eventNumber();
141 int runNO = eventHeader->runNumber();
142
143
144
146 evt_count++;
148 {
150 evt_count = 0;
151 }
152
153
154
155 log << MSG::INFO << "now begin the event: runNO= " << runNO << " eventNO= " << eventNO
156 << endmsg;
157
158 int runFlag = 0;
159 if ( runNO <
RUN0 ) runFlag = 0;
160 if ( runNO >=
RUN1 && runNO <
RUN2 ) runFlag = 1;
161 if ( runNO >=
RUN2 && runNO <
RUN3 ) runFlag = 2;
162 if ( runNO >=
RUN3 && runNO <
RUN4 ) runFlag = 3;
163 if ( runNO >=
RUN4 && runNO <
RUN5 ) runFlag = 4;
164 if ( runNO >=
RUN5 && runNO <
RUN6 ) runFlag = 5;
165 if ( runNO >=
RUN6 && runNO <
RUN7 ) runFlag = 6;
166 if ( runNO >=
RUN7 ) runFlag = 7;
167
168 double tes = -999.0;
169 int esTimeflag = -1;
170 SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
171 if ( ( !aevtimeCol ) || ( aevtimeCol->size() == 0 ) )
172 {
173 tes = -9999.0;
174 esTimeflag = -1;
175 }
176 else
177 {
178 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
179 for ( ; iter_evt != aevtimeCol->end(); iter_evt++ )
180 {
181 tes = ( *iter_evt )->getTest();
182 esTimeflag = ( *iter_evt )->getStat();
183 }
184 }
185 if ( runFlag == 2 )
186 {
187 if ( tes > 1000 ) return;
188 }
189 else if ( runFlag == 3 )
190 {
191 if ( tes > 700 ) return;
192 }
193 else
194 {
195 if ( tes > 1400 ) return;
196 }
197
198 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), "/Event/EvtRec/EvtRecEvent" );
199 if ( !evtRecEvent )
200 {
201 log << MSG::ERROR << "EvtRecEvent not found" << endmsg;
202 return;
203 }
204
205 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), "/Event/EvtRec/EvtRecTrackCol" );
206 if ( !evtRecTrkCol )
207 {
208 log << MSG::ERROR << "EvtRecTrackCol" << endmsg;
209 return;
210 }
211
213 iGood.clear();
214 int nCharge = 0;
215 double db = 0, dz = 0, pt0 = 0, charge0 = 0;
216 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
217 {
219 RecMdcDedx* it = ( *itTrk )->mdcDedx();
220 if ( it == 0 ) continue;
221 HepVector a;
222
223
224
225 if ( ( *itTrk )->isMdcKalTrackValid() )
226 {
229 {
230 DstMdcKalTrack* dstTrk = trk;
232 }
234 }
235 else if ( ( *itTrk )->isMdcTrackValid() )
236 {
239 }
240 else continue;
241 db = a[0];
242 dz = a[3];
243 pt0 = fabs( 1.0 / a[2] );
244 charge0 = ( a[2] > 0 ) ? 1 : -1;
245
246
247 if ( fabs( dz ) >=
VZ0CUT )
continue;
248 if ( db >=
VR0CUT )
continue;
251 iGood.push_back( it->
trackId() );
252 nCharge += charge0;
253 }
254
255 double poca_x = 0, poca_y = 0, poca_z = 0;
256 float sintheta = 0, costheta = 0, ptrk = 0, ptrk_t = 0, charge = 0, trackId = 0;
257 int Nhit = 0, usedhit = 0, Nhits = 0, Nphlisthit = 0;
258 double dEdx_meas_hit = 0, dEdx_meas = 0, dEdx_meas_esat = 0, dEdx_meas_norun = 0;
259 double dEdx_hit[100] = { 0 }, pathlength_hit[100] = { 0 }, wid_hit[100] = { 0 },
260 layid_hit[100] = { 0 }, dd_in_hit[100] = { 0 }, eangle_hit[100] = { 0 },
261 zhit_hit[100] = { 0 };
262 int trk_recalg = -1;
263 Identifier mdcid;
264
265 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
266 {
268 RecMdcDedx* it = ( *itTrk )->mdcDedx();
269 if ( it == 0 ) continue;
271 for ( unsigned int i = 0; i < iGood.size(); i++ )
272 {
274 }
275 if (
flag ==
false )
continue;
276
277 HepVector a;
278 HepSymMatrix tkErrM;
279 if ( ( *itTrk )->isMdcKalTrackValid() )
280 {
282
285
288 {
289 DstMdcKalTrack* dstTrk = trk;
292 }
293 else
294 {
297 }
298 }
299 else if ( ( *itTrk )->isMdcTrackValid() )
300 {
304
308 }
309 else continue;
310
311 sintheta =
sin( M_PI_2 - atan( a[4] ) );
312 costheta =
cos( M_PI_2 - atan( a[4] ) );
313 ptrk_t = 1.0 / fabs( a[2] );
314 ptrk = ptrk_t * sqrt( 1 + a[4] * a[4] );
315 charge = ( a[2] > 0 ) ? 1 : -1;
316
319
321 trk_recalg = it->
status();
323
325 {
326 if ( runFlag == 3 && ( ptrk > 1.88 || ptrk < 1.80 ) ) continue;
327 if ( runFlag == 4 && ( ptrk > 1.72 || ptrk < 1.45 ) ) continue;
328 if ( runFlag == 5 && ( ptrk > 2.00 || ptrk < 1.70 ) ) continue;
329 if ( runFlag == 6 && ( ptrk > 1.90 || ptrk < 1.00 ) ) continue;
330 if ( runFlag == 7 && ( ptrk > 1.90 || ptrk < 0.90 ) ) continue;
331 if (
abs( costheta ) > 0.9 )
continue;
332
333 if ( Nhit < 20 ) continue;
334 if ( usedhit < 6 ) continue;
335 }
336
337 int layid = 0, localwid = 0, w0id = 0, wid = 0, lr = 0;
338 double p_hit = 0, adc_raw = 0, tdc_raw = 0, zhit = 0, driftd = 0, driftD = 0, driftT = 0,
339 dd_in = 0, dd_ex = 0, eangle = 0;
340 double ph = 0, pathlength = 0, rphi_path = 0;
341 long k = 0;
342
344 DedxHitRefVec::iterator it_gothit = gothits.begin();
345 for ( ; it_gothit != gothits.end(); it_gothit++ )
346 {
347 if ( ( *it_gothit )->isMdcHitValid() )
348 {
349 RecMdcHit* itor = ( *it_gothit )->getMdcHit();
353 w0id =
geosvc->Layer( layid )->Wirst();
354 wid = w0id + localwid;
358
360 if ( lr == 2 ) cout << "lr = " << lr << endl;
363 driftd =
abs( driftD );
366 if ( lr == 0 || lr == 2 )
367 {
368 dd_in = -
abs( dd_in );
369 dd_ex = -
abs( dd_ex );
370 }
371 if ( lr == 1 )
372 {
373 dd_in =
abs( dd_in );
374 dd_ex =
abs( dd_ex );
375 }
378 eangle = eangle /
M_PI;
379 }
380 else if ( ( *it_gothit )->isMdcKalHelixSegValid() )
381 {
382 RecMdcKalHelixSeg* itor = ( *it_gothit )->getMdcKalHelixSeg();
384 p_hit = 1.0 / fabs( a_hit_in( 3 ) ) * sqrt( 1 + a_hit_in( 5 ) * a_hit_in( 5 ) );
385
389 w0id =
geosvc->Layer( layid )->Wirst();
390 wid = w0id + localwid;
394
396 if ( lr == 2 ) cout << "lr = " << lr << endl;
397 driftD = itor->
getDD();
398 driftd =
abs( driftD );
399 driftT = itor->
getDT();
402 if ( lr == -1 || lr == 2 )
403 {
404 dd_in = dd_in;
405 dd_ex = dd_ex;
406 }
407 else if ( lr == 1 )
408 {
409 dd_in = -dd_in;
410 dd_ex = -dd_ex;
411 }
413 eangle = eangle /
M_PI;
414 }
415 else continue;
416
417 pathlength = ( *it_gothit )->getPathLength();
418 rphi_path = pathlength * sintheta;
419 ph = ( *it_gothit )->getDedx();
420 if ( layid > 3 )
421 {
422 dEdx_hit[k] = adc_raw;
423 pathlength_hit[k] = pathlength;
424 wid_hit[k] = wid;
425 layid_hit[k] = layid;
426 dd_in_hit[k] = dd_in;
427 eangle_hit[k] = eangle;
428 zhit_hit[k] = zhit;
429
430 k++;
431 }
432
433
434 m_charge1 = charge;
435 m_t01 = tes;
436 m_driftT = driftT;
437 m_tdc_raw = tdc_raw;
438 m_phraw = adc_raw;
439 m_exraw = ph;
440 m_localwid = localwid;
441 m_wire = wid;
442 m_runNO1 = runNO;
443 m_evtNO1 = eventNO;
444 m_runFlag1 = runFlag;
445 m_doca_in = dd_in;
446 m_doca_ex = dd_ex;
447 m_driftdist = driftD;
448 m_eangle = eangle;
449 m_zhit = zhit;
450 m_costheta1 = costheta;
451 m_pathL = pathlength;
452 m_layer = layid;
453 m_ptrk1 = ptrk;
454 m_ptrk_hit = p_hit;
455 m_trackId1 = trackId;
456 StatusCode status;
458 {
459 if (
cut_wire == m_wire ) status = m_nt2->write();
460 }
461 else status = m_nt2->write();
462 if ( status.isFailure() ) log << MSG::ERROR << "Cannot fill Ntuple n102" << endmsg;
463 }
464
465 Nphlisthit = k;
470
471
472 m_poca_x = poca_x;
473 m_poca_y = poca_y;
474 m_poca_z = poca_z;
475 m_ptrk_t = ptrk_t;
476 m_ptrk = ptrk;
477 m_sintheta = sintheta;
478 m_costheta = costheta;
479 m_charge = charge;
480 m_runNO = runNO;
481 m_runFlag = runFlag;
482 m_evtNO = eventNO;
483 m_t0 = tes;
484 m_trackId = trackId;
485 m_recalg = trk_recalg;
486
487 m_nhit = Nhit;
488 m_nhits = Nhits;
489 m_nphlisthit = Nphlisthit;
490 m_usedhit = usedhit;
491 for ( int j = 0; j < Nphlisthit; j++ )
492 {
493 m_dEdx_hit[j] = dEdx_hit[j];
494 m_pathlength_hit[j] = pathlength_hit[j];
495 m_wid_hit[j] = wid_hit[j];
496 m_layid_hit[j] = layid_hit[j];
497 m_dd_in_hit[j] = dd_in_hit[j];
498 m_eangle_hit[j] = eangle_hit[j];
499 m_zhit_hit[j] = zhit_hit[j];
500 }
501
502
503 m_dEdx_meas = dEdx_meas;
504
505
506
508 m_chidedxe = it->
chiE();
509 m_chidedxmu = it->
chiMu();
510 m_chidedxpi = it->
chiPi();
511 m_chidedxk = it->
chiK();
512 m_chidedxp = it->
chiP();
513 for ( int i = 0; i < 5; i++ )
514 {
518 }
519 StatusCode status;
520 if (
cut_wire < 0 ) status = m_nt1->write();
521 if ( status.isFailure() ) { log << MSG::ERROR << "Cannot fill Ntuple n103" << endmsg; }
522 }
523
524}
static int evt_threshold(0)
EvtRecTrackCol::iterator EvtRecTrackIterator
double sin(const BesAngle a)
double cos(const BesAngle a)
SmartRefVector< RecMdcDedxHit > DedxHitRefVec
const HepVector & getZHelix(const int pid) const
const HepSymMatrix & getFError(const int pid) const
const HepVector & getFHelix(const int pid) const
const HepSymMatrix err() const
const HepVector helix() const
......
static int layer(const Identifier &id)
Values of different levels (failure returns 0).
static int wire(const Identifier &id)
RecMdcKalTrack * getMdcKalTrack(void)
double getDedxNoRun(void)
DedxHitRefVec getVecDedxHits() const
double getPidProb(int pid) const
RecMdcTrack * getMdcTrack(void)
double getSigmaDedx(int pid) const
double getDedxExpect(int pid) const
const int getFlagLR(void) const
const double getZhit(void) const
const double getAdc(void) const
const Identifier getMdcId(void) const
const double getTdc(void) const
const double getDriftDistRight(void) const
const double getDriftT(void) const
const double getEntra(void) const
const double getDriftDistLeft(void) const
const double getDoca(void) const
Identifier getMdcId(void) const
int getFlagLR(void) const
HepVector & getHelixIncl(void)
double getEntra(void) const
double getDocaIncl(void) const
double getDocaExcl(void) const
double getTdc(void) const
double getZhit(void) const
double getAdc(void) const
const HepVector & getZHelix() const
const HepSymMatrix & getFError() const
const HepVector & getFHelix() const