136{
137 StatusCode sc=StatusCode::SUCCESS;
138 m_event++;
139
140 MsgStream log(
msgSvc(),name());
141 log<<MSG::INFO<<"in execute()"<<endreq;
142 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
145
146 m_run = eventHeader->runNumber();
147 m_rec = eventHeader->eventNumber();
149 int event = m_rec;
150 int time = eventHeader->time();
151 if(m_event > 1 &&
runNo != m_runNo)
return sc;
154
155 if(m_rec%1000==0)cout<<"Run "<<m_run<<" Event "<<m_rec<<endl;
156
157 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.005, m_ecms);
158 double psipBeta = (p4psip.vect()).mag()/(p4psip.e());
159
160 m_pass[0]+=1;
161
162 int NCharge = evtRecEvent->totalCharged();
163 if(NCharge != 0)return sc;
164
165 m_pass[1]+=1;
166
167 double Emax1=0.0;
168 double Emax2=0.0;
169 int Imax[2];
170
171 if(((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )<2) || ((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )>15)) return sc;
172 m_pass[2]+=1;
173
174 HepLorentzVector p4Gam_1;
175 HepLorentzVector p4Gam_2;
176
177 for(int i=evtRecEvent->totalCharged(); i<evtRecEvent->totalTracks(); i++)
178 {
180 if(!(*itTrk)->isEmcShowerValid()) continue;
182
183 HepLorentzVector p4Gam;
184 double Phi = emcTrk->
phi();
185 double Theta = emcTrk->
theta();
186 double Ener=emcTrk->
energy();
187
188 p4Gam.setPx(Ener *
sin(Theta) *
cos(
Phi));
189 p4Gam.setPy(Ener *
sin(Theta) *
sin(
Phi));
190 p4Gam.setPz(Ener *
cos(Theta));
191 p4Gam.setE(Ener);
192 p4Gam.boost(-0.011, 0 ,-0.005 * 1.0 / 3.686);
193
194 if(Ener>Emax2)
195 {
196 Emax2=Ener;
197 Imax[1]=i;
198 p4Gam_2 = p4Gam;
199 }
200 if(Ener>Emax1)
201 {
202 Emax2=Emax1;
203 p4Gam_2 = p4Gam_1;
204 Imax[1]=Imax[0];
205 Emax1=Ener;
206 p4Gam_1 = p4Gam;
207 Imax[0]=i;
208 }
209 }
210 m_e1_lab = Emax1;
211 m_e2_lab = Emax2;
212 m_e_lab = Emax1 + Emax2;
213 log << MSG::INFO << "Emax1 = " << Emax1 <<"Emax2= "<<Emax2<< endreq;
214 if(Emax1 < m_max1 || Emax2 < m_max2) return sc;
215 m_pass[3]+=1;
216
217
218 double emcphi[2],emctht[2];
219 for(int i = 0;i < 2; i++)
220 {
221 if(i>=evtRecTrkCol->size()) break;
223 if(!(*itTrk)->isEmcShowerValid()) continue;
225 emcphi[i]=emcTrk->
phi();
226 emctht[i]=emcTrk->
theta();
227 }
228 double dltphi_lab = (fabs(emcphi[0] - emcphi[1]) -
pai) * 180.0 /
pai;
229 double dlttheta_lab = (fabs(emctht[0] + emctht[1]) -
pai) * 180.0 /
pai;
230 m_costheta1_lab =
cos(emctht[0]);
231 m_costheta2_lab =
cos(emctht[1]);
232 m_phi1_lab = emcphi[0] * 180.0 /
pai;
233 m_phi2_lab = emcphi[1] * 180.0 /
pai;
234 m_dlttheta_lab = dlttheta_lab;
235 m_dltphi_lab = dltphi_lab;
236
237 if(fabs(m_costheta1_lab) > m_costheta || fabs(m_costheta2_lab) > m_costheta) return sc;
238 m_pass[4]+=1;
239
240
241
242 double px1 = p4Gam_1.px();
243 double py1 = p4Gam_1.py();
244 double pz1 = p4Gam_1.pz();
245 double pxy1 = sqrt(px1 * px1 + py1 * py1);
246 double e1 = p4Gam_1.e();
247
248 double px2 = p4Gam_2.px();
249 double py2 = p4Gam_2.py();
250 double pz2 = p4Gam_2.pz();
251 double pxy2 = sqrt(px2 * px2 + py2 * py2);
252 double e2 = p4Gam_2.e();
253
254
257 if(atan(py1 * 1.0 / px1) > 0){
258 if(px1 > 0)
phi1 = atan(py1 * 1.0 / px1);
259 else phi1 = -(
pai - atan(py1 * 1.0 / px1));
260 }
261 else{
262 if(px1 > 0)
phi1 = atan(py1 * 1.0 / px1);
263 else phi1 =
pai + atan(py1 * 1.0 / px1);
264 }
265
266 if(atan(py2 * 1.0 / px2) > 0){
267 if(px2 > 0)
phi2 = atan(py2 * 1.0 / px2);
268 else phi2 = -(
pai - atan(py2 * 1.0 / px2));
269 }
270 else{
271 if(px2 > 0)
phi2 = atan(py2 * 1.0 / px2);
272 else phi2 =
pai + atan(py2 * 1.0 / px2);
273 }
274
276
279
280 if(pz1 > 0)
theta1 = atan(pxy1 * 1.0 / pz1);
281 else theta1 = (
pai + atan(pxy1 * 1.0 / pz1));
282
283 if(pz2 > 0)
theta2 = atan(pxy2 * 1.0 / pz2);
284 else theta2 = (
pai + atan(pxy2 * 1.0 / pz2));
285
287
288 double xBoost = p4Gam_1.px() + p4Gam_2.px();
289 double yBoost = p4Gam_1.py() + p4Gam_2.py();
290 double zBoost = p4Gam_1.pz() + p4Gam_2.pz();
291 m_xBoost = xBoost;
292 m_yBoost = yBoost;
293 m_zBoost = zBoost;
294
299 m_dlttheta = dlttheta;
304
305
306 if(fabs(m_dltphi) < m_dphi1) Ndata1++;
307 if(fabs(m_dltphi) > m_dphi1 && fabs(m_dltphi) < m_dphi2) Ndata2++;
308
309
310 m_tuple1->write();
311
312
313 return StatusCode::SUCCESS;
314}
double Phi(RecMdcKalTrack *trk)
EvtRecTrackCol::iterator EvtRecTrackIterator
double sin(const BesAngle a)
double cos(const BesAngle a)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol