33 {
34
35 int nGood = iGood.size();
36 int nGam = iGam.size();
37
38 iGoodtag.clear();
39 iGamtag.clear();
40
41 double mass_bcgg, delE_tag_temp;
42 int m_chargetag, m_chargepi1, m_chargepi2;
43 int ipi1_temp, ipi2_temp, iGam1_temp, iGam2_temp;
44 HepLorentzVector pddd, pddd_temp;
45
46 IDataProviderSvc* eventSvc = NULL;
47 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
49 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
50
51 int runNo = eventHeader->runNumber();
52 int rec = eventHeader->eventNumber();
53
54 double xecm = 2 * Ebeam;
55
56 pipipi0md = false;
57 double tagmode = 0;
58
59 if ( ( evtRecEvent->totalCharged() < 2 || nGam < 2 ) ) { return; }
60
62
63 ISimplePIDSvc* simple_pid;
64 Gaudi::svcLocator()->service( "SimplePIDSvc", simple_pid );
65
66 double deltaE_tem = 0.20;
67 int ncount1 = 0;
68
69 Hep3Vector xorigin( 0, 0, 0 );
70 IVertexDbSvc* vtxsvc;
71 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
73 {
76 xorigin.setX( dbv[0] );
77 xorigin.setY( dbv[1] );
78 xorigin.setZ( dbv[2] );
79 }
80
81 double xv = xorigin.x();
82 double yv = xorigin.y();
83 double zv = xorigin.z();
84
86 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
87
88
89 HepLorentzVector p2gfit;
90 HepLorentzVector p2gg;
91
92 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
93 {
95
96 int ipi1 = ( *itTrk )->trackId();
97
98 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
99 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
101
102 m_chargepi1 = mdcKalTrk1->
charge();
103 if (
abs( m_chargepi1 ) != 1 )
continue;
104
106 HepSymMatrix Ea1 = mdcKalTrk1->
getZError();
107 VFHelix helixip3_1( point0, a1, Ea1 );
108 helixip3_1.pivot( IP );
109 HepVector vecipa1 = helixip3_1.a();
110
111 double dr1 = fabs( vecipa1[0] );
112 double dz1 = fabs( vecipa1[3] );
113 double costheta1 =
cos( mdcKalTrk1->
theta() );
114
115 if ( dr1 >= 1.0 ) continue;
116 if ( dz1 >= 10.0 ) continue;
117 if ( fabs( costheta1 ) >= 0.93 ) continue;
118
119 if ( PID_flag == 5 )
120 {
123 continue;
124 }
125
127
128
129
130
131 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
132 {
134
135 int ipi2 = ( *itTrk )->trackId();
136 if ( ipi1 == ipi2 ) continue;
137
138 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
139 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk )->mdcKalTrack();
141
142 m_chargepi2 = mdcKalTrk2->
charge();
143 if ( ( m_chargepi1 + m_chargepi2 ) != 0 ) continue;
144
145
147 HepSymMatrix Ea2 = mdcKalTrk2->
getZError();
148 VFHelix helixip3_2( point0, a2, Ea2 );
149 helixip3_2.pivot( IP );
150 HepVector vecipa2 = helixip3_2.a();
151
152 double dr2 = fabs( vecipa2[0] );
153 double dz2 = fabs( vecipa2[3] );
154 double costheta2 =
cos( mdcKalTrk2->
theta() );
155 if ( dr2 >= 1.0 ) continue;
156 if ( dz2 >= 10.0 ) continue;
157 if ( fabs( costheta2 ) >= 0.93 ) continue;
158
159
160 if ( PID_flag == 5 )
161 {
164 continue;
165 }
166
168
169 for ( int m = 0; m < nGam - 1; m++ )
170 {
171 if ( iGam[m] == -1 ) continue;
172 int iGam1 = iGam[m];
173 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[m] ) )->emcShower();
174 double eraw1 = g1Trk->
energy();
176 double the1 = g1Trk->
theta();
177 HepLorentzVector ptrkg1, ptrkg10, ptrkg12;
178 ptrkg1.setPx( eraw1 *
sin( the1 ) *
cos(
phi1 ) );
179 ptrkg1.setPy( eraw1 *
sin( the1 ) *
sin(
phi1 ) );
180 ptrkg1.setPz( eraw1 *
cos( the1 ) );
181 ptrkg1.setE( eraw1 );
182 ptrkg10 = ptrkg1;
183 ptrkg12 = ptrkg1.boost( -0.011, 0, 0 );
184
185 for (
int n = m + 1;
n < nGam;
n++ )
186 {
187 if ( iGam[
n] == -1 )
continue;
188 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[
n] ) )->emcShower();
190 double eraw2 = g2Trk->
energy();
192 double the2 = g2Trk->
theta();
193 HepLorentzVector ptrkg2, ptrkg20, ptrkg22;
194 ptrkg2.setPx( eraw2 *
sin( the2 ) *
cos(
phi2 ) );
195 ptrkg2.setPy( eraw2 *
sin( the2 ) *
sin(
phi2 ) );
196 ptrkg2.setPz( eraw2 *
cos( the2 ) );
197 ptrkg2.setE( eraw2 );
198 ptrkg20 = ptrkg2;
199 ptrkg22 = ptrkg2.boost( -0.011, 0, 0 );
200
201
202 HepLorentzVector ptrkpi0;
203 ptrkpi0 = ptrkg12 + ptrkg22;
204 double m_xmpi0_tem = ptrkpi0.m();
205 if ( m_xmpi0_tem > 0.150 || m_xmpi0_tem < 0.115 ) continue;
206
207 bool IsEndcap1 = false;
208 bool IsEndcap2 = false;
209 if ( fabs(
cos( the1 ) ) > 0.86 && fabs(
cos( the1 ) ) < 0.92 ) IsEndcap1 =
true;
210 if ( fabs(
cos( the2 ) ) > 0.86 && fabs(
cos( the2 ) ) < 0.92 ) IsEndcap2 =
true;
211 if ( IsEndcap1 && IsEndcap2 ) continue;
212
219
222
223 double pi0_chisq = kmfit->
chisq( 0 );
224 if ( pi0_chisq >= 2500 ) continue;
225 HepLorentzVector p2gfit = kmfit->
pfit( 0 ) + kmfit->
pfit( 1 );
226 p2gfit.boost( -0.011, 0, 0 );
227
228
229 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
230 HepSymMatrix Evx( 3, 0 );
231 double bx = 1E+6;
232 Evx[0][0] = bx * bx;
233 double by = 1E+6;
234 Evx[1][1] = by * by;
235 double bz = 1E+6;
236 Evx[2][2] = bz * bz;
237 VertexParameter vxpar;
240
241
247 if ( !vtxfit->
Fit( 0 ) )
continue;
249
250 WTrackParameter wpip = vtxfit->
wtrk( 0 );
251 WTrackParameter wpim = vtxfit->
wtrk( 1 );
252
253 HepVector pip_val = HepVector( 7, 0 );
254 HepVector pim_val = HepVector( 7, 0 );
257
258 HepLorentzVector P_PIP( pip_val[0], pip_val[1], pip_val[2], pip_val[3] );
259 HepLorentzVector P_PIM( pim_val[0], pim_val[1], pim_val[2], pim_val[3] );
260
261 P_PIM.boost( -0.011, 0, 0 );
262 P_PIP.boost( -0.011, 0, 0 );
263 HepLorentzVector ppipi = P_PIM + P_PIP;
264 pddd = P_PIM + P_PIP + p2gfit;
265
266 double mpipi = ppipi.m();
267
268
269 double ppipipi0 = pddd.rho();
270
271 double temp1 = (
ecms / 2 ) * (
ecms / 2 ) - ppipipi0 * ppipipi0;
272 if ( temp1 < 0 ) temp1 = 0;
273 double mass_bc_tem = sqrt( temp1 );
274 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
275
276 double delE_tag_tag =
ecms / 2 - pddd.e();
277
278 if ( fabs( delE_tag_tag ) < deltaE_tem )
279 {
280 deltaE_tem = fabs( delE_tag_tag );
281 delE_tag_temp = delE_tag_tag;
282 mass_bcgg = mass_bc_tem;
283
284 pddd_temp = pddd;
285
286 ipi1_temp = ipi1;
287 ipi2_temp = ipi2;
288 iGam1_temp = iGam1;
289 iGam2_temp = iGam2;
290
291 ncount1 = 1;
292 }
293 }
294 }
295 }
296 }
297
298 if ( ncount1 == 1 )
299 {
300 tagmode = 16;
301 if ( m_chargetag < 0 ) tagmode = -16;
302 tagmd = tagmode;
303 mass_bc = mass_bcgg;
304 delE_tag = delE_tag_temp;
305 cqtm = 0;
306
307 iGoodtag.push_back( ipi1_temp );
308 iGoodtag.push_back( ipi2_temp );
309 iGamtag.push_back( iGam1_temp );
310 iGamtag.push_back( iGam2_temp );
311 iGamtag.push_back( 9999 );
312 iGamtag.push_back( 9999 );
313
314 ptag = pddd_temp;
315
316 pipipi0md = true;
317 }
318}
HepGeom::Point3D< double > HepPoint3D
EvtRecTrackCol::iterator EvtRecTrackIterator
double sin(const BesAngle a)
double cos(const BesAngle a)
const double theta() const
static void setPidType(PidType pidType)
virtual double probKaon()=0
virtual void preparePID(EvtRecTrack *track)=0
virtual double probPion()=0
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
void setChisqCut(const double chicut=200, const double chiter=0.05)
HepLorentzVector pfit(int n) const
void BuildVirtualParticle(int number)
void AddResonance(int number, double mres, std::vector< int > tlis)
static KalmanKinematicFit * instance()
const HepVector & getZHelix() const
const HepSymMatrix & getZError() const
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
WTrackParameter wtrk(int n) const
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
static VertexFit * instance()
void setEvx(const HepSymMatrix &eVx)
void setVx(const HepPoint3D &vx)
_EXTERN_ std::string EvtRecEvent
1.16.1