38{
39
40 int nGood=iGood.size();
41 int nGam=iGam.size();
42
43 iGoodtag.clear();
44 iGamtag.clear();
45
46 double mass_bcgg,delE_tag_temp;
47 int m_chargetag, m_chargek,m_chargepi;
48 int ika_temp,ipi_temp, iGam1_temp, iGam2_temp, iGam3_temp, iGam4_temp;
49 HepLorentzVector pddd;
50 HepLorentzVector pddd_temp;
51
52 int cqtm_temp;
53 IDataProviderSvc* eventSvc = NULL;
54 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
56 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
57
58 int runNo=eventHeader->runNumber();
59 int rec=eventHeader->eventNumber();
60
61 double xecm=2*Ebeam;
62
63 kpipi0pi0md=false;
64 double tagmode=0;
65
66 if((evtRecEvent->totalCharged() < 2||nGam<4)){ return; }
67
69
70
72 Gaudi::svcLocator()->service("SimplePIDSvc", simple_pid);
73
74 double deltaE_tem = 0.20;
75 int ncount1 = 0;
76
77 Hep3Vector xorigin(0,0,0);
79 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
81 {
84 xorigin.setX(dbv[0]);
85 xorigin.setY(dbv[1]);
86 xorigin.setZ(dbv[2]);
87 }
88
89 double xv=xorigin.x();
90 double yv=xorigin.y();
91 double zv=xorigin.z();
92
94 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
95
96 HepLorentzVector p2gfit1, p2gfit2;
97 HepLorentzVector p2gg;
98
99 HepLorentzVector ptrk1_temp, ptrk2_temp, ptrk3_temp, ptrk4_temp, ptrk5_temp, ptrk6_temp, ptrk7_temp, ptrk8_temp;
100 for(int i = 0; i < evtRecEvent->totalCharged(); i++) {
102
103 int ika= (*itTrk)->trackId();
104
105 if(!(*itTrk)->isMdcKalTrackValid()) continue;
108
109 m_chargek=mdcKalTrk1->
charge();
110 if(Charge_candidate_D != 0) {
111 if(m_chargek != -Charge_candidate_D) continue;
112 }
113 if(Charge_candidate_D == 0) {
114 if(
abs(m_chargek) != 1)
continue;
115 }
116
119 VFHelix helixip3_1(point0,a1,Ea1);
120 helixip3_1.pivot(IP);
121 HepVector vecipa1 = helixip3_1.a();
122
123 double dr1 = fabs(vecipa1[0]);
124 double dz1 = fabs(vecipa1[3]);
125 double costheta1 =
cos(mdcKalTrk1->
theta());
126 if ( dr1 >= 1.0) continue;
127 if ( dz1 >= 10.0) continue;
128 if ( fabs(costheta1) >= 0.93) continue;
129
130 if(PID_flag == 5) {
133 }
134
135
136
138
139
140
141
142 for(int j = 0; j< evtRecEvent->totalCharged(); j++) {
144
145 int ipi= (*itTrk)->trackId();
146 if(ipi==ika) continue;
147
148 if(!(*itTrk)->isMdcKalTrackValid()) continue;
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) {
171 }
172
173
175
176 for(int m = 0; m < nGam-1; m++) {
177 if(iGam[m]==-1) continue;
178 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[m]))->emcShower();
179 double eraw1 = g1Trk->
energy();
180 double phi1 = g1Trk->
phi();
181 double the1 = g1Trk->
theta();
182 HepLorentzVector ptrkg1,ptrkg10,ptrkg12;
183 ptrkg1.setPx(eraw1*
sin(the1)*
cos(phi1));
184 ptrkg1.setPy(eraw1*
sin(the1)*
sin(phi1));
185 ptrkg1.setPz(eraw1*
cos(the1));
186 ptrkg1.setE(eraw1);
187 ptrkg10 = ptrkg1;
188 ptrkg12 = ptrkg1.boost(-0.011,0,0);
189
190 for(int n = m+1; n < nGam; n++) {
191 if(iGam[n]==-1) continue;
192 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[n]))->emcShower();
193 double eraw2 = g2Trk->
energy();
194 double phi2 = g2Trk->
phi();
195 double the2 = g2Trk->
theta();
196 HepLorentzVector ptrkg2,ptrkg20,ptrkg22;
197 ptrkg2.setPx(eraw2*
sin(the2)*
cos(phi2));
198 ptrkg2.setPy(eraw2*
sin(the2)*
sin(phi2));
199 ptrkg2.setPz(eraw2*
cos(the2));
200 ptrkg2.setE(eraw2);
201 ptrkg20 = ptrkg2;
202 ptrkg22 = ptrkg2.boost(-0.011,0,0);
203
204
205 HepLorentzVector ptrkpi0;
206 ptrkpi0 = ptrkg12+ptrkg22;
207 double m_xmpi0_tem = ptrkpi0.mag();
208 if(m_xmpi0_tem>0.150||m_xmpi0_tem<0.115) continue;
209
210 bool IsEndcap1 = false; bool IsEndcap2 = false;
211 if(fabs(
cos(the1)) > 0.86 && fabs(
cos(the1)) < 0.92) IsEndcap1 =
true;
212 if(fabs(
cos(the2)) > 0.86 && fabs(
cos(the2)) < 0.92) IsEndcap2 =
true;
213 if(IsEndcap1 && IsEndcap2) continue;
214
215
224
225 double pi0_chisq = kmfit->
chisq(0);
226 if ( pi0_chisq >= 2500) continue;
227 HepLorentzVector p2gfit1 = kmfit->
pfit(0) + kmfit->
pfit(1);
228 p2gfit1.boost(-0.011,0,0);
229
230 for(
int s = 0;
s < nGam-1;
s++) {
231 if(iGam[
s]==-1)
continue;
232 RecEmcShower *g3Trk = (*(evtRecTrkCol->begin()+iGam[
s]))->emcShower();
233 if(iGam[
s] == iGam[m] || iGam[
s] == iGam[n] )
continue;
234 double eraw3 = g3Trk->
energy();
235 double phi3 = g3Trk->
phi();
236 double the3 = g3Trk->
theta();
237 HepLorentzVector ptrkg3,ptrkg30,ptrkg32;
238 ptrkg3.setPx(eraw3*
sin(the3)*
cos(phi3));
239 ptrkg3.setPy(eraw3*
sin(the3)*
sin(phi3));
240 ptrkg3.setPz(eraw3*
cos(the3));
241 ptrkg3.setE(eraw3);
242 ptrkg30 = ptrkg3;
243 ptrkg32 = ptrkg3.boost(-0.011,0,0);
244
245 for(
int r =
s+1; r < nGam; r++) {
246 if(iGam[r]==-1) continue;
247 RecEmcShower *g4Trk = (*(evtRecTrkCol->begin()+iGam[r]))->emcShower();
248 if(iGam[r] == iGam[m] || iGam[r] == iGam[n] ) continue;
249 double eraw4 = g4Trk->
energy();
250 double phi4 = g4Trk->
phi();
251 double the4 = g4Trk->
theta();
252 HepLorentzVector ptrkg4,ptrkg40,ptrkg42;
253 ptrkg4.setPx(eraw4*
sin(the4)*
cos(phi4));
254 ptrkg4.setPy(eraw4*
sin(the4)*
sin(phi4));
255 ptrkg4.setPz(eraw4*
cos(the4));
256 ptrkg4.setE(eraw4);
257 ptrkg40 = ptrkg4;
258 ptrkg42 = ptrkg4.boost(-0.011,0,0);
259
260
261 HepLorentzVector ptrkpi0_2;
262 ptrkpi0_2 = ptrkg32+ptrkg42;
263 double m_xmpi0_2_tem = ptrkpi0_2.mag();
264 if(m_xmpi0_2_tem>0.150||m_xmpi0_2_tem<0.115) continue;
265
266 bool IsEndcap3 = false; bool IsEndcap4 = false;
267 if(fabs(
cos(the3)) > 0.86 && fabs(
cos(the3)) < 0.92) IsEndcap3 =
true;
268 if(fabs(
cos(the4)) > 0.86 && fabs(
cos(the4)) < 0.92) IsEndcap4 =
true;
269 if(IsEndcap3 && IsEndcap4) continue;
270
271
278
281
282 double pi0_2_chisq = kmfit2->
chisq(0);
283 if ( pi0_2_chisq >= 2500) continue;
284 HepLorentzVector p2gfit2 = kmfit2->
pfit(0) + kmfit2->
pfit(1);
285 p2gfit2.boost(-0.011,0,0);
286
287
288
289 HepPoint3D vx(xorigin.x(), xorigin.y(), xorigin.z());
290 HepSymMatrix Evx(3, 0);
291 double bx = 1E+6; Evx[0][0] = bx*bx;
292 double by = 1E+6; Evx[1][1] = by*by;
293 double bz = 1E+6; Evx[2][2] = bz*bz;
295
296
302 if(!vtxfit->
Fit(0))
continue;
304
307
308 HepVector kam_val = HepVector(7,0);
310 HepVector pip_val = HepVector(7,0);
312
313 HepLorentzVector P_KAM(kam_val[0],kam_val[1],kam_val[2],kam_val[3]);
314 HepLorentzVector P_PIP(pip_val[0],pip_val[1],pip_val[2],pip_val[3]);
315
316 P_KAM.boost(-0.011,0,0);
317 P_PIP.boost(-0.011,0,0);
318 pddd = P_KAM + P_PIP + p2gfit1 + p2gfit2;
319
320 double pkpipi0pi0=pddd.rho();
321
322 double temp1 = (
ecms/2)*(
ecms/2)-pkpipi0pi0*pkpipi0pi0 ;
323 if(temp1<0) temp1 =0;
324 double mass_bc_tem = sqrt(temp1);
325 if(mass_bc_tem < 1.82 || mass_bc_tem > 1.89) continue;
326
327 double delE_tag_tag =
ecms/2-pddd.e();
328
329 if(fabs(delE_tag_tag)<deltaE_tem) {
330 deltaE_tem = fabs(delE_tag_tag);
331 delE_tag_temp = delE_tag_tag;
332 mass_bcgg = mass_bc_tem;
333
334 pddd_temp = pddd;
335 cqtm_temp = m_chargek;
336
337 ika_temp=ika;
338 ipi_temp=ipi;
339
340 iGam1_temp = iGam[m];
341 iGam2_temp = iGam[n];
342 iGam3_temp = iGam[
s];
343 iGam4_temp = iGam[r];
344 ncount1 = 1;
345
346 }
347 }
348 }
349 }
350 }
351 }
352 }
353
354 if(ncount1 == 1){
355 tagmode=25;
356 if(cqtm_temp<0) tagmode=-25;
357 tagmd=tagmode;
358 mass_bc = mass_bcgg;
359 delE_tag = delE_tag_temp;
360
361 cqtm = -1.0*cqtm_temp;
362
363 iGoodtag.push_back(ipi_temp);
364 iGoodtag.push_back(ika_temp);
365
366 iGamtag.push_back(iGam1_temp);
367 iGamtag.push_back(iGam2_temp);
368 iGamtag.push_back(iGam3_temp);
369 iGamtag.push_back(iGam4_temp);
370
371 ptag = pddd_temp;
372
373 kpipi0pi0md = true;
374
375 }
376
377
378}
double sin(const BesAngle a)
double cos(const BesAngle a)
EvtRecTrackCol::iterator EvtRecTrackIterator
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.9.6