BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
rhopi Class Reference

#include <rhopi.h>

+ Inheritance diagram for rhopi:

Public Member Functions

 rhopi (EvtVector4R pd1, EvtVector4R pd2, EvtVector4R pd3)
 
double F00 (double s)
 
double F10 (double s)
 
double amps1 (double s, int i, int j)
 
double amps ()
 
 rhopi (const std::string &name, ISvcLocator *pSvcLocator)
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 

Detailed Description

Definition at line 60 of file UserDIY.cc.

Constructor & Destructor Documentation

◆ rhopi() [1/2]

rhopi::rhopi ( EvtVector4R  pd1,
EvtVector4R  pd2,
EvtVector4R  pd3 
)
inline

Definition at line 62 of file UserDIY.cc.

62 {
63 _pd[0]=pd1;
64 _pd[1]=pd2;
65 _pd[2]=pd3;
66 }

◆ rhopi() [2/2]

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

Definition at line 58 of file rhopi.cxx.

58 :
59 Algorithm(name, pSvcLocator) {
60
61 m_tuple4 = 0;
62 m_tuple5 = 0;
63
64 for(int i = 0; i < 12; i++) m_pass[i] = 0;
65
66//Declare the properties
67 declareProperty("Vr0cut", m_vr0cut=5.0);
68 declareProperty("Vz0cut", m_vz0cut=15.0);
69 declareProperty("EnergyThreshold", m_energyThreshold=0.03);
70 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
71 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
72}

Member Function Documentation

◆ amps()

double rhopi::amps ( )

Definition at line 104 of file UserDIY.cc.

104 {
105 double temp,s12,s13,s23;
106 s12=(_pd[0]+_pd[1]).mass2();
107 s13=(_pd[0]+_pd[2]).mass2();
108 s23=(_pd[1]+_pd[2]).mass2();
109 temp=amps1(s12,0,1)+amps1(s13,0,2)+amps1(s23,1,2);
110 return temp;
111 }
double amps1(double s, int i, int j)
Definition: UserDIY.cc:87
double double double double * s12
Definition: qcdloop1.h:77
double double double double double * s23
Definition: qcdloop1.h:77

◆ amps1()

double rhopi::amps1 ( double  s,
int  i,
int  j 
)

Definition at line 87 of file UserDIY.cc.

87 {
88 double mrho=0.771,wrho=0.1492,dpro;
89 EvtComplex img(0.0,1.0);
90 dpro=pow(abs(s-mrho*mrho+img*sqrt(s)*wrho),2.);
91 EvtVector4R prho;
92 prho=_pd[i]+_pd[j];
93 EvtHelSys angles(prho,_pd[i]),labAngles;
94 double theta,phi,ct1,st1,phi1,st,ct,temp;
95 theta=angles.getHelAng(1);
96 phi =angles.getHelAng(2);
97 ct1 =labAngles.Angles(prho,1);
98 phi1=labAngles.Angles(prho,2);
99 st=sin(theta);ct=cos(theta);
100 temp=pow(F00(s),2.)*pow(F10(s),2.)*pow(st,2.)/dpro; // *(1+pow(ct1,2.)+pow(st1,2.)*cos(2*(phi1+phi)));
101 return temp;
102 }
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
Double_t phi1
XmlRpcServer s
Definition: HelloServer.cpp:11
double Angles(EvtVector4R, int)
Definition: EvtHelSys.cc:108
double F10(double s)
Definition: UserDIY.cc:81
double F00(double s)
Definition: UserDIY.cc:76

Referenced by amps().

◆ execute()

StatusCode rhopi::execute ( )

Definition at line 202 of file rhopi.cxx.

