111 MsgStream log(
msgSvc(), name() );
113 log << MSG::INFO <<
"in initialize()" << endmsg;
115 status = service(
"THistSvc", m_thistsvc );
116 if ( status.isFailure() )
118 log << MSG::INFO <<
"Unable to retrieve pointer to THistSvc" << endmsg;
122 e_z0 =
new TH1F(
"e_z0",
"e_z0", 100, 0, 50 );
123 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/e_z0", e_z0 );
124 e_r0 =
new TH1F(
"e_r0",
"e_r0", 100, 0, 10 );
125 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/e_r0", e_r0 );
127 m_ee_mass =
new TH1F(
"ee_mass",
"ee_mass", 500, m_ecms - 0.5, m_ecms + 0.5 );
128 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_mass", m_ee_mass );
129 m_ee_acoll =
new TH1F(
"ee_acoll",
"ee_acoll", 60, 0, 6 );
130 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_acoll", m_ee_acoll );
131 m_ee_eop_ep =
new TH1F(
"ee_eop_ep",
"ee_eop_ep", 100, 0.4, 1.4 );
132 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eop_ep", m_ee_eop_ep );
133 m_ee_eop_em =
new TH1F(
"ee_eop_em",
"ee_eop_em", 100, 0.4, 1.4 );
134 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eop_em", m_ee_eop_em );
135 m_ee_costheta_ep =
new TH1F(
"ee_costheta_ep",
"ee_costheta_ep", 100, -1, 1 );
136 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_costheta_ep", m_ee_costheta_ep );
137 m_ee_costheta_em =
new TH1F(
"ee_costheta_em",
"ee_costheta_em", 100, -1, 1 );
138 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_costheta_em", m_ee_costheta_em );
140 m_ee_phi_ep =
new TH1F(
"ee_phi_ep",
"ee_phi_ep", 120, -3.2, 3.2 );
141 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_phi_ep", m_ee_phi_ep );
142 m_ee_phi_em =
new TH1F(
"ee_phi_em",
"ee_phi_em", 120, -3.2, 3.2 );
143 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_phi_em", m_ee_phi_em );
145 m_ee_nneu =
new TH1I(
"ee_nneu",
"ee_nneu", 5, 0, 5 );
146 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_nneu", m_ee_nneu );
149 new TH1F(
"ee_eemc_ep",
"ee_eemc_ep", 100, m_ecms / 2. - 0.8, m_ecms / 2. + 0.3 );
150 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eemc_ep", m_ee_eemc_ep );
152 new TH1F(
"ee_eemc_em",
"ee_eemc_em", 100, m_ecms / 2. - 0.8, m_ecms / 2. + 0.3 );
153 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eemc_em", m_ee_eemc_em );
154 m_ee_x_ep =
new TH1F(
"ee_x_ep",
"ee_x_ep", 100, -1.0, 1.0 );
155 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_x_ep", m_ee_x_ep );
156 m_ee_y_ep =
new TH1F(
"ee_y_ep",
"ee_y_ep", 100, -1.0, 1.0 );
157 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_y_ep", m_ee_y_ep );
158 m_ee_z_ep =
new TH1F(
"ee_z_ep",
"ee_z_ep", 100, -10.0, 10.0 );
159 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_z_ep", m_ee_z_ep );
160 m_ee_x_em =
new TH1F(
"ee_x_em",
"ee_x_em", 100, -1.0, 1.0 );
161 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_x_em", m_ee_x_em );
162 m_ee_y_em =
new TH1F(
"ee_y_em",
"ee_y_em", 100, -1.0, 1.0 );
163 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_y_em", m_ee_y_em );
164 m_ee_z_em =
new TH1F(
"ee_z_em",
"ee_z_em", 100, -10.0, 10.0 );
165 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_z_em", m_ee_z_em );
167 m_ee_px_ep =
new TH1F(
"ee_px_ep",
"ee_px_ep", 200, -2.2, 2.2 );
168 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_px_ep", m_ee_px_ep );
169 m_ee_py_ep =
new TH1F(
"ee_py_ep",
"ee_py_ep", 200, -2.2, 2.2 );
170 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_py_ep", m_ee_py_ep );
171 m_ee_pz_ep =
new TH1F(
"ee_pz_ep",
"ee_pz_ep", 200, -2., 2.5 );
172 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pz_ep", m_ee_pz_ep );
173 m_ee_p_ep =
new TH1F(
"ee_p_ep",
"ee_p_ep", 100, m_ecms / 2. - 0.8, m_ecms / 2. + 0.5 );
174 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_p_ep", m_ee_p_ep );
175 m_ee_px_em =
new TH1F(
"ee_px_em",
"ee_px_em", 100, -2.2, 2.2 );
176 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_px_em", m_ee_px_em );
177 m_ee_py_em =
new TH1F(
"ee_py_em",
"ee_py_em", 100, -2.2, 2.2 );
178 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_py_em", m_ee_py_em );
179 m_ee_pz_em =
new TH1F(
"ee_pz_em",
"ee_pz_em", 100, -2.5, 2. );
180 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pz_em", m_ee_pz_em );
181 m_ee_p_em =
new TH1F(
"ee_p_em",
"ee_p_em", 100, m_ecms / 2. - 0.8, m_ecms / 2. + 0.5 );
182 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_p_em", m_ee_p_em );
183 m_ee_deltatof =
new TH1F(
"ee_deltatof",
"ee_deltatof", 50, 0.0, 10.0 );
184 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_deltatof", m_ee_deltatof );
186 m_ee_pidchidedx_ep =
new TH1F(
"ee_pidchidedx_ep",
"ee_pidchidedx_ep", 160, -4, 4 );
187 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchidedx_ep", m_ee_pidchidedx_ep );
188 m_ee_pidchidedx_em =
new TH1F(
"ee_pidchidedx_em",
"ee_pidchidedx_em", 160, -4, 4 );
189 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchidedx_em", m_ee_pidchidedx_em );
190 m_ee_pidchitof1_ep =
new TH1F(
"ee_pidchitof1_ep",
"ee_pidchitof1_ep", 160, -4, 4 );
191 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof1_ep", m_ee_pidchitof1_ep );
192 m_ee_pidchitof1_em =
new TH1F(
"ee_pidchitof1_em",
"ee_pidchitof1_em", 160, -4, 4 );
193 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof1_em", m_ee_pidchitof1_em );
194 m_ee_pidchitof2_ep =
new TH1F(
"ee_pidchitof2_ep",
"ee_pidchitof2_ep", 160, -4, 4 );
195 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof2_ep", m_ee_pidchitof2_ep );
196 m_ee_pidchitof2_em =
new TH1F(
"ee_pidchitof2_em",
"ee_pidchitof2_em", 160, -4, 4 );
197 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof2_em", m_ee_pidchitof2_em );
199 NTuplePtr nt1(
ntupleSvc(),
"DQAFILE/Bhabha" );
200 if ( nt1 ) m_tuple1 = nt1;
203 m_tuple1 =
ntupleSvc()->book(
"DQAFILE/Bhabha", CLID_ColumnWiseTuple,
"N-Tuple" );
206 status = m_tuple1->addItem(
"run", m_run );
207 status = m_tuple1->addItem(
"rec", m_rec );
208 status = m_tuple1->addItem(
"Nchrg", m_ncharg );
209 status = m_tuple1->addItem(
"Nneu", m_nneu, 0, 40 );
210 status = m_tuple1->addItem(
"NGch", m_ngch, 0, 40 );
211 status = m_tuple1->addItem(
"NGam", m_nGam );
213 status = m_tuple1->addItem(
"bhabhatag", m_bhabhatag );
215 status = m_tuple1->addItem(
"acoll", m_acoll );
216 status = m_tuple1->addItem(
"acopl", m_acopl );
217 status = m_tuple1->addItem(
"deltatof", m_deltatof );
218 status = m_tuple1->addItem(
"eop1", m_eop1 );
219 status = m_tuple1->addItem(
"eop2", m_eop2 );
220 status = m_tuple1->addItem(
"eoeb1", m_eoeb1 );
221 status = m_tuple1->addItem(
"eoeb2", m_eoeb2 );
222 status = m_tuple1->addItem(
"poeb1", m_poeb1 );
223 status = m_tuple1->addItem(
"poeb2", m_poeb2 );
224 status = m_tuple1->addItem(
"etoeb1", m_etoeb1 );
225 status = m_tuple1->addItem(
"etoeb2", m_etoeb2 );
226 status = m_tuple1->addItem(
"mucinfo1", m_mucinfo1 );
227 status = m_tuple1->addItem(
"mucinfo2", m_mucinfo2 );
229 status = m_tuple1->addIndexedItem(
"delang", m_nneu, m_delang );
230 status = m_tuple1->addIndexedItem(
"delphi", m_nneu, m_delphi );
231 status = m_tuple1->addIndexedItem(
"delthe", m_nneu, m_delthe );
232 status = m_tuple1->addIndexedItem(
"npart", m_nneu, m_npart );
233 status = m_tuple1->addIndexedItem(
"nemchits", m_nneu, m_nemchits );
234 status = m_tuple1->addIndexedItem(
"module", m_nneu, m_module );
235 status = m_tuple1->addIndexedItem(
"x", m_nneu, m_x );
236 status = m_tuple1->addIndexedItem(
"y", m_nneu, m_y );
237 status = m_tuple1->addIndexedItem(
"z", m_nneu, m_z );
238 status = m_tuple1->addIndexedItem(
"px", m_nneu, m_px );
239 status = m_tuple1->addIndexedItem(
"py", m_nneu, m_py );
240 status = m_tuple1->addIndexedItem(
"pz", m_nneu, m_pz );
241 status = m_tuple1->addIndexedItem(
"theta", m_nneu, m_theta );
242 status = m_tuple1->addIndexedItem(
"phi", m_nneu, m_phi );
243 status = m_tuple1->addIndexedItem(
"dx", m_nneu, m_dx );
244 status = m_tuple1->addIndexedItem(
"dy", m_nneu, m_dy );
245 status = m_tuple1->addIndexedItem(
"dz", m_nneu, m_dz );
246 status = m_tuple1->addIndexedItem(
"dtheta", m_nneu, m_dtheta );
247 status = m_tuple1->addIndexedItem(
"dphi", m_nneu, m_dphi );
248 status = m_tuple1->addIndexedItem(
"energy", m_nneu, m_energy );
249 status = m_tuple1->addIndexedItem(
"dE", m_nneu, m_dE );
250 status = m_tuple1->addIndexedItem(
"eSeed", m_nneu, m_eSeed );
251 status = m_tuple1->addIndexedItem(
"nSeed", m_nneu, m_nSeed );
252 status = m_tuple1->addIndexedItem(
"e3x3", m_nneu, m_e3x3 );
253 status = m_tuple1->addIndexedItem(
"e5x5", m_nneu, m_e5x5 );
254 status = m_tuple1->addIndexedItem(
"secondMoment", m_nneu, m_secondMoment );
255 status = m_tuple1->addIndexedItem(
"latMoment", m_nneu, m_latMoment );
256 status = m_tuple1->addIndexedItem(
"a20Moment", m_nneu, m_a20Moment );
257 status = m_tuple1->addIndexedItem(
"a42Moment", m_nneu, m_a42Moment );
258 status = m_tuple1->addIndexedItem(
"getTime", m_nneu, m_getTime );
259 status = m_tuple1->addIndexedItem(
"getEAll", m_nneu, m_getEAll );
261 status = m_tuple1->addIndexedItem(
"charge", m_ngch, m_charge );
262 status = m_tuple1->addIndexedItem(
"vx", m_ngch, m_vx0 );
263 status = m_tuple1->addIndexedItem(
"vy", m_ngch, m_vy0 );
264 status = m_tuple1->addIndexedItem(
"vz", m_ngch, m_vz0 );
265 status = m_tuple1->addIndexedItem(
"r0", m_ngch, m_vr0 );
267 status = m_tuple1->addIndexedItem(
"px", m_ngch, m_px );
268 status = m_tuple1->addIndexedItem(
"py", m_ngch, m_py );
269 status = m_tuple1->addIndexedItem(
"pz", m_ngch, m_pz );
270 status = m_tuple1->addIndexedItem(
"p", m_ngch, m_p );
272 status = m_tuple1->addIndexedItem(
"kal_vx", m_ngch, m_kal_vx0 );
273 status = m_tuple1->addIndexedItem(
"kal_vy", m_ngch, m_kal_vy0 );
274 status = m_tuple1->addIndexedItem(
"kal_vz", m_ngch, m_kal_vz0 );
276 status = m_tuple1->addIndexedItem(
"kal_px", m_ngch, m_kal_px );
277 status = m_tuple1->addIndexedItem(
"kal_py", m_ngch, m_kal_py );
278 status = m_tuple1->addIndexedItem(
"kal_pz", m_ngch, m_kal_pz );
279 status = m_tuple1->addIndexedItem(
"kal_p", m_ngch, m_kal_p );
281 status = m_tuple1->addIndexedItem(
"probPH", m_ngch, m_probPH );
282 status = m_tuple1->addIndexedItem(
"normPH", m_ngch, m_normPH );
283 status = m_tuple1->addIndexedItem(
"chie", m_ngch, m_chie );
284 status = m_tuple1->addIndexedItem(
"chimu", m_ngch, m_chimu );
285 status = m_tuple1->addIndexedItem(
"chipi", m_ngch, m_chipi );
286 status = m_tuple1->addIndexedItem(
"chik", m_ngch, m_chik );
287 status = m_tuple1->addIndexedItem(
"chip", m_ngch, m_chip );
288 status = m_tuple1->addIndexedItem(
"ghit", m_ngch, m_ghit );
289 status = m_tuple1->addIndexedItem(
"thit", m_ngch, m_thit );
291 status = m_tuple1->addIndexedItem(
"e_emc", m_ngch, m_e_emc );
292 status = m_tuple1->addIndexedItem(
"phi_emc", m_ngch, m_phi_emc );
293 status = m_tuple1->addIndexedItem(
"theta_emc", m_ngch, m_theta_emc );
295 status = m_tuple1->addIndexedItem(
"nhit_muc", m_ngch, m_nhit_muc );
296 status = m_tuple1->addIndexedItem(
"nlay_muc", m_ngch, m_nlay_muc );
297 status = m_tuple1->addIndexedItem(
"t_btof", m_ngch, m_t_btof );
298 status = m_tuple1->addIndexedItem(
"t_etof", m_ngch, m_t_etof );
299 status = m_tuple1->addIndexedItem(
"qual_etof", m_ngch, m_qual_etof );
300 status = m_tuple1->addIndexedItem(
"tof_etof", m_ngch, m_tof_etof );
301 status = m_tuple1->addIndexedItem(
"te_etof", m_ngch, m_te_etof );
302 status = m_tuple1->addIndexedItem(
"tmu_etof", m_ngch, m_tmu_etof );
303 status = m_tuple1->addIndexedItem(
"tpi_etof", m_ngch, m_tpi_etof );
304 status = m_tuple1->addIndexedItem(
"tk_etof", m_ngch, m_tk_etof );
305 status = m_tuple1->addIndexedItem(
"tp_etof", m_ngch, m_tp_etof );
307 status = m_tuple1->addIndexedItem(
"qual_btof1", m_ngch, m_qual_btof1 );
308 status = m_tuple1->addIndexedItem(
"tof_btof1", m_ngch, m_tof_btof1 );
309 status = m_tuple1->addIndexedItem(
"te_btof1", m_ngch, m_te_btof1 );
310 status = m_tuple1->addIndexedItem(
"tmu_btof1", m_ngch, m_tmu_btof1 );
311 status = m_tuple1->addIndexedItem(
"tpi_btof1", m_ngch, m_tpi_btof1 );
312 status = m_tuple1->addIndexedItem(
"tk_btof1", m_ngch, m_tk_btof1 );
313 status = m_tuple1->addIndexedItem(
"tp_btof1", m_ngch, m_tp_btof1 );
315 status = m_tuple1->addIndexedItem(
"qual_btof2", m_ngch, m_qual_btof2 );
316 status = m_tuple1->addIndexedItem(
"tof_btof2", m_ngch, m_tof_btof2 );
317 status = m_tuple1->addIndexedItem(
"te_btof2", m_ngch, m_te_btof2 );
318 status = m_tuple1->addIndexedItem(
"tmu_btof2", m_ngch, m_tmu_btof2 );
319 status = m_tuple1->addIndexedItem(
"tpi_btof2", m_ngch, m_tpi_btof2 );
320 status = m_tuple1->addIndexedItem(
"tk_btof2", m_ngch, m_tk_btof2 );
321 status = m_tuple1->addIndexedItem(
"tp_btof2", m_ngch, m_tp_btof2 );
322 status = m_tuple1->addIndexedItem(
"pidcode", m_ngch, m_pidcode );
323 status = m_tuple1->addIndexedItem(
"pidprob", m_ngch, m_pidprob );
324 status = m_tuple1->addIndexedItem(
"pidchiDedx", m_ngch, m_pidchiDedx );
325 status = m_tuple1->addIndexedItem(
"pidchiTof1", m_ngch, m_pidchiTof1 );
326 status = m_tuple1->addIndexedItem(
"pidchiTof2", m_ngch, m_pidchiTof2 );
328 status = m_tuple1->addItem(
"dedx_GoodHits_ep", m_dedx_goodhits_ep );
329 status = m_tuple1->addItem(
"dedx_chi_ep", m_dedx_chiep );
330 status = m_tuple1->addItem(
"dedx_GoodHits_em", m_dedx_goodhits_em );
331 status = m_tuple1->addItem(
"dedx_chi_em", m_dedx_chiem );
333 status = m_tuple1->addItem(
"px_cms_ep", m_px_cms_ep );
334 status = m_tuple1->addItem(
"py_cms_ep", m_py_cms_ep );
335 status = m_tuple1->addItem(
"pz_cms_ep", m_pz_cms_ep );
336 status = m_tuple1->addItem(
"e_cms_ep", m_e_cms_ep );
337 status = m_tuple1->addItem(
"p_cms_ep", m_p_cms_ep );
339 status = m_tuple1->addItem(
"cos_ep", m_cos_ep );
340 status = m_tuple1->addItem(
"kal_p_ep", m_kal_p_ep );
341 status = m_tuple1->addItem(
"kal_px_ep", m_kal_px_ep );
342 status = m_tuple1->addItem(
"kal_py_ep", m_kal_py_ep );
343 status = m_tuple1->addItem(
"kal_pz_ep", m_kal_pz_ep );
345 status = m_tuple1->addItem(
"px_cms_em", m_px_cms_em );
346 status = m_tuple1->addItem(
"py_cms_em", m_py_cms_em );
347 status = m_tuple1->addItem(
"pz_cms_em", m_pz_cms_em );
348 status = m_tuple1->addItem(
"e_cms_em", m_e_cms_em );
349 status = m_tuple1->addItem(
"p_cms_em", m_p_cms_em );
351 status = m_tuple1->addItem(
"cos_em", m_cos_em );
352 status = m_tuple1->addItem(
"kal_p_em", m_kal_p_em );
353 status = m_tuple1->addItem(
"kal_px_em", m_kal_px_em );
354 status = m_tuple1->addItem(
"kal_py_em", m_kal_py_em );
355 status = m_tuple1->addItem(
"kal_pz_em", m_kal_pz_em );
357 status = m_tuple1->addItem(
"mass_ee", m_mass_ee );
358 status = m_tuple1->addItem(
"px_ee", m_px_ee );
359 status = m_tuple1->addItem(
"py_ee", m_py_ee );
360 status = m_tuple1->addItem(
"pz_ee", m_pz_ee );
361 status = m_tuple1->addItem(
"e_ee", m_e_ee );
362 status = m_tuple1->addItem(
"p_ee", m_p_ee );
364 status = m_tuple1->addItem(
"nep", m_nep );
365 status = m_tuple1->addItem(
"nem", m_nem );
369 log << MSG::ERROR <<
" Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
370 return StatusCode::FAILURE;
378 log << MSG::INFO <<
"successfully return from initialize()" << endmsg;
379 return StatusCode::SUCCESS;
384 const double beamEnergy = m_ecms / 2.;
385 const HepLorentzVector
p_cms( m_ecms *
sin( m_beamangle * 0.5 ), 0.0, 0.0, m_ecms );
386 const Hep3Vector
u_cms = -
p_cms.boostVector();
388 MsgStream log(
msgSvc(), name() );
389 log << MSG::INFO <<
"in execute()" << endmsg;
391 setFilterPassed(
false );
393 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(),
"/Event/EventHeader" );
396 log << MSG::FATAL <<
"Could not find Event Header" << endmsg;
397 return StatusCode::SUCCESS;
400 m_run = eventHeader->runNumber();
401 m_rec = eventHeader->eventNumber();
406 log << MSG::FATAL <<
"Could not find EvtRecEvent" << endmsg;
407 return StatusCode::SUCCESS;
409 log << MSG::INFO <<
"ncharg, nneu, tottks = " << evtRecEvent->totalCharged() <<
" , "
410 << evtRecEvent->totalNeutral() <<
" , " << evtRecEvent->totalTracks() << endmsg;
412 m_ncharg = evtRecEvent->totalCharged();
413 m_nneu = evtRecEvent->totalNeutral();
416 HepSymMatrix Evx( 3, 0 );
420 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc ).ignore();
428 Evx[0][0] = vv[0] * vv[0];
429 Evx[0][1] = vv[0] * vv[1];
430 Evx[1][1] = vv[1] * vv[1];
431 Evx[1][2] = vv[1] * vv[2];
432 Evx[2][2] = vv[2] * vv[2];
439 log << MSG::FATAL <<
"Could not find EvtRecTrackCol" << endmsg;
440 return StatusCode::SUCCESS;
447 for (
int i = 0; i < evtRecEvent->totalCharged(); i++ )
450 if ( !( *itTrk )->isMdcTrackValid() )
continue;
451 if ( !( *itTrk )->isMdcKalTrackValid() )
continue;
454 double pch = mdcTrk->
p();
455 double x0 = mdcTrk->
x();
456 double y0 = mdcTrk->
y();
457 double z0 = mdcTrk->
z();
458 double phi0 = mdcTrk->
helix( 1 );
462 double Rxy = ( x0 - xv ) *
cos( phi0 ) + ( y0 - yv ) *
sin( phi0 );
469 if ( fabs( z0 ) >= m_vz0cut )
continue;
471 if ( fabs( Rxy ) >= m_vr0cut )
continue;
473 iGood.push_back( i );
474 nCharge += mdcTrk->
charge();
477 int nGood = iGood.size();
479 log << MSG::DEBUG <<
"ngood, totcharge = " << nGood <<
" , " << nCharge << endmsg;
481 if ( ( nGood != 2 ) || ( nCharge != 0 ) ) {
return StatusCode::SUCCESS; }
492 Vint ipip, ipim, iep, iem, imup, imum;
502 for (
int i = 0; i < m_ngch; i++ )
506 m_pidchiDedx[i] = -99;
507 m_pidchiTof1[i] = -99;
508 m_pidchiTof2[i] = -99;
531 HepLorentzVector ptrk;
532 ptrk.setPx( mdcTrk->
px() );
533 ptrk.setPy( mdcTrk->
py() );
534 ptrk.setPz( mdcTrk->
pz() );
535 double p3 = ptrk.mag();
538 m_pidprob[i] = pid->
prob( 0 );
539 m_pidchiDedx[i] = pid->
chiDedx( 0 );
540 m_pidchiTof1[i] = pid->
chiTof1( 0 );
541 m_pidchiTof2[i] = pid->
chiTof2( 0 );
542 if ( mdcTrk->
charge() > 0 ) { iep.push_back( iGood[i] ); }
543 if ( mdcTrk->
charge() < 0 ) { iem.push_back( iGood[i] ); }
557 for (
int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
559 if ( i >= evtRecTrkCol->size() )
break;
561 if ( !( *itTrk )->isEmcShowerValid() )
continue;
563 Hep3Vector emcpos( emcTrk->
x(), emcTrk->
y(), emcTrk->
z() );
568 int module = emcTrk->
module();
569 double x = emcTrk->
x();
570 double y = emcTrk->
y();
571 double z = emcTrk->
z();
572 double dx = emcTrk->
dx();
573 double dy = emcTrk->
dy();
574 double dth = emcTrk->
dtheta();
575 double dph = emcTrk->
dphi();
576 double dz = emcTrk->
dz();
578 double dE = emcTrk->
dE();
579 double eSeed = emcTrk->
eSeed();
580 double e3x3 = emcTrk->
e3x3();
581 double e5x5 = emcTrk->
e5x5();
584 double getTime = emcTrk->
time();
585 double getEAll = emcTrk->
getEAll();
591 m_nemchits[iphoton] =
n;
592 m_npart[iphoton] = npart;
593 m_module[iphoton] =
module;
594 m_theta[iphoton] = EmcPos.theta();
595 m_phi[iphoton] = EmcPos.phi();
602 m_dtheta[iphoton] = dth;
603 m_dphi[iphoton] = dph;
604 m_energy[iphoton] =
energy;
606 m_eSeed[iphoton] = eSeed;
607 m_nSeed[iphoton] = nseed;
608 m_e3x3[iphoton] = e3x3;
609 m_e5x5[iphoton] = e5x5;
610 m_secondMoment[iphoton] = secondMoment;
611 m_latMoment[iphoton] = latMoment;
612 m_getTime[iphoton] = getTime;
613 m_getEAll[iphoton] = getEAll;
614 m_a20Moment[iphoton] = a20Moment;
615 m_a42Moment[iphoton] = a42Moment;
622 for (
int j = 0; j < nGood; j++ )
625 if ( !( *jtTrk )->isMdcTrackValid() )
continue;
627 double jtcharge = jtmdcTrk->
charge();
628 if ( !( *jtTrk )->isExtTrackValid() )
continue;
633 double angd = extpos.angle( emcpos );
634 double thed = extpos.theta() - emcpos.theta();
635 double phid = extpos.deltaPhi( emcpos );
636 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
637 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
639 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
640 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
641 if ( angd < dang ) dang = angd;
647 dthe = dthe * 180 / ( CLHEP::pi );
648 dphi = dphi * 180 / ( CLHEP::pi );
649 dang = dang * 180 / ( CLHEP::pi );
650 double eraw = emcTrk->
energy();
652 double phi = emcTrk->
phi();
653 double the = emcTrk->
theta();
655 m_delphi[iphoton] = dphi;
656 m_delthe[iphoton] = dthe;
657 m_delang[iphoton] = dang;
660 if ( fc_theta < 0.80 )
662 if ( eraw <= m_energyThreshold / 2 )
continue;
664 else if ( fc_theta > 0.86 && fc_theta < 0.92 )
666 if ( eraw <= m_energyThreshold )
continue;
670 if ( getTime > m_gammathCut || getTime < m_gammatlCut )
continue;
673 if ( dang < m_gammaTrkCut )
continue;
677 if ( iphoton >= 40 )
return StatusCode::SUCCESS;
680 if ( iphoton > 0 ) counter[4]++;
681 int nGam = iGam.size();
692 for (
int i = 0; i < m_nGam; i++ )
695 if ( !( *itTrk )->isEmcShowerValid() )
continue;
697 double eraw = emcTrk->
energy();
698 double phi = emcTrk->
phi();
699 double the = emcTrk->
theta();
700 HepLorentzVector ptrk;
701 ex_gam += eraw *
sin( the ) *
cos( phi );
702 ey_gam += eraw *
sin( the ) *
sin( phi );
703 ez_gam += eraw *
cos( the );
704 et_gam += eraw *
sin( the );
706 if ( eraw >= egam_ext ) { egam_ext = eraw; }
723 for (
int i = 0; i < m_ngch; i++ )
728 if ( !( *itTrk )->isMdcTrackValid() )
continue;
729 if ( !( *itTrk )->isMdcKalTrackValid() )
continue;
734 m_charge[i] = mdcTrk->
charge();
735 m_vx0[i] = mdcTrk->
x();
736 m_vy0[i] = mdcTrk->
y();
737 m_vz0[i] = mdcTrk->
z();
739 m_px[i] = mdcTrk->
px();
740 m_py[i] = mdcTrk->
py();
741 m_pz[i] = mdcTrk->
pz();
742 m_p[i] = mdcTrk->
p();
753 m_kal_vx0[i] = mdcKalTrk->
dr() *
cos( mdcKalTrk->
fi0() );
755 m_kal_vy0[i] = mdcKalTrk->
dr() *
sin( mdcKalTrk->
fi0() );
757 m_kal_vz0[i] = mdcKalTrk->
dz();
759 m_kal_px[i] = mdcKalTrk->
px();
760 m_kal_py[i] = mdcKalTrk->
py();
761 m_kal_pz[i] = mdcKalTrk->
pz();
763 t_pxy2 = m_kal_px[i] * m_kal_px[i] + m_kal_py[i] * m_kal_py[i];
764 m_kal_p[i] = sqrt( t_pxy2 + m_kal_pz[i] * m_kal_pz[i] );
765 double ptrk = m_kal_p[i];
766 px_had += m_kal_px[i];
767 py_had += m_kal_py[i];
768 pz_had += m_kal_pz[i];
769 pt_had += sqrt( t_pxy2 );
771 e_had += sqrt( m_kal_p[i] * m_kal_p[i] + mdcKalTrk->
mass() * mdcKalTrk->
mass() );
773 if ( ( *itTrk )->isMdcDedxValid() )
777 m_probPH[i] = dedxTrk->
probPH();
778 m_normPH[i] = dedxTrk->
normPH();
780 m_chie[i] = dedxTrk->
chiE();
781 m_chimu[i] = dedxTrk->
chiMu();
782 m_chipi[i] = dedxTrk->
chiPi();
783 m_chik[i] = dedxTrk->
chiK();
784 m_chip[i] = dedxTrk->
chiP();
789 if ( ( *itTrk )->isEmcShowerValid() )
792 m_e_emc[i] = emcTrk->
energy();
793 m_phi_emc[i] = emcTrk->
phi();
794 m_theta_emc[i] = emcTrk->
theta();
800 if ( ( *itTrk )->isMucTrackValid() )
803 m_nhit_muc[i] = mucTrk->
numHits();
807 if ( ( *itTrk )->isTofTrackValid() )
810 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
811 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
813 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
818 status->
setStatus( ( *iter_tof )->status() );
821 if ( ( status->
is_cluster() ) ) m_t_etof[i] = ( *iter_tof )->tof();
824 if ( status )
delete status;
827 if ( status->
layer() != 0 )
829 if ( status )
delete status;
833 double path = ( *iter_tof )->path();
834 double tof = ( *iter_tof )->tof();
835 double ph = ( *iter_tof )->ph();
836 double rhit = ( *iter_tof )->zrhit();
837 double qual = 0.0 + ( *iter_tof )->quality();
838 double cntr = 0.0 + ( *iter_tof )->tofID();
840 for (
int j = 0; j < 5; j++ )
842 double gb = ptrk /
xmass[j];
843 double beta = gb / sqrt( 1 + gb * gb );
844 texp[j] = path / beta /
velc;
846 m_qual_etof[i] = qual;
851 if ( ( status->
is_cluster() ) ) m_t_btof[i] = ( *iter_tof )->tof();
854 if ( status )
delete status;
857 if ( status->
layer() == 1 )
859 double path = ( *iter_tof )->path();
860 double tof = ( *iter_tof )->tof();
861 double ph = ( *iter_tof )->ph();
862 double rhit = ( *iter_tof )->zrhit();
863 double qual = 0.0 + ( *iter_tof )->quality();
864 double cntr = 0.0 + ( *iter_tof )->tofID();
866 for (
int j = 0; j < 5; j++ )
868 double gb = ptrk /
xmass[j];
869 double beta = gb / sqrt( 1 + gb * gb );
870 texp[j] = path / beta /
velc;
872 m_qual_btof1[i] = qual;
873 m_tof_btof1[i] = tof;
876 if ( status->
layer() == 2 )
878 double path = ( *iter_tof )->path();
879 double tof = ( *iter_tof )->tof();
880 double ph = ( *iter_tof )->ph();
881 double rhit = ( *iter_tof )->zrhit();
882 double qual = 0.0 + ( *iter_tof )->quality();
883 double cntr = 0.0 + ( *iter_tof )->tofID();
885 for (
int j = 0; j < 5; j++ )
887 double gb = ptrk /
xmass[j];
888 double beta = gb / sqrt( 1 + gb * gb );
889 texp[j] = path / beta /
velc;
891 m_qual_btof2[i] = qual;
892 m_tof_btof2[i] = tof;
895 if ( status )
delete status;
903 if ( m_ngch != 2 || nCharge != 0 )
return StatusCode::SUCCESS;
909 HepLorentzVector p41e, p42e, p4le;
910 Hep3Vector p31e, p32e, p3le;
911 HepLorentzVector p41m, p42m, p4lm;
912 Hep3Vector p31m, p32m, p3lm;
913 HepLorentzVector p41h, p42h, p4lh;
914 Hep3Vector p31h, p32h, p3lh;
930 double ex1{ 0.0 }, ey1{ 0.0 }, ez1{ 0.0 }, epx1{ 0.0 }, epy1{ 0.0 }, epz1{ 0.0 },
931 epp1{ 0.0 }, pidchidedx1{ 0.0 }, pidchitof11{ 0.0 }, pidchitof21{ 0.0 }, eoeb1{ 0.0 },
932 exoeb1{ 0.0 }, eyoeb1{ 0.0 }, ezoeb1{ 0.0 }, etoeb1{ 0.0 }, kalpp{ 0.0 }, cmsp{ 0.0 };
933 double ex2{ 0.0 }, ey2{ 0.0 }, ez2{ 0.0 }, epx2{ 0.0 }, epy2{ 0.0 }, epz2{ 0.0 },
934 epp2{ 0.0 }, pidchidedx2{ 0.0 }, pidchitof12{ 0.0 }, pidchitof22{ 0.0 }, eoeb2{ 0.0 },
935 exoeb2{ 0.0 }, eyoeb2{ 0.0 }, ezoeb2{ 0.0 }, etoeb2{ 0.0 }, kalpm{ 0.0 }, cmsm{ 0.0 };
937 for (
int i = 0; i < m_ngch; i++ )
939 if ( m_charge[i] > 0 )
941 itTrk1 = evtRecTrkCol->begin() + iGood[i];
942 mdcKalTrk1 = ( *itTrk1 )->mdcKalTrack();
944 if ( !( *itTrk1 )->isMdcDedxValid() )
continue;
945 RecMdcDedx* dedxTrk1 = ( *itTrk1 )->mdcDedx();
948 m_dedx_chiep = dedxTrk1->
chiE();
958 m_px_cms_ep = p41e.px();
959 m_py_cms_ep = p41e.py();
960 m_pz_cms_ep = p41e.pz();
961 m_e_cms_ep = p41e.e();
963 sqrt( p41e.px() * p41e.px() + p41e.py() * p41e.py() + p41e.pz() * p41e.pz() );
966 m_kal_p_ep = m_kal_p[i];
982 pidchidedx1 = m_pidchiDedx[i];
983 pidchitof11 = m_pidchiTof1[i];
984 pidchitof21 = m_pidchiTof2[i];
986 eoeb1 = m_e_emc[i] / beamEnergy;
988 if ( p41e.rho() > 0 ) eope1 = m_e_emc[i] / p41e.rho();
990 exoeb1 = m_e_emc[i] *
sin( m_theta_emc[i] ) *
cos( m_phi_emc[i] ) / beamEnergy;
991 eyoeb1 = m_e_emc[i] *
sin( m_theta_emc[i] ) *
sin( m_phi_emc[i] ) / beamEnergy;
992 ezoeb1 = m_e_emc[i] *
cos( m_theta_emc[i] ) / beamEnergy;
993 etoeb1 = m_e_emc[i] *
sin( m_theta_emc[i] ) / beamEnergy;
995 if ( m_nhit_muc[i] >= 2 && m_nlay_muc[i] >= 2 ) mucinfo1 = 1;
998 if ( m_charge[i] < 0 )
1000 itTrk2 = evtRecTrkCol->begin() + iGood[i];
1001 mdcKalTrk2 = ( *itTrk2 )->mdcKalTrack();
1004 if ( !( *itTrk2 )->isMdcDedxValid() )
continue;
1005 RecMdcDedx* dedxTrk2 = ( *itTrk2 )->mdcDedx();
1008 m_dedx_chiem = dedxTrk2->
chiE();
1013 p42e.boost(
u_cms );
1016 m_px_cms_em = p42e.px();
1017 m_py_cms_em = p42e.py();
1018 m_pz_cms_em = p42e.pz();
1019 m_e_cms_em = p42e.e();
1021 sqrt( p42e.px() * p42e.px() + p42e.py() * p42e.py() + p42e.pz() * p42e.pz() );
1023 m_kal_p_em = m_kal_p[i];
1039 pidchidedx2 = m_pidchiDedx[i];
1040 pidchitof12 = m_pidchiTof1[i];
1041 pidchitof22 = m_pidchiTof2[i];
1043 eoeb2 = m_e_emc[i] / beamEnergy;
1045 if ( p42e.rho() > 0 ) eope2 = m_e_emc[i] / p42e.rho();
1047 exoeb2 = m_e_emc[i] *
sin( m_theta_emc[i] ) *
cos( m_phi_emc[i] ) / beamEnergy;
1048 eyoeb2 = m_e_emc[i] *
sin( m_theta_emc[i] ) *
sin( m_phi_emc[i] ) / beamEnergy;
1049 ezoeb2 = m_e_emc[i] *
cos( m_theta_emc[i] ) / beamEnergy;
1050 etoeb2 = m_e_emc[i] *
sin( m_theta_emc[i] ) / beamEnergy;
1052 if ( m_nhit_muc[i] >= 2 && m_nlay_muc[i] >= 2 ) mucinfo2 = 1;
1056 p4le = ( e01 > e02 ) ? p41e : p42e;
1057 p4lm = ( e01 > e02 ) ? p41m : p42m;
1058 p3le = ( e01 > e02 ) ? p31e : p32e;
1059 p3lm = ( e01 > e02 ) ? p31m : p32m;
1061 double acolle = 180. - p31e.angle( p32e ) * 180.0 / CLHEP::pi;
1062 double acople = 180. - ( p31e.perpPart() ).angle( p32e.perpPart() ) * 180.0 / CLHEP::pi;
1063 double poeb1e = p41e.rho() / beamEnergy;
1064 double poeb2e = p42e.rho() / beamEnergy;
1066 int ilarge = ( e01 > e02 ) ? iip : iim;
1068 double eoebl = m_e_emc[ilarge] / beamEnergy;
1069 if ( p4le.rho() > 0 ) eopl = m_e_emc[ilarge] / p4le.rho();
1071 m_e_emc[ilarge] *
sin( m_theta_emc[ilarge] ) *
cos( m_phi_emc[ilarge] ) / beamEnergy;
1073 m_e_emc[ilarge] *
sin( m_theta_emc[ilarge] ) *
sin( m_phi_emc[ilarge] ) / beamEnergy;
1074 double ezoebl = m_e_emc[ilarge] *
cos( m_theta_emc[ilarge] ) / beamEnergy;
1075 double etoebl = m_e_emc[ilarge] *
sin( m_theta_emc[ilarge] ) / beamEnergy;
1076 int mucinfol = ( m_nhit_muc[ilarge] >= 2 && m_nlay_muc[ilarge] >= 2 ) ? 1 : 0;
1078 int pidel = ( e01 > e02 ) ? m_nep : m_nem;
1080 if ( m_t_btof[iip] * m_t_btof[iim] != 0 ) deltatof = fabs( m_t_btof[iip] - m_t_btof[iim] );
1081 if ( m_t_etof[iip] * m_t_etof[iim] != 0 ) deltatof = fabs( m_t_etof[iip] - m_t_etof[iim] );
1084 if ( m_use_acolle && acolle > m_acoll_e_cut )
return StatusCode::SUCCESS;
1088 if ( m_use_acople && acople > m_acopl_e_cut )
return StatusCode::SUCCESS;
1092 if ( m_use_eoeb && sqrt( ( eoeb1 - 1 ) * ( eoeb1 - 1 ) + ( eoeb2 - 1 ) * ( eoeb2 - 1 ) ) >
1094 return StatusCode::SUCCESS;
1098 if ( m_use_deltatof && m_useTOF && ( deltatof > m_dtof_e_cut ) )
return StatusCode::SUCCESS;
1102 if ( m_use_poeb && poeb1e < m_tpoeb_e_cut && poeb2e < m_tpoeb_e_cut &&
1103 ( sqrt( ( poeb1e - 1 ) * ( poeb1e - 1 ) + ( poeb2e - 1 ) * ( poeb2e - 1 ) ) >
1105 return StatusCode::SUCCESS;
1109 if ( m_use_mucinfo && m_useMUC && ( mucinfo1 != 0 && mucinfo2 != 0 ) )
1110 return StatusCode::SUCCESS;
1114 if ( m_use_ne && m_usePID && ( m_nep != 1 || m_nem != 1 ) )
return StatusCode::SUCCESS;
1123 m_cos_ep = p41e.cosTheta();
1124 m_cos_em = p42e.cosTheta();
1125 m_mass_ee = ( p41e + p42e ).m();
1126 m_px_ee = ( p41e + p42e ).px();
1127 m_py_ee = ( p41e + p42e ).py();
1128 m_pz_ee = ( p41e + p42e ).pz();
1129 m_e_ee = ( p41e + p42e ).e();
1130 m_p_ee = ( p41e + p42e ).rho();
1132 m_deltatof = deltatof;
1139 m_mucinfo1 = mucinfo1;
1140 m_mucinfo2 = mucinfo2;
1150 ( *itTrk1 )->tagElectron();
1151 ( *itTrk2 )->tagElectron();
1157 ( *itTrk1 )->setQuality( 0 );
1158 ( *itTrk2 )->setQuality( 0 );
1163 setFilterPassed(
true );
1165 m_ee_mass->Fill( ( p41e + p42e ).m() );
1166 m_ee_acoll->Fill( acolle );
1167 m_ee_eop_ep->Fill( eope1 );
1168 m_ee_eop_em->Fill( eope2 );
1169 m_ee_costheta_ep->Fill( p41e.cosTheta() );
1170 m_ee_costheta_em->Fill( p42e.cosTheta() );
1172 m_ee_phi_ep->Fill( p41e.phi() );
1173 m_ee_phi_em->Fill( p42e.phi() );
1175 m_ee_nneu->Fill( m_nGam );
1177 m_ee_eemc_ep->Fill( e01 );
1178 m_ee_eemc_em->Fill( e02 );
1180 m_ee_x_ep->Fill( ex1 );
1181 m_ee_y_ep->Fill( ey1 );
1182 m_ee_z_ep->Fill( ez1 );
1184 m_ee_x_em->Fill( ex2 );
1185 m_ee_y_em->Fill( ey2 );
1186 m_ee_z_em->Fill( ez2 );
1188 m_ee_px_ep->Fill( epx1 );
1189 m_ee_py_ep->Fill( epy1 );
1190 m_ee_pz_ep->Fill( epz1 );
1191 m_ee_p_ep->Fill( epp1 );
1193 m_ee_px_em->Fill( epx2 );
1194 m_ee_py_em->Fill( epy2 );
1195 m_ee_pz_em->Fill( epz2 );
1196 m_ee_p_em->Fill( epp2 );
1198 m_ee_deltatof->Fill( deltatof );
1200 m_ee_pidchidedx_ep->Fill( pidchidedx1 );
1201 m_ee_pidchidedx_em->Fill( pidchidedx2 );
1202 m_ee_pidchitof1_ep->Fill( pidchitof11 );
1203 m_ee_pidchitof1_em->Fill( pidchitof12 );
1204 m_ee_pidchitof2_ep->Fill( pidchitof21 );
1205 m_ee_pidchitof2_em->Fill( pidchitof22 );
1207 if ( m_writentuple ) m_tuple1->write().ignore();
1209 return StatusCode::SUCCESS;