BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
DQASelBhabha.cxx
Go to the documentation of this file.
1#include <vector>
2// #include <iostream.h>
3
4#include "GaudiKernel/MsgStream.h"
5#include "GaudiKernel/AlgFactory.h"
6#include "GaudiKernel/ISvcLocator.h"
7#include "GaudiKernel/SmartDataPtr.h"
8#include "GaudiKernel/IDataProviderSvc.h"
9#include "GaudiKernel/PropertyMgr.h"
11#include "GaudiKernel/Bootstrap.h"
12#include "GaudiKernel/ISvcLocator.h"
13
15#include "EventModel/Event.h"
16
21
22#include "TMath.h"
23#include "GaudiKernel/INTupleSvc.h"
24#include "GaudiKernel/NTuple.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/IHistogramSvc.h"
27#include "CLHEP/Vector/ThreeVector.h"
28#include "CLHEP/Vector/LorentzVector.h"
29#include "CLHEP/Vector/TwoVector.h"
30
31using CLHEP::Hep3Vector;
32using CLHEP::Hep2Vector;
33using CLHEP::HepLorentzVector;
34#include "CLHEP/Geometry/Point3D.h"
35
37#include "VertexFit/VertexFit.h"
40
41
43
44#ifndef ENABLE_BACKWARDS_COMPATIBILITY
46#endif
47
51
52
53
54#include "Identifier/TofID.h"
55
57
58
59using CLHEP::HepLorentzVector;
60const double mpsi2s=3.68609;
61const double mjpsi =3.09691;
62const double mpi = 0.13957;
63const double mk = 0.493677;
64const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
65const double velc = 299.792458; // tof path unit in mm
66typedef std::vector<int> Vint;
67typedef std::vector<HepLorentzVector> Vp4;
68//declare one counter
69static int counter[13]={0,0,0,0,0,0,0,0,0,0,0,0,0};
70static int nbhabha=0;
71
72/////////////////////////////////////////////////////////////////////////////
73
74DQASelBhabha::DQASelBhabha(const std::string& name, ISvcLocator* pSvcLocator) :
75 Algorithm(name, pSvcLocator) {
76
77 //Declare the properties
78 declareProperty("writentuple",m_writentuple = true );
79 declareProperty("useVertexDB", m_useVertexDB = false );
80 declareProperty("ecms",m_ecms = 3.097);
81 declareProperty("beamangle",m_beamangle = 0.022);
82 declareProperty("Vr0cut", m_vr0cut=1.0);
83 declareProperty("Vz0cut", m_vz0cut=10.0);
84 declareProperty("Coscut", m_coscut=0.93);
85
86 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
87 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
88 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
89 declareProperty("GammaTrkCut", m_gammaTrkCut=20.0);
90 declareProperty("GammaTLCut", m_gammatlCut=0);
91 declareProperty("GammaTHCut", m_gammathCut=14);
92
93
94 declareProperty ("acoll_e_cut", m_acoll_e_cut=6.);
95 declareProperty ("acopl_e_cut", m_acopl_e_cut=6.);
96 declareProperty ("poeb_e_cut", m_poeb_e_cut=0.5);
97 declareProperty ("dtof_e_cut", m_dtof_e_cut=4.);
98 declareProperty ("eoeb_e_cut", m_eoeb_e_cut=0.4);
99 declareProperty ("etotal_e_cut", m_etotal_e_cut=0.8);
100 declareProperty ("tpoeb_e_cut", m_tpoeb_e_cut=0.95);
101 declareProperty ("tptotal_e_cut", m_tptotal_e_cut=0.16);
102 declareProperty ("tetotal_e_cut", m_tetotal_e_cut=0.65);
103
104
105 //normally, MDC+EMC, otherwise EMC only
106 declareProperty ("m_useEMConly", m_useEMConly= false);
107 declareProperty ("m_usePID", m_usePID= false);// sub-system is under study
108 declareProperty ("m_useMDC", m_useMDC= true);
109 declareProperty ("m_useDEDX", m_useDEDX= false);// not used
110 declareProperty ("m_useTOF", m_useTOF= false);//sub-system is under study
111 declareProperty ("m_useEMC", m_useEMC= true);
112 declareProperty ("m_useMUC", m_useMUC= false);// efficiency
113
114 declareProperty ("m_use_acolle", m_use_acolle = false);
115 declareProperty ("m_use_acople",m_use_acople = false);
116 declareProperty ("m_use_eoeb",m_use_eoeb = false);
117 declareProperty ("m_use_deltatof", m_use_deltatof = false);
118 declareProperty ("m_use_poeb", m_use_poeb= false);
119 declareProperty ("m_use_mucinfo",m_use_mucinfo= false);
120 declareProperty ("m_use_ne",m_use_ne = false);
121}
122
123
124// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
126
127 MsgStream log(msgSvc(), name());
128
129 log << MSG::INFO << "in initialize()" << endmsg;
130 StatusCode status;
131 status = service("THistSvc", m_thistsvc);
132 if(status.isFailure() ){
133 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
134 return status;
135 }
136
137
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);
142
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);
155
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);
160
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);
163
164
165
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);
182
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);
201
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);
214
215
216
217 NTuplePtr nt1(ntupleSvc(), "DQAFILE/Bhabha");
218 if ( nt1 ) m_tuple1 = nt1;
219 else {
220 m_tuple1 = ntupleSvc()->book ("DQAFILE/Bhabha", CLID_ColumnWiseTuple, "N-Tuple");
221 if ( m_tuple1 ) {
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);
228
229
230 status = m_tuple1->addItem ("bhabhatag", m_bhabhatag);
231
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);
245
246
247
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);
279
280
281
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);
287
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) ;
292
293
294
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);
298
299
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) ;
304
305
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) ;
315
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) ;
319
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 );
331
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 );
339
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);
352
353
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);
358
359 status = m_tuple1->addItem ("px_cms_ep", m_px_cms_ep); //momentum of electron+
360 status = m_tuple1->addItem ("py_cms_ep", m_py_cms_ep); //momentum of electron+
361 status = m_tuple1->addItem ("pz_cms_ep", m_pz_cms_ep); //momentum of electron+
362 status = m_tuple1->addItem ("e_cms_ep", m_e_cms_ep); //momentum of electron+
363 status = m_tuple1->addItem ("p_cms_ep", m_p_cms_ep); //momentum of electron+
364
365 status = m_tuple1->addItem ("cos_ep", m_cos_ep); //momentum of electron+
366 status = m_tuple1->addItem ("kal_p_ep", m_kal_p_ep); //momentum of electron+
367 status = m_tuple1->addItem ("kal_px_ep", m_kal_px_ep); //momentum of electron+
368 status = m_tuple1->addItem ("kal_py_ep", m_kal_py_ep); //momentum of electron+
369 status = m_tuple1->addItem ("kal_pz_ep", m_kal_pz_ep); //momentum of electron+
370
371 status = m_tuple1->addItem ("px_cms_em", m_px_cms_em); //momentum of electron-
372 status = m_tuple1->addItem ("py_cms_em", m_py_cms_em); //momentum of electron-
373 status = m_tuple1->addItem ("pz_cms_em", m_pz_cms_em); //momentum of electron-
374 status = m_tuple1->addItem ("e_cms_em", m_e_cms_em); //momentum of electron-
375 status = m_tuple1->addItem ("p_cms_em", m_p_cms_em); //momentum of electron-
376
377 status = m_tuple1->addItem ("cos_em", m_cos_em); //momentum of electron-
378 status = m_tuple1->addItem ("kal_p_em", m_kal_p_em); //momentum of electron-
379 status = m_tuple1->addItem ("kal_px_em", m_kal_px_em); //momentum of electron-
380 status = m_tuple1->addItem ("kal_py_em", m_kal_py_em); //momentum of electron-
381 status = m_tuple1->addItem ("kal_pz_em", m_kal_pz_em); //momentum of electron-
382
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); //
389
390 status = m_tuple1->addItem ( "nep", m_nep );
391 status = m_tuple1->addItem ( "nem", m_nem );
392
393
394 }
395 else {
396 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
397 return StatusCode::FAILURE;
398 }
399 }
400
401 //
402 //--------end of book--------
403 //
404
405 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
406 return StatusCode::SUCCESS;
407
408}
409
410// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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();
415 counter[0]++;
416 MsgStream log(msgSvc(), name());
417 log << MSG::INFO << "in execute()" << endreq;
418
419 setFilterPassed(false);
420
421 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
422 if (!eventHeader)
423 {
424 log << MSG::FATAL << "Could not find Event Header" << endreq;
425 return StatusCode::SUCCESS;
426 }
427
428 m_run = eventHeader->runNumber();
429 m_rec = eventHeader->eventNumber();
430
431
432
433 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
434 if (!evtRecEvent)
435 {
436 log << MSG::FATAL << "Could not find EvtRecEvent" << endreq;
437 return StatusCode::SUCCESS;
438 }
439 log << MSG::INFO <<"ncharg, nneu, tottks = "
440 << evtRecEvent->totalCharged() << " , "
441 << evtRecEvent->totalNeutral() << " , "
442 << evtRecEvent->totalTracks() <<endreq;
443
444 m_ncharg = evtRecEvent->totalCharged();
445 m_nneu = evtRecEvent->totalNeutral();
446
447
448 HepPoint3D vx(0., 0., 0.);
449 HepSymMatrix Evx(3, 0);
450 if (m_useVertexDB) {
451 IVertexDbSvc* vtxsvc;
452 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
453 if(vtxsvc->isVertexValid()){
454 double* dbv = vtxsvc->PrimaryVertex();
455 double* vv = vtxsvc->SigmaPrimaryVertex();
456 vx.setX(dbv[0]);
457 vx.setY(dbv[1]);
458 vx.setZ(dbv[2]);
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];
464 }
465 }
466
467 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
468 if (!evtRecTrkCol)
469 {
470 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endreq;
471 return StatusCode::SUCCESS;
472 }
473 Vint iGood;
474 iGood.clear();
475
476 int nCharge = 0;
477
478 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
479 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
480 if(!(*itTrk)->isMdcTrackValid()) continue;
481 if(!(*itTrk)->isMdcKalTrackValid()) continue;
482
483 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
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);
489 double xv=vx.x();
490 double yv=vx.y();
491 double zv=vx.z();
492 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
493 double m_vx0 = x0;
494 double m_vy0 = y0;
495 double m_vz0 = z0;
496 double m_vr0 = Rxy;
497
498 e_z0->Fill(z0);
499 if(fabs(z0) >= m_vz0cut) continue;
500 e_r0->Fill(Rxy);
501 if(fabs(Rxy) >= m_vr0cut) continue;
502
503 iGood.push_back(i);
504 nCharge += mdcTrk->charge();
505
506 }
507
508 int nGood = iGood.size();
509 m_ngch=nGood;
510 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
511
512 if((nGood != 2)||(nCharge!=0)){
513 return StatusCode::SUCCESS;
514 }
515
516
517
518 counter[1]++;
519
520 //
521 // Finish Good Charged Track Selection
522 //
523
524
525
526
527 //
528 // Particle ID
529 //
530 Vint ipip, ipim, iep,iem,imup,imum;
531 ipip.clear();
532 ipim.clear();
533 iep.clear();
534 iem.clear();
535 imup.clear();
536 imum.clear();
537
538
540
541 for(int i = 0; i < m_ngch; i++) {
542 m_pidcode[i]=-99;
543 m_pidprob[i]=-99;
544 m_pidchiDedx[i]=-99;
545 m_pidchiTof1[i]=-99;
546 m_pidchiTof2[i]=-99;
547
548 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
549
550 pid->init();
551 pid->setMethod(pid->methodProbability());
552 pid->setChiMinCut(4);
553 pid->setRecTrack(*itTrk);
554 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());//|pid->useEmc()|pid->useMuc()); // use PID sub-system
555 pid->identify(pid->onlyElectron()|pid->onlyMuon()|pid->onlyPion()); // seperater Pion/Kaon/Proton
556 pid->calculate();
557 if(!(pid->IsPidInfoValid())) continue;
558 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
559
560 double prob_pi = pid->probPion();
561 double prob_K = pid->probKaon();
562 double prob_p = pid->probProton();
563 double prob_e = pid->probElectron();
564 double prob_mu = pid->probMuon();
565
566 HepLorentzVector ptrk;
567 ptrk.setPx(mdcTrk->px()) ;
568 ptrk.setPy(mdcTrk->py()) ;
569 ptrk.setPz(mdcTrk->pz()) ;
570 double p3 = ptrk.mag() ;
571
572 m_pidcode[i]=0;
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]);
579 }
580 if (mdcTrk->charge() < 0) {
581 iem.push_back(iGood[i]);
582 }
583 }
584
585 m_nep = iep.size() ;
586 m_nem = iem.size() ;
587
588 counter[2]++;
589
590
591 //
592 // Good neutral track selection
593 //
594 Vint iGam;
595 iGam.clear();
596 int iphoton=0;
597 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
598 if(i>=evtRecTrkCol->size())break;
599 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
600 if(!(*itTrk)->isEmcShowerValid()) continue;
601 RecEmcShower *emcTrk = (*itTrk)->emcShower();
602 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
603
604 RecEmcID showerId = emcTrk->getShowerId();
605 unsigned int npart = EmcID::barrel_ec(showerId);
606 int n = emcTrk->numHits();
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();
616 double energy = emcTrk->energy();
617 double dE = emcTrk->dE();
618 double eSeed = emcTrk->eSeed();
619 double e3x3 = emcTrk->e3x3();
620 double e5x5 = emcTrk->e5x5();
621 double secondMoment = emcTrk->secondMoment();
622 double latMoment = emcTrk->latMoment();
623 double getTime = emcTrk->time();
624 double getEAll = emcTrk->getEAll();
625 double a20Moment = emcTrk->a20Moment();
626 double a42Moment = emcTrk->a42Moment();
627 double nseed=0;
628
629 HepPoint3D EmcPos(x,y,z);
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();
635 m_x[iphoton]=x;
636 m_y[iphoton]=y;
637 m_z[iphoton]=z;
638 m_dx[iphoton]=dx;
639 m_dy[iphoton]=dy;
640 m_dz[iphoton]=dz;
641 m_dtheta[iphoton]=dth;
642 m_dphi[iphoton]=dph;
643 m_energy[iphoton]=energy;
644 m_dE[iphoton]=dE;
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;
655
656 double dthe = 200.;
657 double dphi = 200.;
658 double dang = 200.;
659
660 // find the nearest charged track
661 for(int j = 0; j < nGood; j++) {
662 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() +iGood[j];
663 if (!(*jtTrk)->isMdcTrackValid()) continue;
664 RecMdcTrack *jtmdcTrk = (*jtTrk)->mdcTrack();
665 double jtcharge = jtmdcTrk->charge();
666 if(!(*jtTrk)->isExtTrackValid()) continue;
667 RecExtTrack *extTrk = (*jtTrk)->extTrack();
668 if(extTrk->emcVolumeNumber() == -1) continue;
669 Hep3Vector extpos = extTrk->emcPosition();
670 // double ctht = extpos.cosTheta(emcpos);
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;
676
677 if(fabs(thed) < fabs(dthe)) dthe = thed;
678 if(fabs(phid) < fabs(dphi)) dphi = phid;
679 if(angd < dang) dang = angd;
680
681 }
682
683
684
685 //
686 // good photon cut will be set here
687 //
688 dthe = dthe * 180 / (CLHEP::pi);
689 dphi = dphi * 180 / (CLHEP::pi);
690 dang = dang * 180 / (CLHEP::pi);
691 double eraw = emcTrk->energy();
692 double theta1 = emcTrk->theta();
693 double phi = emcTrk->phi();
694 double the = emcTrk->theta();
695
696 m_delphi[iphoton]=dphi;
697 m_delthe[iphoton]=dthe;
698 m_delang[iphoton]=dang;
699// if(energy < m_energyThreshold) continue;
700 double fc_theta = fabs(cos(theta1));
701 if (fc_theta < 0.80 ) { // Barrel EMC
702 if ( eraw <= m_energyThreshold/2) continue;
703 } else if (fc_theta > 0.86 && fc_theta < 0.92 ) { // Endcap EMC
704 if (eraw <= m_energyThreshold ) continue;
705 }
706 else continue;
707
708 if(getTime>m_gammathCut||getTime<m_gammatlCut)continue;
709
710 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
711 if(dang< m_gammaTrkCut) continue;
712
713 iphoton++;
714 iGam.push_back(i);
715 if(iphoton>=40)return StatusCode::SUCCESS;
716 }
717
718 if(iphoton>0) counter[4]++;
719 int nGam = iGam.size();
720 m_nGam=nGam;
721
722 counter[3]++;
723
724 double egam_ext=0;
725 double ex_gam=0;
726 double ey_gam=0;
727 double ez_gam=0;
728 double et_gam=0;
729 double e_gam=0;
730 for(int i = 0; i < m_nGam; i++) {
731 EvtRecTrackIterator itTrk = evtRecTrkCol->begin()+ iGam[i];
732 if(!(*itTrk)->isEmcShowerValid()) continue;
733 RecEmcShower* emcTrk = (*itTrk)->emcShower();
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);
742 e_gam+=eraw ;
743 if(eraw>=egam_ext)
744 {
745 egam_ext=eraw;
746 }
747
748 }
749
750
751
752
753
754 double px_had=0;
755 double py_had=0;
756 double pz_had=0;
757 double t_pxy2 = 0;
758 double pt_had=0;
759 double p_had=0;
760 double e_had=0;
761
762 //
763 // check good charged track's infomation
764 //
765 int ii ;
766 m_e_emc[0]=-0.1;
767 m_e_emc[1]=-0.1;
768 for(int i = 0; i < m_ngch; i++ ){
769
770 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
771
772 if(!(*itTrk)->isMdcTrackValid()) continue; // MDC information
773 if(!(*itTrk)->isMdcKalTrackValid()) continue;
774 // if(!(*itTrk)->isEmcShowerValid()) return StatusCode::SUCCESS;///dbg
775 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
776 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
777
778 m_charge[i] = mdcTrk->charge();
779 m_vx0[i] = mdcTrk->x();
780 m_vy0[i] = mdcTrk->y();
781 m_vz0[i] = mdcTrk->z();
782
783 m_px[i] = mdcTrk->px();
784 m_py[i] = mdcTrk->py();
785 m_pz[i] = mdcTrk->pz();
786 m_p[i] = mdcTrk->p();
787
788
790
791
792 /// if(m_pidcode[i]==3)mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
793 /// if(m_pidcode[i]==4)mdcKalTrk->setPidType(RecMdcKalTrack::proton);
794
795 //m_kal_vx0[i] = mdcKalTrk->x();
796 //m_kal_vy0[i] = mdcKalTrk->y();
797 //m_kal_vz0[i] = mdcKalTrk->z();
798
799 m_kal_vx0[i]= mdcKalTrk->dr()*cos(mdcKalTrk->fi0());
800
801 m_kal_vy0[i] = mdcKalTrk->dr()*sin(mdcKalTrk->fi0());
802
803 m_kal_vz0[i]= mdcKalTrk->dz();
804
805
806 m_kal_px[i] = mdcKalTrk->px();
807 m_kal_py[i] = mdcKalTrk->py();
808 m_kal_pz[i] = mdcKalTrk->pz();
809// pxy() and p() are not filled in the reconstruction algorithm
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];
813 px_had+=m_kal_px[i];
814 py_had+=m_kal_py[i];
815 pz_had+=m_kal_pz[i];
816 pt_had+=sqrt(t_pxy2);
817 p_had+=m_kal_p[i];
818 e_had+=sqrt(m_kal_p[i]*m_kal_p[i]+mdcKalTrk->mass()*mdcKalTrk->mass());
819
820
821 if((*itTrk)->isMdcDedxValid()) { // DEDX information
822
823 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
824 m_probPH[i]= dedxTrk->probPH();
825 m_normPH[i]= dedxTrk->normPH();
826
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();
832 m_ghit[i] = dedxTrk->numGoodHits();
833 m_thit[i] = dedxTrk->numTotalHits();
834 }
835
836 if((*itTrk)->isEmcShowerValid()) {
837 RecEmcShower *emcTrk = (*itTrk)->emcShower();
838 m_e_emc[i] = emcTrk->energy();
839 m_phi_emc[i] = emcTrk->phi();
840 m_theta_emc[i] = emcTrk->theta();
841 }
842
843 m_nhit_muc[i] = 0;
844 m_nlay_muc[i] = 0;
845
846 if((*itTrk)->isMucTrackValid()){
847 RecMucTrack* mucTrk = (*itTrk)->mucTrack() ;
848 m_nhit_muc[i] = mucTrk->numHits() ;
849 m_nlay_muc[i] = mucTrk->numLayers() ;
850 }
851
852 if((*itTrk)->isTofTrackValid()) { //TOF information
853
854 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
855 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
856
857 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
858 m_t_etof[i]=0;
859 m_t_btof[i]=0;
860 TofHitStatus *status = new TofHitStatus;
861 status->setStatus((*iter_tof)->status());
862 if(!(status->is_barrel())){//endcap
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;}//layer1
866
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();
873 double texp[5];
874 for(int j = 0; j < 5; j++) {
875 double gb = ptrk/xmass[j];
876 double beta = gb/sqrt(1+gb*gb);
877 texp[j] = path /beta/velc;
878 }
879 m_qual_etof[i] = qual;
880 m_tof_etof[i] = tof ;
881 }
882 else {//barrel
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){ //layer1
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();
892 double texp[5];
893 for(int j = 0; j < 5; j++) {
894 double gb = ptrk/xmass[j];
895 double beta = gb/sqrt(1+gb*gb);
896 texp[j] = path /beta/velc;
897 }
898 m_qual_btof1[i] = qual;
899 m_tof_btof1[i] = tof ;
900 }
901
902 if(status->layer()==2){//layer2
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();
909 double texp[5];
910 for(int j = 0; j < 5; j++) {
911 double gb = ptrk/xmass[j];
912 double beta = gb/sqrt(1+gb*gb);
913 texp[j] = path /beta/velc;
914 }
915 m_qual_btof2[i] = qual;
916 m_tof_btof2[i] = tof ;
917 }
918 }
919 if(status) delete status;
920 }
921 }
922
923 }
924 counter[5]++;
925
926
927
928 m_bhabhatag=0;
929
930 if(m_ngch != 2 || nCharge != 0 ) return StatusCode::SUCCESS;
931 EvtRecTrackIterator itTrk1;
932 EvtRecTrackIterator itTrk2;
933 RecMdcKalTrack *mdcKalTrk1;
934 RecMdcKalTrack *mdcKalTrk2;
935
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;
942 WTrackParameter w1_ini;
943 WTrackParameter w2_ini;
944
945 int iip=-1;
946 int iim=-1;
947 int mucinfo1=0;
948 int mucinfo2=0;
949
950 double e01=0;
951 double e02=0;
952 double eope1=0;
953 double eope2=0;
954 double eopl=0;
955 double deltatof=0;
956
957
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;
960
961
962 for(int i = 0; i < m_ngch; i++ ){
963 if(m_charge[i]>0){
964 itTrk1= evtRecTrkCol->begin() + iGood[i];
965 mdcKalTrk1 = (*itTrk1)->mdcKalTrack();
966
967 if(!(*itTrk1)->isMdcDedxValid())continue;
968 RecMdcDedx* dedxTrk1 = (*itTrk1)->mdcDedx();
969
970 m_dedx_goodhits_ep=dedxTrk1->numGoodHits();
971 m_dedx_chiep=dedxTrk1->chiE();
972
973 iip=i;
974
975 w1_ini=WTrackParameter (xmass[0],mdcKalTrk1->getZHelixE(),mdcKalTrk1->getZErrorE());
976 p41e =w1_ini.p();
977
978 p41e.boost(u_cms);
979 p31e = p41e.vect();
980
981 m_px_cms_ep=p41e.px();
982 m_py_cms_ep=p41e.py();
983 m_pz_cms_ep=p41e.pz();
984 m_e_cms_ep=p41e.e();
985 m_p_cms_ep=sqrt(p41e.px()*p41e.px()+p41e.py()*p41e.py()+p41e.pz()*p41e.pz());
986 cmsp=m_p_cms_ep;
987
988 m_kal_p_ep=m_kal_p[i];
989 kalpp=m_kal_p_ep;
990 e01=m_e_emc[i];
991
992 ex1=m_kal_vx0[i];
993 ey1=m_kal_vy0[i];
994 ez1=m_kal_vz0[i];
995 epx1=m_kal_px[i];
996 epy1=m_kal_py[i];
997 epz1=m_kal_pz[i];
998 epp1=m_kal_p[i];
999
1000 m_kal_px_ep=epx1;
1001 m_kal_py_ep=epy1;
1002 m_kal_pz_ep=epz1;
1003
1004 pidchidedx1=m_pidchiDedx[i];
1005 pidchitof11=m_pidchiTof1[i];
1006 pidchitof21=m_pidchiTof2[i];
1007
1008 eoeb1=m_e_emc[i]/beamEnergy;
1009
1010 if(p41e.rho()>0) eope1=m_e_emc[i]/p41e.rho();
1011
1012
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;
1017
1018 if(m_nhit_muc[i]>=2&&m_nlay_muc[i]>=2) mucinfo1=1;
1019 }
1020
1021
1022
1023
1024 if(m_charge[i]<0){
1025 itTrk2= evtRecTrkCol->begin() + iGood[i];
1026 mdcKalTrk2 = (*itTrk2)->mdcKalTrack();
1027 iim=i;
1028
1029 if(!(*itTrk2)->isMdcDedxValid())continue;
1030 RecMdcDedx* dedxTrk2 = (*itTrk2)->mdcDedx();
1031
1032 m_dedx_goodhits_em=dedxTrk2->numGoodHits();
1033 m_dedx_chiem=dedxTrk2->chiE();
1034
1035
1036 w2_ini=WTrackParameter (xmass[0],mdcKalTrk2 ->getZHelixE(),mdcKalTrk2 ->getZErrorE());
1037 p42e =w2_ini.p();
1038
1039 p42e.boost(u_cms);
1040 p32e = p42e.vect();
1041
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());
1047 cmsm=m_p_cms_em;
1048 m_kal_p_em=m_kal_p[i];
1049 kalpm=m_kal_p_em;
1050 e02=m_e_emc[i];
1051
1052 ex2=m_kal_vx0[i];
1053 ey2=m_kal_vy0[i];
1054 ez2=m_kal_vz0[i];
1055 epx2=m_kal_px[i];
1056 epy2=m_kal_py[i];
1057 epz2=m_kal_pz[i];
1058 epp2=m_kal_p[i];
1059
1060 m_kal_px_em=epx2;
1061 m_kal_py_em=epy2;
1062 m_kal_pz_em=epz2;
1063
1064 pidchidedx2=m_pidchiDedx[i];
1065 pidchitof12=m_pidchiTof1[i];
1066 pidchitof22=m_pidchiTof2[i];
1067
1068 eoeb2=m_e_emc[i]/beamEnergy;
1069
1070
1071 if(p42e.rho()>0) eope2=m_e_emc[i]/p42e.rho();
1072
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;
1077
1078 if(m_nhit_muc[i]>=2&&m_nlay_muc[i]>=2) mucinfo2=1;
1079
1080
1081 }
1082 }
1083
1084
1085 p4le=( e01 > e02 ) ?p41e:p42e;
1086 p4lm=( e01 > e02 ) ?p41m:p42m;
1087 p3le=( e01 > e02 ) ?p31e:p32e;
1088 p3lm=( e01 > e02 ) ?p31m:p32m;
1089
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;
1094
1095 int ilarge=( e01 > e02 ) ?iip:iim;
1096
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;
1104
1105 int pidel=( e01 > e02 ) ? m_nep : m_nem;
1106
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]);
1109
1110
1111 //acolle cut
1112 if (m_use_acolle&&acolle>m_acoll_e_cut) return StatusCode::SUCCESS;
1113 counter[6]++;
1114
1115
1116 //acople cut
1117 if (m_use_acople&&acople>m_acopl_e_cut) return StatusCode::SUCCESS;
1118 counter[7]++;
1119
1120 //eoeb cut
1121 if (m_use_eoeb&&sqrt((eoeb1-1)*(eoeb1-1)+(eoeb2-1)*(eoeb2-1))>m_tetotal_e_cut) return StatusCode::SUCCESS;
1122 counter[8]++;
1123
1124
1125
1126 //deltatof cut
1127 if (m_use_deltatof&&m_useTOF&&(deltatof>m_dtof_e_cut)) return StatusCode::SUCCESS;
1128 counter[9]++;
1129
1130
1131 //poeb cut
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;
1133 counter[10]++;
1134
1135
1136 //mucinfo cut
1137 if (m_use_mucinfo&&m_useMUC&&(mucinfo1!=0&&mucinfo2!=0)) return StatusCode::SUCCESS;
1138 counter[11]++;
1139
1140
1141 //ne cut
1142 if (m_use_ne&&m_usePID&&(m_nep!=1||m_nem!=1)) return StatusCode::SUCCESS;
1143 counter[12]++;
1144
1145
1146
1147
1148
1149 m_acoll=acolle;
1150 m_acopl=acople;
1151 m_poeb1=poeb1e;
1152 m_poeb2=poeb2e;
1153 m_eop1=eope1;
1154 m_eop2=eope2;
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();
1163
1164 m_deltatof=deltatof;
1165
1166 m_eoeb1=eoeb1;
1167 m_eoeb2=eoeb2;
1168
1169 m_etoeb1=etoeb1;
1170 m_etoeb2=etoeb2;
1171 m_mucinfo1=mucinfo1;
1172 m_mucinfo2=mucinfo2;
1173
1174
1175 nbhabha++;
1176
1177
1178 ////////////////////////////////////////////////////////////
1179 // DQA
1180 // set tag and quality
1181
1182
1183 // evtRecTrk->tagElectron(), tagMuon(), tagPion(), tagKaon(), tagProton()
1184
1185 (*itTrk1)->tagElectron();
1186 (*itTrk2)->tagElectron();
1187 // Quality: defined by whether dE/dx or TOF is used to identify particle
1188 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
1189 // 1 - only dE/dx (can be used for TOF calibration)
1190 // 2 - only TOF (can be used for dE/dx calibration)
1191 // 3 - Both dE/dx and TOF
1192 (*itTrk1)->setQuality(0);
1193 (*itTrk2)->setQuality(0);
1194
1195 // DQA
1196 // Add the line below at the end of execute(), (before return)
1197 //
1198 setFilterPassed(true);
1199 ////////////////////////////////////////////////////////////
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 ());
1206
1207 m_ee_phi_ep->Fill(p41e.phi ());
1208 m_ee_phi_em->Fill(p42e.phi ());
1209
1210 m_ee_nneu->Fill(m_nGam);
1211
1212
1213 m_ee_eemc_ep->Fill(e01);
1214 m_ee_eemc_em->Fill(e02);
1215
1216 m_ee_x_ep->Fill(ex1);
1217 m_ee_y_ep->Fill(ey1);
1218 m_ee_z_ep->Fill(ez1);
1219
1220 m_ee_x_em->Fill(ex2);
1221 m_ee_y_em->Fill(ey2);
1222 m_ee_z_em->Fill(ez2);
1223
1224
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);
1229
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);
1234
1235 m_ee_deltatof->Fill(deltatof);
1236
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);
1243
1244 if(m_writentuple) m_tuple1 -> write();
1245
1246 return StatusCode::SUCCESS;
1247
1248}
1249
1250// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1252 MsgStream log(msgSvc(), name());
1253 log << MSG::INFO << "in finalize()" << endmsg;
1254
1255 double p[13];
1256
1257 for(int i=0;i<13;i++){
1258 p[i]=(counter[i]*0.1/(counter[0]*0.1))*100;
1259 }
1260
1261 cout<<"**************************************************************************************"<<endl<<endl;
1262 cout<<"Cuts of Bhabha Selection"<<'\t'<<'\t'<<'\t'<<"Numbers"<<'\t'<<'\t'<<'\t'<<"Ratio"<<endl<<endl;
1263
1264 cout<<"Total Number Before All Cuts"<<'\t'<<'\t'<<'\t'<<counter[0]<<'\t'<<'\t'<<'\t'<<p[0]<<"%"<<endl;
1265 cout<<"Finish Good Charged Track Selection"<<'\t'<<'\t'<<counter[1]<<'\t'<<'\t'<<'\t'<<p[1]<<"%"<<endl;
1266 cout<<"Finish Particle ID"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[2]<<'\t'<<'\t'<<'\t'<<p[2]<<"%"<<endl;
1267 cout<<"Finish Good Neutral Track Selection"<<'\t'<<'\t'<<counter[3]<<'\t'<<'\t'<<'\t'<<p[3]<<"%"<<endl;
1268 cout<<"Good Neutral Track Not Zero"<<'\t'<<'\t'<<'\t'<<counter[4]<<endl;
1269 cout<<"Finish Check Good Charged Track's Info."<<'\t'<<'\t'<<counter[5]<<'\t'<<'\t'<<'\t'<<p[5]<<"%"<<endl;
1270 if(m_use_acolle) cout<<"After Acolle Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[6]<<'\t'<<'\t'<<'\t'<<p[6]<<"%"<<endl;
1271 else cout<<"No Acolle Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1272 if(m_use_acople) cout<<"After Acople Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[7]<<'\t'<<'\t'<<'\t'<<p[7]<<"%"<<endl;
1273 else cout<<"No Acople Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1274 if(m_use_eoeb) cout<<"After Eoeb Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<counter[8]<<'\t'<<'\t'<<'\t'<<p[8]<<"%"<<endl;
1275 else cout<<"No Eoeb Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1276 if(m_use_deltatof) cout<<"After Deltatof Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[9]<<'\t'<<'\t'<<'\t'<<p[9]<<"%"<<endl;
1277 else cout<<"No Deltatof Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1278 if(m_use_poeb) cout<<"After Poeb Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<counter[10]<<'\t'<<'\t'<<'\t'<<p[10]<<"%"<<endl;
1279 else cout<<"No Poeb Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1280 if(m_use_mucinfo) cout<<"After Mucinfo Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[11]<<'\t'<<'\t'<<'\t'<<p[11]<<"%"<<endl;
1281 else cout<<"No Mucinfo Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1282 if(m_use_ne) cout<<"After PID_ne Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<counter[12]<<'\t'<<'\t'<<'\t'<<p[12]<<"%"<<endl<<endl;
1283 else cout<<"No PID_ne Cut"<<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"NULL"<<'\t'<<'\t'<<'\t'<<"NULL"<<endl;
1284
1285 cout<<"**************************************************************************************"<<endl<<endl<<endl;
1286
1287 cout<<endl;
1288 return StatusCode::SUCCESS;
1289}
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
HepGeom::Point3D< double > HepPoint3D
const double mjpsi
std::vector< HepLorentzVector > Vp4
const double xmass[5]
const double mpsi2s
const double velc
const double mk
const double mpi
std::vector< int > Vint
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
const double xmass[5]
Definition: Gam4pikp.cxx:50
const double velc
Definition: Gam4pikp.cxx:51
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
************Class m_ypar INTEGER m_KeyWgt INTEGER m_KeyIHVP INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition: KK2f.h:50
INTupleSvc * ntupleSvc()
IMessageSvc * msgSvc()
@ theta1
Definition: TrkKalDeriv.h:24
StatusCode finalize()
StatusCode initialize()
DQASelBhabha(const std::string &name, ISvcLocator *pSvcLocator)
StatusCode execute()
double dy() const
double latMoment() const
Definition: DstEmcShower.h:52
double a42Moment() const
Definition: DstEmcShower.h:54
double eSeed() const
Definition: DstEmcShower.h:47
double dphi() const
Definition: DstEmcShower.h:44
double theta() const
Definition: DstEmcShower.h:38
int module() const
Definition: DstEmcShower.h:33
double e3x3() const
Definition: DstEmcShower.h:48
double dz() const
double phi() const
Definition: DstEmcShower.h:39
double dx() const
Definition: DstEmcShower.cxx:3
double secondMoment() const
Definition: DstEmcShower.h:51
double x() const
Definition: DstEmcShower.h:35
double e5x5() const
Definition: DstEmcShower.h:49
double time() const
Definition: DstEmcShower.h:50
double z() const
Definition: DstEmcShower.h:37
int numHits() const
Definition: DstEmcShower.h:30
double a20Moment() const
Definition: DstEmcShower.h:53
double energy() const
Definition: DstEmcShower.h:45
double dE() const
Definition: DstEmcShower.h:46
double dtheta() const
Definition: DstEmcShower.h:43
double y() const
Definition: DstEmcShower.h:36
const Hep3Vector emcPosition() const
Definition: DstExtTrack.h:126
const int emcVolumeNumber() const
Definition: DstExtTrack.h:132
double probPH() const
Definition: DstMdcDedx.h:66
double chiE() const
Definition: DstMdcDedx.h:59
int numTotalHits() const
Definition: DstMdcDedx.h:65
int numGoodHits() const
Definition: DstMdcDedx.h:64
double normPH() const
Definition: DstMdcDedx.h:67
double chiPi() const
Definition: DstMdcDedx.h:61
double chiK() const
Definition: DstMdcDedx.h:62
double chiMu() const
Definition: DstMdcDedx.h:60
double chiP() const
Definition: DstMdcDedx.h:63
const double dz(void) const
const double px() const
const double dr(void) const
static void setPidType(PidType pidType)
const double pz() const
const double py() const
const double fi0(void) const
const double mass() const
const double py() const
Definition: DstMdcTrack.h:56
const int charge() const
Definition: DstMdcTrack.h:53
const double px() const
Definition: DstMdcTrack.h:55
const double pz() const
Definition: DstMdcTrack.h:57
const HepVector helix() const
......
const double z() const
Definition: DstMdcTrack.h:63
const double p() const
Definition: DstMdcTrack.h:58
const double y() const
Definition: DstMdcTrack.h:62
const double x() const
Definition: DstMdcTrack.h:61
int numHits() const
Definition: DstMucTrack.h:41
int numLayers() const
Definition: DstMucTrack.h:42
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
int useTof2() const
int methodProbability() const
int useDedx() const
int onlyMuon() const
int onlyElectron() const
int onlyPion() const
int useTof1() const
void setChiMinCut(const double chi=4)
void setRecTrack(EvtRecTrack *trk)
double probKaon() const
Definition: ParticleID.h:124
void setMethod(const int method)
Definition: ParticleID.h:94
double prob(int n) const
Definition: ParticleID.h:114
double chiTof2(int n) const
void identify(const int pidcase)
Definition: ParticleID.h:103
double probMuon() const
Definition: ParticleID.h:122
double probElectron() const
Definition: ParticleID.h:121
void usePidSys(const int pidsys)
Definition: ParticleID.h:97
static ParticleID * instance()
Definition: ParticleID.cxx:22
bool IsPidInfoValid() const
double probPion() const
Definition: ParticleID.h:123
double chiTof1(int n) const
void calculate()
Definition: ParticleID.cxx:97
void init()
Definition: ParticleID.cxx:27
double probProton() const
Definition: ParticleID.h:125
double chiDedx(int n) const
RecEmcID getShowerId() const
Definition: RecEmcShower.h:55
RecEmcEnergy getEAll() const
Definition: RecEmcShower.h:92
HepSymMatrix & getZErrorE()
HepVector & getZHelixE()
bool is_barrel() const
Definition: TofHitStatus.h:26
unsigned int layer() const
Definition: TofHitStatus.h:28
bool is_cluster() const
Definition: TofHitStatus.h:25
void setStatus(unsigned int status)
bool is_counter() const
Definition: TofHitStatus.h:24
HepLorentzVector p() const
double y[1000]
double x[1000]
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:117
float ptrk
const float pi
Definition: vector3.h:133