BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
DQASelBhabha.cxx
Go to the documentation of this file.
1#include <vector>
2
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/AlgFactory.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include "GaudiKernel/IDataProviderSvc.h"
8#include "GaudiKernel/PropertyMgr.h"
9#include "VertexFit/IVertexDbSvc.h"
10#include "GaudiKernel/Bootstrap.h"
11#include "GaudiKernel/ISvcLocator.h"
12
13#include "EventModel/EventModel.h"
14#include "EventModel/Event.h"
15
16#include "EvtRecEvent/EvtRecEvent.h"
17#include "EvtRecEvent/EvtRecTrack.h"
18#include "DstEvent/TofHitStatus.h"
19#include "EventModel/EventHeader.h"
20
21#include "TMath.h"
22#include "GaudiKernel/INTupleSvc.h"
23#include "GaudiKernel/NTuple.h"
24#include "GaudiKernel/Bootstrap.h"
25#include "GaudiKernel/IHistogramSvc.h"
26#include "CLHEP/Vector/ThreeVector.h"
27#include "CLHEP/Vector/LorentzVector.h"
28#include "CLHEP/Vector/TwoVector.h"
29
30using CLHEP::Hep3Vector;
31using CLHEP::Hep2Vector;
32using CLHEP::HepLorentzVector;
33#include "CLHEP/Geometry/Point3D.h"
34
35#include "VertexFit/KinematicFit.h"
36#include "VertexFit/VertexFit.h"
37#include "VertexFit/IVertexDbSvc.h"
38#include "ParticleID/ParticleID.h"
39
40//
41#include "DQASelBhabha/DQASelBhabha.h"
42
43#ifndef ENABLE_BACKWARDS_COMPATIBILITY
44 typedef HepGeom::Point3D<double> HepPoint3D;
45#endif
46using CLHEP::HepLorentzVector;
47const double mpsi2s=3.68609;
48const double mpi = 0.13957;
49const double mk = 0.493677;
50const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
51const double velc = 299.792458; // tof path unit in mm
52typedef std::vector<int> Vint;
53typedef std::vector<HepLorentzVector> Vp4;
54//declare one counter
55static int counter[10]={0,0,0,0,0,0,0,0,0,0};
56static int nbhabha=0;
57/////////////////////////////////////////////////////////////////////////////
58
59DQASelBhabha::DQASelBhabha(const std::string& name, ISvcLocator* pSvcLocator) :
60 Algorithm(name, pSvcLocator) {
61
62 //Declare the properties
63 declareProperty("writentuple",m_writentuple = false);
64 declareProperty("ecms",m_ecms = 3.097);
65 declareProperty("beamangle",m_beamangle = 0.022);
66 declareProperty("Vr0cut", m_vr0cut=1.0);
67 declareProperty("Vz0cut", m_vz0cut=8.0);
68 declareProperty("Coscut", m_coscut=0.93);
69
70 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
71 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
72 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
73 declareProperty("GammaTrkCut", m_gammaTrkCut=20.0);
74 declareProperty("GammaTLCut", m_gammatlCut=0);
75 declareProperty("GammaTHCut", m_gammathCut=60);
76
77
78 declareProperty ("acoll_e_cut", m_acoll_e_cut=6.);
79 declareProperty ("acopl_e_cut", m_acopl_e_cut=6.);
80 declareProperty ("poeb_e_cut", m_poeb_e_cut=0.5);
81 declareProperty ("dtof_e_cut", m_dtof_e_cut=4.);
82 declareProperty ("eoeb_e_cut", m_eoeb_e_cut=0.4);
83 declareProperty ("etotal_e_cut", m_etotal_e_cut=0.8);
84 declareProperty ("tpoeb_e_cut", m_tpoeb_e_cut=0.95);
85 declareProperty ("tptotal_e_cut", m_tptotal_e_cut=0.16);
86 declareProperty ("tetotal_e_cut", m_tetotal_e_cut=0.65);
87
88
89 //normally, MDC+EMC, otherwise EMC only
90 declareProperty ("m_useEMConly", m_useEMConly= false);
91 declareProperty ("m_usePID", m_usePID= false);// sub-system is under study
92 declareProperty ("m_useMDC", m_useMDC= true);
93 declareProperty ("m_useDEDX", m_useDEDX= false);// not used
94 declareProperty ("m_useTOF", m_useTOF= false);//sub-system is under study
95 declareProperty ("m_useEMC", m_useEMC= true);
96 declareProperty ("m_useMUC", m_useMUC= false);// efficiency
97
98}
99
100
101// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
102StatusCode DQASelBhabha::initialize(){
103 MsgStream log(msgSvc(), name());
104
105 log << MSG::INFO << "in initialize()" << endmsg;
106 StatusCode status;
107 status = service("THistSvc", m_thistsvc);
108 if(status.isFailure() ){
109 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
110 return status;
111 }
112
113
114 m_ee_mass = new TH1F( "ee_mass", "ee_mass", 80, m_ecms-0.3, m_ecms+0.5 );
115 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_mass", m_ee_mass);
116 m_ee_acoll = new TH1F( "ee_acoll", "ee_acoll", 60, 0, 6 );
117 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_acoll", m_ee_acoll);
118 m_ee_eop_ep = new TH1F( "ee_eop_ep", "ee_eop_ep", 100,0.4,1.4 );
119 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eop_ep", m_ee_eop_ep);
120 m_ee_eop_em = new TH1F( "ee_eop_em", "ee_eop_em", 100,0.4,1.4 );
121 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eop_em", m_ee_eop_em);
122 m_ee_costheta_ep = new TH1F( "ee_costheta_ep", "ee_costheta_ep", 100,-1,1 );
123 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_costheta_ep", m_ee_costheta_ep);
124 m_ee_costheta_em = new TH1F( "ee_costheta_em", "ee_costheta_em", 100,-1,1 );
125 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_costheta_em", m_ee_costheta_em);
126
127 m_ee_phi_ep = new TH1F( "ee_phi_ep", "ee_phi_ep", 120,-3.2,3.2 );
128 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_phi_ep", m_ee_phi_ep);
129 m_ee_phi_em = new TH1F( "ee_phi_em", "ee_phi_em", 120,-3.2,3.2 );
130 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_phi_em", m_ee_phi_em);
131
132 m_ee_nneu = new TH1I( "ee_nneu", "ee_nneu", 5,0,5 );
133 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_nneu", m_ee_nneu);
134
135 m_ee_eemc_ep=new TH1F("ee_eemc_ep","ee_eemc_ep",100,1.0,2.0);
136 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eemc_ep", m_ee_eemc_ep);
137 m_ee_eemc_em=new TH1F("ee_eemc_em","ee_eemc_em",100,1.0,2.0);
138 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eemc_em", m_ee_eemc_em);
139 m_ee_x_ep=new TH1F("ee_x_ep","ee_x_ep",100,-1.0,1.0);
140 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_x_ep", m_ee_x_ep);
141 m_ee_y_ep=new TH1F("ee_y_ep","ee_y_ep",100,-1.0,1.0);
142 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_y_ep", m_ee_y_ep);
143 m_ee_z_ep=new TH1F("ee_z_ep","ee_z_ep",100,-10.0,10.0);
144 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_z_ep", m_ee_z_ep);
145 m_ee_x_em=new TH1F("ee_x_em","ee_x_em",100,-1.0,1.0);
146 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_x_em", m_ee_x_em);
147 m_ee_y_em=new TH1F("ee_y_em","ee_y_em",100,-1.0,1.0);
148 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_y_em", m_ee_y_em);
149 m_ee_z_em=new TH1F("ee_z_em","ee_z_em",100,-10.0,10.0);
150 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_z_em", m_ee_z_em);
151
152 m_ee_px_ep=new TH1F("ee_px_ep","ee_px_ep",200,-2.0,2.0);
153 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_px_ep", m_ee_px_ep);
154 m_ee_py_ep=new TH1F("ee_py_ep","ee_py_ep",200,-2.0,2.0);
155 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_py_ep", m_ee_py_ep);
156 m_ee_pz_ep=new TH1F("ee_pz_ep","ee_pz_ep",200,-2.0,2.0);
157 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pz_ep", m_ee_pz_ep);
158 m_ee_p_ep=new TH1F("ee_p_ep","ee_p_ep",100,1.0,2.0);
159 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_p_ep", m_ee_p_ep);
160 m_ee_px_em=new TH1F("ee_px_em","ee_px_em",100,-2.0,2.0);
161 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_px_em", m_ee_px_em);
162 m_ee_py_em=new TH1F("ee_py_em","ee_py_em",100,-2.0,2.0);
163 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_py_em", m_ee_py_em);
164 m_ee_pz_em=new TH1F("ee_pz_em","ee_pz_em",100,-2.0,2.0);
165 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pz_em", m_ee_pz_em);
166 m_ee_p_em=new TH1F("ee_p_em","ee_p_em",100,1.0,2.0);
167 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_p_em", m_ee_p_em);
168 m_ee_deltatof=new TH1F("ee_deltatof","ee_deltatof",50,0.0,10.0);
169 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_deltatof", m_ee_deltatof);
170
171 m_ee_pidchidedx_ep=new TH1F("ee_pidchidedx_ep","ee_pidchidedx_ep",160,-4,4);
172 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchidedx_ep", m_ee_pidchidedx_ep);
173 m_ee_pidchidedx_em=new TH1F("ee_pidchidedx_em","ee_pidchidedx_em",160,-4,4);
174 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchidedx_em", m_ee_pidchidedx_em);
175 m_ee_pidchitof1_ep=new TH1F("ee_pidchitof1_ep","ee_pidchitof1_ep",160,-4,4);
176 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof1_ep", m_ee_pidchitof1_ep);
177 m_ee_pidchitof1_em=new TH1F("ee_pidchitof1_em","ee_pidchitof1_em",160,-4,4);
178 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof1_em", m_ee_pidchitof1_em);
179 m_ee_pidchitof2_ep=new TH1F("ee_pidchitof2_ep","ee_pidchitof2_ep",160,-4,4);
180 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof2_ep", m_ee_pidchitof2_ep);
181 m_ee_pidchitof2_em=new TH1F("ee_pidchitof2_em","ee_pidchitof2_em",160,-4,4);
182 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof2_em", m_ee_pidchitof2_em);
183
184
185
186
187
188
189
190 NTuplePtr nt1(ntupleSvc(), "DQAFILE/Bhabha");
191 if ( nt1 ) m_tuple1 = nt1;
192 else {
193 m_tuple1 = ntupleSvc()->book ("DQAFILE/Bhabha", CLID_ColumnWiseTuple, "N-Tuple");
194 if ( m_tuple1 ) {
195 status = m_tuple1->addItem ("run", m_run);
196 status = m_tuple1->addItem ("rec", m_rec);
197 status = m_tuple1->addItem ("Nchrg", m_ncharg);
198 status = m_tuple1->addItem ("Nneu", m_nneu,0,40);
199 status = m_tuple1->addItem ("NGch", m_ngch, 0, 40);
200 status = m_tuple1->addItem ("NGam", m_nGam);
201
202
203 status = m_tuple1->addItem ("bhabhatag", m_bhabhatag);
204
205 status = m_tuple1->addItem ("acoll", m_acoll);
206 status = m_tuple1->addItem ("acopl", m_acopl);
207 status = m_tuple1->addItem ("deltatof", m_deltatof);
208 status = m_tuple1->addItem ("eop1", m_eop1);
209 status = m_tuple1->addItem ("eop2", m_eop2);
210 status = m_tuple1->addItem ("eoeb1", m_eoeb1);
211 status = m_tuple1->addItem ("eoeb2", m_eoeb2);
212 status = m_tuple1->addItem ("poeb1", m_poeb1);
213 status = m_tuple1->addItem ("poeb2", m_poeb2);
214 status = m_tuple1->addItem ("etoeb1", m_etoeb1);
215 status = m_tuple1->addItem ("etoeb2", m_etoeb2);
216 status = m_tuple1->addItem ("mucinfo1", m_mucinfo1);
217 status = m_tuple1->addItem ("mucinfo2", m_mucinfo2);
218
219
220 status = m_tuple1->addIndexedItem ("delang",m_nneu, m_delang);
221 status = m_tuple1->addIndexedItem ("delphi",m_nneu, m_delphi);
222 status = m_tuple1->addIndexedItem ("delthe",m_nneu, m_delthe);
223 status = m_tuple1->addIndexedItem ("npart",m_nneu, m_npart);
224 status = m_tuple1->addIndexedItem ("nemchits",m_nneu, m_nemchits);
225 status = m_tuple1->addIndexedItem ("module",m_nneu, m_module);
226 status = m_tuple1->addIndexedItem ("x",m_nneu, m_x);
227 status = m_tuple1->addIndexedItem ("y",m_nneu, m_y);
228 status = m_tuple1->addIndexedItem ("z",m_nneu, m_z);
229 status = m_tuple1->addIndexedItem ("px",m_nneu, m_px);
230 status = m_tuple1->addIndexedItem ("py",m_nneu, m_py);
231 status = m_tuple1->addIndexedItem ("pz",m_nneu, m_pz);
232 status = m_tuple1->addIndexedItem ("theta",m_nneu, m_theta);
233 status = m_tuple1->addIndexedItem ("phi",m_nneu, m_phi);
234 status = m_tuple1->addIndexedItem ("dx",m_nneu, m_dx);
235 status = m_tuple1->addIndexedItem ("dy",m_nneu, m_dy);
236 status = m_tuple1->addIndexedItem ("dz",m_nneu, m_dz);
237 status = m_tuple1->addIndexedItem ("dtheta",m_nneu, m_dtheta);
238 status = m_tuple1->addIndexedItem ("dphi",m_nneu, m_dphi);
239 status = m_tuple1->addIndexedItem ("energy",m_nneu, m_energy);
240 status = m_tuple1->addIndexedItem ("dE",m_nneu, m_dE);
241 status = m_tuple1->addIndexedItem ("eSeed",m_nneu, m_eSeed);
242 status = m_tuple1->addIndexedItem ("nSeed",m_nneu, m_nSeed);
243 status = m_tuple1->addIndexedItem ("e3x3",m_nneu, m_e3x3);
244 status = m_tuple1->addIndexedItem ("e5x5",m_nneu, m_e5x5);
245 status = m_tuple1->addIndexedItem ("secondMoment",m_nneu, m_secondMoment);
246 status = m_tuple1->addIndexedItem ("latMoment",m_nneu, m_latMoment);
247 status = m_tuple1->addIndexedItem ("a20Moment",m_nneu, m_a20Moment);
248 status = m_tuple1->addIndexedItem ("a42Moment",m_nneu, m_a42Moment);
249 status = m_tuple1->addIndexedItem ("getTime",m_nneu, m_getTime);
250 status = m_tuple1->addIndexedItem ("getEAll",m_nneu, m_getEAll);
251
252
253
254 status = m_tuple1->addIndexedItem("charge", m_ngch, m_charge);
255 status = m_tuple1->addIndexedItem ("vx", m_ngch, m_vx0);
256 status = m_tuple1->addIndexedItem ("vy", m_ngch, m_vy0);
257 status = m_tuple1->addIndexedItem ("vz", m_ngch, m_vz0);
258
259
260 status = m_tuple1->addIndexedItem ("px", m_ngch, m_px) ;
261 status = m_tuple1->addIndexedItem ("py", m_ngch, m_py) ;
262 status = m_tuple1->addIndexedItem ("pz", m_ngch, m_pz) ;
263 status = m_tuple1->addIndexedItem ("p", m_ngch, m_p) ;
264
265
266
267 status = m_tuple1->addIndexedItem ("kal_vx", m_ngch, m_kal_vx0);
268 status = m_tuple1->addIndexedItem ("kal_vy", m_ngch, m_kal_vy0);
269 status = m_tuple1->addIndexedItem ("kal_vz", m_ngch, m_kal_vz0);
270
271
272 status = m_tuple1->addIndexedItem ("kal_px", m_ngch, m_kal_px) ;
273 status = m_tuple1->addIndexedItem ("kal_py", m_ngch, m_kal_py) ;
274 status = m_tuple1->addIndexedItem ("kal_pz", m_ngch, m_kal_pz) ;
275 status = m_tuple1->addIndexedItem ("kal_p", m_ngch, m_kal_p) ;
276
277
278 status = m_tuple1->addIndexedItem ("probPH" , m_ngch, m_probPH) ;
279 status = m_tuple1->addIndexedItem ("normPH" , m_ngch, m_normPH) ;
280 status = m_tuple1->addIndexedItem ("chie" , m_ngch, m_chie) ;
281 status = m_tuple1->addIndexedItem ("chimu" , m_ngch, m_chimu) ;
282 status = m_tuple1->addIndexedItem ("chipi" , m_ngch, m_chipi) ;
283 status = m_tuple1->addIndexedItem ("chik" , m_ngch, m_chik) ;
284 status = m_tuple1->addIndexedItem ("chip" , m_ngch, m_chip) ;
285 status = m_tuple1->addIndexedItem ("ghit" , m_ngch, m_ghit) ;
286 status = m_tuple1->addIndexedItem ("thit" , m_ngch, m_thit) ;
287
288 status = m_tuple1->addIndexedItem ("e_emc" , m_ngch, m_e_emc) ;
289 status = m_tuple1->addIndexedItem ("phi_emc" , m_ngch, m_phi_emc) ;
290 status = m_tuple1->addIndexedItem ("theta_emc" , m_ngch, m_theta_emc) ;
291
292 status = m_tuple1->addIndexedItem ("nhit_muc" , m_ngch, m_nhit_muc) ;
293 status = m_tuple1->addIndexedItem ("nlay_muc" , m_ngch, m_nlay_muc) ;
294 status = m_tuple1->addIndexedItem ("t_btof" , m_ngch, m_t_btof );
295 status = m_tuple1->addIndexedItem ("t_etof" , m_ngch, m_t_etof );
296 status = m_tuple1->addIndexedItem ("qual_etof" , m_ngch, m_qual_etof );
297 status = m_tuple1->addIndexedItem ("tof_etof" , m_ngch, m_tof_etof );
298 status = m_tuple1->addIndexedItem ("te_etof" , m_ngch, m_te_etof );
299 status = m_tuple1->addIndexedItem ("tmu_etof" , m_ngch, m_tmu_etof );
300 status = m_tuple1->addIndexedItem ("tpi_etof" , m_ngch, m_tpi_etof );
301 status = m_tuple1->addIndexedItem ("tk_etof" , m_ngch, m_tk_etof );
302 status = m_tuple1->addIndexedItem ("tp_etof" , m_ngch, m_tp_etof );
303
304 status = m_tuple1->addIndexedItem ("qual_btof1", m_ngch, m_qual_btof1 );
305 status = m_tuple1->addIndexedItem ("tof_btof1" , m_ngch, m_tof_btof1 );
306 status = m_tuple1->addIndexedItem ("te_btof1" , m_ngch, m_te_btof1 );
307 status = m_tuple1->addIndexedItem ("tmu_btof1" , m_ngch, m_tmu_btof1 );
308 status = m_tuple1->addIndexedItem ("tpi_btof1" , m_ngch, m_tpi_btof1 );
309 status = m_tuple1->addIndexedItem ("tk_btof1" , m_ngch, m_tk_btof1 );
310 status = m_tuple1->addIndexedItem ("tp_btof1" , m_ngch, m_tp_btof1 );
311
312 status = m_tuple1->addIndexedItem ("qual_btof2", m_ngch, m_qual_btof2 );
313 status = m_tuple1->addIndexedItem ("tof_btof2" , m_ngch, m_tof_btof2 );
314 status = m_tuple1->addIndexedItem ("te_btof2" , m_ngch, m_te_btof2 );
315 status = m_tuple1->addIndexedItem ("tmu_btof2" , m_ngch, m_tmu_btof2 );
316 status = m_tuple1->addIndexedItem ("tpi_btof2" , m_ngch, m_tpi_btof2 );
317 status = m_tuple1->addIndexedItem ("tk_btof2" , m_ngch, m_tk_btof2 );
318 status = m_tuple1->addIndexedItem ("tp_btof2" , m_ngch, m_tp_btof2 );
319 status = m_tuple1->addIndexedItem ("pidcode" , m_ngch, m_pidcode);
320 status = m_tuple1->addIndexedItem ("pidprob" , m_ngch, m_pidprob);
321 status = m_tuple1->addIndexedItem ("pidchiDedx" , m_ngch, m_pidchiDedx);
322 status = m_tuple1->addIndexedItem ("pidchiTof1" , m_ngch, m_pidchiTof1);
323 status = m_tuple1->addIndexedItem ("pidchiTof2" , m_ngch, m_pidchiTof2);
324
325 status = m_tuple1->addItem ("px_cms_ep", m_px_cms_ep); //momentum of electron+
326 status = m_tuple1->addItem ("py_cms_ep", m_py_cms_ep); //momentum of electron+
327 status = m_tuple1->addItem ("pz_cms_ep", m_pz_cms_ep); //momentum of electron+
328 status = m_tuple1->addItem ("e_cms_ep", m_e_cms_ep); //momentum of electron+
329 status = m_tuple1->addItem ("cos_ep", m_cos_ep); //momentum of electron+
330 status = m_tuple1->addItem ("px_cms_em", m_px_cms_em); //momentum of electron-
331 status = m_tuple1->addItem ("py_cms_em", m_py_cms_em); //momentum of electron-
332 status = m_tuple1->addItem ("pz_cms_em", m_pz_cms_em); //momentum of electron-
333 status = m_tuple1->addItem ("e_cms_em", m_e_cms_em); //momentum of electron-
334 status = m_tuple1->addItem ("cos_em", m_cos_em); //momentum of electron-
335 status = m_tuple1->addItem ("mass_ee", m_mass_ee); //
336 status = m_tuple1->addItem ("emax", m_emax); //
337 status = m_tuple1->addItem ("esum", m_esum); //
338 status = m_tuple1->addItem ( "npip", m_npip );
339 status = m_tuple1->addItem ( "npim", m_npim );
340 status = m_tuple1->addItem ( "nkp", m_nkp );
341 status = m_tuple1->addItem ( "nkm", m_nkm );
342 status = m_tuple1->addItem ( "np", m_np );
343 status = m_tuple1->addItem ( "npb", m_npb );
344
345 status = m_tuple1->addItem ( "nep", m_nep );
346 status = m_tuple1->addItem ( "nem", m_nem );
347 status = m_tuple1->addItem ( "nmup", m_nmup );
348 status = m_tuple1->addItem ( "nmum", m_nmum );
349
350 }
351 else {
352 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
353 return StatusCode::FAILURE;
354 }
355 }
356
357 //
358 //--------end of book--------
359 //
360
361 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
362 return StatusCode::SUCCESS;
363
364
365
366
367}
368
369// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
370StatusCode DQASelBhabha::execute() {
371
372 const double beamEnergy = m_ecms/2.;
373 const HepLorentzVector p_cms(m_ecms*sin(m_beamangle*0.5),0.0,0.0,m_ecms);
374 const Hep3Vector u_cms = -p_cms.boostVector();
375 MsgStream log(msgSvc(), name());
376 log << MSG::INFO << "in execute()" << endreq;
377
378 setFilterPassed(false);
379
380 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
381 if (!eventHeader)
382 {
383 log << MSG::FATAL << "Could not find Event Header" << endreq;
384 return StatusCode::SUCCESS;
385 }
386
387 m_run = eventHeader->runNumber();
388 m_rec = eventHeader->eventNumber();
389
390
391
392
393 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
394 if (!evtRecEvent)
395 {
396 log << MSG::FATAL << "Could not find EvtRecEvent" << endreq;
397 return StatusCode::SUCCESS;
398 }
399 log << MSG::INFO <<"ncharg, nneu, tottks = "
400 << evtRecEvent->totalCharged() << " , "
401 << evtRecEvent->totalNeutral() << " , "
402 << evtRecEvent->totalTracks() <<endreq;
403 // if(evtRecEvent->totalNeutral()>30)return sc;
404 m_ncharg = evtRecEvent->totalCharged();
405
406 m_nneu = evtRecEvent->totalNeutral();
407
408
409
410 HepPoint3D vx(0., 0., 0.);
411 HepSymMatrix Evx(3, 0);
412 IVertexDbSvc* vtxsvc;
413 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
414 if(vtxsvc->isVertexValid()){
415 double* dbv = vtxsvc->PrimaryVertex();
416 double* vv = vtxsvc->SigmaPrimaryVertex();
417 // if (m_reader.isRunNumberValid( m_run)) {
418 // HepVector dbv = m_reader.PrimaryVertex( m_run);
419 // HepVector vv = m_reader.SigmaPrimaryVertex( m_run);
420 vx.setX(dbv[0]);
421 vx.setY(dbv[1]);
422 vx.setZ(dbv[2]);
423 Evx[0][0]=vv[0]*vv[0];
424 Evx[0][1]=vv[0]*vv[1];
425 Evx[1][1]=vv[1]*vv[1];
426 Evx[1][2]=vv[1]*vv[2];
427 Evx[2][2]=vv[2]*vv[2];
428 }
429
430 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
431 if (!evtRecTrkCol)
432 {
433 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endreq;
434 return StatusCode::SUCCESS;
435 }
436 Vint iGood;
437 iGood.clear();
438
439 int nCharge = 0;
440
441 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
442 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
443 if(!(*itTrk)->isMdcTrackValid()) continue;
444 if(!(*itTrk)->isMdcKalTrackValid()) continue;
445
446 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
447 double pch=mdcTrk->p();
448 double x0=mdcTrk->x();
449 double y0=mdcTrk->y();
450 double z0=mdcTrk->z();
451 double phi0=mdcTrk->helix(1);
452 double xv=vx.x();
453 double yv=vx.y();
454 double zv=vx.z();
455 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
456 double m_vx0 = x0;
457 double m_vy0 = y0;
458 double m_vz0 = z0;
459 double m_vr0 = Rxy;
460 if(fabs(z0) >= m_vz0cut) continue;
461 if(fabs(Rxy) >= m_vr0cut) continue;
462
463
464 if(fabs(m_vz0) >= m_vz0cut) continue;
465 if(m_vr0 >= m_vr0cut) continue;
466
467 // double cost = cos(mdcTrk->theta());
468 // if(fabs(cost) >= m_coscut ) continue;
469// iGood.push_back((*itTrk)->trackId());
470 iGood.push_back(i);
471 nCharge += mdcTrk->charge();
472
473 }
474
475
476
477
478
479 //
480 // Finish Good Charged Track Selection
481 //
482 int nGood = iGood.size();
483 m_ngch=nGood;
484 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
485
486 if((nGood != 2)||(nCharge!=0)){
487 return StatusCode::SUCCESS;
488 }
489
490 counter[1]++;
491
492 //
493 // Particle ID
494 //
495 Vint ipip, ipim, iep,iem,imup,imum;
496 ipip.clear();
497 ipim.clear();
498 iep.clear();
499 iem.clear();
500 imup.clear();
501 imum.clear();
502
503
504 ParticleID *pid = ParticleID::instance();
505 for(int i = 0; i < m_ngch; i++) {
506 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
507 // if(pid) delete pid;
508 pid->init();
509 pid->setMethod(pid->methodProbability());
510 pid->setChiMinCut(4);
511 pid->setRecTrack(*itTrk);
512 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());//|pid->useEmc()|pid->useMuc()); // use PID sub-system
513 pid->identify(pid->onlyElectron()|pid->onlyMuon()|pid->onlyPion()); // seperater Pion/Kaon/Proton
514 pid->calculate();
515 if(!(pid->IsPidInfoValid())) continue;
516 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
517 /// RecMdcKalTrack* mdcKalTrk = 0 ;
518 /// if((*itTrk)->isMdcKalTrackValid()) mdcKalTrk = (*itTrk)->mdcKalTrack();
519 double prob_pi = pid->probPion();
520 double prob_K = pid->probKaon();
521 double prob_p = pid->probProton();
522 double prob_e = pid->probElectron();
523 double prob_mu = pid->probMuon();
524 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
525 HepLorentzVector ptrk;
526 ptrk.setPx(mdcTrk->px()) ;
527 ptrk.setPy(mdcTrk->py()) ;
528 ptrk.setPz(mdcTrk->pz()) ;
529 double p3 = ptrk.mag() ;
530
531 m_pidcode[i]=0;
532 m_pidprob[i]=pid->prob(0);
533 m_pidchiDedx[i]=pid->chiDedx(0);
534 m_pidchiTof1[i]=pid->chiTof1(0);
535 m_pidchiTof2[i]=pid->chiTof2(0);
536 if(mdcTrk->charge() > 0) {
537 iep.push_back(iGood[i]);
538
539 }
540 if (mdcTrk->charge() < 0) {
541 iem.push_back(iGood[i]);
542
543 }
544
545
546
547
548 }
549
550 m_nep = iep.size() ;
551 m_nem = iem.size() ;
552 m_nmup = imup.size() ;
553 m_nmum = imum.size() ;
554
555 counter[2]++;
556
557 //
558 // Good neutral track selection
559 //
560 Vint iGam;
561 iGam.clear();
562 int iphoton=0;
563 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
564 if(i>=evtRecTrkCol->size())break;
565 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
566 if(!(*itTrk)->isEmcShowerValid()) continue;
567 RecEmcShower *emcTrk = (*itTrk)->emcShower();
568 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
569
570 RecEmcID showerId = emcTrk->getShowerId();
571 unsigned int npart = EmcID::barrel_ec(showerId);
572 int n = emcTrk->numHits();
573 int module=emcTrk->module();
574 double x = emcTrk->x();
575 double y = emcTrk->y();
576 double z = emcTrk->z();
577 double dx = emcTrk->dx();
578 double dy = emcTrk->dy();
579 double dth = emcTrk->dtheta();
580 double dph = emcTrk->dphi();
581 double dz = emcTrk->dz();
582 double energy = emcTrk->energy();
583 double dE = emcTrk->dE();
584 double eSeed = emcTrk->eSeed();
585 double e3x3 = emcTrk->e3x3();
586 double e5x5 = emcTrk->e5x5();
587 double secondMoment = emcTrk->secondMoment();
588 double latMoment = emcTrk->latMoment();
589 double getTime = emcTrk->time();
590 double getEAll = emcTrk->getEAll();
591 double a20Moment = emcTrk->a20Moment();
592 double a42Moment = emcTrk->a42Moment();
593 // int phigap=emcTrk->PhiGap();
594 // int thetagap=emcTrk->ThetaGap();
595 // double getETof2x1 = emcTrk->getETof2x1();
596 // double getETof2x3 = emcTrk->getETof2x3();
597 // double getELepton = emcTrk->getELepton();
598 double nseed=0;//(emcTrk->getCluster() )->getSeedSize() ;
599 HepPoint3D EmcPos(x,y,z);
600 m_nemchits[iphoton]=n;
601 m_npart[iphoton]=npart;
602 m_module[iphoton]=module;
603 m_theta[iphoton]=EmcPos.theta();
604 m_phi[iphoton]=EmcPos.phi();
605 m_x[iphoton]=x;
606 m_y[iphoton]=y;
607 m_z[iphoton]=z;
608 m_dx[iphoton]=dx;
609 m_dy[iphoton]=dy;
610 m_dz[iphoton]=dz;
611 m_dtheta[iphoton]=dth;
612 m_dphi[iphoton]=dph;
613 m_energy[iphoton]=energy;
614 m_dE[iphoton]=dE;
615 m_eSeed[iphoton]=eSeed;
616 m_nSeed[iphoton]=nseed;
617 m_e3x3[iphoton]=e3x3;
618 m_e5x5[iphoton]=e5x5;
619 m_secondMoment[iphoton]=secondMoment;
620 m_latMoment[iphoton]=latMoment;
621 m_getTime[iphoton]=getTime;
622 m_getEAll[iphoton]=getEAll;
623 m_a20Moment[iphoton]=a20Moment;
624 m_a42Moment[iphoton]=a42Moment;
625
626 // m_getELepton[iphoton]=getELepton;
627 // m_getETof2x1[iphoton]=getETof2x1;
628 // m_getETof2x3[iphoton]=getETof2x3;
629 // m_PhiGap[iphoton]=phigap;
630 // m_ThetaGap[iphoton]=thetagap;
631 double dthe = 200.;
632 double dphi = 200.;
633 double dang = 200.;
634
635 // find the nearest charged track
636 for(int j = 0; j < nGood; j++) {
637
638
639 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() +iGood[j];
640 if (!(*jtTrk)->isMdcTrackValid()) continue;
641 RecMdcTrack *jtmdcTrk = (*jtTrk)->mdcTrack();
642 double jtcharge = jtmdcTrk->charge();
643 if(!(*jtTrk)->isExtTrackValid()) continue;
644 RecExtTrack *extTrk = (*jtTrk)->extTrack();
645 if(extTrk->emcVolumeNumber() == -1) continue;
646 Hep3Vector extpos = extTrk->emcPosition();
647 // double ctht = extpos.cosTheta(emcpos);
648 double angd = extpos.angle(emcpos);
649 double thed = extpos.theta() - emcpos.theta();
650 double phid = extpos.deltaPhi(emcpos);
651 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
652 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
653
654 if(fabs(thed) < fabs(dthe)) dthe = thed;
655 if(fabs(phid) < fabs(dphi)) dphi = phid;
656 if(angd < dang) dang = angd;
657
658 }
659
660
661
662 //
663 // good photon cut will be set here
664 //
665
666 dthe = dthe * 180 / (CLHEP::pi);
667 dphi = dphi * 180 / (CLHEP::pi);
668 dang = dang * 180 / (CLHEP::pi);
669 double eraw = emcTrk->energy();
670 double phi = emcTrk->phi();
671 double the = emcTrk->theta();
672
673 m_delphi[iphoton]=dphi;
674 m_delthe[iphoton]=dthe;
675 m_delang[iphoton]=dang;
676 if(energy < m_energyThreshold) continue;
677 if(getTime>m_gammathCut||getTime<m_gammatlCut)continue;
678 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
679 if(dang< m_gammaTrkCut) continue;
680 iphoton++;
681 iGam.push_back(i);
682 if(iphoton>=40)return StatusCode::SUCCESS;
683 }
684
685 int nGam = iGam.size();
686 m_nGam=nGam;
687 // std::cout << "num Good Photon " << m_nGam << " , " <<evtRecEvent->totalNeutral()<<std::endl;
688 //std::cout<<"dbg_4"<<std::endl;
689 counter[3]++;
690
691 double egam_ext=0;
692 double ex_gam=0;
693 double ey_gam=0;
694 double ez_gam=0;
695 double et_gam=0;
696 double e_gam=0;
697 for(int i = 0; i < m_nGam; i++) {
698 EvtRecTrackIterator itTrk = evtRecTrkCol->begin()+ iGam[i];
699 if(!(*itTrk)->isEmcShowerValid()) continue;
700 RecEmcShower* emcTrk = (*itTrk)->emcShower();
701 double eraw = emcTrk->energy();
702 double phi = emcTrk->phi();
703 double the = emcTrk->theta();
704 HepLorentzVector ptrk;
705 ex_gam+=eraw*sin(the)*cos(phi);
706 ey_gam+=eraw*sin(the)*sin(phi);
707 ez_gam+=eraw*cos(the);
708 et_gam+=eraw*sin(the);
709 e_gam+=eraw ;
710 if(eraw>=egam_ext)
711 {
712 egam_ext=eraw;
713 }
714
715 }
716
717
718
719
720
721 double px_had=0;
722 double py_had=0;
723 double pz_had=0;
724 double pt_had=0;
725 double p_had=0;
726 double e_had=0;
727 //
728 // check good charged track's infomation
729 //
730 int ii ;
731 m_e_emc[0]=-0.1;
732 m_e_emc[1]=-0.1;
733 for(int i = 0; i < m_ngch; i++ ){
734
735 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
736
737 if(!(*itTrk)->isMdcTrackValid()) continue; // MDC information
738 if(!(*itTrk)->isMdcKalTrackValid()) continue;
739 // if(!(*itTrk)->isEmcShowerValid()) return StatusCode::SUCCESS;///dbg
740 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
741 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
742 ii=i;
743
744
745 m_charge[ii] = mdcTrk->charge();
746 m_vx0[ii] = mdcTrk->x();
747 m_vy0[ii] = mdcTrk->y();
748 m_vz0[ii] = mdcTrk->z();
749
750
751 m_px[ii] = mdcTrk->px();
752 m_py[ii] = mdcTrk->py();
753 m_pz[ii] = mdcTrk->pz();
754 m_p[ii] = mdcTrk->p();
755
756
757 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
758
759
760 /// if(m_pidcode[ii]==3)mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
761 /// if(m_pidcode[ii]==4)mdcKalTrk->setPidType(RecMdcKalTrack::proton);
762 m_kal_vx0[ii] = mdcKalTrk->x();
763 m_kal_vy0[ii] = mdcKalTrk->y();
764 m_kal_vz0[ii] = mdcKalTrk->z();
765
766
767 m_kal_px[ii] = mdcKalTrk->px();
768 m_kal_py[ii] = mdcKalTrk->py();
769 m_kal_pz[ii] = mdcKalTrk->pz();
770 m_kal_p[ii] = mdcKalTrk->p();
771
772
773 px_had+=mdcKalTrk->px();
774 py_had+=mdcKalTrk->py();
775 pz_had+=mdcKalTrk->pz();
776 pt_had+=mdcKalTrk->pxy();
777 p_had+=mdcKalTrk->p();
778 e_had+=sqrt(mdcKalTrk->p()*mdcKalTrk->p()+mdcKalTrk->mass()*mdcKalTrk->mass());
779
780 double ptrk = mdcKalTrk->p() ;
781
782
783 if((*itTrk)->isMdcDedxValid()) { // DEDX information
784
785 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
786 m_probPH[ii]= dedxTrk->probPH();
787 m_normPH[ii]= dedxTrk->normPH();
788
789 m_chie[ii] = dedxTrk->chiE();
790 m_chimu[ii] = dedxTrk->chiMu();
791 m_chipi[ii] = dedxTrk->chiPi();
792 m_chik[ii] = dedxTrk->chiK();
793 m_chip[ii] = dedxTrk->chiP();
794 m_ghit[ii] = dedxTrk->numGoodHits();
795 m_thit[ii] = dedxTrk->numTotalHits();
796 }
797
798 if((*itTrk)->isEmcShowerValid()) {
799
800 RecEmcShower *emcTrk = (*itTrk)->emcShower();
801 m_e_emc[ii] = emcTrk->energy();
802 m_phi_emc[ii] = emcTrk->phi();
803 m_theta_emc[ii] = emcTrk->theta();
804 }
805
806
807 if((*itTrk)->isMucTrackValid()){
808
809 RecMucTrack* mucTrk = (*itTrk)->mucTrack() ;
810 m_nhit_muc[ii] = mucTrk->numHits() ;
811 m_nlay_muc[ii] = mucTrk->numLayers() ;
812
813 }
814
815 if((*itTrk)->isTofTrackValid()) { //TOF information
816
817 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
818
819 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
820
821 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
822 TofHitStatus *status = new TofHitStatus;
823 status->setStatus((*iter_tof)->status());
824
825 if(!(status->is_barrel())){//endcap
826 if( (status->is_cluster()) ) m_t_etof[ii] = (*iter_tof)->tof();
827 if( !(status->is_counter()) ){if(status) delete status; continue; }// ?
828 if( status->layer()!=0 ) {if(status) delete status;continue;}//layer1
829 double path=(*iter_tof)->path(); // ?
830 double tof = (*iter_tof)->tof();
831 double ph = (*iter_tof)->ph();
832 double rhit = (*iter_tof)->zrhit();
833 double qual = 0.0 + (*iter_tof)->quality();
834 double cntr = 0.0 + (*iter_tof)->tofID();
835 double texp[5];
836 for(int j = 0; j < 5; j++) {
837 double gb = ptrk/xmass[j];
838 double beta = gb/sqrt(1+gb*gb);
839 texp[j] = path /beta/velc;
840 }
841
842 m_qual_etof[ii] = qual;
843 m_tof_etof[ii] = tof ;
844 }
845 else {//barrel
846 if( (status->is_cluster()) ) m_t_btof[ii] = (*iter_tof)->tof();
847 if( !(status->is_counter()) ){if(status) delete status; continue;} // ?
848 if(status->layer()==1){ //layer1
849 double path=(*iter_tof)->path(); // ?
850 double tof = (*iter_tof)->tof();
851 double ph = (*iter_tof)->ph();
852 double rhit = (*iter_tof)->zrhit();
853 double qual = 0.0 + (*iter_tof)->quality();
854 double cntr = 0.0 + (*iter_tof)->tofID();
855 double texp[5];
856 for(int j = 0; j < 5; j++) {
857 double gb = ptrk/xmass[j];
858 double beta = gb/sqrt(1+gb*gb);
859 texp[j] = path /beta/velc;
860 }
861
862 m_qual_btof1[ii] = qual;
863 m_tof_btof1[ii] = tof ;
864 }
865
866 if(status->layer()==2){//layer2
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
880 m_qual_btof2[ii] = qual;
881 m_tof_btof2[ii] = tof ;
882 }
883 }
884 if(status) delete status;
885 }
886 }
887
888 }
889 counter[4]++;
890
891 //std::cout<<"dbg_5"<<std::endl;
892 //tag
893
894 m_bhabhatag=0;
895
896 if(m_ngch != 2 || nCharge != 0 ) return StatusCode::SUCCESS;
897 EvtRecTrackIterator itTrk1;
898
899 EvtRecTrackIterator itTrk2;
900
901 RecMdcKalTrack *mdcKalTrk1;
902
903 RecMdcKalTrack *mdcKalTrk2;
904
905 HepLorentzVector p41e,p42e,p4le;
906 Hep3Vector p31e,p32e,p3le;
907 HepLorentzVector p41m,p42m,p4lm;
908 Hep3Vector p31m,p32m,p3lm;
909 HepLorentzVector p41h,p42h,p4lh;
910 Hep3Vector p31h,p32h,p3lh;
911 WTrackParameter w1_ini;
912 WTrackParameter w2_ini;
913 int iip=-1;
914 int iim=-1;
915 for(int i = 0; i < m_ngch; i++ ){
916 if(m_charge[i]>0)itTrk1= evtRecTrkCol->begin() + iGood[i];
917 if(m_charge[i]<0) itTrk2= evtRecTrkCol->begin() + iGood[i];
918 if(m_charge[i]>0) mdcKalTrk1 = (*itTrk1)->mdcKalTrack();
919 if(m_charge[i]<0) mdcKalTrk2 = (*itTrk2)->mdcKalTrack();
920 if(m_charge[i]>0)iip=i;
921 if(m_charge[i]<0)iim=i;
922
923
924 if(m_charge[i]>0) w1_ini=WTrackParameter (xmass[0],mdcKalTrk1->getZHelixE(),mdcKalTrk1->getZErrorE());
925 if(m_charge[i]<0) w2_ini=WTrackParameter (xmass[0],mdcKalTrk2 ->getZHelixE(),mdcKalTrk2 ->getZErrorE());
926 if(m_charge[i]>0) p41e =w1_ini.p();
927 if(m_charge[i]<0) p42e =w2_ini.p();
928 if(m_charge[i]>0) p41e.boost(u_cms);
929 if(m_charge[i]<0) p42e.boost(u_cms);
930 if(m_charge[i]>0) p31e = p41e.vect();
931 if(m_charge[i]<0) p32e = p42e.vect();
932
933
934
935
936 if(m_charge[i]>0){
937 m_px_cms_ep=p41e.px();
938 m_py_cms_ep=p41e.py();
939 m_pz_cms_ep=p41e.pz();
940 m_e_cms_ep=p41e.e();
941 }
942 if(m_charge[i]<0){
943 m_px_cms_em=p42e.px();
944 m_py_cms_em=p42e.py();
945 m_pz_cms_em=p42e.pz();
946 m_e_cms_em=p42e.e();
947 }
948
949
950 }
951
952 double e01=(iip!=-1)?m_e_emc[iip]:0;//m_e_cms_ep;
953 double e02=(iim!=-1)?m_e_emc[iim]:0;//m_e_cms_em;
954 int ilarge=( e01 > e02 ) ?iip:iim;
955 p4le=( e01 > e02 ) ?p41e:p42e;
956 p4lm=( e01 > e02 ) ?p41m:p42m;
957
958 p3le=( e01 > e02 ) ?p31e:p32e;
959 p3lm=( e01 > e02 ) ?p31m:p32m;
960
961 double acolle= 180.-p31e.angle(p32e)* 180.0 / CLHEP::pi;
962 double acople= 180.- (p31e.perpPart()).angle(p32e.perpPart ())* 180.0 / CLHEP::pi;
963 double poeb1e=p41e.rho()/beamEnergy;
964 double poeb2e=p42e.rho()/beamEnergy;
965 double poeble=p4le.rho()/beamEnergy ;
966
967
968 double eemc1=m_e_emc[iip];
969 double eemc2=m_e_emc[iim];
970
971
972 double ex1=m_kal_vx0[iip];
973 double ey1=m_kal_vy0[iip];
974 double ez1=m_kal_vz0[iip];
975 double epx1=m_kal_px[iip];
976 double epy1=m_kal_py[iip];
977 double epz1=m_kal_pz[iip];
978 double epp1=m_kal_p[iip];
979 double ex2=m_kal_vx0[iim];
980 double ey2=m_kal_vy0[iim];
981 double ez2=m_kal_vz0[iim];
982 double epx2=m_kal_px[iim];
983 double epy2=m_kal_py[iim];
984 double epz2=m_kal_pz[iim];
985 double epp2=m_kal_p[iim];
986
987 double pidchidedx1=m_pidchiDedx[iip];
988 double pidchitof11= m_pidchiTof1[iip];
989 double pidchitof21= m_pidchiTof2[iip];
990 double pidchidedx2=m_pidchiDedx[iim];
991 double pidchitof12= m_pidchiTof1[iim];
992 double pidchitof22= m_pidchiTof2[iim];
993
994 double eoeb1=m_e_emc[iip]/beamEnergy;
995 double eoeb2=m_e_emc[iim]/beamEnergy;
996
997 double eope1=0;
998 if(p41e.rho()>0)eope1=m_e_emc[iip]/p41e.rho();
999 double eope2=0;
1000 if(p42e.rho()>0)eope2=m_e_emc[iim]/p42e.rho();
1001
1002
1003
1004 double exoeb1= m_e_emc[iip]*sin(m_theta_emc[iip])*cos(m_phi_emc[iip])/beamEnergy;
1005 double eyoeb1= m_e_emc[iip]*sin(m_theta_emc[iip])*sin(m_phi_emc[iip])/beamEnergy;
1006 double ezoeb1=m_e_emc[iip]*cos(m_theta_emc[iip])/beamEnergy;
1007 double etoeb1=m_e_emc[iip]*sin(m_theta_emc[iip])/beamEnergy;
1008
1009 double exoeb2= m_e_emc[iim]*sin(m_theta_emc[iim])*cos(m_phi_emc[iim])/beamEnergy;
1010 double eyoeb2= m_e_emc[iim]*sin(m_theta_emc[iim])*sin(m_phi_emc[iim])/beamEnergy;
1011 double ezoeb2=m_e_emc[iim]*cos(m_theta_emc[iim])/beamEnergy;
1012 double etoeb2=m_e_emc[iim]*sin(m_theta_emc[iim])/beamEnergy;
1013
1014 double eoebl=m_e_emc[ilarge]/beamEnergy;
1015
1016 double eopl=0;
1017 if(p4le.rho()>0)eopl=m_e_emc[ilarge]/p4le.rho();
1018
1019 double exoebl= m_e_emc[ilarge]*sin(m_theta_emc[ilarge])*cos(m_phi_emc[ilarge])/beamEnergy;
1020 double eyoebl= m_e_emc[ilarge]*sin(m_theta_emc[ilarge])*sin(m_phi_emc[ilarge])/beamEnergy;
1021 double ezoebl=m_e_emc[ilarge]*cos(m_theta_emc[ilarge])/beamEnergy;
1022 double etoebl=m_e_emc[ilarge]*sin(m_theta_emc[ilarge])/beamEnergy;
1023
1024 int mucinfo1=(m_nhit_muc[iip]>=2&&m_nlay_muc[iip]>=2 ) ? 1 : 0;
1025 int mucinfo2=(m_nhit_muc[iim]>=2&&m_nlay_muc[iim]>=2) ? 1 : 0;
1026 int mucinfol=(m_nhit_muc[ilarge]>=2&&m_nlay_muc[ilarge]>=2) ? 1 : 0;
1027 int pidel=( e01 > e02 ) ? m_nep : m_nem;
1028 int pidmul=( e01 > e02 ) ? m_nmup : m_nmum;
1029 double deltatof=0;
1030
1031 if(m_t_btof[iip]*m_t_btof[iim]!=0) deltatof+=fabs(m_t_btof[iip]-m_t_btof[iim]);
1032 if(m_t_etof[iip]*m_t_etof[iim]!=0) deltatof+=fabs(m_t_etof[iip]-m_t_etof[iim]);
1033
1034
1035
1036
1037
1038
1039
1040 if (acolle<m_acoll_e_cut) m_bhabhatag+=1;
1041 if (acople<m_acopl_e_cut)m_bhabhatag+=10;
1042 if (fabs(m_ecms-mpsi2s)<0.001){
1043 if (sqrt((eoeb1-1)*(eoeb1-1)+(eoeb2-1)*(eoeb2-1))<m_tetotal_e_cut)m_bhabhatag+=100;
1044 }
1045 else{
1046 if ((fabs(poeb1e-1)<m_poeb_e_cut)&&(fabs(poeb2e-1)<m_poeb_e_cut))m_bhabhatag+=100;
1047 }
1048 if (!m_useTOF||(deltatof<m_dtof_e_cut))m_bhabhatag+=1000;
1049 if (poeb1e>m_tpoeb_e_cut||poeb2e>m_tpoeb_e_cut||(sqrt((poeb1e-1)*(poeb1e-1)+(poeb2e-1)*(poeb2e-1))<m_tptotal_e_cut))m_bhabhatag+=10000;
1050 if (!m_useMUC||(mucinfo1==0||mucinfo2==0))m_bhabhatag+=100000;
1051 if (!m_usePID||(m_nep==1&&m_nem==1))m_bhabhatag+=1000000;
1052
1053
1054
1055
1056
1057 m_acoll=acolle;
1058 m_acopl=acople;
1059 m_poeb1=poeb1e;
1060 m_poeb2=poeb2e;
1061 m_eop1=eope1;
1062 m_eop2=eope2;
1063 m_cos_ep=p41e.cosTheta ();
1064 m_cos_em=p42e.cosTheta ();
1065 m_mass_ee=(p41e+p42e).m();
1066 m_deltatof=deltatof;
1067
1068 m_eoeb1=eoeb1;
1069 m_eoeb2=eoeb2;
1070
1071 m_etoeb1=etoeb1;
1072 m_etoeb2=etoeb2;
1073 m_mucinfo1=mucinfo1;
1074 m_mucinfo2=mucinfo2;
1075
1076 if(m_bhabhatag==1111111){
1077 nbhabha++;
1078 if(m_writentuple)m_tuple1 -> write();
1079 ////////////////////////////////////////////////////////////
1080 // DQA
1081 // set tag and quality
1082
1083 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
1084
1085 (*itTrk1)->setPartId(1);
1086 (*itTrk2)->setPartId(1);
1087 // Quality: defined by whether dE/dx or TOF is used to identify particle
1088 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
1089 // 1 - only dE/dx (can be used for TOF calibration)
1090 // 2 - only TOF (can be used for dE/dx calibration)
1091 // 3 - Both dE/dx and TOF
1092 (*itTrk1)->setQuality(0);
1093 (*itTrk2)->setQuality(0);
1094
1095 // DQA
1096 // Add the line below at the end of execute(), (before return)
1097 //
1098 setFilterPassed(true);
1099 ////////////////////////////////////////////////////////////
1100 m_ee_mass->Fill((p41e+p42e).m());
1101 m_ee_acoll->Fill(acolle);
1102 m_ee_eop_ep->Fill(eope1);
1103 m_ee_eop_em->Fill(eope2);
1104 m_ee_costheta_ep->Fill(p41e.cosTheta ());
1105 m_ee_costheta_em->Fill(p42e.cosTheta ());
1106 m_ee_phi_ep->Fill(p41e.phi ());
1107 m_ee_phi_em->Fill(p42e.phi ());
1108 m_ee_nneu->Fill(m_nGam);
1109
1110
1111
1112 m_ee_eemc_ep->Fill(eemc1);
1113 m_ee_eemc_em->Fill(eemc2);
1114
1115 m_ee_x_ep->Fill(ex1);
1116 m_ee_y_ep->Fill(ey1);
1117 m_ee_z_ep->Fill(ez1);
1118
1119 m_ee_x_em->Fill(ex2);
1120 m_ee_y_em->Fill(ey2);
1121 m_ee_z_em->Fill(ez2);
1122
1123
1124 m_ee_px_ep->Fill(epx1);
1125 m_ee_py_ep->Fill(epy1);
1126 m_ee_pz_ep->Fill(epz1);
1127 m_ee_p_ep->Fill(epp1);
1128
1129 m_ee_px_em->Fill(epx2);
1130 m_ee_py_em->Fill(epy2);
1131 m_ee_pz_em->Fill(epz2);
1132 m_ee_p_em->Fill(epp2);
1133
1134 m_ee_deltatof->Fill(deltatof);
1135
1136 m_ee_pidchidedx_ep->Fill(pidchidedx1);
1137 m_ee_pidchidedx_em->Fill(pidchidedx2);
1138 m_ee_pidchitof1_ep->Fill(pidchitof11);
1139 m_ee_pidchitof1_em->Fill(pidchitof12);
1140 m_ee_pidchitof2_ep->Fill(pidchitof21);
1141 m_ee_pidchitof2_em->Fill(pidchitof22);
1142
1143
1144
1145
1146 }
1147
1148
1149
1150
1151
1152
1153
1154
1155 return StatusCode::SUCCESS;
1156
1157
1158
1159
1160
1161
1162}
1163
1164// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1165StatusCode DQASelBhabha::finalize() {
1166
1167 MsgStream log(msgSvc(), name());
1168 log << MSG::INFO << "in finalize()" << endmsg;
1169 return StatusCode::SUCCESS;
1170}
1171
1172
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:44
Vp4
std::vector< HepLorentzVector > Vp4
Definition: DQASelBhabha.cxx:53
xmass
const double xmass[5]
Definition: DQASelBhabha.cxx:50
mpsi2s
const double mpsi2s
Definition: DQASelBhabha.cxx:47
velc
const double velc
Definition: DQASelBhabha.cxx:51
mk
const double mk
Definition: DQASelBhabha.cxx:49
mpi
const double mpi
Definition: DQASelBhabha.cxx:48
Vint
std::vector< int > Vint
Definition: DQASelBhabha.cxx:52
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_nphot 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:1165
DQASelBhabha::initialize
StatusCode initialize()
Definition: DQASelBhabha.cxx:102
DQASelBhabha::DQASelBhabha
DQASelBhabha(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DQASelBhabha.cxx:59
DQASelBhabha::execute
StatusCode execute()
Definition: DQASelBhabha.cxx:370
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