BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQAPi2p2 Class Reference

#include <DQAPi2p2.h>

Inheritance diagram for DQAPi2p2:

Public Member Functions

 DQAPi2p2 (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Detailed Description

Definition at line 7 of file DQAPi2p2.h.

Constructor & Destructor Documentation

◆ DQAPi2p2()

DQAPi2p2::DQAPi2p2 ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 51 of file DQAPi2p2.cxx.

52 : Algorithm( name, pSvcLocator ) {
53
54 // Declare the properties
55 declareProperty( "saventuple", m_saventuple = false );
56
57 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
58 declareProperty( "Vz0cut", m_vz0cut = 5.0 );
59 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
60 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
61 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
62 declareProperty( "GammaAngleCut", m_gammaAngleCut = 20.0 );
63 declareProperty( "Test4C", m_test4C = 1 );
64 declareProperty( "Test5C", m_test5C = 1 );
65 declareProperty( "CheckDedx", m_checkDedx = 1 );
66 declareProperty( "CheckTof", m_checkTof = 1 );
67}

Referenced by DQAPi2p2().

Member Function Documentation

◆ execute()

StatusCode DQAPi2p2::execute ( )

sigma;

sigma;

sigma;

sigma;

sigma;

Definition at line 161 of file DQAPi2p2.cxx.

161 {
162
163 // std::cout << "execute()" << std::endl;
164
165 MsgStream log( msgSvc(), name() );
166 log << MSG::INFO << "in execute()" << endmsg;
167
168 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
169 int runNo = eventHeader->runNumber();
170 int event = eventHeader->eventNumber();
171 log << MSG::DEBUG << "run, evtnum = " << runNo << " , " << event << endmsg;
172 m_nrun = runNo;
173 m_nrec = event;
174
175 setFilterPassed( false );
176
177 if ( ( event % 100000 ) == 0 ) { cout << "event: " << event << endl; }
178 Ncut0++;
179
180 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
181 // log << MSG::INFO << "get event tag OK" << endmsg;
182 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
183 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
184
185 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
186
187 SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
188 double temp_p_pp, temp_p_pm;
189
190 int t_j = 0;
191
192 //
193 // check x0, y0, z0, r0
194 // suggest cut: |z0|<5 && r0<1
195 //
196 Vint iGood, ipip, ipim, ikaonp, ikaonm, iprotonp, iprotonm;
197 iGood.clear();
198 ipip.clear();
199 ipim.clear();
200 ikaonp.clear();
201 ikaonm.clear();
202 iprotonp.clear();
203 iprotonm.clear();
204 Vp4 ppip, ppim;
205 ppip.clear();
206 ppim.clear();
207
208 int nCharge = 0;
209
210 Hep3Vector xorigin( 0, 0, 0 );
211
212 // if (m_reader.isRunNumberValid(runNo))
213 IVertexDbSvc* vtxsvc;
214 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
215 if ( vtxsvc->isVertexValid() )
216 {
217 double* dbv = vtxsvc->PrimaryVertex();
218 double* vv = vtxsvc->SigmaPrimaryVertex();
219 // HepVector dbv = m_reader.PrimaryVertex(runNo);
220 // HepVector vv = m_reader.SigmaPrimaryVertex(runNo);
221 xorigin.setX( dbv[0] );
222 xorigin.setY( dbv[1] );
223 xorigin.setZ( dbv[2] );
224 }
225
226 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
227 {
228 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
229 if ( !( *itTrk )->isMdcTrackValid() ) continue;
230 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
231 double pch = mdcTrk->p();
232 double x0 = mdcTrk->x();
233 double y0 = mdcTrk->y();
234 double z0 = mdcTrk->z();
235 double phi0 = mdcTrk->helix( 1 );
236 double xv = xorigin.x();
237 double yv = xorigin.y();
238 double Rxy = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
239 // m_vx0 = x0;
240 // m_vy0 = y0;
241 // m_vz0 = z0;
242 // m_vr0 = Rxy;
243
244 HepVector a = mdcTrk->helix();
245 HepSymMatrix Ea = mdcTrk->err();
246 HepPoint3D point0( 0., 0., 0. ); // the initial point for MDC recosntruction
247 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
248 VFHelix helixip( point0, a, Ea );
249 helixip.pivot( IP );
250 HepVector vecipa = helixip.a();
251 double Rvxy0 = fabs( vecipa[0] ); // the nearest distance to IP in xy plane
252 double Rvz0 = vecipa[3]; // the nearest distance to IP in z direction
253 double Rvphi0 = vecipa[1];
254
255 if ( fabs( Rvz0 ) >= m_vz0cut ) continue;
256 if ( fabs( Rvxy0 ) >= m_vr0cut ) continue;
257 if ( cos( mdcTrk->theta() ) > 0.93 ) continue;
258 if ( pch > 5 ) continue;
259
260 iGood.push_back( i );
261 nCharge += mdcTrk->charge();
262 if ( t_j < 4 )
263 {
264 if ( ( *itTrk )->isMdcDedxValid() )
265 {
266 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
267
268 m_dedxchi_e[t_j] = dedxTrk->chiE();
269 m_dedxchi_mu[t_j] = dedxTrk->chiMu();
270 m_dedxchi_pi[t_j] = dedxTrk->chiPi();
271 m_dedxchi_kaon[t_j] = dedxTrk->chiK();
272 m_dedxchi_proton[t_j] = dedxTrk->chiP();
273
274 m_dedxngoodhit[t_j] = dedxTrk->numGoodHits();
275 }
276 if ( ( *itTrk )->isTofTrackValid() )
277 {
278 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
279 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
280 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
281 {
282 TofHitStatus* status = new TofHitStatus;
283 status->setStatus( ( *iter_tof )->status() );
284 if ( status->is_cluster() )
285 {
286 double time = ( *iter_tof )->tof();
287 double sigma = 1.1 * ( *iter_tof )->sigma( 0 ) / 1000;
288 double expE_tof = ( *iter_tof )->texpElectron();
289 double expMu_tof = ( *iter_tof )->texpMuon();
290 double expPi_tof = ( *iter_tof )->texpPion();
291 double expK_tof = ( *iter_tof )->texpKaon();
292 double expP_tof = ( *iter_tof )->texpProton();
293
294 if ( sigma != 0 )
295 {
296
297 m_tofchi_e[t_j] = ( time - expE_tof ); /// sigma;
298 m_tofchi_mu[t_j] = ( time - expMu_tof ); /// sigma;
299 m_tofchi_pi[t_j] = ( time - expPi_tof ); /// sigma;
300 m_tofchi_kaon[t_j] = ( time - expK_tof ); /// sigma;
301 m_tofchi_proton[t_j] = ( time - expP_tof ); /// sigma;
302 }
303 }
304 delete status;
305 }
306 }
307
308 t_j++;
309 }
310 }
311
312 //
313 // Finish Good Charged Track Selection
314 //
315 int nGood = iGood.size();
316 m_nGood = nGood;
317 m_nCharge = nCharge;
318 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
319 if ( ( nGood != 4 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
320 Ncut1++;
321
322 for ( int i = 0; i < nGood; i++ )
323 {
324 m_eop[i] = 5.;
325 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
326 if ( !( *itTrk )->isMdcTrackValid() ) continue;
327 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
328 double p = mdcTrk->p();
329 m_eop[i] = 6.;
330 if ( !( *itTrk )->isEmcShowerValid() ) continue;
331 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
332 double eraw = emcTrk->energy();
333 m_eop[i] = eraw / p;
334 }
335
336 Vint iGam;
337 iGam.clear();
338 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
339 {
340 if ( i >= evtRecTrkCol->size() ) break;
341 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
342 if ( !( *itTrk )->isEmcShowerValid() ) continue;
343 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
344 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
345 // find the nearest charged track
346 double dthe = 200.;
347 double dphi = 200.;
348 double dang = 200.;
349 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
350 {
351 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
352 if ( !( *jtTrk )->isExtTrackValid() ) continue;
353 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
354 if ( extTrk->emcVolumeNumber() == -1 ) continue;
355 Hep3Vector extpos = extTrk->emcPosition();
356 // double ctht = extpos.cosTheta(emcpos);
357 double angd = extpos.angle( emcpos );
358 double thed = extpos.theta() - emcpos.theta();
359 double phid = extpos.deltaPhi( emcpos );
360 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
361 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
362 if ( angd < dang )
363 {
364 dang = angd;
365 dthe = thed;
366 dphi = phid;
367 }
368 }
369 if ( dang >= 200 ) continue;
370 double eraw = emcTrk->energy();
371 dthe = dthe * 180 / ( CLHEP::pi );
372 dphi = dphi * 180 / ( CLHEP::pi );
373 dang = dang * 180 / ( CLHEP::pi );
374 if ( eraw < m_energyThreshold ) continue;
375 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
376 if ( fabs( dang ) < m_gammaAngleCut ) continue;
377 double dtime = emcTrk->time();
378
379 if ( dtime < 0 ) continue;
380 if ( dtime > 14 ) continue;
381
382 //
383 // good photon cut will be set here
384 //
385 iGam.push_back( ( *itTrk )->trackId() );
386 }
387
388 //
389 // Finish Good Photon Selection
390 //
391 int nGam = iGam.size();
392 m_nGam = nGam;
393
394 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
395 << endmsg;
396 // if(nGam<1){
397 // return StatusCode::SUCCESS;
398 // }
399 Ncut2++;
400
401 //
402 // Assign 4-momentum to each photon
403 //
404
405 //
406 // Assign 4-momentum to each charged track
407 //
408 Vp4 pCh;
409 pCh.clear();
410 Vint iPid, iCh;
411 iPid.clear();
412 iCh.clear();
413 int npi = 0, npip = 0, npim = 0;
414 int nkaon = 0, nkaonp = 0, nkaonm = 0;
415 int nproton = 0, nprotonp = 0, nprotonm = 0;
416 ParticleID* pid = ParticleID::instance();
417
418 for ( int i = 0; i < nGood; i++ )
419 {
420 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
421 // if(pid) delete pid;
422 pid->init();
423 pid->setMethod( pid->methodProbability() );
424 // pid->setMethod(pid->methodLikelihood()); //for Likelihood Method
425
426 pid->setChiMinCut( 4 );
427 pid->setRecTrack( *itTrk );
428 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
429 // pid->useDedx()
430 pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() ); // seperater
431 // Pion/Kaon
432 // pid->identify(pid->onlyPion());
433 // pid->identify(pid->onlyKaon());
434 pid->calculate();
435 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
436 // m_ptrk_pid = mdcTrk->p();
437 // m_cost_pid = cos(mdcTrk->theta());
438 iCh.push_back( mdcTrk->charge() );
439 if ( !( pid->IsPidInfoValid() ) )
440 {
441 iPid.push_back( 0 );
442 HepLorentzVector ptrk;
443 ptrk.setPx( mdcTrk->px() );
444 ptrk.setPy( mdcTrk->py() );
445 ptrk.setPz( mdcTrk->pz() );
446 double p3 = ptrk.vect().mag();
447 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
448 pCh.push_back( ptrk );
449 }
450 if ( !( pid->IsPidInfoValid() ) ) continue;
451
452 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
453 // if(pid->probPion() < 0.001) continue;
454 // if(pid->pdf(2)<pid->pdf(3)) continue; // for Likelihood Method(0=electron 1=muon
455 // 2=pion 3=kaon 4=proton)
456
457 RecMdcKalTrack* mdcKalTrk =
458 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
459 // RecMdcTrack
460 // if (!(*itTrk)->isMdcKalTrackValid()) continue;
461 // RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon,
462 // pion, kaon and proton;The default setting is pion
463
464 HepLorentzVector ptrk;
465 ptrk.setPx( mdcTrk->px() );
466 ptrk.setPy( mdcTrk->py() );
467 ptrk.setPz( mdcTrk->pz() );
468
469 if ( ( pid->probPion() >= pid->probKaon() ) && ( pid->probPion() >= pid->probProton() ) )
470 {
471 npi = npi + 1;
472 iPid.push_back( 2 );
473 if ( ( *itTrk )->isMdcKalTrackValid() )
474 {
476 ptrk.setPx( mdcKalTrk->px() );
477 ptrk.setPy( mdcKalTrk->py() );
478 ptrk.setPz( mdcKalTrk->pz() );
479 }
480 double p3 = ptrk.vect().mag();
481 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
482 pCh.push_back( ptrk );
483 if ( mdcTrk->charge() > 0 ) npip++;
484 if ( mdcTrk->charge() < 0 ) npim++;
485 }
486 else if ( ( pid->probKaon() >= pid->probPion() ) &&
487 ( pid->probKaon() >= pid->probProton() ) )
488 {
489 nkaon = nkaon + 1;
490 iPid.push_back( 3 );
491 if ( ( *itTrk )->isMdcKalTrackValid() )
492 {
494 ptrk.setPx( mdcKalTrk->px() );
495 ptrk.setPy( mdcKalTrk->py() );
496 ptrk.setPz( mdcKalTrk->pz() );
497 }
498 double p3 = ptrk.vect().mag();
499 ptrk.setE( sqrt( p3 * p3 + mkaon * mkaon ) );
500 pCh.push_back( ptrk );
501 }
502 else if ( ( pid->probProton() >= pid->probPion() ) &&
503 ( pid->probProton() >= pid->probKaon() ) )
504 {
505 nproton = nproton + 1;
506 iPid.push_back( 4 );
507 if ( ( *itTrk )->isMdcKalTrackValid() )
508 {
510 ptrk.setPx( mdcKalTrk->px() );
511 ptrk.setPy( mdcKalTrk->py() );
512 ptrk.setPz( mdcKalTrk->pz() );
513 }
514 double p3 = ptrk.vect().mag();
515 ptrk.setE( sqrt( p3 * p3 + mproton * mproton ) );
516 pCh.push_back( ptrk );
517 if ( mdcTrk->charge() > 0 ) nprotonp++;
518 if ( mdcTrk->charge() < 0 ) nprotonm++;
519 }
520 else
521 {
522 iPid.push_back( 0 );
523 HepLorentzVector ptrk;
524 ptrk.setPx( mdcTrk->px() );
525 ptrk.setPy( mdcTrk->py() );
526 ptrk.setPz( mdcTrk->pz() );
527 double p3 = ptrk.vect().mag();
528 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
529 pCh.push_back( ptrk );
530 // cout<<"errrrrrrrrrrrrrrr"<<endl;
531 }
532 }
533 m_npi = npi;
534 m_nkaon = nkaon;
535 m_nproton = nproton;
536 // int npip = ipip.size();
537 // int npim = ipim.size();
538 // if(npip*npim != 1) return SUCCESS;
539 m_istat_pmiss = 0;
540 m_istat_pbmiss = 0;
541 m_istat_pipmiss = 0;
542 m_istat_pimmiss = 0;
543 for ( int i = 0; i < 4; i++ )
544 {
545 m_ptrk_pmiss[i] = 5;
546 m_ptrk_pbmiss[i] = 5;
547 m_ptrk_pipmiss[i] = 5;
548 m_ptrk_pimmiss[i] = 5;
549 }
550 for ( int i = 0; i < 3; i++ )
551 {
552 m_id_pmiss[i] = 0;
553 m_id_pbmiss[i] = 0;
554 m_id_pipmiss[i] = 0;
555 m_id_pimmiss[i] = 0;
556 }
557 ////////////////////////////////////////////////////////////////////////
558 // classify
559 /////////////////////////////////////////////////////
560 HepLorentzVector ecms( 0.03406837, 0, 0, 3.0971873 );
561 HepLorentzVector pmiss;
562 m_mpmiss = 5.;
563 m_mpbmiss = 5.;
564 m_mpipmiss = 5.;
565 m_mpimmiss = 5.;
566 m_ppmiss = 5.;
567 m_ppbmiss = 5.;
568 m_ppipmiss = 5.;
569 m_ppimmiss = 5.;
570
571 int istat_nhit;
572
573 if ( ( npip == 1 ) && ( npim == 1 ) && ( nprotonm == 1 ) )
574 {
575 pmiss.setPx( 0 );
576 pmiss.setPy( 0 );
577 pmiss.setPz( 0 );
578 pmiss.setE( 0 );
579 istat_nhit = 0;
580
581 int j = 0;
582 for ( int i = 0; i < nGood; i++ )
583 {
584 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
585 // if(!(*itTrk)->isMdcTrackValid()) continue;
586 RecMdcKalTrack* mdcKalTrk =
587 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
588 // RecMdcTrack
589 // RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon,
590 // pion, kaon and proton;The default setting is pion
591 double p;
592 double eraw;
593 if ( ( *itTrk )->isEmcShowerValid() )
594 {
595 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
596 eraw = emcTrk->energy();
597 }
598
599 // pi+ pi- anti-proton and miss proton
600 if ( ( iPid[i] * iCh[i] ) == 2 )
601 {
602 m_index_pmiss[j] = i;
603 pmiss = pmiss + pCh[i];
605 p = mdcKalTrk->p();
606 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
607 j++;
608 }
609 else if ( ( iPid[i] * iCh[i] ) == -2 )
610 {
611 m_index_pmiss[j] = i;
612 pmiss = pmiss + pCh[i];
614 p = mdcKalTrk->p();
615 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
616
617 j++;
618 }
619 else if ( ( iPid[i] * iCh[i] ) == -4 )
620 {
621 m_index_pmiss[j] = i;
622 pmiss = pmiss + pCh[i];
624 p = mdcKalTrk->p();
625 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
626
627 j++;
628 }
629 else
630 {
631 if ( nGood == 4 )
632 {
633 m_index_pmiss[3] = i;
635 p = mdcKalTrk->p();
636 HepLorentzVector ptrk;
637 ptrk.setPx( mdcKalTrk->px() );
638 ptrk.setPy( mdcKalTrk->py() );
639 ptrk.setPz( mdcKalTrk->pz() );
640 double p3 = ptrk.vect().mag();
641 ptrk.setE( sqrt( p3 * p3 + mproton * mproton ) );
642 m_ptrk_pmiss[0] = ptrk.px();
643 m_ptrk_pmiss[1] = ptrk.py();
644 m_ptrk_pmiss[2] = ptrk.pz();
645 m_ptrk_pmiss[3] = ptrk.e();
646 }
647 }
648 } // end loop of good charge particle
649
650 for ( int i = 0; i < nGood; i++ )
651 {
652 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
653 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
654 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
655 if ( iCh[i] < 0 ) continue;
656 // if(nGood==6) protonpTrk = (*(evtRecTrkCol->begin()+iGood[i]))->mdcKalTrack();
657 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
658 for ( int ii = 0; ii < 3; ii++ )
659 {
660 pid->init();
661 pid->setMethod( pid->methodProbability() );
662 pid->setChiMinCut( 4 );
663 pid->setRecTrack( *itTrk );
664 if ( ii == 0 ) { pid->usePidSys( pid->useDedx() ); } // use dedx only
665 else if ( ii == 1 )
666 { pid->usePidSys( pid->useTof1() | pid->useTof2() ); } // use tof only
667 else if ( ii == 2 )
668 {
669 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() );
670 } // use both dedx and tof
671
672 pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() ); // seperater
673 // Pion/Kaon
674 pid->calculate();
675 if ( !( pid->IsPidInfoValid() ) ) continue;
676 if ( ( pid->probProton() >= pid->probPion() ) &&
677 ( pid->probProton() >= pid->probKaon() ) )
678 m_id_pmiss[ii]++;
679 }
680 }
681
682 pmiss = ecms - pmiss;
683 m_mpmiss = pmiss.m();
684 m_ppmiss = pmiss.rho();
685
686 if ( fabs( m_mpmiss - mproton ) < 0.02 && istat_nhit == 0 )
687 {
688 m_istat_pmiss = 1;
689 ( *( evtRecTrkCol->begin() + iGood[m_index_pmiss[3]] ) )->setPartId( 3 );
690 ( *( evtRecTrkCol->begin() + iGood[m_index_pmiss[3]] ) )->setQuality( 0 );
691 }
692
693 Ncut3++;
694 } // end miss proton
695 if ( ( npip == 1 ) && ( npim == 1 ) && ( nprotonp == 1 ) )
696 {
697 pmiss.setPx( 0 );
698 pmiss.setPy( 0 );
699 pmiss.setPz( 0 );
700 pmiss.setE( 0 );
701 istat_nhit = 0;
702
703 int j = 0;
704 for ( int i = 0; i < nGood; i++ )
705 {
706 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
707 // if(!(*itTrk)->isMdcTrackValid()) continue;
708 RecMdcKalTrack* mdcKalTrk =
709 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
710 // RecMdcTrack
711 // RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon,
712 // pion, kaon and proton;The default setting is pion
713 double p;
714 double eraw;
715 if ( ( *itTrk )->isEmcShowerValid() )
716 {
717 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
718 eraw = emcTrk->energy();
719 }
720 // m_eop[i]=eraw/p;
721
722 // pi+ pi- proton and miss anti-proton
723 if ( ( iPid[i] * iCh[i] ) == 2 )
724 {
725 m_index_pbmiss[j] = i;
726 pmiss = pmiss + pCh[i];
728 p = mdcKalTrk->p();
729 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
730
731 j++;
732 }
733 else if ( ( iPid[i] * iCh[i] ) == -2 )
734 {
735 m_index_pbmiss[j] = i;
736
737 pmiss = pmiss + pCh[i];
739 p = mdcKalTrk->p();
740 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
741
742 j++;
743 }
744 else if ( ( iPid[i] * iCh[i] ) == 4 )
745 {
746 m_index_pbmiss[j] = i;
747 pmiss = pmiss + pCh[i];
749 p = mdcKalTrk->p();
750 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
751
752 j++;
753 }
754 else
755 {
756 if ( nGood == 4 )
757 {
758 m_index_pbmiss[3] = i;
760 p = mdcKalTrk->p();
761 HepLorentzVector ptrk;
762 ptrk.setPx( mdcKalTrk->px() );
763 ptrk.setPy( mdcKalTrk->py() );
764 ptrk.setPz( mdcKalTrk->pz() );
765 double p3 = ptrk.vect().mag();
766 ptrk.setE( sqrt( p3 * p3 + mproton * mproton ) );
767 m_ptrk_pbmiss[0] = ptrk.px();
768 m_ptrk_pbmiss[1] = ptrk.py();
769 m_ptrk_pbmiss[2] = ptrk.pz();
770 m_ptrk_pbmiss[3] = ptrk.e();
771 }
772 }
773 } // end loop of good charge particle
774
775 for ( int i = 0; i < nGood; i++ )
776 {
777 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
778 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
779 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
780 if ( iCh[i] > 0 ) continue;
781 // if(nGood==6) protonpTrk = (*(evtRecTrkCol->begin()+iGood[i]))->mdcKalTrack();
782 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
783 for ( int ii = 0; ii < 3; ii++ )
784 {
785 pid->init();
786 pid->setMethod( pid->methodProbability() );
787 pid->setChiMinCut( 4 );
788 pid->setRecTrack( *itTrk );
789 if ( ii == 0 ) { pid->usePidSys( pid->useDedx() ); } // use dedx only
790 else if ( ii == 1 )
791 { pid->usePidSys( pid->useTof1() | pid->useTof2() ); } // use tof only
792 else if ( ii == 2 )
793 {
794 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() );
795 } // use both dedx and tof
796
797 pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() ); // seperater
798 // Pion/Kaon
799 pid->calculate();
800 if ( !( pid->IsPidInfoValid() ) ) continue;
801 if ( ( pid->probProton() >= pid->probPion() ) &&
802 ( pid->probProton() >= pid->probKaon() ) )
803 m_id_pbmiss[ii]++;
804 }
805 }
806 pmiss = ecms - pmiss;
807 m_mpbmiss = pmiss.m();
808 m_ppbmiss = pmiss.rho();
809 if ( fabs( m_mpbmiss - mproton ) < 0.02 && istat_nhit == 0 )
810 {
811 m_istat_pbmiss = 1;
812 ( *( evtRecTrkCol->begin() + iGood[m_index_pbmiss[3]] ) )->setPartId( 3 );
813 ( *( evtRecTrkCol->begin() + iGood[m_index_pbmiss[3]] ) )->setQuality( 0 );
814 }
815 Ncut4++;
816 } // end miss anti-proton
817 int initial_pip, initial_pim;
818 //***************************************
819 //// test the istat of jpsi mass
820 //****************************************
821
822 if ( ( npim == 1 ) && ( nprotonp == 1 ) && ( nprotonm == 1 ) )
823 {
824 pmiss.setPx( 0 );
825 pmiss.setPy( 0 );
826 pmiss.setPz( 0 );
827 pmiss.setE( 0 );
828 istat_nhit = 0;
829
830 int j = 0;
831 for ( int i = 0; i < nGood; i++ )
832 {
833 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
834 // if(!(*itTrk)->isMdcTrackValid()) continue;
835 RecMdcKalTrack* mdcKalTrk =
836 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
837 // RecMdcTrack
838 // RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon,
839 // pion, kaon and proton;The default setting is pion
840 double p;
841 double eraw;
842 if ( ( *itTrk )->isEmcShowerValid() )
843 {
844 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
845 eraw = emcTrk->energy();
846 }
847 // m_eop[i]=eraw/p;
848
849 // pi+ pi- proton and miss anti-proton
850 if ( ( iPid[i] * iCh[i] ) == 4 )
851 {
852 m_index_pipmiss[j] = i;
853 pmiss = pmiss + pCh[i];
855 p = mdcKalTrk->p();
856 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
857
858 j++;
859 }
860 else if ( ( iPid[i] * iCh[i] ) == -4 )
861 {
862 m_index_pipmiss[j] = i;
863
864 pmiss = pmiss + pCh[i];
866 p = mdcKalTrk->p();
867 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
868
869 j++;
870 }
871 else if ( ( iPid[i] * iCh[i] ) == -2 )
872 {
873 m_index_pipmiss[j] = i;
874
875 pmiss = pmiss + pCh[i];
877 p = mdcKalTrk->p();
878 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
879
880 j++;
881 }
882 else
883 {
884 if ( nGood == 4 )
885 {
886 m_index_pipmiss[3] = i;
887
889 p = mdcKalTrk->p();
890 HepLorentzVector ptrk;
891 ptrk.setPx( mdcKalTrk->px() );
892 ptrk.setPy( mdcKalTrk->py() );
893 ptrk.setPz( mdcKalTrk->pz() );
894 double p3 = ptrk.vect().mag();
895 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
896 m_ptrk_pipmiss[0] = ptrk.px();
897 m_ptrk_pipmiss[1] = ptrk.py();
898 m_ptrk_pipmiss[2] = ptrk.pz();
899 m_ptrk_pipmiss[3] = ptrk.e();
900 }
901 }
902 } // end loop of good charge particle
903
904 for ( int i = 0; i < nGood; i++ )
905 {
906 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
907 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
908 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
909 if ( iCh[i] < 0 ) continue;
910 // if(nGood==6) protonpTrk = (*(evtRecTrkCol->begin()+iGood[i]))->mdcKalTrack();
911 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
912 for ( int ii = 0; ii < 3; ii++ )
913 {
914 pid->init();
915 pid->setMethod( pid->methodProbability() );
916 pid->setChiMinCut( 4 );
917 pid->setRecTrack( *itTrk );
918 if ( ii == 0 ) { pid->usePidSys( pid->useDedx() ); } // use dedx only
919 else if ( ii == 1 )
920 { pid->usePidSys( pid->useTof1() | pid->useTof2() ); } // use tof only
921 else if ( ii == 2 )
922 {
923 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() );
924 } // use both dedx and tof
925
926 pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() ); // seperater
927 // Pion/Kaon
928 pid->calculate();
929 if ( !( pid->IsPidInfoValid() ) ) continue;
930 if ( ( pid->probPion() >= pid->probProton() ) &&
931 ( pid->probPion() >= pid->probKaon() ) )
932 m_id_pipmiss[ii]++;
933 }
934 }
935 pmiss = ecms - pmiss;
936 m_mpipmiss = pmiss.m();
937 m_ppipmiss = pmiss.rho();
938 if ( fabs( m_mpipmiss - mpi ) < 0.05 && istat_nhit == 0 )
939 {
940 m_istat_pipmiss = 1;
941 ( *( evtRecTrkCol->begin() + iGood[m_index_pipmiss[3]] ) )->setPartId( 2 );
942 ( *( evtRecTrkCol->begin() + iGood[m_index_pipmiss[3]] ) )->setQuality( 0 );
943 }
944 Ncut5++;
945 } // end miss pip
946 if ( ( npip == 1 ) && ( nprotonp == 1 ) && ( nprotonm == 1 ) )
947 {
948 pmiss.setPx( 0 );
949 pmiss.setPy( 0 );
950 pmiss.setPz( 0 );
951 pmiss.setE( 0 );
952 istat_nhit = 0;
953
954 int j = 0;
955 for ( int i = 0; i < nGood; i++ )
956 {
957 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
958 // if(!(*itTrk)->isMdcTrackValid()) continue;
959 RecMdcKalTrack* mdcKalTrk =
960 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
961 // RecMdcTrack
962 // RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon,
963 // pion, kaon and proton;The default setting is pion
964 double p;
965 double eraw;
966 if ( ( *itTrk )->isEmcShowerValid() )
967 {
968 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
969 eraw = emcTrk->energy();
970 }
971 // m_eop[i]=eraw/p;
972
973 // pi+ pi- proton and miss anti-proton
974
975 if ( ( iPid[i] * iCh[i] ) == 4 )
976 {
977 m_index_pimmiss[j] = i;
978 pmiss = pmiss + pCh[i];
980 p = mdcKalTrk->p();
981 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
982
983 j++;
984 }
985 else if ( ( iPid[i] * iCh[i] ) == -4 )
986 {
987 m_index_pimmiss[j] = i;
988
989 pmiss = pmiss + pCh[i];
991 p = mdcKalTrk->p();
992 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
993
994 j++;
995 }
996 else if ( ( iPid[i] * iCh[i] ) == 2 )
997 {
998 m_index_pimmiss[j] = i;
999 pmiss = pmiss + pCh[i];
1001 p = mdcKalTrk->p();
1002 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
1003
1004 j++;
1005 }
1006 else
1007 {
1008 if ( nGood == 4 )
1009 {
1010 m_index_pimmiss[3] = i;
1012 p = mdcKalTrk->p();
1013 HepLorentzVector ptrk;
1014 ptrk.setPx( mdcKalTrk->px() );
1015 ptrk.setPy( mdcKalTrk->py() );
1016 ptrk.setPz( mdcKalTrk->pz() );
1017 double p3 = ptrk.vect().mag();
1018 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
1019 m_ptrk_pimmiss[0] = ptrk.px();
1020 m_ptrk_pimmiss[1] = ptrk.py();
1021 m_ptrk_pimmiss[2] = ptrk.pz();
1022 m_ptrk_pimmiss[3] = ptrk.e();
1023 }
1024 }
1025 } // end loop of good charge particle
1026
1027 for ( int i = 0; i < nGood; i++ )
1028 {
1029 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
1030 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
1031 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
1032 if ( iCh[i] > 0 ) continue;
1033 // if(nGood==6) protonpTrk = (*(evtRecTrkCol->begin()+iGood[i]))->mdcKalTrack();
1034 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
1035 for ( int ii = 0; ii < 3; ii++ )
1036 {
1037 pid->init();
1038 pid->setMethod( pid->methodProbability() );
1039 pid->setChiMinCut( 4 );
1040 pid->setRecTrack( *itTrk );
1041 if ( ii == 0 ) { pid->usePidSys( pid->useDedx() ); } // use dedx only
1042 else if ( ii == 1 )
1043 { pid->usePidSys( pid->useTof1() | pid->useTof2() ); } // use tof only
1044 else if ( ii == 2 )
1045 {
1046 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() );
1047 } // use both dedx and tof
1048
1049 pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() ); // seperater
1050 // Pion/Kaon
1051 pid->calculate();
1052 if ( !( pid->IsPidInfoValid() ) ) continue;
1053 if ( ( pid->probPion() >= pid->probProton() ) &&
1054 ( pid->probPion() >= pid->probKaon() ) )
1055 m_id_pimmiss[ii]++;
1056 }
1057 }
1058 pmiss = ecms - pmiss;
1059 m_mpimmiss = pmiss.m();
1060 m_ppimmiss = pmiss.rho();
1061 if ( fabs( m_mpimmiss - mpi ) < 0.05 && istat_nhit == 0 )
1062 {
1063 m_istat_pimmiss = 1;
1064 ( *( evtRecTrkCol->begin() + iGood[m_index_pimmiss[3]] ) )->setPartId( 2 );
1065 ( *( evtRecTrkCol->begin() + iGood[m_index_pimmiss[3]] ) )->setQuality( 0 );
1066 }
1067 Ncut6++;
1068 } // end miss pip
1069
1070 if ( m_istat_pmiss != 1 && m_istat_pbmiss != 1 && m_istat_pipmiss != 1 &&
1071 m_istat_pimmiss != 1 )
1072 { return StatusCode::SUCCESS; }
1073
1074 // Ncut3++;
1075 if ( m_saventuple ) m_tuple0->write().ignore();
1076
1077 setFilterPassed( true );
1078
1079 return StatusCode::SUCCESS;
1080}
HepGeom::Point3D< double > HepPoint3D
int runNo
Definition DQA_TO_DB.cxx:13
const double mkaon
Double_t time
EvtRecTrackCol::iterator EvtRecTrackIterator
double mpi
double pi
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
std::vector< int > Vint
Definition Gam4pikp.cxx:37
const double mproton
Definition PipiJpsi.cxx:49
IMessageSvc * msgSvc()
const HepSymMatrix err() const
const HepVector helix() const
......
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
static ParticleID * instance()
bool IsPidInfoValid() const
void calculate()
void init()
void setStatus(unsigned int status)
const double ecms
Definition inclkstar.cxx:26

