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, iGam3_temp, iGam4_temp;
44 HepLorentzVector pddd;
45 HepLorentzVector pddd_temp;
46
47 int cqtm_temp;
48 IDataProviderSvc* eventSvc = NULL;
49 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
51 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
52
53 int runNo = eventHeader->runNumber();
54 int rec = eventHeader->eventNumber();
55
56 double xecm = 2 * Ebeam;
57
58 kpipi0pi0md = false;
59 double tagmode = 0;
60
61 if ( ( evtRecEvent->totalCharged() < 2 || nGam < 4 ) ) { return; }
62
64
65 ISimplePIDSvc* simple_pid;
66 Gaudi::svcLocator()->service( "SimplePIDSvc", simple_pid );
67
68 double deltaE_tem = 0.20;
69 int ncount1 = 0;
70
71 Hep3Vector xorigin( 0, 0, 0 );
72 IVertexDbSvc* vtxsvc;
73 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
75 {
78 xorigin.setX( dbv[0] );
79 xorigin.setY( dbv[1] );
80 xorigin.setZ( dbv[2] );
81 }
82
83 double xv = xorigin.x();
84 double yv = xorigin.y();
85 double zv = xorigin.z();
86
88 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
89
90 HepLorentzVector p2gfit1, p2gfit2;
91 HepLorentzVector p2gg;
92
93 HepLorentzVector ptrk1_temp, ptrk2_temp, ptrk3_temp, ptrk4_temp, ptrk5_temp, ptrk6_temp,
94 ptrk7_temp, ptrk8_temp;
95 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
96 {
98
99 int ika = ( *itTrk )->trackId();
100
101 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
102 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
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
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 {
132 continue;
133 }
134
135
137
138
139
140
141 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
142 {
144
145 int ipi = ( *itTrk )->trackId();
146 if ( ipi == ika ) continue;
147
148 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
149 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk )->mdcKalTrack();
151
152 m_chargepi = mdcKalTrk2->
charge();
153 if ( ( m_chargek + m_chargepi ) != 0 ) continue;
154
156 HepSymMatrix Ea2 = mdcKalTrk2->
getZError();
157 VFHelix helixip3_2( point0, a2, Ea2 );
158 helixip3_2.pivot( IP );
159 HepVector vecipa2 = helixip3_2.a();
160
161 double dr2 = fabs( vecipa2[0] );
162 double dz2 = fabs( vecipa2[3] );
163 double costheta2 =
cos( mdcKalTrk2->
theta() );
164 if ( dr2 >= 1.0 ) continue;
165 if ( dz2 >= 10.0 ) continue;
166 if ( fabs( costheta2 ) >= 0.93 ) continue;
167
168 if ( PID_flag == 5 )
169 {
172 continue;
173 }
174
175
177
178 for ( int m = 0; m < nGam - 1; m++ )
179 {
180 if ( iGam[m] == -1 ) continue;
181 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[m] ) )->emcShower();
182 double eraw1 = g1Trk->
energy();
184 double the1 = g1Trk->
theta();
185 HepLorentzVector ptrkg1, ptrkg10, ptrkg12;
186 ptrkg1.setPx( eraw1 *
sin( the1 ) *
cos(
phi1 ) );
187 ptrkg1.setPy( eraw1 *
sin( the1 ) *
sin(
phi1 ) );
188 ptrkg1.setPz( eraw1 *
cos( the1 ) );
189 ptrkg1.setE( eraw1 );
190 ptrkg10 = ptrkg1;
191 ptrkg12 = ptrkg1.boost( -0.011, 0, 0 );
192
193 for (
int n = m + 1;
n < nGam;
n++ )
194 {
195 if ( iGam[
n] == -1 )
continue;
196 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[
n] ) )->emcShower();
197 double eraw2 = g2Trk->
energy();
199 double the2 = g2Trk->
theta();
200 HepLorentzVector ptrkg2, ptrkg20, ptrkg22;
201 ptrkg2.setPx( eraw2 *
sin( the2 ) *
cos(
phi2 ) );
202 ptrkg2.setPy( eraw2 *
sin( the2 ) *
sin(
phi2 ) );
203 ptrkg2.setPz( eraw2 *
cos( the2 ) );
204 ptrkg2.setE( eraw2 );
205 ptrkg20 = ptrkg2;
206 ptrkg22 = ptrkg2.boost( -0.011, 0, 0 );
207
208
209 HepLorentzVector ptrkpi0;
210 ptrkpi0 = ptrkg12 + ptrkg22;
211 double m_xmpi0_tem = ptrkpi0.mag();
212 if ( m_xmpi0_tem > 0.150 || m_xmpi0_tem < 0.115 ) continue;
213
214 bool IsEndcap1 = false;
215 bool IsEndcap2 = false;
216 if ( fabs(
cos( the1 ) ) > 0.86 && fabs(
cos( the1 ) ) < 0.92 ) IsEndcap1 =
true;
217 if ( fabs(
cos( the2 ) ) > 0.86 && fabs(
cos( the2 ) ) < 0.92 ) IsEndcap2 =
true;
218 if ( IsEndcap1 && IsEndcap2 ) continue;
219
220
229
230 double pi0_chisq = kmfit->
chisq( 0 );
231 if ( pi0_chisq >= 2500 ) continue;
232 HepLorentzVector p2gfit1 = kmfit->
pfit( 0 ) + kmfit->
pfit( 1 );
233 p2gfit1.boost( -0.011, 0, 0 );
234
235 for (
int s = 0;
s < nGam - 1;
s++ )
236 {
237 if ( iGam[
s] == -1 )
continue;
238 RecEmcShower* g3Trk = ( *( evtRecTrkCol->begin() + iGam[
s] ) )->emcShower();
239 if ( iGam[
s] == iGam[m] || iGam[
s] == iGam[
n] )
continue;
240 double eraw3 = g3Trk->
energy();
241 double phi3 = g3Trk->
phi();
242 double the3 = g3Trk->
theta();
243 HepLorentzVector ptrkg3, ptrkg30, ptrkg32;
244 ptrkg3.setPx( eraw3 *
sin( the3 ) *
cos( phi3 ) );
245 ptrkg3.setPy( eraw3 *
sin( the3 ) *
sin( phi3 ) );
246 ptrkg3.setPz( eraw3 *
cos( the3 ) );
247 ptrkg3.setE( eraw3 );
248 ptrkg30 = ptrkg3;
249 ptrkg32 = ptrkg3.boost( -0.011, 0, 0 );
250
251 for (
int r =
s + 1; r < nGam; r++ )
252 {
253 if ( iGam[r] == -1 ) continue;
254 RecEmcShower* g4Trk = ( *( evtRecTrkCol->begin() + iGam[r] ) )->emcShower();
255 if ( iGam[r] == iGam[m] || iGam[r] == iGam[
n] )
continue;
256 double eraw4 = g4Trk->
energy();
257 double phi4 = g4Trk->
phi();
258 double the4 = g4Trk->
theta();
259 HepLorentzVector ptrkg4, ptrkg40, ptrkg42;
260 ptrkg4.setPx( eraw4 *
sin( the4 ) *
cos( phi4 ) );
261 ptrkg4.setPy( eraw4 *
sin( the4 ) *
sin( phi4 ) );
262 ptrkg4.setPz( eraw4 *
cos( the4 ) );
263 ptrkg4.setE( eraw4 );
264 ptrkg40 = ptrkg4;
265 ptrkg42 = ptrkg4.boost( -0.011, 0, 0 );
266
267
268 HepLorentzVector ptrkpi0_2;
269 ptrkpi0_2 = ptrkg32 + ptrkg42;
270 double m_xmpi0_2_tem = ptrkpi0_2.mag();
271 if ( m_xmpi0_2_tem > 0.150 || m_xmpi0_2_tem < 0.115 ) continue;
272
273 bool IsEndcap3 = false;
274 bool IsEndcap4 = false;
275 if ( fabs(
cos( the3 ) ) > 0.86 && fabs(
cos( the3 ) ) < 0.92 ) IsEndcap3 =
true;
276 if ( fabs(
cos( the4 ) ) > 0.86 && fabs(
cos( the4 ) ) < 0.92 ) IsEndcap4 =
true;
277 if ( IsEndcap3 && IsEndcap4 ) continue;
278
279
286
289
290 double pi0_2_chisq = kmfit2->
chisq( 0 );
291 if ( pi0_2_chisq >= 2500 ) continue;
292 HepLorentzVector p2gfit2 = kmfit2->
pfit( 0 ) + kmfit2->
pfit( 1 );
293 p2gfit2.boost( -0.011, 0, 0 );
294
295
296 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
297 HepSymMatrix Evx( 3, 0 );
298 double bx = 1E+6;
299 Evx[0][0] = bx * bx;
300 double by = 1E+6;
301 Evx[1][1] = by * by;
302 double bz = 1E+6;
303 Evx[2][2] = bz * bz;
304 VertexParameter vxpar;
307
308
314 if ( !vtxfit->
Fit( 0 ) )
continue;
316
317 WTrackParameter wkam = vtxfit->
wtrk( 0 );
318 WTrackParameter wpip = vtxfit->
wtrk( 1 );
319
320 HepVector kam_val = HepVector( 7, 0 );
322 HepVector pip_val = HepVector( 7, 0 );
324
325 HepLorentzVector P_KAM( kam_val[0], kam_val[1], kam_val[2], kam_val[3] );
326 HepLorentzVector P_PIP( pip_val[0], pip_val[1], pip_val[2], pip_val[3] );
327
328 P_KAM.boost( -0.011, 0, 0 );
329 P_PIP.boost( -0.011, 0, 0 );
330 pddd = P_KAM + P_PIP + p2gfit1 + p2gfit2;
331
332 double pkpipi0pi0 = pddd.rho();
333
334 double temp1 = (
ecms / 2 ) * (
ecms / 2 ) - pkpipi0pi0 * pkpipi0pi0;
335 if ( temp1 < 0 ) temp1 = 0;
336 double mass_bc_tem = sqrt( temp1 );
337 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
338
339 double delE_tag_tag =
ecms / 2 - pddd.e();
340
341 if ( fabs( delE_tag_tag ) < deltaE_tem )
342 {
343 deltaE_tem = fabs( delE_tag_tag );
344 delE_tag_temp = delE_tag_tag;
345 mass_bcgg = mass_bc_tem;
346
347 pddd_temp = pddd;
348 cqtm_temp = m_chargek;
349
350 ika_temp = ika;
351 ipi_temp = ipi;
352
353 iGam1_temp = iGam[m];
354 iGam2_temp = iGam[
n];
355 iGam3_temp = iGam[
s];
356 iGam4_temp = iGam[r];
357 ncount1 = 1;
358 }
359 }
360 }
361 }
362 }
363 }
364 }
365
366 if ( ncount1 == 1 )
367 {
368 tagmode = 25;
369 if ( cqtm_temp < 0 ) tagmode = -25;
370 tagmd = tagmode;
371 mass_bc = mass_bcgg;
372 delE_tag = delE_tag_temp;
373
374 cqtm = -1.0 * cqtm_temp;
375
376 iGoodtag.push_back( ipi_temp );
377 iGoodtag.push_back( ika_temp );
378
379 iGamtag.push_back( iGam1_temp );
380 iGamtag.push_back( iGam2_temp );
381 iGamtag.push_back( iGam3_temp );
382 iGamtag.push_back( iGam4_temp );
383
384 ptag = pddd_temp;
385
386 kpipi0pi0md = true;
387 }
388}
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