202 {
203
204 StatusCode sc = StatusCode::SUCCESS;
205
206 MsgStream log(msgSvc(), name());
207// log << MSG::INFO << "in execute()" << endreq;
208
209 m_pass[0] += 1;
210
211 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
212 //SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::Recon::EvtRecEvent);
213 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
214
215// SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::Recon::EvtRecTrackCol);
216 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
217 //
218 // check x0, y0, z0, r0
219 // suggest cut: |z0|<10 && r0<5
220 //
221
222 Vint ipip, ipim, iGood;
223 iGood.clear();
224 ipip.clear();
225 ipim.clear();
226
227 Vp4 ppip, ppim;
228 ppip.clear();
229 ppim.clear();
230
231 int TotCharge = 0;
232 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
233 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
234 if(!(*itTrk)->isMdcTrackValid()) continue;
235 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
236
237 m_pi_vx->Fill(mdcTrk->x());
238 m_pi_vy->Fill(mdcTrk->y());
239 m_pi_vz->Fill(mdcTrk->z());
240
241 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
242 if(mdcTrk->r() >= m_vr0cut) continue;
243 iGood.push_back((*itTrk)->trackId());
244 TotCharge += mdcTrk->charge();
245 }
246 //
247 // Finish Good Charged Track Selection
248 //
249 int nGood = iGood.size();
250 m_ngoodch = nGood;
251
252 Vint iGam;
253 iGam.clear();
254 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
255 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
256 if(!(*itTrk)->isEmcShowerValid()) continue;
257 RecEmcShower *emcTrk = (*itTrk)->emcShower();
258 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
259 // find the nearest charged track
260 double dthe = 200.;
261 double dphi = 200.;
262 double dang = 200.;
263 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
264 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
265 if(!(*jtTrk)->isExtTrackValid()) continue;
266 RecExtTrack *extTrk = (*jtTrk)->extTrack();
267 if(extTrk->emcVolumeNumber() == -1) continue;
268 Hep3Vector extpos = extTrk->emcPosition();
269 // double ctht = extpos.cosTheta(emcpos);
270 double angd = extpos.angle(emcpos);
271 double thed = extpos.theta() - emcpos.theta();
272 double phid = extpos.deltaPhi(emcpos);
273 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
274 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
275 angd = fmod(angd, CLHEP::twopi);
276
277 if(angd < dang){
278 dang = angd;
279 dthe = thed;
280 dphi = phid;
281 }
282 }
283
284 if(dang >= 200.) continue;
285
286 double eraw = emcTrk->energy();
287 dthe = dthe * 180 / (CLHEP::pi);
288 dphi = dphi * 180 / (CLHEP::pi);
289 dang = dang * 180 / (CLHEP::pi);
290 if(eraw < m_energyThreshold) continue;
291// if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
292 //
293 // good photon cut will be set here
294 //
295 iGam.push_back((*itTrk)->trackId());
296 }
297 //
298 // Finish Good Photon Selection
299 //
300 int nGam = iGam.size();
301 m_ngoodneu = nGam;
302 //
303 // Assign 4-momentum to each photon
304 //
305
306 HepLorentzVector ptcharg, ptneu, ptchargp, ptchargm;
307 Vp4 pGam;
308 pGam.clear();
309 for(int i = 0; i < nGam; i++) {
310 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
311 RecEmcShower* emcTrk = (*itTrk)->emcShower();
312 double eraw = emcTrk->energy();
313 double phi = emcTrk->phi();
314 double the = emcTrk->theta();
315 HepLorentzVector ptrk;
316 ptrk.setPx(eraw*sin(the)*cos(phi));
317 ptrk.setPy(eraw*sin(the)*sin(phi));
318 ptrk.setPz(eraw*cos(the));
319 ptrk.setE(eraw);
320 pGam.push_back(ptrk);
321 }
322
323//
324// Dedx Check ================
325//
326 for(int i = 0; i < nGood; i++) {
327 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
328 if(!(*itTrk)->isMdcTrackValid()) continue;
329 if(!(*itTrk)->isMdcDedxValid())continue;
330 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
331 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
332 m_chie_dedx -> Fill(dedxTrk->chiE());
333 m_chimu_dedx -> Fill(dedxTrk->chiMu());
334 m_chipi_dedx -> Fill(dedxTrk->chiPi());
335 m_chik_dedx -> Fill(dedxTrk->chiK());
336 m_chip_dedx -> Fill(dedxTrk->chiP());
337 }
338
339 //
340 // Assign 4-momentum to each charged track
341 //
342
343 int TotQ = 0;
344
346 for(int i = 0; i < nGood; i++) {
347 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
348 // if(pid) delete pid;
349 pid->init();
350 pid->setMethod(pid->methodProbability());
351 pid->setChiMinCut(4);
352 pid->setRecTrack(*itTrk);
353 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID sub-system
354 pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
355 pid->calculate();
356 if(!(pid->IsPidInfoValid())) continue;
357
358 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
359 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
360 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
361
362 TotQ += mdcKalTrk->charge();
363
364// RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
365
366// if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
367 if(pid->probPion() < pid->probKaon()) continue;
369 HepLorentzVector ptrk;
370 ptrk.setPx(mdcKalTrk->px());
371 ptrk.setPy(mdcKalTrk->py());
372 ptrk.setPz(mdcKalTrk->pz());
373 double p3 = ptrk.mag();
374 ptrk.setE(sqrt(p3*p3+mpi*mpi));
375 if(mdcTrk->charge() >0) {
376 ipip.push_back(iGood[i]);
377 ppip.push_back(ptrk);
378 } else {
379 ipim.push_back(iGood[i]);
380 ppim.push_back(ptrk);
381 }
382 }
383
384 int npip = ipip.size();
385 int npim = ipim.size();
386 m_npip = npip;
387 m_npim = npim;
388
389 if(nGood != 2 || TotCharge!=0) return sc;
390 m_pass[1] += 1;
391 if(nGam<2 ) return sc;
392 m_pass[2] += 1;
393
394// if(npip != 1) return sc;
395// if(npip*npim != 1) return sc;
396
397 HepLorentzVector BoostP(0.0,0.0,0.0,ecms);
398 BoostP[0] = 2*sin(0.011)*(ecms/2);
399 Hep3Vector u_Boost = -BoostP.boostVector();
400
401 if(npip == 1) {
402//
403//---- Find the Best Pi0 -----
404//
405 HepLorentzVector p2g, ppi0;
406 double delmpi = 999.;
407 int ig11=-1, ig21=-1;
408
409 for(int i = 0; i < nGam - 1; i++){
410 for(int j = i+1; j < nGam; j++) {
411 HepLorentzVector ppi0 = pGam[i] + pGam[j];
412 if(fabs(ppi0.m()-mpi0c) > delmpi) continue;
413 if(ppi0.m() > 0.2) continue;
414 delmpi = fabs(ppi0.m()-mpi0c);
415 ig11 = i;
416 ig21 = j;
417 }
418 }
419 if(ig11==-1 || ig21== -1) return sc;
420 p2g = pGam[ig11] + pGam[ig21];
421 m_m2graw = p2g.m();
422
423 HepLorentzVector Ppip1 = ppip[0];
424 Ppip1.boost(u_Boost);
425 p2g.boost(u_Boost);
426//
427// ----- Pion Efficiency -----
428//
429//
430//------ Missing pi-
431//
432 HepLorentzVector Pmiss = -Ppip1 - p2g;
433 m_pmiss = Pmiss.rho();
434 double emiss = ecms - Ppip1.e() - p2g.e();
435 m_emiss = emiss;
436 m_mmiss = sqrt(fabs(emiss*emiss - Pmiss.rho()*Pmiss.rho()));
437
438 HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
439 pmisspi.setPx(Pmiss.px());
440 pmisspi.setPy(Pmiss.py());
441 pmisspi.setPz(Pmiss.pz());
442 pmisspi.setE(emiss);
443
444 Prho0 = Ppip1 + pmisspi;
445 Prhop = Ppip1 + p2g;
446 Prhom = pmisspi + p2g;
447 ptot = Ppip1 + pmisspi + p2g;
448
449 m_pmrho0 = Prho0.rho();
450 m_mmrho0 = Prho0.m();
451// m_prho0 = (Ppip1 + Ppim1).rho();
452// m_mrho0 = (Ppip1 + Ppim1).m();
453 m_prhop = Prhop.rho();
454 m_mrhop = Prhop.m();
455// m_prhom = (p2g + Ppim1).rho();
456// m_mrhom = (p2g + Ppim1).m();
457 m_pmrhom = Prhom.rho();
458 m_mmrhom = Prhom.m();
459 m_ppipraw = Ppip1.rho();
460
461 m_tuple5->write();
462 }
463
464 if(npip*npim != 1) return sc;
465 m_pass[3] += 1;
466
467 HepLorentzVector Ppim1 = ppim[0];
468 Ppim1.boost(u_Boost);
469 HepLorentzVector Ppip1 = ppip[0];
470 Ppip1.boost(u_Boost);
471
472 HepLorentzVector Pneumiss = -(Ppim1 + Ppip1);
473 m_pneumiss = Pneumiss.rho();
474 double eneumiss = ecms - Ppim1.e() - Ppip1.e();
475 m_eneumiss = eneumiss;
476 m_mneumiss = sqrt(fabs(eneumiss*eneumiss - Pneumiss.rho()*Pneumiss.rho()));
477
478//
479// Let's play vertex fit here
480//
481 HepPoint3D vx(0., 0., 0.);
482 HepSymMatrix Evx(3, 0);
483 double bx = 1E+6;
484 double by = 1E+6;
485 double bz = 1E+6;
486 Evx[0][0] = bx*bx;
487 Evx[1][1] = by*by;
488 Evx[2][2] = bz*bz;
489
490 VertexParameter vxpar;
491 vxpar.setVx(vx);
492 vxpar.setEvx(Evx);
493
494 VertexFit* vtxfit = VertexFit::instance();
496
497 RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
499 WTrackParameter wpip0(mpi, pipTrk->helix(), pipTrk->err());
500
502 WTrackParameter wkp0(mk, pipTrk->helix(), pipTrk->err());
503
504 RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
506 WTrackParameter wpim0(mpi, pimTrk->helix(), pimTrk->err());
507
509 WTrackParameter wkm0(mk, pimTrk->helix(), pimTrk->err());
510
511
512 vtxfit->init();
513 vtxfit->AddTrack(0, wpip0);
514 vtxfit->AddTrack(1, wpim0);
515 vtxfit->AddVertex(0, vxpar,0, 1);
516 if(!vtxfit->Fit(0)) return sc;
517
518
519
520 m_pass[4] += 1;
521
522 WTrackParameter wpip = vtxfit->wtrk(0);
523 WTrackParameter wpim = vtxfit->wtrk(1);
524
525//
526// Apply Kinematic 4C fit
527//
528
529 double chisq4c = 9999.;
530
531 int ig1 = -1;
532 int ig2 = -1;
533
534 for(int i = 0; i < nGam-1; i++) {
535 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
536 for(int j = i+1; j < nGam; j++) {
537 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
538
539 kmfit->init();
540 kmfit->AddTrack(0, wpip);
541 kmfit->AddTrack(1, wpim);
542 kmfit->AddTrack(2, 0.0, g1Trk);
543 kmfit->AddTrack(3, 0.0, g2Trk);
544 kmfit->AddFourMomentum(0, BoostP);
545
546 if(!kmfit->Fit()) continue;
547 double chi2 = kmfit->chisq();
548 if(chi2 > chisq4c) continue;
549 chisq4c = chi2;
550 ig1 = i;
551 ig2 = j;
552 } // gamma1
553 } // gamma2
554
555
556 if(chisq4c > 500 || ig1 < 0) return sc;
557 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[ig1]))->emcShower();
558 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[ig2]))->emcShower();
559
560 kmfit->init();
561 kmfit->AddTrack(0, wpip);
562 kmfit->AddTrack(1, wpim);
563 kmfit->AddTrack(2, 0.0, g1Trk);
564 kmfit->AddTrack(3, 0.0, g2Trk);
565 kmfit->AddFourMomentum(0, BoostP);
566 if(!kmfit->Fit()) return sc;
567 double chi = kmfit->chisq();
568 m_chisq4c = chi;
569 m_chisq_4c -> Fill(chi);
570
571 HepLorentzVector p2gam = kmfit->pfit(2) + kmfit->pfit(3);
572 m_ppi0=p2gam.rho();
573 m_mpi0=p2gam.m();
574
575 m_pi0_mass -> Fill(p2gam.m());
576
577 m_pass[5] += 1;
578
579//
580// Apply Kinematic 5C fit
581//
582//
583// 5c for J/psi-->pi+pi-pi0
584//
585 kmfit->init();
586 kmfit->AddTrack(0, wpip);
587 kmfit->AddTrack(1, wpim);
588 kmfit->AddTrack(2, 0.0, g1Trk);
589 kmfit->AddTrack(3, 0.0, g2Trk);
590 kmfit->AddResonance(0,mpi0c,2,3);
591 kmfit->AddFourMomentum(1, BoostP);
592 if(!kmfit->Fit()) return sc;
593 double chis = kmfit->chisq();
594
595 m_pass[6] += 1;
596//
597// 5c for J/psi-->K+K-pi0
598//
599 double chisk5c = 9999.;
600
601 vtxfit->init();
602 vtxfit->AddTrack(0, wkp0);
603 vtxfit->AddTrack(1, wkm0);
604 vtxfit->AddVertex(0, vxpar,0, 1);
605
606 if(vtxfit->Fit(0)){
607
608 WTrackParameter wkp = vtxfit->wtrk(0);
609 WTrackParameter wkm = vtxfit->wtrk(1);
610
611 kmfit->init();
612 kmfit->AddTrack(0, wkp);
613 kmfit->AddTrack(1, wkm);
614 kmfit->AddTrack(2, 0.0, g1Trk);
615 kmfit->AddTrack(3, 0.0, g2Trk);
616 kmfit->AddResonance(0,mpi0c,2,3);
617 kmfit->AddFourMomentum(1, BoostP);
618 if(kmfit->Fit()) chisk5c = kmfit->chisq();
619 }
620 if(chisk5c<chis) return sc;
621
622 m_pass[7] += 1;
623 //
624 // 5c Fit again for J/psi-->pi+pi-pi0
625 //
626 kmfit->init();
627 kmfit->AddTrack(0, wpip);
628 kmfit->AddTrack(1, wpim);
629 kmfit->AddTrack(2, 0.0, g1Trk);
630 kmfit->AddTrack(3, 0.0, g2Trk);
631 kmfit->AddResonance(0,mpi0c,2,3);
632 kmfit->AddFourMomentum(1, BoostP);
633 if(!kmfit->Fit(0)) return sc;
634 if(!kmfit->Fit()) return sc;
635
636 m_pass[8] += 1;
637 m_chisq5c = kmfit->chisq();
638
639 m_chisq_5c -> Fill(kmfit->chisq());
640
641 HepLorentzVector ppi1 = kmfit->pfit(0);
642 HepLorentzVector ppi2 = kmfit->pfit(1);
643 HepLorentzVector pgam1 = kmfit->pfit(2);
644 HepLorentzVector pgam2 = kmfit->pfit(3);
645 HepLorentzVector p2gam1 = kmfit->pfit(2) + kmfit->pfit(3);
646 HepLorentzVector p2pi = kmfit->pfit(0) + kmfit->pfit(1);
647 m_ppi0fit = p2gam1.rho();
648 m_mpi0fit = p2gam1.m();
649 m_ppip = ppi1.rho();
650 m_ppim = ppi2.rho();
651 m_p2pi = p2pi.rho();
652 m_m2pi = p2pi.m();
653 m_ppip0 = (ppi1+p2gam1).rho();
654 m_mpip0 = (ppi1+p2gam1).m();
655 m_ppim0 = (ppi2+p2gam1).rho();
656 m_mpim0 = (ppi2+p2gam1).m();
657
658 m_pip_mom -> Fill(ppi1.rho());
659 m_pim_mom -> Fill(ppi2.rho());
660 m_rho0_mass -> Fill(p2pi.m());
661 m_rhop_mass -> Fill((ppi1+p2gam1).m());
662 m_rhom_mass -> Fill((ppi2+p2gam1).m());
663
664 //--------------Angle between two charged pions ---------
665 double theta2pi= ppi1.px()*ppi2.px()+ppi1.py()*ppi2.py()+ppi1.pz()*ppi2.pz();
666 double time2pi = ppi1.rho()*ppi2.rho();
667 if(time2pi == 0.) return sc;
668
669 m_pass[9] += 1;
670
671 theta2pi = theta2pi/time2pi;
672
673 // cout <<"time2pi=====" << time2pi <<endl;
674
675 if(theta2pi <= -1.) theta2pi=180.;
676 else if(theta2pi >= 1.) theta2pi=0.;
677 else if(fabs(theta2pi) < 1.) theta2pi = acos(theta2pi)*180./3.14159;
678 m_theta2pi=theta2pi;
679
680 //----- Theta angle of gamma in pi0 system --------
681 HepLorentzVector pg1inpi0, pg2inpi0;
682
683 // HepLorentzVector betapi0 = p2gam1/(p2gam1.e());
684 double betapi0 = p2gam1.beta();
685 // ==================================
686 Hep3Vector twogam(p2gam.px(),p2gam.py(),p2gam.pz());
687 // CLHEP::boostOf(pg1inpi0,gam1,betapi0);
688 m_g1inpi0the = ((boostOf(pgam1,twogam/p2gam.rho(),betapi0)).theta())*180./3.14159;
689 m_g2inpi0the = ((boostOf(pgam2,twogam/p2gam.rho(),betapi0)).theta())*180./3.14159;
690 /*
691 //
692 // ----- Pion Efficiency -----
693 //
694 HepLorentzVector p2g = pGam[ig1] + pGam[ig2];
695 m_m2graw = p2g.m();
696 p2g.boost(u_Boost);
697 //
698 //------ Missing pi-
699 //
700 HepLorentzVector Pmiss = -Ppip1 - p2g;
701 m_pmiss = Pmiss.rho();
702 double emiss = ecms - Ppip1.e() - p2g.e();
703 m_emiss = emiss;
704 m_mmiss = sqrt(fabs(emiss*emiss - Pmiss.rho()*Pmiss.rho()));
705
706 HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
707 pmisspi.setPx(Pmiss.px());
708 pmisspi.setPy(Pmiss.py());
709 pmisspi.setPz(Pmiss.pz());
710 pmisspi.setE(emiss);
711
712 Prho0 = Ppip1 + pmisspi;
713 Prhop = Ppip1 + p2g;
714 Prhom = pmisspi + p2g;
715 ptot = Ppip1 + pmisspi + p2g;
716
717 m_pmrho0 = Prho0.rho();
718 m_mmrho0 = Prho0.m();
719 m_prho0 = (Ppip1 + Ppim1).rho();
720 m_mrho0 = (Ppip1 + Ppim1).m();
721 m_prhop = Prhop.rho();
722 m_mrhop = Prhop.m();
723 m_prhom = (p2g + Ppim1).rho();
724 m_mrhom = (p2g + Ppim1).m();
725 m_pmrhom = Prhom.rho();
726 m_mmrhom = Prhom.m();
727 */
728 m_tuple4->write();
729
730 return StatusCode::SUCCESS;
731}
curve Fill()
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
const double mk
Definition: Gam4pikp.cxx:48
const double mpi
Definition: Gam4pikp.cxx:47
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
IMessageSvc * msgSvc()
double theta() const
Definition: DstEmcShower.h:38
double phi() const
Definition: DstEmcShower.h:39
double x() const
Definition: DstEmcShower.h:35
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
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 HepVector & helix() const
const double px() const
static void setPidType(PidType pidType)
const double pz() const
const double py() const
const HepSymMatrix & err() const
const int charge() const
const double r() const
Definition: DstMdcTrack.h:64
const int charge() const
Definition: DstMdcTrack.h:53
const double z() const
Definition: DstMdcTrack.h:63
const double y() const
Definition: DstMdcTrack.h:62
const double x() const
Definition: DstMdcTrack.h:61
static KinematicFit * instance()
void AddResonance(int number, double mres, std::vector< int > tlis)
double chisq() const
Definition: KinematicFit.h:150
void AddFourMomentum(int number, HepLorentzVector p4)
HepLorentzVector pfit(int n) const
Definition: KinematicFit.h:154
int useTof2() 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
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
WTrackParameter wtrk(int n) const
Definition: VertexFit.h:79
void init()
Definition: VertexFit.cxx:29
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition: VertexFit.cxx:89
static VertexFit * instance()
Definition: VertexFit.cxx:15
bool Fit()
Definition: VertexFit.cxx:301
void setEvx(const HepSymMatrix &eVx)
void setVx(const HepPoint3D &vx)
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 double mpi0c
Definition: rhopi.cxx:40
const float pi
Definition: vector3.h:133

