BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
K3pipi0.cxx
Go to the documentation of this file.
1//
2// K3pipi0.cxx is the single D0 tag code to reconstruct D0 or anti-D0 through the final states
3// of K3pipi0 from D0 decays. K3pipi0.cxx was transfered from the Fortran routine "K3pipi0.f"
4// which was orignally used for study of the D0D0-bar production and D0 decays at the BES-II
5// experiment during the time period from 2002 to 2008.
6//
7// The orignal Fortran routine "K3pipi0.f" used at the BES-II experiment was coded by G. Rong
8// in 2001.
9//
10// K3pipi0.cxx was transfered by G. Rong and J. Liu in December, 2005.
11//
12// Since 2008, G. Rong and L.L. Jiang have been working on developing this code to analyze of
13// the data taken at 3.773 GeV with the BES-III detector at the BEPC-II collider.
14//
15// During developing this code, many People made significant contributions to this code. These
16// are
17// G. Rong, L.L. Jiang, J. Liu, H.L. Ma, J.C. Chen, D.H. Zhang,
18// M.G. Zhao, B. Zheng, L. Li, Y. Fang, Z.Y. Yi, H.H. Liu, Z.Q. Liu et al.
19//
20// By G. Rong and L.L. Jiang
21// March, 2009
22//
23// ==========================================================================================
24//
25#include "SD0TagAlg/K3pipi0.h"
26#include "SD0TagAlg/SingleBase.h"
27
28K3pipi0::K3pipi0() {}
29
30K3pipi0::~K3pipi0() {}
31
32void K3pipi0::MTotal( double event, SmartDataPtr<EvtRecTrackCol> evtRecTrkCol, Vint iGood,
33 Vint iGam, double Ebeam, int PID_flag, int Charge_candidate_D ) {
34
35 int nGood = iGood.size();
36 int nGam = iGam.size();
37 iGoodtag.clear();
38 iGamtag.clear();
39
40 double mass_bcgg, delE_tag_temp;
41 int m_chargetag, m_chargek, m_chargepi1, m_chargepi2, m_chargepi3;
42 int ika_temp, ipi1_temp, ipi2_temp, ipi3_temp, ipi4_temp, iGam1_temp, iGam2_temp;
43 HepLorentzVector kmfit1, kmfit2, kmfit3, kmfit4, pddd;
44
45 int cqtm_temp;
46 HepLorentzVector pddd_temp;
47 IDataProviderSvc* eventSvc = NULL;
48 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
49 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc, EventModel::EvtRec::EvtRecEvent );
50 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
51
52 int runNo = eventHeader->runNumber();
53 int rec = eventHeader->eventNumber();
54
55 double xecm = 2 * Ebeam;
56
57 k3pipi0md = false;
58 double tagmode = 0;
59
60 if ( ( evtRecEvent->totalCharged() < 4 || nGam < 2 ) ) { return; }
61
62 double ecms = xecm;
63
64 ISimplePIDSvc* simple_pid;
65 Gaudi::svcLocator()->service( "SimplePIDSvc", simple_pid );
66
67 double deltaE_tem = 0.20;
68 int ncount1 = 0;
69
70 HepLorentzVector p2gfit;
71 HepLorentzVector p2gg;
72
73 Hep3Vector xorigin( 0, 0, 0 );
74 IVertexDbSvc* vtxsvc;
75 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
76 if ( vtxsvc->isVertexValid() )
77 {
78 double* dbv = vtxsvc->PrimaryVertex();
79 double* vv = vtxsvc->SigmaPrimaryVertex();
80 xorigin.setX( dbv[0] );
81 xorigin.setY( dbv[1] );
82 xorigin.setZ( dbv[2] );
83 }
84
85 double xv = xorigin.x();
86 double yv = xorigin.y();
87 double zv = xorigin.z();
88
89 HepPoint3D point0( 0., 0., 0. );
90 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
91
92 HepLorentzVector ptrk1_temp, ptrk2_temp, ptrk3_temp, ptrk4_temp, ptrk5_temp, ptrk6_temp,
93 ptrk7_temp;
94 //////////////////////////////////////////////////////////////////
95 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
96 {
97 EvtRecTrackIterator itTrk1 = evtRecTrkCol->begin() + i;
98
99 int ika = ( *itTrk1 )->trackId();
100
101 if ( !( *itTrk1 )->isMdcKalTrackValid() ) continue;
102 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk1 )->mdcKalTrack();
103 RecMdcKalTrack::setPidType( RecMdcKalTrack::kaon );
104 /////////////////////////////////////////
105 m_chargek = mdcKalTrk1->charge();
106 if ( Charge_candidate_D != 0 )
107 {
108 if ( m_chargek != -Charge_candidate_D ) continue;
109 }
110 if ( Charge_candidate_D == 0 )
111 {
112 if ( abs( m_chargek ) != 1 ) continue;
113 }
114 /////////////////////////////////////////
115 HepVector a1 = mdcKalTrk1->getZHelixK();
116 HepSymMatrix Ea1 = mdcKalTrk1->getZErrorK();
117 VFHelix helixip3_1( point0, a1, Ea1 );
118 helixip3_1.pivot( IP );
119 HepVector vecipa1 = helixip3_1.a();
120
121 double dr1 = fabs( vecipa1[0] );
122 double dz1 = fabs( vecipa1[3] );
123 double costheta1 = cos( mdcKalTrk1->theta() );
124 if ( dr1 >= 1.0 ) continue;
125 if ( dz1 >= 10.0 ) continue;
126 if ( fabs( costheta1 ) >= 0.93 ) continue;
127 /////////////////////////////////////////
128 if ( PID_flag == 5 )
129 {
130 simple_pid->preparePID( *itTrk1 );
131 if ( simple_pid->probKaon() < 0.0 || simple_pid->probKaon() < simple_pid->probPion() )
132 continue;
133 }
134
135 /////////////////////////////////////////
136
137 WTrackParameter kam( xmass[3], mdcKalTrk1->getZHelixK(), mdcKalTrk1->getZErrorK() );
138
139 //
140 // select pi1
141 //
142 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
143 {
144 EvtRecTrackIterator itTrk2 = evtRecTrkCol->begin() + j;
145
146 int ipi1 = ( *itTrk2 )->trackId();
147 if ( ipi1 == ika ) continue;
148
149 if ( !( *itTrk2 )->isMdcKalTrackValid() ) continue;
150 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk2 )->mdcKalTrack();
151 RecMdcKalTrack::setPidType( RecMdcKalTrack::pion );
152 /////////////////////////////////////////
153 m_chargepi1 = mdcKalTrk2->charge();
154 if ( ( m_chargek + m_chargepi1 ) != 0 ) continue;
155 /////////////////////////////////////////
156 HepVector a2 = mdcKalTrk2->getZHelix();
157 HepSymMatrix Ea2 = mdcKalTrk2->getZError();
158 VFHelix helixip3_2( point0, a2, Ea2 );
159 helixip3_2.pivot( IP );
160 HepVector vecipa2 = helixip3_2.a();
161
162 double dr2 = fabs( vecipa2[0] );
163 double dz2 = fabs( vecipa2[3] );
164 double costheta2 = cos( mdcKalTrk2->theta() );
165 if ( dr2 >= 1.0 ) continue;
166 if ( dz2 >= 10.0 ) continue;
167 if ( fabs( costheta2 ) >= 0.93 ) continue;
168 /////////////////////////////////////////
169 if ( PID_flag == 5 )
170 {
171 simple_pid->preparePID( *itTrk2 );
172 if ( simple_pid->probPion() < 0.0 || simple_pid->probPion() < simple_pid->probKaon() )
173 continue;
174 }
175 /////////////////////////////////////////
176 WTrackParameter pip1( xmass[2], mdcKalTrk2->getZHelix(), mdcKalTrk2->getZError() );
177
178 //
179 // select pi2
180 //
181 for ( int k = 0; k < evtRecEvent->totalCharged(); k++ )
182 {
183 EvtRecTrackIterator itTrk3 = evtRecTrkCol->begin() + k;
184
185 int ipi2 = ( *itTrk3 )->trackId();
186 if ( ipi2 == ika || ipi2 == ipi1 ) continue;
187
188 if ( !( *itTrk3 )->isMdcKalTrackValid() ) continue;
189 RecMdcKalTrack* mdcKalTrk3 = ( *itTrk3 )->mdcKalTrack();
190 RecMdcKalTrack::setPidType( RecMdcKalTrack::pion );
191 /////////////////////////////////////////
192 m_chargepi2 = mdcKalTrk3->charge();
193 if ( ( m_chargek + m_chargepi2 ) != 0 ) continue;
194 /////////////////////////////////////////
195 HepVector a3 = mdcKalTrk3->getZHelix();
196 HepSymMatrix Ea3 = mdcKalTrk3->getZError();
197 VFHelix helixip3_3( point0, a3, Ea3 );
198 helixip3_3.pivot( IP );
199 HepVector vecipa3 = helixip3_3.a();
200
201 double dr3 = fabs( vecipa3[0] );
202 double dz3 = fabs( vecipa3[3] );
203 double costheta3 = cos( mdcKalTrk3->theta() );
204 if ( dr3 >= 1.0 ) continue;
205 if ( dz3 >= 10.0 ) continue;
206 if ( fabs( costheta3 ) >= 0.93 ) continue;
207 /////////////////////////////////////////
208 if ( PID_flag == 5 )
209 {
210 simple_pid->preparePID( *itTrk3 );
211 if ( simple_pid->probPion() < 0.0 ||
212 simple_pid->probPion() < simple_pid->probKaon() )
213 continue;
214 }
215 /////////////////////////////////////////
216 WTrackParameter pip2( xmass[2], mdcKalTrk3->getZHelix(), mdcKalTrk3->getZError() );
217
218 //
219 // select pi3
220 //
221 for ( int l = 0; l < evtRecEvent->totalCharged(); l++ )
222 {
223 EvtRecTrackIterator itTrk4 = evtRecTrkCol->begin() + l;
224
225 int ipi3 = ( *itTrk4 )->trackId();
226 if ( ipi3 == ika || ipi3 == ipi1 || ipi3 == ipi2 ) continue;
227
228 if ( !( *itTrk4 )->isMdcKalTrackValid() ) continue;
229 RecMdcKalTrack* mdcKalTrk4 = ( *itTrk4 )->mdcKalTrack();
230 RecMdcKalTrack::setPidType( RecMdcKalTrack::pion );
231 /////////////////////////////////////////
232 m_chargepi3 = mdcKalTrk4->charge();
233 if ( ( m_chargepi2 + m_chargepi3 ) != 0 ) continue;
234 /////////////////////////////////////////
235 HepVector a4 = mdcKalTrk4->getZHelix();
236 HepSymMatrix Ea4 = mdcKalTrk4->getZError();
237 VFHelix helixip3_4( point0, a4, Ea4 );
238 helixip3_4.pivot( IP );
239 HepVector vecipa4 = helixip3_4.a();
240
241 double dr4 = fabs( vecipa4[0] );
242 double dz4 = fabs( vecipa4[3] );
243 double costheta4 = cos( mdcKalTrk4->theta() );
244 if ( dr4 >= 1.0 ) continue;
245 if ( dz4 >= 10.0 ) continue;
246 if ( fabs( costheta4 ) >= 0.93 ) continue;
247 /////////////////////////////////////////
248 if ( PID_flag == 5 )
249 {
250 simple_pid->preparePID( *itTrk4 );
251 if ( simple_pid->probPion() < 0.0 ||
252 simple_pid->probPion() < simple_pid->probKaon() )
253 continue;
254 }
255 /////////////////////////////////////////
256 WTrackParameter pip3( xmass[2], mdcKalTrk4->getZHelix(), mdcKalTrk4->getZError() );
257
258 for ( int m = 0; m < nGam - 1; m++ )
259 {
260 if ( iGam[m] == -1 ) continue;
261 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[m] ) )->emcShower();
262 double eraw1 = g1Trk->energy();
263 double phi1 = g1Trk->phi();
264 double the1 = g1Trk->theta();
265 HepLorentzVector ptrkg1, ptrkg10, ptrkg12;
266 ptrkg1.setPx( eraw1 * sin( the1 ) * cos( phi1 ) );
267 ptrkg1.setPy( eraw1 * sin( the1 ) * sin( phi1 ) );
268 ptrkg1.setPz( eraw1 * cos( the1 ) );
269 ptrkg1.setE( eraw1 );
270 ptrkg10 = ptrkg1;
271 ptrkg12 = ptrkg1.boost( -0.011, 0, 0 );
272
273 for ( int n = m + 1; n < nGam; n++ )
274 {
275 if ( iGam[n] == -1 ) continue;
276 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[n] ) )->emcShower();
277 double eraw2 = g2Trk->energy();
278 double phi2 = g2Trk->phi();
279 double the2 = g2Trk->theta();
280 HepLorentzVector ptrkg2, ptrkg20, ptrkg22;
281 ptrkg2.setPx( eraw2 * sin( the2 ) * cos( phi2 ) );
282 ptrkg2.setPy( eraw2 * sin( the2 ) * sin( phi2 ) );
283 ptrkg2.setPz( eraw2 * cos( the2 ) );
284 ptrkg2.setE( eraw2 );
285 ptrkg20 = ptrkg2;
286 ptrkg22 = ptrkg2.boost( -0.011, 0, 0 );
287
288 /////////////////////////////////////////////////////////////
289 HepLorentzVector ptrkpi0;
290 ptrkpi0 = ptrkg12 + ptrkg22;
291 double m_xmpi0_tem = ptrkpi0.mag();
292 if ( m_xmpi0_tem > 0.150 || m_xmpi0_tem < 0.115 ) continue;
293 /////////////////////////////////////////////////////////////
294 bool IsEndcap1 = false;
295 bool IsEndcap2 = false;
296 if ( fabs( cos( the1 ) ) > 0.86 && fabs( cos( the1 ) ) < 0.92 ) IsEndcap1 = true;
297 if ( fabs( cos( the2 ) ) > 0.86 && fabs( cos( the2 ) ) < 0.92 ) IsEndcap2 = true;
298 if ( IsEndcap1 && IsEndcap2 ) continue;
299 /////////////////////////////////////////////////////////////
300
301 KalmanKinematicFit* kmfit = KalmanKinematicFit::instance();
302 kmfit->init();
303 kmfit->setChisqCut( 2500 );
304 kmfit->AddTrack( 0, 0.0, g1Trk );
305 kmfit->AddTrack( 1, 0.0, g2Trk );
306 kmfit->AddResonance( 0, mpi0, 0, 1 );
307
308 kmfit->Fit( 0 ); // Perform fit
309 kmfit->BuildVirtualParticle( 0 );
310
311 double pi0_chisq = kmfit->chisq( 0 );
312 if ( pi0_chisq >= 2500 ) continue;
313 HepLorentzVector p2gfit = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
314 p2gfit.boost( -0.011, 0, 0 );
315
316 //////////////////////////////////////////////////////////////
317 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
318 HepSymMatrix Evx( 3, 0 );
319 double bx = 1E+6;
320 Evx[0][0] = bx * bx;
321 double by = 1E+6;
322 Evx[1][1] = by * by;
323 double bz = 1E+6;
324 Evx[2][2] = bz * bz;
325 VertexParameter vxpar;
326 vxpar.setVx( vx );
327 vxpar.setEvx( Evx );
328 //////////////////////////////////////////////////////////////
329
330 VertexFit* vtxfit = VertexFit::instance();
331 vtxfit->init();
332 vtxfit->AddTrack( 0, kam );
333 vtxfit->AddTrack( 1, pip1 );
334 vtxfit->AddTrack( 2, pip2 );
335 vtxfit->AddTrack( 3, pip3 );
336 vtxfit->AddVertex( 0, vxpar, 0, 1, 2, 3 );
337 if ( !vtxfit->Fit( 0 ) ) continue;
338 vtxfit->Swim( 0 );
339
340 WTrackParameter wkam = vtxfit->wtrk( 0 );
341 WTrackParameter wpip1 = vtxfit->wtrk( 1 );
342 WTrackParameter wpip2 = vtxfit->wtrk( 2 );
343 WTrackParameter wpip3 = vtxfit->wtrk( 3 );
344
345 HepVector kam_val = HepVector( 7, 0 );
346 HepVector pip1_val = HepVector( 7, 0 );
347 HepVector pip2_val = HepVector( 7, 0 );
348 HepVector pip3_val = HepVector( 7, 0 );
349 kam_val = wkam.w();
350 pip1_val = wpip1.w();
351 pip2_val = wpip2.w();
352 pip3_val = wpip3.w();
353
354 HepLorentzVector P_KAM( kam_val[0], kam_val[1], kam_val[2], kam_val[3] );
355 HepLorentzVector P_PIP1( pip1_val[0], pip1_val[1], pip1_val[2], pip1_val[3] );
356 HepLorentzVector P_PIP2( pip2_val[0], pip2_val[1], pip2_val[2], pip2_val[3] );
357 HepLorentzVector P_PIP3( pip3_val[0], pip3_val[1], pip3_val[2], pip3_val[3] );
358
359 P_KAM.boost( -0.011, 0, 0 );
360 P_PIP1.boost( -0.011, 0, 0 );
361 P_PIP2.boost( -0.011, 0, 0 );
362 P_PIP3.boost( -0.011, 0, 0 );
363 pddd = P_KAM + P_PIP1 + P_PIP2 + P_PIP3 + p2gfit;
364
365 double pk3pipi0 = pddd.rho();
366
367 double temp1 = ( ecms / 2 ) * ( ecms / 2 ) - pk3pipi0 * pk3pipi0;
368 if ( temp1 < 0 ) temp1 = 0;
369 double mass_bc_tem = sqrt( temp1 );
370 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
371
372 double delE_tag_tag = ecms / 2 - pddd.e();
373
374 if ( fabs( delE_tag_tag ) < deltaE_tem )
375 {
376 deltaE_tem = fabs( delE_tag_tag );
377 delE_tag_temp = delE_tag_tag;
378 mass_bcgg = mass_bc_tem;
379
380 pddd_temp = pddd;
381 cqtm_temp = m_chargek;
382
383 ika_temp = ika;
384 ipi1_temp = ipi1;
385 ipi2_temp = ipi2;
386 ipi3_temp = ipi3;
387 iGam1_temp = iGam[m];
388 iGam2_temp = iGam[n];
389
390 ncount1 = 1;
391 }
392 }
393 }
394 }
395 }
396 }
397 }
398 if ( ncount1 == 1 )
399 {
400 tagmode = 24;
401 if ( cqtm_temp < 0 ) tagmode = -24;
402 tagmd = tagmode;
403 mass_bc = mass_bcgg;
404 delE_tag = delE_tag_temp;
405 cqtm = -1.0 * cqtm_temp;
406
407 iGoodtag.push_back( ika_temp );
408 iGoodtag.push_back( ipi1_temp );
409 iGoodtag.push_back( ipi2_temp );
410 iGoodtag.push_back( ipi3_temp );
411
412 iGamtag.push_back( iGam1_temp );
413 iGamtag.push_back( iGam2_temp );
414 iGamtag.push_back( 9999 );
415 iGamtag.push_back( 9999 );
416
417 ptag = pddd_temp;
418
419 k3pipi0md = true;
420 }
421}
runNo
int runNo
Definition DQA_TO_DB.cxx:13
n
const Int_t n
Definition DataBase/tau_mode.c:57
phi2
Double_t phi2
Definition DataBase/tau_mode.c:8
phi1
Double_t phi1
Definition DataBase/tau_mode.c:8
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi0
double mpi0
Definition EvtPhokharaDef.hh:17
xmass
const double xmass[5]
Definition Gam4pikp.cxx:35
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
Vint
std::vector< int > Vint
Definition K3pipi0.h:16
ISimplePIDSvc::probKaon
virtual double probKaon()=0
ISimplePIDSvc::preparePID
virtual void preparePID(EvtRecTrack *track)=0
ISimplePIDSvc::probPion
virtual double probPion()=0
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
K3pipi0::K3pipi0
K3pipi0()
Definition K3pipi0.cxx:28
K3pipi0::~K3pipi0
~K3pipi0()
Definition K3pipi0.cxx:30
K3pipi0::MTotal
void MTotal(double event, SmartDataPtr< EvtRecTrackCol > evtRecTrkCol, Vint iGood, Vint iGam, double Ebeam, int PID_flag, int Charge_candidate_D)
Definition K3pipi0.cxx:32
KalmanKinematicFit::setChisqCut
void setChisqCut(const double chicut=200, const double chiter=0.05)
Definition Analysis/VertexFit/include/VertexFit/KalmanKinematicFit.h:122
KalmanKinematicFit::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition KalmanKinematicFit.cxx:900
KalmanKinematicFit::AddResonance
void AddResonance(int number, double mres, std::vector< int > tlis)
Definition KalmanKinematicFit.cxx:142
KalmanKinematicFit::instance
static KalmanKinematicFit * instance()
Definition KalmanKinematicFit.cxx:21
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
VFHelix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Analysis/VertexFit/include/VertexFit/Helix.h:177
VFHelix::a
const HepVector & a(void) const
returns helix parameters.
Definition Analysis/VertexFit/include/VertexFit/Helix.h:195
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:91
VertexFit::init
void init()
Definition VertexFit.cxx:27
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288
ecms
const double ecms
Definition inclkstar.cxx:26