BOSS 7.0.5
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"
10#include "VertexFit/IVertexDbSvc.h"
11#include "GaudiKernel/Bootstrap.h"
12#include "GaudiKernel/ISvcLocator.h"
13
14#include "EventModel/EventModel.h"
15#include "EventModel/Event.h"
16
17#include "EvtRecEvent/EvtRecEvent.h"
18#include "EvtRecEvent/EvtRecTrack.h"
19#include "DstEvent/TofHitStatus.h"
20#include "EventModel/EventHeader.h"
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
36#include "VertexFit/KinematicFit.h"
37#include "VertexFit/VertexFit.h"
38#include "VertexFit/IVertexDbSvc.h"
39#include "ParticleID/ParticleID.h"
40
41
42#include "DQASelBhabha/DQASelBhabha.h"
43
44#ifndef ENABLE_BACKWARDS_COMPATIBILITY
45 typedef HepGeom::Point3D<double> HepPoint3D;
46#endif
47
48#include "EmcRecEventModel/RecEmcEventModel.h"
49#include "EmcRecGeoSvc/EmcRecBarrelGeo.h"
50#include "EmcRecGeoSvc/EmcRecGeoSvc.h"
51
52
53
54#include "Identifier/TofID.h"
55
56#include "EvTimeEvent/RecEsTime.h"
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// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
125StatusCode DQASelBhabha::initialize(){
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// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
411StatusCode DQASelBhabha::execute() {
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
539 ParticleID *pid = ParticleID::instance();
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
789 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
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// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1251StatusCode DQASelBhabha::finalize() {
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}
u_cms
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: DQASelBhabha.cxx:45
mjpsi
const double mjpsi
Definition: DQASelBhabha.cxx:61
Vp4
std::vector< HepLorentzVector > Vp4
Definition: DQASelBhabha.cxx:67
xmass
const double xmass[5]
Definition: DQASelBhabha.cxx:64
mpsi2s
const double mpsi2s
Definition: DQASelBhabha.cxx:60
velc
const double velc
Definition: DQASelBhabha.cxx:65
mk
const double mk
Definition: DQASelBhabha.cxx:63
mpi
const double mpi
Definition: DQASelBhabha.cxx:62
Vint
std::vector< int > Vint
Definition: DQASelBhabha.cxx:66
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
xmass
const double xmass[5]
Definition: Gam4pikp.cxx:50
velc
const double velc
Definition: Gam4pikp.cxx:51
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
energy
************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
DQASelBhabha::finalize
StatusCode finalize()
Definition: DQASelBhabha.cxx:1251
DQASelBhabha::initialize
StatusCode initialize()
Definition: DQASelBhabha.cxx:125
DQASelBhabha::DQASelBhabha
DQASelBhabha(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DQASelBhabha.cxx:74
DQASelBhabha::execute
StatusCode execute()
Definition: DQASelBhabha.cxx:411
DstEmcShower::dy
double dy() const
Definition: DstEmcShower.cxx:11
DstEmcShower::dz
double dz() const
Definition: DstEmcShower.cxx:19
DstEmcShower::dx
double dx() const
Definition: DstEmcShower.cxx:3
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::instance
static ParticleID * instance()
Definition: ParticleID.cxx:22
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition: ParticleID.cxx:97
ParticleID::init
void init()
Definition: ParticleID.cxx:27
ParticleID::chiDedx
double chiDedx(int n) const
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
ptrk
float ptrk
Definition: prepare_data.cxx:28