◆ F00()

double rhopi::F00 ( double  s)

Definition at line 76 of file UserDIY.cc.

76 {
77 double mpi=0.1395;
78 return sqrt(s-4*mpi*mpi)/sqrt(s);
79 }

Referenced by amps1().

◆ F10()

double rhopi::F10 ( double  s)

Definition at line 81 of file UserDIY.cc.

81 {
82 double mpi=0.1395,mpsi=3.096916;
83 double tep=sqrt((mpsi*mpsi-pow(mpi+sqrt(s),2.))*(mpsi*mpsi-pow(mpi-sqrt(s),2.)));
84 return tep;
85 }

Referenced by amps1().

◆ finalize()

StatusCode rhopi::finalize ( )

Definition at line 734 of file rhopi.cxx.

734 {
735
736 MsgStream log(msgSvc(), name());
737 log << MSG::ALWAYS << "Total Entries : " << m_pass[0] << endreq;
738 log << MSG::ALWAYS << "Ncharge=2 Cut : " << m_pass[1] << endreq;
739 log << MSG::ALWAYS << "Ngam<2 cut : " << m_pass[2] << endreq;
740 log << MSG::ALWAYS << "Nch+=1,Nch-=1 cut : " << m_pass[3] << endreq;
741 log << MSG::ALWAYS << "Vertex Fit : " << m_pass[4] << endreq;
742 log << MSG::ALWAYS << "4c-Fit : " << m_pass[5] << endreq;
743 log << MSG::ALWAYS << "5c-Fit : " << m_pass[6] << endreq;
744 log << MSG::ALWAYS << "chi3pi<chikkpi cut : " << m_pass[7] << endreq;
745 log << MSG::ALWAYS << "5c-Fit Again: " << m_pass[8] << endreq;
746 log << MSG::ALWAYS << "Final : " << m_pass[9] << endreq;
747 return StatusCode::SUCCESS;
748}

