127 MsgStream log(
msgSvc(), name());
129 log << MSG::INFO <<
"in initialize()" << endmsg;
131 status = service(
"THistSvc", m_thistsvc);
132 if(status.isFailure() ){
133 log << MSG::INFO <<
"Unable to retrieve pointer to THistSvc" << endreq;
138 e_z0 =
new TH1F(
"e_z0",
"e_z0",100,0,50);
139 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/e_z0", e_z0);
140 e_r0 =
new TH1F(
"e_r0",
"e_r0",100,0,10);
141 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/e_r0", e_r0);
143 m_ee_mass =
new TH1F(
"ee_mass",
"ee_mass", 500, m_ecms-0.5, m_ecms+0.5 );
144 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_mass", m_ee_mass);
145 m_ee_acoll =
new TH1F(
"ee_acoll",
"ee_acoll", 60, 0, 6 );
146 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_acoll", m_ee_acoll);
147 m_ee_eop_ep =
new TH1F(
"ee_eop_ep",
"ee_eop_ep", 100,0.4,1.4 );
148 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eop_ep", m_ee_eop_ep);
149 m_ee_eop_em =
new TH1F(
"ee_eop_em",
"ee_eop_em", 100,0.4,1.4 );
150 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eop_em", m_ee_eop_em);
151 m_ee_costheta_ep =
new TH1F(
"ee_costheta_ep",
"ee_costheta_ep", 100,-1,1 );
152 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_costheta_ep", m_ee_costheta_ep);
153 m_ee_costheta_em =
new TH1F(
"ee_costheta_em",
"ee_costheta_em", 100,-1,1 );
154 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_costheta_em", m_ee_costheta_em);
156 m_ee_phi_ep =
new TH1F(
"ee_phi_ep",
"ee_phi_ep", 120,-3.2,3.2 );
157 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_phi_ep", m_ee_phi_ep);
158 m_ee_phi_em =
new TH1F(
"ee_phi_em",
"ee_phi_em", 120,-3.2,3.2 );
159 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_phi_em", m_ee_phi_em);
161 m_ee_nneu =
new TH1I(
"ee_nneu",
"ee_nneu", 5,0,5 );
162 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_nneu", m_ee_nneu);
166 m_ee_eemc_ep=
new TH1F(
"ee_eemc_ep",
"ee_eemc_ep",100,m_ecms/2.-0.8,m_ecms/2.+0.3);
167 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eemc_ep", m_ee_eemc_ep);
168 m_ee_eemc_em=
new TH1F(
"ee_eemc_em",
"ee_eemc_em",100,m_ecms/2.-0.8,m_ecms/2.+0.3);
169 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_eemc_em", m_ee_eemc_em);
170 m_ee_x_ep=
new TH1F(
"ee_x_ep",
"ee_x_ep",100,-1.0,1.0);
171 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_x_ep", m_ee_x_ep);
172 m_ee_y_ep=
new TH1F(
"ee_y_ep",
"ee_y_ep",100,-1.0,1.0);
173 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_y_ep", m_ee_y_ep);
174 m_ee_z_ep=
new TH1F(
"ee_z_ep",
"ee_z_ep",100,-10.0,10.0);
175 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_z_ep", m_ee_z_ep);
176 m_ee_x_em=
new TH1F(
"ee_x_em",
"ee_x_em",100,-1.0,1.0);
177 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_x_em", m_ee_x_em);
178 m_ee_y_em=
new TH1F(
"ee_y_em",
"ee_y_em",100,-1.0,1.0);
179 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_y_em", m_ee_y_em);
180 m_ee_z_em=
new TH1F(
"ee_z_em",
"ee_z_em",100,-10.0,10.0);
181 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_z_em", m_ee_z_em);
183 m_ee_px_ep=
new TH1F(
"ee_px_ep",
"ee_px_ep",200,-2.2,2.2);
184 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_px_ep", m_ee_px_ep);
185 m_ee_py_ep=
new TH1F(
"ee_py_ep",
"ee_py_ep",200,-2.2,2.2);
186 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_py_ep", m_ee_py_ep);
187 m_ee_pz_ep=
new TH1F(
"ee_pz_ep",
"ee_pz_ep",200,-2.,2.5);
188 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pz_ep", m_ee_pz_ep);
189 m_ee_p_ep=
new TH1F(
"ee_p_ep",
"ee_p_ep",100, m_ecms/2.-0.8 , m_ecms/2.+0.5 );
190 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_p_ep", m_ee_p_ep);
191 m_ee_px_em=
new TH1F(
"ee_px_em",
"ee_px_em",100,-2.2,2.2);
192 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_px_em", m_ee_px_em);
193 m_ee_py_em=
new TH1F(
"ee_py_em",
"ee_py_em",100,-2.2,2.2);
194 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_py_em", m_ee_py_em);
195 m_ee_pz_em=
new TH1F(
"ee_pz_em",
"ee_pz_em",100,-2.5,2.);
196 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pz_em", m_ee_pz_em);
197 m_ee_p_em=
new TH1F(
"ee_p_em",
"ee_p_em",100,m_ecms/2.-0.8,m_ecms/2.+0.5);
198 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_p_em", m_ee_p_em);
199 m_ee_deltatof=
new TH1F(
"ee_deltatof",
"ee_deltatof",50,0.0,10.0);
200 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_deltatof", m_ee_deltatof);
202 m_ee_pidchidedx_ep=
new TH1F(
"ee_pidchidedx_ep",
"ee_pidchidedx_ep",160,-4,4);
203 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchidedx_ep", m_ee_pidchidedx_ep);
204 m_ee_pidchidedx_em=
new TH1F(
"ee_pidchidedx_em",
"ee_pidchidedx_em",160,-4,4);
205 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchidedx_em", m_ee_pidchidedx_em);
206 m_ee_pidchitof1_ep=
new TH1F(
"ee_pidchitof1_ep",
"ee_pidchitof1_ep",160,-4,4);
207 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof1_ep", m_ee_pidchitof1_ep);
208 m_ee_pidchitof1_em=
new TH1F(
"ee_pidchitof1_em",
"ee_pidchitof1_em",160,-4,4);
209 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof1_em", m_ee_pidchitof1_em);
210 m_ee_pidchitof2_ep=
new TH1F(
"ee_pidchitof2_ep",
"ee_pidchitof2_ep",160,-4,4);
211 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof2_ep", m_ee_pidchitof2_ep);
212 m_ee_pidchitof2_em=
new TH1F(
"ee_pidchitof2_em",
"ee_pidchitof2_em",160,-4,4);
213 status = m_thistsvc->regHist(
"/DQAHist/Bhabha/ee_pidchitof2_em", m_ee_pidchitof2_em);
217 NTuplePtr nt1(
ntupleSvc(),
"DQAFILE/Bhabha");
218 if ( nt1 ) m_tuple1 = nt1;
220 m_tuple1 =
ntupleSvc()->book (
"DQAFILE/Bhabha", CLID_ColumnWiseTuple,
"N-Tuple");
222 status = m_tuple1->addItem (
"run", m_run);
223 status = m_tuple1->addItem (
"rec", m_rec);
224 status = m_tuple1->addItem (
"Nchrg", m_ncharg);
225 status = m_tuple1->addItem (
"Nneu", m_nneu,0,40);
226 status = m_tuple1->addItem (
"NGch", m_ngch, 0, 40);
227 status = m_tuple1->addItem (
"NGam", m_nGam);
230 status = m_tuple1->addItem (
"bhabhatag", m_bhabhatag);
232 status = m_tuple1->addItem (
"acoll", m_acoll);
233 status = m_tuple1->addItem (
"acopl", m_acopl);
234 status = m_tuple1->addItem (
"deltatof", m_deltatof);
235 status = m_tuple1->addItem (
"eop1", m_eop1);
236 status = m_tuple1->addItem (
"eop2", m_eop2);
237 status = m_tuple1->addItem (
"eoeb1", m_eoeb1);
238 status = m_tuple1->addItem (
"eoeb2", m_eoeb2);
239 status = m_tuple1->addItem (
"poeb1", m_poeb1);
240 status = m_tuple1->addItem (
"poeb2", m_poeb2);
241 status = m_tuple1->addItem (
"etoeb1", m_etoeb1);
242 status = m_tuple1->addItem (
"etoeb2", m_etoeb2);
243 status = m_tuple1->addItem (
"mucinfo1", m_mucinfo1);
244 status = m_tuple1->addItem (
"mucinfo2", m_mucinfo2);
248 status = m_tuple1->addIndexedItem (
"delang",m_nneu, m_delang);
249 status = m_tuple1->addIndexedItem (
"delphi",m_nneu, m_delphi);
250 status = m_tuple1->addIndexedItem (
"delthe",m_nneu, m_delthe);
251 status = m_tuple1->addIndexedItem (
"npart",m_nneu, m_npart);
252 status = m_tuple1->addIndexedItem (
"nemchits",m_nneu, m_nemchits);
253 status = m_tuple1->addIndexedItem (
"module",m_nneu, m_module);
254 status = m_tuple1->addIndexedItem (
"x",m_nneu, m_x);
255 status = m_tuple1->addIndexedItem (
"y",m_nneu, m_y);
256 status = m_tuple1->addIndexedItem (
"z",m_nneu, m_z);
257 status = m_tuple1->addIndexedItem (
"px",m_nneu, m_px);
258 status = m_tuple1->addIndexedItem (
"py",m_nneu, m_py);
259 status = m_tuple1->addIndexedItem (
"pz",m_nneu, m_pz);
260 status = m_tuple1->addIndexedItem (
"theta",m_nneu, m_theta);
261 status = m_tuple1->addIndexedItem (
"phi",m_nneu, m_phi);
262 status = m_tuple1->addIndexedItem (
"dx",m_nneu, m_dx);
263 status = m_tuple1->addIndexedItem (
"dy",m_nneu, m_dy);
264 status = m_tuple1->addIndexedItem (
"dz",m_nneu, m_dz);
265 status = m_tuple1->addIndexedItem (
"dtheta",m_nneu, m_dtheta);
266 status = m_tuple1->addIndexedItem (
"dphi",m_nneu, m_dphi);
267 status = m_tuple1->addIndexedItem (
"energy",m_nneu, m_energy);
268 status = m_tuple1->addIndexedItem (
"dE",m_nneu, m_dE);
269 status = m_tuple1->addIndexedItem (
"eSeed",m_nneu, m_eSeed);
270 status = m_tuple1->addIndexedItem (
"nSeed",m_nneu, m_nSeed);
271 status = m_tuple1->addIndexedItem (
"e3x3",m_nneu, m_e3x3);
272 status = m_tuple1->addIndexedItem (
"e5x5",m_nneu, m_e5x5);
273 status = m_tuple1->addIndexedItem (
"secondMoment",m_nneu, m_secondMoment);
274 status = m_tuple1->addIndexedItem (
"latMoment",m_nneu, m_latMoment);
275 status = m_tuple1->addIndexedItem (
"a20Moment",m_nneu, m_a20Moment);
276 status = m_tuple1->addIndexedItem (
"a42Moment",m_nneu, m_a42Moment);
277 status = m_tuple1->addIndexedItem (
"getTime",m_nneu, m_getTime);
278 status = m_tuple1->addIndexedItem (
"getEAll",m_nneu, m_getEAll);
282 status = m_tuple1->addIndexedItem(
"charge", m_ngch, m_charge);
283 status = m_tuple1->addIndexedItem (
"vx", m_ngch, m_vx0);
284 status = m_tuple1->addIndexedItem (
"vy", m_ngch, m_vy0);
285 status = m_tuple1->addIndexedItem (
"vz", m_ngch, m_vz0);
286 status = m_tuple1->addIndexedItem (
"r0", m_ngch, m_vr0);
288 status = m_tuple1->addIndexedItem (
"px", m_ngch, m_px) ;
289 status = m_tuple1->addIndexedItem (
"py", m_ngch, m_py) ;
290 status = m_tuple1->addIndexedItem (
"pz", m_ngch, m_pz) ;
291 status = m_tuple1->addIndexedItem (
"p", m_ngch, m_p) ;
295 status = m_tuple1->addIndexedItem (
"kal_vx", m_ngch, m_kal_vx0);
296 status = m_tuple1->addIndexedItem (
"kal_vy", m_ngch, m_kal_vy0);
297 status = m_tuple1->addIndexedItem (
"kal_vz", m_ngch, m_kal_vz0);
300 status = m_tuple1->addIndexedItem (
"kal_px", m_ngch, m_kal_px) ;
301 status = m_tuple1->addIndexedItem (
"kal_py", m_ngch, m_kal_py) ;
302 status = m_tuple1->addIndexedItem (
"kal_pz", m_ngch, m_kal_pz) ;
303 status = m_tuple1->addIndexedItem (
"kal_p", m_ngch, m_kal_p) ;
306 status = m_tuple1->addIndexedItem (
"probPH" , m_ngch, m_probPH) ;
307 status = m_tuple1->addIndexedItem (
"normPH" , m_ngch, m_normPH) ;
308 status = m_tuple1->addIndexedItem (
"chie" , m_ngch, m_chie) ;
309 status = m_tuple1->addIndexedItem (
"chimu" , m_ngch, m_chimu) ;
310 status = m_tuple1->addIndexedItem (
"chipi" , m_ngch, m_chipi) ;
311 status = m_tuple1->addIndexedItem (
"chik" , m_ngch, m_chik) ;
312 status = m_tuple1->addIndexedItem (
"chip" , m_ngch, m_chip) ;
313 status = m_tuple1->addIndexedItem (
"ghit" , m_ngch, m_ghit) ;
314 status = m_tuple1->addIndexedItem (
"thit" , m_ngch, m_thit) ;
316 status = m_tuple1->addIndexedItem (
"e_emc" , m_ngch, m_e_emc) ;
317 status = m_tuple1->addIndexedItem (
"phi_emc" , m_ngch, m_phi_emc) ;
318 status = m_tuple1->addIndexedItem (
"theta_emc" , m_ngch, m_theta_emc) ;
320 status = m_tuple1->addIndexedItem (
"nhit_muc" , m_ngch, m_nhit_muc) ;
321 status = m_tuple1->addIndexedItem (
"nlay_muc" , m_ngch, m_nlay_muc) ;
322 status = m_tuple1->addIndexedItem (
"t_btof" , m_ngch, m_t_btof );
323 status = m_tuple1->addIndexedItem (
"t_etof" , m_ngch, m_t_etof );
324 status = m_tuple1->addIndexedItem (
"qual_etof" , m_ngch, m_qual_etof );
325 status = m_tuple1->addIndexedItem (
"tof_etof" , m_ngch, m_tof_etof );
326 status = m_tuple1->addIndexedItem (
"te_etof" , m_ngch, m_te_etof );
327 status = m_tuple1->addIndexedItem (
"tmu_etof" , m_ngch, m_tmu_etof );
328 status = m_tuple1->addIndexedItem (
"tpi_etof" , m_ngch, m_tpi_etof );
329 status = m_tuple1->addIndexedItem (
"tk_etof" , m_ngch, m_tk_etof );
330 status = m_tuple1->addIndexedItem (
"tp_etof" , m_ngch, m_tp_etof );
332 status = m_tuple1->addIndexedItem (
"qual_btof1", m_ngch, m_qual_btof1 );
333 status = m_tuple1->addIndexedItem (
"tof_btof1" , m_ngch, m_tof_btof1 );
334 status = m_tuple1->addIndexedItem (
"te_btof1" , m_ngch, m_te_btof1 );
335 status = m_tuple1->addIndexedItem (
"tmu_btof1" , m_ngch, m_tmu_btof1 );
336 status = m_tuple1->addIndexedItem (
"tpi_btof1" , m_ngch, m_tpi_btof1 );
337 status = m_tuple1->addIndexedItem (
"tk_btof1" , m_ngch, m_tk_btof1 );
338 status = m_tuple1->addIndexedItem (
"tp_btof1" , m_ngch, m_tp_btof1 );
340 status = m_tuple1->addIndexedItem (
"qual_btof2", m_ngch, m_qual_btof2 );
341 status = m_tuple1->addIndexedItem (
"tof_btof2" , m_ngch, m_tof_btof2 );
342 status = m_tuple1->addIndexedItem (
"te_btof2" , m_ngch, m_te_btof2 );
343 status = m_tuple1->addIndexedItem (
"tmu_btof2" , m_ngch, m_tmu_btof2 );
344 status = m_tuple1->addIndexedItem (
"tpi_btof2" , m_ngch, m_tpi_btof2 );
345 status = m_tuple1->addIndexedItem (
"tk_btof2" , m_ngch, m_tk_btof2 );
346 status = m_tuple1->addIndexedItem (
"tp_btof2" , m_ngch, m_tp_btof2 );
347 status = m_tuple1->addIndexedItem (
"pidcode" , m_ngch, m_pidcode);
348 status = m_tuple1->addIndexedItem (
"pidprob" , m_ngch, m_pidprob);
349 status = m_tuple1->addIndexedItem (
"pidchiDedx" , m_ngch, m_pidchiDedx);
350 status = m_tuple1->addIndexedItem (
"pidchiTof1" , m_ngch, m_pidchiTof1);
351 status = m_tuple1->addIndexedItem (
"pidchiTof2" , m_ngch, m_pidchiTof2);
354 status = m_tuple1->addItem (
"dedx_GoodHits_ep" , m_dedx_goodhits_ep);
355 status = m_tuple1->addItem (
"dedx_chi_ep" , m_dedx_chiep);
356 status = m_tuple1->addItem (
"dedx_GoodHits_em" , m_dedx_goodhits_em);
357 status = m_tuple1->addItem (
"dedx_chi_em" , m_dedx_chiem);
359 status = m_tuple1->addItem (
"px_cms_ep", m_px_cms_ep);
360 status = m_tuple1->addItem (
"py_cms_ep", m_py_cms_ep);
361 status = m_tuple1->addItem (
"pz_cms_ep", m_pz_cms_ep);
362 status = m_tuple1->addItem (
"e_cms_ep", m_e_cms_ep);
363 status = m_tuple1->addItem (
"p_cms_ep", m_p_cms_ep);
365 status = m_tuple1->addItem (
"cos_ep", m_cos_ep);
366 status = m_tuple1->addItem (
"kal_p_ep", m_kal_p_ep);
367 status = m_tuple1->addItem (
"kal_px_ep", m_kal_px_ep);
368 status = m_tuple1->addItem (
"kal_py_ep", m_kal_py_ep);
369 status = m_tuple1->addItem (
"kal_pz_ep", m_kal_pz_ep);
371 status = m_tuple1->addItem (
"px_cms_em", m_px_cms_em);
372 status = m_tuple1->addItem (
"py_cms_em", m_py_cms_em);
373 status = m_tuple1->addItem (
"pz_cms_em", m_pz_cms_em);
374 status = m_tuple1->addItem (
"e_cms_em", m_e_cms_em);
375 status = m_tuple1->addItem (
"p_cms_em", m_p_cms_em);
377 status = m_tuple1->addItem (
"cos_em", m_cos_em);
378 status = m_tuple1->addItem (
"kal_p_em", m_kal_p_em);
379 status = m_tuple1->addItem (
"kal_px_em", m_kal_px_em);
380 status = m_tuple1->addItem (
"kal_py_em", m_kal_py_em);
381 status = m_tuple1->addItem (
"kal_pz_em", m_kal_pz_em);
383 status = m_tuple1->addItem (
"mass_ee", m_mass_ee);
384 status = m_tuple1->addItem (
"px_ee", m_px_ee);
385 status = m_tuple1->addItem (
"py_ee", m_py_ee);
386 status = m_tuple1->addItem (
"pz_ee", m_pz_ee);
387 status = m_tuple1->addItem (
"e_ee", m_e_ee);
388 status = m_tuple1->addItem (
"p_ee", m_p_ee);
390 status = m_tuple1->addItem (
"nep", m_nep );
391 status = m_tuple1->addItem (
"nem", m_nem );
396 log << MSG::ERROR <<
" Cannot book N-tuple:" << long(m_tuple1) << endmsg;
397 return StatusCode::FAILURE;
405 log << MSG::INFO <<
"successfully return from initialize()" <<endmsg;
406 return StatusCode::SUCCESS;
412 const double beamEnergy = m_ecms/2.;
413 const HepLorentzVector
p_cms(m_ecms*
sin(m_beamangle*0.5),0.0,0.0,m_ecms);
414 const Hep3Vector
u_cms = -
p_cms.boostVector();
416 MsgStream log(
msgSvc(), name());
417 log << MSG::INFO <<
"in execute()" << endreq;
419 setFilterPassed(
false);
421 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
424 log << MSG::FATAL <<
"Could not find Event Header" << endreq;
425 return StatusCode::SUCCESS;
428 m_run = eventHeader->runNumber();
429 m_rec = eventHeader->eventNumber();
436 log << MSG::FATAL <<
"Could not find EvtRecEvent" << endreq;
437 return StatusCode::SUCCESS;
439 log << MSG::INFO <<
"ncharg, nneu, tottks = "
440 << evtRecEvent->totalCharged() <<
" , "
441 << evtRecEvent->totalNeutral() <<
" , "
442 << evtRecEvent->totalTracks() <<endreq;
444 m_ncharg = evtRecEvent->totalCharged();
445 m_nneu = evtRecEvent->totalNeutral();
449 HepSymMatrix Evx(3, 0);
452 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
459 Evx[0][0]=vv[0]*vv[0];
460 Evx[0][1]=vv[0]*vv[1];
461 Evx[1][1]=vv[1]*vv[1];
462 Evx[1][2]=vv[1]*vv[2];
463 Evx[2][2]=vv[2]*vv[2];
470 log << MSG::FATAL <<
"Could not find EvtRecTrackCol" << endreq;
471 return StatusCode::SUCCESS;
478 for(
int i = 0; i < evtRecEvent->totalCharged(); i++){
480 if(!(*itTrk)->isMdcTrackValid())
continue;
481 if(!(*itTrk)->isMdcKalTrackValid())
continue;
484 double pch=mdcTrk->
p();
485 double x0=mdcTrk->
x();
486 double y0=mdcTrk->
y();
487 double z0=mdcTrk->
z();
488 double phi0=mdcTrk->
helix(1);
492 double Rxy=(x0-xv)*
cos(phi0)+(y0-yv)*
sin(phi0);
499 if(fabs(z0) >= m_vz0cut)
continue;
501 if(fabs(Rxy) >= m_vr0cut)
continue;
504 nCharge += mdcTrk->
charge();
508 int nGood = iGood.size();
510 log << MSG::DEBUG <<
"ngood, totcharge = " << nGood <<
" , " << nCharge << endreq;
512 if((nGood != 2)||(nCharge!=0)){
513 return StatusCode::SUCCESS;
530 Vint ipip, ipim, iep,iem,imup,imum;
541 for(
int i = 0; i < m_ngch; i++) {
566 HepLorentzVector
ptrk;
567 ptrk.setPx(mdcTrk->
px()) ;
568 ptrk.setPy(mdcTrk->
py()) ;
569 ptrk.setPz(mdcTrk->
pz()) ;
573 m_pidprob[i]=pid->
prob(0);
574 m_pidchiDedx[i]=pid->
chiDedx(0);
575 m_pidchiTof1[i]=pid->
chiTof1(0);
576 m_pidchiTof2[i]=pid->
chiTof2(0);
577 if(mdcTrk->
charge() > 0) {
578 iep.push_back(iGood[i]);
580 if (mdcTrk->
charge() < 0) {
581 iem.push_back(iGood[i]);
597 for(
int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
598 if(i>=evtRecTrkCol->size())
break;
600 if(!(*itTrk)->isEmcShowerValid())
continue;
602 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
607 int module=emcTrk->
module();
608 double x = emcTrk->
x();
609 double y = emcTrk->
y();
610 double z = emcTrk->
z();
611 double dx = emcTrk->
dx();
612 double dy = emcTrk->
dy();
613 double dth = emcTrk->
dtheta();
614 double dph = emcTrk->
dphi();
615 double dz = emcTrk->
dz();
617 double dE = emcTrk->
dE();
618 double eSeed = emcTrk->
eSeed();
619 double e3x3 = emcTrk->
e3x3();
620 double e5x5 = emcTrk->
e5x5();
623 double getTime = emcTrk->
time();
624 double getEAll = emcTrk->
getEAll();
630 m_nemchits[iphoton]=
n;
631 m_npart[iphoton]=npart;
632 m_module[iphoton]=
module;
633 m_theta[iphoton]=EmcPos.theta();
634 m_phi[iphoton]=EmcPos.phi();
641 m_dtheta[iphoton]=dth;
645 m_eSeed[iphoton]=eSeed;
646 m_nSeed[iphoton]=nseed;
647 m_e3x3[iphoton]=e3x3;
648 m_e5x5[iphoton]=e5x5;
649 m_secondMoment[iphoton]=secondMoment;
650 m_latMoment[iphoton]=latMoment;
651 m_getTime[iphoton]=getTime;
652 m_getEAll[iphoton]=getEAll;
653 m_a20Moment[iphoton]=a20Moment;
654 m_a42Moment[iphoton]=a42Moment;
661 for(
int j = 0; j < nGood; j++) {
663 if (!(*jtTrk)->isMdcTrackValid())
continue;
665 double jtcharge = jtmdcTrk->
charge();
666 if(!(*jtTrk)->isExtTrackValid())
continue;
671 double angd = extpos.angle(emcpos);
672 double thed = extpos.theta() - emcpos.theta();
673 double phid = extpos.deltaPhi(emcpos);
674 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
675 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
677 if(fabs(thed) < fabs(dthe)) dthe = thed;
678 if(fabs(phid) < fabs(dphi)) dphi = phid;
679 if(angd < dang) dang = angd;
688 dthe = dthe * 180 / (CLHEP::pi);
689 dphi = dphi * 180 / (CLHEP::pi);
690 dang = dang * 180 / (CLHEP::pi);
691 double eraw = emcTrk->
energy();
693 double phi = emcTrk->
phi();
694 double the = emcTrk->
theta();
696 m_delphi[iphoton]=dphi;
697 m_delthe[iphoton]=dthe;
698 m_delang[iphoton]=dang;
701 if (fc_theta < 0.80 ) {
702 if ( eraw <= m_energyThreshold/2)
continue;
703 }
else if (fc_theta > 0.86 && fc_theta < 0.92 ) {
704 if (eraw <= m_energyThreshold )
continue;
708 if(getTime>m_gammathCut||getTime<m_gammatlCut)
continue;
711 if(dang< m_gammaTrkCut)
continue;
715 if(iphoton>=40)
return StatusCode::SUCCESS;
718 if(iphoton>0) counter[4]++;
719 int nGam = iGam.size();
730 for(
int i = 0; i < m_nGam; i++) {
732 if(!(*itTrk)->isEmcShowerValid())
continue;
734 double eraw = emcTrk->
energy();
735 double phi = emcTrk->
phi();
736 double the = emcTrk->
theta();
737 HepLorentzVector
ptrk;
738 ex_gam+=eraw*
sin(the)*
cos(phi);
739 ey_gam+=eraw*
sin(the)*
sin(phi);
740 ez_gam+=eraw*
cos(the);
741 et_gam+=eraw*
sin(the);
768 for(
int i = 0; i < m_ngch; i++ ){
772 if(!(*itTrk)->isMdcTrackValid())
continue;
773 if(!(*itTrk)->isMdcKalTrackValid())
continue;
778 m_charge[i] = mdcTrk->
charge();
779 m_vx0[i] = mdcTrk->
x();
780 m_vy0[i] = mdcTrk->
y();
781 m_vz0[i] = mdcTrk->
z();
783 m_px[i] = mdcTrk->
px();
784 m_py[i] = mdcTrk->
py();
785 m_pz[i] = mdcTrk->
pz();
786 m_p[i] = mdcTrk->
p();
799 m_kal_vx0[i]= mdcKalTrk->
dr()*
cos(mdcKalTrk->
fi0());
801 m_kal_vy0[i] = mdcKalTrk->
dr()*
sin(mdcKalTrk->
fi0());
803 m_kal_vz0[i]= mdcKalTrk->
dz();
806 m_kal_px[i] = mdcKalTrk->
px();
807 m_kal_py[i] = mdcKalTrk->
py();
808 m_kal_pz[i] = mdcKalTrk->
pz();
810 t_pxy2 = m_kal_px[i]*m_kal_px[i] + m_kal_py[i]*m_kal_py[i];
811 m_kal_p[i] = sqrt(t_pxy2 + m_kal_pz[i]*m_kal_pz[i]);
812 double ptrk = m_kal_p[i];
816 pt_had+=sqrt(t_pxy2);
818 e_had+=sqrt(m_kal_p[i]*m_kal_p[i]+mdcKalTrk->
mass()*mdcKalTrk->
mass());
821 if((*itTrk)->isMdcDedxValid()) {
824 m_probPH[i]= dedxTrk->
probPH();
825 m_normPH[i]= dedxTrk->
normPH();
827 m_chie[i] = dedxTrk->
chiE();
828 m_chimu[i] = dedxTrk->
chiMu();
829 m_chipi[i] = dedxTrk->
chiPi();
830 m_chik[i] = dedxTrk->
chiK();
831 m_chip[i] = dedxTrk->
chiP();
836 if((*itTrk)->isEmcShowerValid()) {
838 m_e_emc[i] = emcTrk->
energy();
839 m_phi_emc[i] = emcTrk->
phi();
840 m_theta_emc[i] = emcTrk->
theta();
846 if((*itTrk)->isMucTrackValid()){
848 m_nhit_muc[i] = mucTrk->
numHits() ;
852 if((*itTrk)->isTofTrackValid()) {
854 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
855 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
857 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
861 status->
setStatus((*iter_tof)->status());
863 if( (status->
is_cluster()) ) m_t_etof[i] = (*iter_tof)->tof();
864 if( !(status->
is_counter()) ){
if(status)
delete status;
continue; }
865 if( status->
layer()!=0 ) {
if(status)
delete status;
continue;}
867 double path=(*iter_tof)->path();
868 double tof = (*iter_tof)->tof();
869 double ph = (*iter_tof)->ph();
870 double rhit = (*iter_tof)->zrhit();
871 double qual = 0.0 + (*iter_tof)->quality();
872 double cntr = 0.0 + (*iter_tof)->tofID();
874 for(
int j = 0; j < 5; j++) {
876 double beta = gb/sqrt(1+gb*gb);
877 texp[j] = path /beta/
velc;
879 m_qual_etof[i] = qual;
880 m_tof_etof[i] = tof ;
883 if( (status->
is_cluster()) ) m_t_btof[i] = (*iter_tof)->tof();
884 if( !(status->
is_counter()) ){
if(status)
delete status;
continue;}
885 if(status->
layer()==1){
886 double path=(*iter_tof)->path();
887 double tof = (*iter_tof)->tof();
888 double ph = (*iter_tof)->ph();
889 double rhit = (*iter_tof)->zrhit();
890 double qual = 0.0 + (*iter_tof)->quality();
891 double cntr = 0.0 + (*iter_tof)->tofID();
893 for(
int j = 0; j < 5; j++) {
895 double beta = gb/sqrt(1+gb*gb);
896 texp[j] = path /beta/
velc;
898 m_qual_btof1[i] = qual;
899 m_tof_btof1[i] = tof ;
902 if(status->
layer()==2){
903 double path=(*iter_tof)->path();
904 double tof = (*iter_tof)->tof();
905 double ph = (*iter_tof)->ph();
906 double rhit = (*iter_tof)->zrhit();
907 double qual = 0.0 + (*iter_tof)->quality();
908 double cntr = 0.0 + (*iter_tof)->tofID();
910 for(
int j = 0; j < 5; j++) {
912 double beta = gb/sqrt(1+gb*gb);
913 texp[j] = path /beta/
velc;
915 m_qual_btof2[i] = qual;
916 m_tof_btof2[i] = tof ;
919 if(status)
delete status;
930 if(m_ngch != 2 || nCharge != 0 )
return StatusCode::SUCCESS;
936 HepLorentzVector p41e,p42e,p4le;
937 Hep3Vector p31e,p32e,p3le;
938 HepLorentzVector p41m,p42m,p4lm;
939 Hep3Vector p31m,p32m,p3lm;
940 HepLorentzVector p41h,p42h,p4lh;
941 Hep3Vector p31h,p32h,p3lh;
958 double ex1, ey1, ez1, epx1,epy1, epz1,epp1,pidchidedx1,pidchitof11,pidchitof21,eoeb1,exoeb1,eyoeb1,ezoeb1,etoeb1,kalpp,cmsp;
959 double ex2, ey2, ez2, epx2,epy2, epz2,epp2,pidchidedx2,pidchitof12,pidchitof22,eoeb2,exoeb2,eyoeb2,ezoeb2,etoeb2,kalpm,cmsm;
962 for(
int i = 0; i < m_ngch; i++ ){
964 itTrk1= evtRecTrkCol->begin() + iGood[i];
965 mdcKalTrk1 = (*itTrk1)->mdcKalTrack();
967 if(!(*itTrk1)->isMdcDedxValid())
continue;
971 m_dedx_chiep=dedxTrk1->
chiE();
981 m_px_cms_ep=p41e.px();
982 m_py_cms_ep=p41e.py();
983 m_pz_cms_ep=p41e.pz();
985 m_p_cms_ep=sqrt(p41e.px()*p41e.px()+p41e.py()*p41e.py()+p41e.pz()*p41e.pz());
988 m_kal_p_ep=m_kal_p[i];
1004 pidchidedx1=m_pidchiDedx[i];
1005 pidchitof11=m_pidchiTof1[i];
1006 pidchitof21=m_pidchiTof2[i];
1008 eoeb1=m_e_emc[i]/beamEnergy;
1010 if(p41e.rho()>0) eope1=m_e_emc[i]/p41e.rho();
1013 exoeb1=m_e_emc[i]*
sin(m_theta_emc[i])*
cos(m_phi_emc[i])/beamEnergy;
1014 eyoeb1=m_e_emc[i]*
sin(m_theta_emc[i])*
sin(m_phi_emc[i])/beamEnergy;
1015 ezoeb1=m_e_emc[i]*
cos(m_theta_emc[i])/beamEnergy;
1016 etoeb1=m_e_emc[i]*
sin(m_theta_emc[i])/beamEnergy;
1018 if(m_nhit_muc[i]>=2&&m_nlay_muc[i]>=2) mucinfo1=1;
1025 itTrk2= evtRecTrkCol->begin() + iGood[i];
1026 mdcKalTrk2 = (*itTrk2)->mdcKalTrack();
1029 if(!(*itTrk2)->isMdcDedxValid())
continue;
1033 m_dedx_chiem=dedxTrk2->
chiE();
1042 m_px_cms_em=p42e.px();
1043 m_py_cms_em=p42e.py();
1044 m_pz_cms_em=p42e.pz();
1045 m_e_cms_em=p42e.e();
1046 m_p_cms_em=sqrt(p42e.px()*p42e.px()+p42e.py()*p42e.py()+p42e.pz()*p42e.pz());
1048 m_kal_p_em=m_kal_p[i];
1064 pidchidedx2=m_pidchiDedx[i];
1065 pidchitof12=m_pidchiTof1[i];
1066 pidchitof22=m_pidchiTof2[i];
1068 eoeb2=m_e_emc[i]/beamEnergy;
1071 if(p42e.rho()>0) eope2=m_e_emc[i]/p42e.rho();
1073 exoeb2=m_e_emc[i]*
sin(m_theta_emc[i])*
cos(m_phi_emc[i])/beamEnergy;
1074 eyoeb2=m_e_emc[i]*
sin(m_theta_emc[i])*
sin(m_phi_emc[i])/beamEnergy;
1075 ezoeb2=m_e_emc[i]*
cos(m_theta_emc[i])/beamEnergy;
1076 etoeb2=m_e_emc[i]*
sin(m_theta_emc[i])/beamEnergy;
1078 if(m_nhit_muc[i]>=2&&m_nlay_muc[i]>=2) mucinfo2=1;
1085 p4le=( e01 > e02 ) ?p41e:p42e;
1086 p4lm=( e01 > e02 ) ?p41m:p42m;
1087 p3le=( e01 > e02 ) ?p31e:p32e;
1088 p3lm=( e01 > e02 ) ?p31m:p32m;
1090 double acolle=180.-p31e.angle(p32e)* 180.0 / CLHEP::pi;
1091 double acople=180.- (p31e.perpPart()).angle(p32e.perpPart ())* 180.0 / CLHEP::pi;
1092 double poeb1e=p41e.rho()/beamEnergy;
1093 double poeb2e=p42e.rho()/beamEnergy;
1095 int ilarge=( e01 > e02 ) ?iip:iim;
1097 double eoebl=m_e_emc[ilarge]/beamEnergy;
1098 if(p4le.rho()>0)eopl=m_e_emc[ilarge]/p4le.rho();
1099 double exoebl= m_e_emc[ilarge]*
sin(m_theta_emc[ilarge])*
cos(m_phi_emc[ilarge])/beamEnergy;
1100 double eyoebl= m_e_emc[ilarge]*
sin(m_theta_emc[ilarge])*
sin(m_phi_emc[ilarge])/beamEnergy;
1101 double ezoebl=m_e_emc[ilarge]*
cos(m_theta_emc[ilarge])/beamEnergy;
1102 double etoebl=m_e_emc[ilarge]*
sin(m_theta_emc[ilarge])/beamEnergy;
1103 int mucinfol=(m_nhit_muc[ilarge]>=2&&m_nlay_muc[ilarge]>=2) ? 1 : 0;
1105 int pidel=( e01 > e02 ) ? m_nep : m_nem;
1107 if(m_t_btof[iip]*m_t_btof[iim]!=0) deltatof=fabs(m_t_btof[iip]-m_t_btof[iim]);
1108 if(m_t_etof[iip]*m_t_etof[iim]!=0) deltatof=fabs(m_t_etof[iip]-m_t_etof[iim]);
1112 if (m_use_acolle&&acolle>m_acoll_e_cut)
return StatusCode::SUCCESS;
1117 if (m_use_acople&&acople>m_acopl_e_cut)
return StatusCode::SUCCESS;
1121 if (m_use_eoeb&&sqrt((eoeb1-1)*(eoeb1-1)+(eoeb2-1)*(eoeb2-1))>m_tetotal_e_cut)
return StatusCode::SUCCESS;
1127 if (m_use_deltatof&&m_useTOF&&(deltatof>m_dtof_e_cut))
return StatusCode::SUCCESS;
1132 if (m_use_poeb&&poeb1e<m_tpoeb_e_cut&&poeb2e<m_tpoeb_e_cut&&(sqrt((poeb1e-1)*(poeb1e-1)+(poeb2e-1)*(poeb2e-1))>m_tptotal_e_cut))
return StatusCode::SUCCESS;
1137 if (m_use_mucinfo&&m_useMUC&&(mucinfo1!=0&&mucinfo2!=0))
return StatusCode::SUCCESS;
1142 if (m_use_ne&&m_usePID&&(m_nep!=1||m_nem!=1))
return StatusCode::SUCCESS;
1155 m_cos_ep=p41e.cosTheta ();
1156 m_cos_em=p42e.cosTheta ();
1157 m_mass_ee=(p41e+p42e).m();
1158 m_px_ee=(p41e+p42e).px();
1159 m_py_ee=(p41e+p42e).py();
1160 m_pz_ee=(p41e+p42e).pz();
1161 m_e_ee=(p41e+p42e).e();
1162 m_p_ee=(p41e+p42e).rho();
1164 m_deltatof=deltatof;
1171 m_mucinfo1=mucinfo1;
1172 m_mucinfo2=mucinfo2;
1185 (*itTrk1)->tagElectron();
1186 (*itTrk2)->tagElectron();
1192 (*itTrk1)->setQuality(0);
1193 (*itTrk2)->setQuality(0);
1198 setFilterPassed(
true);
1200 m_ee_mass->Fill((p41e+p42e).m());
1201 m_ee_acoll->Fill(acolle);
1202 m_ee_eop_ep->Fill(eope1);
1203 m_ee_eop_em->Fill(eope2);
1204 m_ee_costheta_ep->Fill(p41e.cosTheta ());
1205 m_ee_costheta_em->Fill(p42e.cosTheta ());
1207 m_ee_phi_ep->Fill(p41e.phi ());
1208 m_ee_phi_em->Fill(p42e.phi ());
1210 m_ee_nneu->Fill(m_nGam);
1213 m_ee_eemc_ep->Fill(e01);
1214 m_ee_eemc_em->Fill(e02);
1216 m_ee_x_ep->Fill(ex1);
1217 m_ee_y_ep->Fill(ey1);
1218 m_ee_z_ep->Fill(ez1);
1220 m_ee_x_em->Fill(ex2);
1221 m_ee_y_em->Fill(ey2);
1222 m_ee_z_em->Fill(ez2);
1225 m_ee_px_ep->Fill(epx1);
1226 m_ee_py_ep->Fill(epy1);
1227 m_ee_pz_ep->Fill(epz1);
1228 m_ee_p_ep->Fill(epp1);
1230 m_ee_px_em->Fill(epx2);
1231 m_ee_py_em->Fill(epy2);
1232 m_ee_pz_em->Fill(epz2);
1233 m_ee_p_em->Fill(epp2);
1235 m_ee_deltatof->Fill(deltatof);
1237 m_ee_pidchidedx_ep->Fill(pidchidedx1);
1238 m_ee_pidchidedx_em->Fill(pidchidedx2);
1239 m_ee_pidchitof1_ep->Fill(pidchitof11);
1240 m_ee_pidchitof1_em->Fill(pidchitof12);
1241 m_ee_pidchitof2_ep->Fill(pidchitof21);
1242 m_ee_pidchitof2_em->Fill(pidchitof22);
1244 if(m_writentuple) m_tuple1 -> write();
1246 return StatusCode::SUCCESS;