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_chargek, m_chargepi;
43 int ika_temp, ipi_temp, iGam1_temp, iGam2_temp;
44 HepLorentzVector pddd, pddd_temp;
45
46 int cqtm_temp;
47 IDataProviderSvc* eventSvc = NULL;
48 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
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 kpipi0md = false;
58 double tagmode = 0;
59
60 if ( ( evtRecEvent->totalCharged() < 2 || nGam < 2 ) ) { return; }
61
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 Hep3Vector xorigin( 0, 0, 0 );
71 IVertexDbSvc* vtxsvc;
72 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
74 {
77 xorigin.setX( dbv[0] );
78 xorigin.setY( dbv[1] );
79 xorigin.setZ( dbv[2] );
80 }
81
82 double xv = xorigin.x();
83 double yv = xorigin.y();
84 double zv = xorigin.z();
85
87 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
88
89 HepLorentzVector ptrk1_temp, ptrk2_temp, ptrk3_temp, ptrk4_temp, ptrk5_temp;
90
91 HepLorentzVector p2gfit;
92 HepLorentzVector p2gg;
93 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
94 {
96
97 int ika = ( *itTrk )->trackId();
98
99 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
100 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
102
103 m_chargek = mdcKalTrk1->
charge();
104 if ( Charge_candidate_D != 0 )
105 {
106 if ( m_chargek != -Charge_candidate_D ) continue;
107 }
108 if ( Charge_candidate_D == 0 )
109 {
110 if (
abs( m_chargek ) != 1 )
continue;
111 }
112
115 VFHelix helixip3_1( point0, a1, Ea1 );
116 helixip3_1.pivot( IP );
117 HepVector vecipa1 = helixip3_1.a();
118
119 double dr1 = fabs( vecipa1[0] );
120 double dz1 = fabs( vecipa1[3] );
121 double costheta1 =
cos( mdcKalTrk1->
theta() );
122 if ( dr1 >= 1.0 ) continue;
123 if ( dz1 >= 10.0 ) continue;
124 if ( fabs( costheta1 ) >= 0.93 ) continue;
125
126
127 if ( PID_flag == 5 )
128 {
131 continue;
132 }
133
134
136
137
138
139
140 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
141 {
143
144 int ipi = ( *itTrk )->trackId();
145 if ( ipi == ika ) continue;
146
147 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
148 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk )->mdcKalTrack();
150
151 m_chargepi = mdcKalTrk2->
charge();
152 if ( ( m_chargek + m_chargepi ) != 0 ) continue;
153
155 HepSymMatrix Ea2 = mdcKalTrk2->
getZError();
156 VFHelix helixip3_2( point0, a2, Ea2 );
157 helixip3_2.pivot( IP );
158 HepVector vecipa2 = helixip3_2.a();
159
160 double dr2 = fabs( vecipa2[0] );
161 double dz2 = fabs( vecipa2[3] );
162 double costheta2 =
cos( mdcKalTrk2->
theta() );
163 if ( dr2 >= 1.0 ) continue;
164 if ( dz2 >= 10.0 ) continue;
165 if ( fabs( costheta2 ) >= 0.93 ) continue;
166
167 if ( PID_flag == 5 )
168 {
171 continue;
172 }
173
174
176
177 for ( int m = 0; m < nGam - 1; m++ )
178 {
179 if ( iGam[m] == -1 ) continue;
180 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[m] ) )->emcShower();
181 double eraw1 = g1Trk->
energy();
183 double the1 = g1Trk->
theta();
184 HepLorentzVector ptrkg1, ptrkg10, ptrkg12;
185 ptrkg1.setPx( eraw1 *
sin( the1 ) *
cos(
phi1 ) );
186 ptrkg1.setPy( eraw1 *
sin( the1 ) *
sin(
phi1 ) );
187 ptrkg1.setPz( eraw1 *
cos( the1 ) );
188 ptrkg1.setE( eraw1 );
189 ptrkg10 = ptrkg1;
190 ptrkg12 = ptrkg1.boost( -0.011, 0, 0 );
191
192 for (
int n = m + 1;
n < nGam;
n++ )
193 {
194 if ( iGam[
n] == -1 )
continue;
195 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[
n] ) )->emcShower();
196 double eraw2 = g2Trk->
energy();
198 double the2 = g2Trk->
theta();
199 HepLorentzVector ptrkg2, ptrkg20, ptrkg22;
200 ptrkg2.setPx( eraw2 *
sin( the2 ) *
cos(
phi2 ) );
201 ptrkg2.setPy( eraw2 *
sin( the2 ) *
sin(
phi2 ) );
202 ptrkg2.setPz( eraw2 *
cos( the2 ) );
203 ptrkg2.setE( eraw2 );
204 ptrkg20 = ptrkg2;
205 ptrkg22 = ptrkg2.boost( -0.011, 0, 0 );
206
207
208 HepLorentzVector ptrkpi0;
209 ptrkpi0 = ptrkg12 + ptrkg22;
210 double m_xmpi0_tem = ptrkpi0.m();
211 if ( m_xmpi0_tem > 0.150 || m_xmpi0_tem < 0.115 ) continue;
212
213 bool IsEndcap1 = false;
214 bool IsEndcap2 = false;
215 if ( fabs(
cos( the1 ) ) > 0.86 && fabs(
cos( the1 ) ) < 0.92 ) IsEndcap1 =
true;
216 if ( fabs(
cos( the2 ) ) > 0.86 && fabs(
cos( the2 ) ) < 0.92 ) IsEndcap2 =
true;
217 if ( IsEndcap1 && IsEndcap2 ) continue;
218
225
228
229 double pi0_chisq = kmfit->
chisq( 0 );
230 if ( pi0_chisq >= 2500 ) continue;
231 HepLorentzVector p2gfit = kmfit->
pfit( 0 ) + kmfit->
pfit( 1 );
232 p2gfit.boost( -0.011, 0, 0 );
233
234
235 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
236 HepSymMatrix Evx( 3, 0 );
237 double bx = 1E+6;
238 Evx[0][0] = bx * bx;
239 double by = 1E+6;
240 Evx[1][1] = by * by;
241 double bz = 1E+6;
242 Evx[2][2] = bz * bz;
243 VertexParameter vxpar;
246
247
253 if ( !vtxfit->
Fit( 0 ) )
continue;
255
256 WTrackParameter wkam = vtxfit->
wtrk( 0 );
257 WTrackParameter wpip = vtxfit->
wtrk( 1 );
258
259 HepVector kam_val = HepVector( 7, 0 );
261 HepVector pip_val = HepVector( 7, 0 );
263
264 HepLorentzVector P_KAM( kam_val[0], kam_val[1], kam_val[2], kam_val[3] );
265 HepLorentzVector P_PIP( pip_val[0], pip_val[1], pip_val[2], pip_val[3] );
266
267 P_KAM.boost( -0.011, 0, 0 );
268 P_PIP.boost( -0.011, 0, 0 );
269 pddd = P_KAM + P_PIP + p2gfit;
270
271 double pkpipi0 = pddd.rho();
272
273 double temp1 = (
ecms / 2 ) * (
ecms / 2 ) - pkpipi0 * pkpipi0;
274 if ( temp1 < 0 ) temp1 = 0;
275 double mass_bc_tem = sqrt( temp1 );
276 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
277
278 double delE_tag_tag =
ecms / 2 - pddd.e();
279
280 if ( fabs( delE_tag_tag ) < deltaE_tem )
281 {
282 deltaE_tem = fabs( delE_tag_tag );
283 delE_tag_temp = delE_tag_tag;
284 mass_bcgg = mass_bc_tem;
285
286 pddd_temp = pddd;
287 cqtm_temp = m_chargek;
288
289 ika_temp = ika;
290 ipi_temp = ipi;
291
292 iGam1_temp = iGam[m];
293 iGam2_temp = iGam[
n];
294 ncount1 = 1;
295 }
296 }
297 }
298 }
299 }
300
301 if ( ncount1 == 1 )
302 {
303 tagmode = 12;
304 if ( cqtm_temp < 0 ) tagmode = -12;
305 tagmd = tagmode;
306 mass_bc = mass_bcgg;
307 delE_tag = delE_tag_temp;
308 cqtm = -1.0 * cqtm_temp;
309
310 iGoodtag.push_back( ipi_temp );
311 iGoodtag.push_back( ika_temp );
312
313 iGamtag.push_back( iGam1_temp );
314 iGamtag.push_back( iGam2_temp );
315 iGamtag.push_back( 9999 );
316 iGamtag.push_back( 9999 );
317
318 ptag = pddd_temp;
319
320 kpipi0md = true;
321 }
322}
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
HepSymMatrix & getZErrorK()
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