◆ initialize()

StatusCode rhopi::initialize ( )

Definition at line 75 of file rhopi.cxx.

75 {
76 MsgStream log(msgSvc(), name());
77
78 log << MSG::INFO << "in initialize()" << endmsg;
79
80 StatusCode status;
81
82 status = service("THistSvc", m_thistsvc);
83 if(status.isFailure() ){
84 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
85 return status;
86 }
87 m_pi_vx = new TH1F( "pi_vx", "pi_vx", 100,-1.0,1.0 );
88 status = m_thistsvc->regHist("/VAL/PHY/pi_vx", m_pi_vx);
89 m_pi_vy = new TH1F( "pi_vy", "pi_vy", 100,-1.0,1.0 );
90 status = m_thistsvc->regHist("/VAL/PHY/pi_vy", m_pi_vy);
91 m_pi_vz = new TH1F( "pi_vz", "pi_vz", 100,-6.0,6.0 );
92 status = m_thistsvc->regHist("/VAL/PHY/pi_vz", m_pi_vz);
93 m_pip_mom = new TH1F( "pip_mom", "pip_moment", 100,0.0,1.6 );
94 status = m_thistsvc->regHist("/VAL/PHY/pip_mom", m_pip_mom);
95 m_pim_mom = new TH1F( "pim_mom", "pim_moment", 100,0.0,1.6 );
96 status = m_thistsvc->regHist("/VAL/PHY/pim_mom", m_pim_mom);
97 m_rhop_mass = new TH1F( "rhop_mass", "rhop_mass", 100,0.0,1.5 );
98 status = m_thistsvc->regHist("/VAL/PHY/rhop_mass", m_rhop_mass);
99 m_rhom_mass = new TH1F( "rhom_mass", "rhom_mass", 100,0.0,1.5 );
100 status = m_thistsvc->regHist("/VAL/PHY/rhom_mass", m_rhom_mass);
101 m_rho0_mass = new TH1F( "rho0_mass", "rho0_mass", 100,0.0,1.5 );
102 status = m_thistsvc->regHist("/VAL/PHY/rho0_mass", m_rho0_mass);
103 m_pi0_mass = new TH1F( "pi0_mass", "pi0_mass", 100,0.0,0.5 );
104 status = m_thistsvc->regHist("/VAL/PHY/pi0_mass", m_pi0_mass);
105 m_chisq_4c = new TH1F( "chisq_4c", "chisq_4c", 100,0.0,150. );
106 status = m_thistsvc->regHist("/VAL/PHY/chisq_4c", m_chisq_4c);
107 m_chisq_5c = new TH1F( "chisq_5c", "chisq_5c", 100,0.0,150. );
108 status = m_thistsvc->regHist("/VAL/PHY/chisq_5c", m_chisq_5c);
109
110 m_cos_pip = new TH1F( "cos_pip", "cos_pip", 100, -1.0,1.0);
111 status = m_thistsvc->regHist("/VAL/PHY/cos_pip", m_cos_pip);
112 m_cos_pim = new TH1F( "cos_pim", "cos_pim", 100, -1.0,1.0);
113 status = m_thistsvc->regHist("/VAL/PHY/cos_pim", m_cos_pim);
114
115 m_chipi_dedx = new TH1F( "chipi_dedx", "chipi_dedx", 60,-5.0,10. );
116 status = m_thistsvc->regHist("/VAL/PHY/chipi_dedx", m_chipi_dedx);
117 m_chie_dedx = new TH1F( "chie_dedx", "chie_dedx", 60,-15.0,5. );
118 status = m_thistsvc->regHist("/VAL/PHY/chie_dedx", m_chie_dedx);
119 m_chimu_dedx = new TH1F( "chimu_dedx", "chimu_dedx", 60,-5.0,10. );
120 status = m_thistsvc->regHist("/VAL/PHY/chimu_dedx", m_chimu_dedx);
121 m_chik_dedx = new TH1F( "chik_dedx", "chik_dedx", 100,-20.0,10. );
122 status = m_thistsvc->regHist("/VAL/PHY/chik_dedx", m_chik_dedx);
123 m_chip_dedx = new TH1F( "chip_dedx", "chip_dedx", 100,-20.0,5. );
124 status = m_thistsvc->regHist("/VAL/PHY/chip_dedx", m_chip_dedx);
125
126 NTuplePtr nt4(ntupleSvc(), "FILE1/h4");
127 if ( nt4 ) m_tuple4 = nt4;
128 else {
129 m_tuple4 = ntupleSvc()->book ("FILE1/h4", CLID_ColumnWiseTuple, "4gam6pi Ntuple");
130 if ( m_tuple4 ) {
131 status = m_tuple4->addItem ("ngam", m_ngoodneu);
132 status = m_tuple4->addItem ("npip", m_npip);
133 status = m_tuple4->addItem ("npim", m_npim);
134 status = m_tuple4->addItem ("ngoodch", m_ngoodch);
135 status = m_tuple4->addItem ("chisq4c", m_chisq4c);
136 status = m_tuple4->addItem ("ppi04c", m_ppi0);
137 status = m_tuple4->addItem ("mpi04c", m_mpi0);
138 status = m_tuple4->addItem ("chisq5c", m_chisq5c);
139 status = m_tuple4->addItem ("ppi05c", m_ppi0fit);
140 status = m_tuple4->addItem ("mpi05c", m_mpi0fit);
141 status = m_tuple4->addItem ("g1inpi0the", m_g1inpi0the);
142 status = m_tuple4->addItem ("g2inpi0the", m_g2inpi0the);
143 status = m_tuple4->addItem ("theta2pi", m_theta2pi);
144 status = m_tuple4->addItem ("ppip", m_ppip);
145 status = m_tuple4->addItem ("ppim", m_ppim);
146 status = m_tuple4->addItem ("p2pi", m_p2pi);
147 status = m_tuple4->addItem ("m2pi", m_m2pi);
148 status = m_tuple4->addItem ("ppip0", m_ppip0);
149 status = m_tuple4->addItem ("mpip0", m_mpip0);
150 status = m_tuple4->addItem ("ppim0", m_ppim0);
151 status = m_tuple4->addItem ("mpim0", m_mpim0);
152 status = m_tuple4->addItem ("eneumiss", m_eneumiss);
153 status = m_tuple4->addItem ("pneumiss", m_pneumiss);
154 status = m_tuple4->addItem ("mneumiss", m_mneumiss);
155// status = m_tuple4->addIndexedItem ("kal_cos", m_ngoodch , m_kal_cos);
156 }
157 else {
158 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
159 return StatusCode::FAILURE;
160 }
161 }
162
163 NTuplePtr nt5(ntupleSvc(), "FILE1/h5");
164 if ( nt5 ) m_tuple5 = nt5;
165 else {
166 m_tuple5 = ntupleSvc()->book ("FILE1/h5", CLID_ColumnWiseTuple, "4gam6pi Ntuple");
167 if ( m_tuple5 ) {
168 status = m_tuple5->addItem ("m2graw", m_m2graw);
169 status = m_tuple5->addItem ("emiss", m_emiss);
170 status = m_tuple5->addItem ("pmiss", m_pmiss);
171 status = m_tuple5->addItem ("mmiss", m_mmiss);
172 status = m_tuple5->addItem ("prho0", m_prho0);
173 status = m_tuple5->addItem ("mrho0", m_mrho0);
174 status = m_tuple5->addItem ("pmrho0", m_pmrho0);
175 status = m_tuple5->addItem ("mmrho0", m_mmrho0);
176 status = m_tuple5->addItem ("prhop", m_prhop);
177 status = m_tuple5->addItem ("mrhop", m_mrhop);
178 status = m_tuple5->addItem ("pmrhom", m_pmrhom);
179 status = m_tuple5->addItem ("mmrhom", m_mmrhom);
180 status = m_tuple5->addItem ("prhom", m_prhom);
181 status = m_tuple5->addItem ("ppipraw", m_ppipraw);
182 status = m_tuple5->addItem ("mrhom", m_mrhom);
183
184 }
185 else {
186 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
187 return StatusCode::FAILURE;
188 }
189 }
190
191
192 //
193 //--------end of book--------
194 //
195
196 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
197 return StatusCode::SUCCESS;
198
199}
INTupleSvc * ntupleSvc()

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