BOSS 7.0.8
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 56 of file rhopi.cxx.

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

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

◆ 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 200 of file rhopi.cxx.

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

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

◆ initialize()

StatusCode rhopi::initialize ( )

Definition at line 73 of file rhopi.cxx.

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

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