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