◆ finalize()

StatusCode DQAPi2p2::finalize ( )

Definition at line 1083 of file DQAPi2p2.cxx.

1083 {
1084 cout << "total number: " << Ncut0 << endl;
1085 cout << "nGood==2, nCharge==0: " << Ncut1 << endl;
1086 cout << "nGam>=2: " << Ncut2 << endl;
1087 cout << "Pass Pid: " << Ncut3 << endl;
1088 cout << "Pass 4C: " << Ncut4 << endl;
1089 cout << "Pass 5C: " << Ncut5 << endl;
1090 cout << "J/psi->pi+ pi- p andti-p: " << Ncut6 << endl;
1091 MsgStream log( msgSvc(), name() );
1092 log << MSG::INFO << "in finalize()" << endmsg;
1093 return StatusCode::SUCCESS;
1094}

◆ initialize()

StatusCode DQAPi2p2::initialize ( )

Definition at line 70 of file DQAPi2p2.cxx.

70 {
71 MsgStream log( msgSvc(), name() );
72
73 log << MSG::INFO << "in initialize()" << endmsg;
74 Ncut0 = Ncut1 = Ncut2 = Ncut3 = Ncut4 = Ncut5 = Ncut6 = 0;
75
76 StatusCode status;
77
78 NTuplePtr nt0( ntupleSvc(), "DQAFILE/pi2p2" );
79 if ( nt0 ) m_tuple0 = nt0;
80 else
81 {
82 m_tuple0 =
83 ntupleSvc()->book( "DQAFILE/pi2p2", CLID_ColumnWiseTuple, "ks N-Tuple example" );
84 if ( m_tuple0 )
85 {
86 status = m_tuple0->addItem( "nrun", m_nrun );
87 status = m_tuple0->addItem( "nrec", m_nrec );
88
89 status = m_tuple0->addItem( "nGam", m_nGam );
90 status = m_tuple0->addItem( "nGood", m_nGood );
91 status = m_tuple0->addItem( "nCharge", m_nCharge );
92 status = m_tuple0->addItem( "npi", m_npi );
93 status = m_tuple0->addItem( "nkaon", m_nkaon );
94 status = m_tuple0->addItem( "nproton", m_nproton );
95
96 status = m_tuple0->addItem( "dedxchi_e", 4, m_dedxchi_e );
97 status = m_tuple0->addItem( "dedxchi_mu", 4, m_dedxchi_mu );
98 status = m_tuple0->addItem( "dedxchi_pi", 4, m_dedxchi_pi );
99 status = m_tuple0->addItem( "dedxchi_kaon", 4, m_dedxchi_kaon );
100 status = m_tuple0->addItem( "dedxchi_proton", 4, m_dedxchi_proton );
101
102 status = m_tuple0->addItem( "tofchi_e", 4, m_tofchi_e );
103 status = m_tuple0->addItem( "tofchi_mu", 4, m_tofchi_mu );
104 status = m_tuple0->addItem( "tofchi_pi", 4, m_tofchi_pi );
105 status = m_tuple0->addItem( "tofchi_kaon", 4, m_tofchi_kaon );
106 status = m_tuple0->addItem( "tofchi_proton", 4, m_tofchi_proton );
107
108 status = m_tuple0->addItem( "trackfitchi", 4, m_trackfitchi );
109 status = m_tuple0->addItem( "dedxngoodhit", 4, m_dedxngoodhit );
110 status = m_tuple0->addItem( "trackfitndof", 4, m_trackfitndof );
111
112 status = m_tuple0->addItem( "index_pmiss", 4, m_index_pmiss );
113 status = m_tuple0->addItem( "index_pbmiss", 4, m_index_pbmiss );
114 status = m_tuple0->addItem( "index_pipmiss", 4, m_index_pipmiss );
115 status = m_tuple0->addItem( "index_pimmiss", 4, m_index_pimmiss );
116
117 status = m_tuple0->addItem( "istat_pmiss", m_istat_pmiss );
118 status = m_tuple0->addItem( "istat_pbmiss", m_istat_pbmiss );
119 status = m_tuple0->addItem( "istat_pipmiss", m_istat_pipmiss );
120 status = m_tuple0->addItem( "istat_pimmiss", m_istat_pimmiss );
121
122 status = m_tuple0->addItem( "mpmiss", m_mpmiss );
123 status = m_tuple0->addItem( "mpbmiss", m_mpbmiss );
124 status = m_tuple0->addItem( "mpipmiss", m_mpipmiss );
125 status = m_tuple0->addItem( "mpimmiss", m_mpimmiss );
126
127 status = m_tuple0->addItem( "ppmiss", m_ppmiss );
128 status = m_tuple0->addItem( "ppbmiss", m_ppbmiss );
129 status = m_tuple0->addItem( "ppipmiss", m_ppipmiss );
130 status = m_tuple0->addItem( "ppimmiss", m_ppimmiss );
131
132 status = m_tuple0->addItem( "ptrk_pmiss", 4, m_ptrk_pmiss );
133 status = m_tuple0->addItem( "ptrk_pbmiss", 4, m_ptrk_pbmiss );
134 status = m_tuple0->addItem( "ptrk_pipmiss", 4, m_ptrk_pipmiss );
135 status = m_tuple0->addItem( "ptrk_pimmiss", 4, m_ptrk_pimmiss );
136
137 status = m_tuple0->addItem( "id_pmiss", 3, m_id_pmiss );
138 status = m_tuple0->addItem( "id_pbmiss", 3, m_id_pbmiss );
139 status = m_tuple0->addItem( "id_pipmiss", 3, m_id_pipmiss );
140 status = m_tuple0->addItem( "id_pimmiss", 3, m_id_pimmiss );
141
142 status = m_tuple0->addItem( "eop", 4, m_eop );
143 }
144 else
145 {
146 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple0 ) << endmsg;
147 return StatusCode::FAILURE;
148 }
149 }
150
151 StatusCode sc;
152 //
153 //--------end of book--------
154 //
155
156 log << MSG::INFO << "successfully return from initialize()" << endmsg;
157 return StatusCode::SUCCESS;
158}
INTupleSvc * ntupleSvc()

The documentation for this class was generated from the following files: