BOSS 7.1.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 10 of file DQAPi2p2.h.

Constructor & Destructor Documentation

◆ DQAPi2p2()

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

Definition at line 61 of file DQAPi2p2.cxx.

61 :
62 Algorithm(name, pSvcLocator) {
63
64 //Declare the properties
65 declareProperty("saventuple", m_saventuple=false);
66
67 declareProperty("Vr0cut", m_vr0cut=1.0);
68 declareProperty("Vz0cut", m_vz0cut=5.0);
69 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
70 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
71 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
72 declareProperty("GammaAngleCut", m_gammaAngleCut=20.0);
73 declareProperty("Test4C", m_test4C = 1);
74 declareProperty("Test5C", m_test5C = 1);
75 declareProperty("CheckDedx", m_checkDedx = 1);
76 declareProperty("CheckTof", m_checkTof = 1);
77}

Member Function Documentation

◆ execute()

StatusCode DQAPi2p2::execute ( )

sigma;

sigma;

sigma;

sigma;

sigma;

Definition at line 176 of file DQAPi2p2.cxx.

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

◆ finalize()

StatusCode DQAPi2p2::finalize ( )

Definition at line 1068 of file DQAPi2p2.cxx.

1068 {
1069 cout<<"total number: "<<Ncut0<<endl;
1070 cout<<"nGood==2, nCharge==0: "<<Ncut1<<endl;
1071 cout<<"nGam>=2: "<<Ncut2<<endl;
1072 cout<<"Pass Pid: "<<Ncut3<<endl;
1073 cout<<"Pass 4C: "<<Ncut4<<endl;
1074 cout<<"Pass 5C: "<<Ncut5<<endl;
1075 cout<<"J/psi->pi+ pi- p andti-p: "<<Ncut6<<endl;
1076 MsgStream log(msgSvc(), name());
1077 log << MSG::INFO << "in finalize()" << endmsg;
1078 return StatusCode::SUCCESS;
1079}

◆ initialize()

StatusCode DQAPi2p2::initialize ( )

Definition at line 80 of file DQAPi2p2.cxx.

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

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