BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ValidPhyJPsill.cxx
Go to the documentation of this file.
1#include "GaudiKernel/IDataProviderSvc.h"
2#include "GaudiKernel/ISvcLocator.h"
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/PropertyMgr.h"
5#include "GaudiKernel/SmartDataPtr.h"
6
7#include "EventModel/Event.h"
8#include "EventModel/EventModel.h"
9
10#include "DstEvent/TofHitStatus.h"
11#include "EventModel/EventHeader.h"
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14#include "McTruth/McParticle.h"
15
16#include "GaudiKernel/Bootstrap.h"
17#include "GaudiKernel/INTupleSvc.h"
18#include "GaudiKernel/NTuple.h"
19#include "TMath.h"
20// #include "GaudiKernel/IHistogramSvc.h"
21#include "CLHEP/Vector/LorentzVector.h"
22#include "CLHEP/Vector/ThreeVector.h"
23#include "CLHEP/Vector/TwoVector.h"
24using CLHEP::Hep2Vector;
25using CLHEP::Hep3Vector;
26using CLHEP::HepLorentzVector;
27#include "CLHEP/Geometry/Point3D.h"
28#ifndef ENABLE_BACKWARDS_COMPATIBILITY
29typedef HepGeom::Point3D<double> HepPoint3D;
30#endif
31#include "ParticleID/ParticleID.h"
32#include "ValidPhyJPsill.h"
33#include "VertexFit/KinematicFit.h"
34#include "VertexFit/VertexFit.h"
35
36#include <vector>
37// const double twopi = 6.2831853;
38// const double pi = 3.1415927;
39const double mpi = 0.13957;
40const double mks0 = 0.497648;
41const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
42// const double velc = 29.9792458; tof_path unit in cm.
43const double velc = 299.792458; // tof path unit in mm
44// const double beamAngle = 0.022; // the angle between the beam directions.
45// const double beamEnergy = 3.686/2.;
46typedef std::vector<int> Vint;
47typedef std::vector<HepLorentzVector> Vp4;
48// const HepLorentzVector p_cms(0.034067,0.0,0.0,3.097);
49// const Hep3Vector u_cms = -p_cms.boostVector();
50
51// declare one counter
52static int counter[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
53/////////////////////////////////////////////////////////////////////////////
54
55static int nbhabha = 0;
56static int ndimu = 0;
57static int nhadron = 0;
58static int n0prong = 0;
59static int n2prong = 0;
60static int n4prong = 0;
61static int ng4prong = 0;
62DECLARE_COMPONENT( ValidPhyJPsill )
63ValidPhyJPsill::ValidPhyJPsill( const std::string& name, ISvcLocator* pSvcLocator )
64 : Algorithm( name, pSvcLocator ) {
65
66 // Declare the properties
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 = 8.0 );
71 declareProperty( "Coscut", m_coscut = 0.93 );
72
73 declareProperty( "EnergyThreshold", m_energyThreshold = 0.01 );
74 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
75 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
76 declareProperty( "GammaTrkCut", m_gammaTrkCut = 20.0 );
77
78 declareProperty( "acoll_e_cut", m_acoll_e_cut = 2. );
79 declareProperty( "acopl_e_cut", m_acopl_e_cut = 2. );
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 declareProperty( "acoll_mu_cut", m_acoll_mu_cut = 2. );
89 declareProperty( "acopl_mu_cut", m_acopl_mu_cut = 2. );
90 declareProperty( "poeb_mu_cut", m_poeb_mu_cut = 0.3 );
91 declareProperty( "dtof_mu_cut", m_dtof_mu_cut = 4. );
92 declareProperty( "eoeb_mu_cut", m_eoeb_mu_cut = 0.35 );
93 declareProperty( "etotal_mu_cut", m_etotal_mu_cut = 0.6 );
94 declareProperty( "tpoebh_mu_cut", m_tpoebh_mu_cut = 0.9 );
95 declareProperty( "tpoebl_mu_cut", m_tpoebl_mu_cut = 0.7 );
96 declareProperty( "tptotal_mu_cut", m_tptotal_mu_cut = 1.0 );
97
98 declareProperty( "TagBhabha", m_tagBhabha = true );
99 declareProperty( "endcap", m_endcap = false );
100 declareProperty( "TagDimu", m_tagDimu = true );
101 declareProperty( "m_useTOF", m_useTOF = true );
102 declareProperty( "m_usePID", m_usePID = true );
103 declareProperty( "m_useEMC", m_useEMC = true );
104 declareProperty( "m_useMUC", m_useMUC = true );
105 declareProperty( "MCdump", m_mcdump = false );
106 declareProperty( "StudyMode", m_studymode = false );
107}
108
109// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
110// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
111StatusCode ValidPhyJPsill::initialize() {
112 MsgStream log( msgSvc(), name() );
113
114 log << MSG::INFO << "in initialize()" << endmsg;
115
116 StatusCode status;
117
118 status = service( "THistSvc", m_thistsvc );
119 if ( status.isFailure() )
120 {
121 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
122 return status;
123 }
124
125 if ( m_tagBhabha )
126 {
127
128 m_ee_mass = new TH1F( "ee_mass", "ee_mass", 80, 2.8, 3.2 );
129 status = m_thistsvc->regHist( "/VAL/PHY/ee_mass", m_ee_mass );
130 m_ee_acoll = new TH1F( "ee_acoll", "ee_acoll", 60, 0, 6 );
131 status = m_thistsvc->regHist( "/VAL/PHY/ee_acoll", m_ee_acoll );
132 m_ee_eop_ep = new TH1F( "ee_eop_ep", "ee_eop_ep", 100, 0.4, 1.4 );
133 status = m_thistsvc->regHist( "/VAL/PHY/ee_eop_ep", m_ee_eop_ep );
134 m_ee_eop_em = new TH1F( "ee_eop_em", "ee_eop_em", 100, 0.4, 1.4 );
135 status = m_thistsvc->regHist( "/VAL/PHY/ee_eop_em", m_ee_eop_em );
136 m_ee_costheta_ep = new TH1F( "ee_costheta_ep", "ee_costheta_ep", 100, -1, 1 );
137 status = m_thistsvc->regHist( "/VAL/PHY/ee_costheta_ep", m_ee_costheta_ep );
138 m_ee_costheta_em = new TH1F( "ee_costheta_em", "ee_costheta_em", 100, -1, 1 );
139 status = m_thistsvc->regHist( "/VAL/PHY/ee_costheta_em", m_ee_costheta_em );
140
141 m_ee_phi_ep = new TH1F( "ee_phi_ep", "ee_phi_ep", 120, -3.2, 3.2 );
142 status = m_thistsvc->regHist( "/VAL/PHY/ee_phi_ep", m_ee_phi_ep );
143 m_ee_phi_em = new TH1F( "ee_phi_em", "ee_phi_em", 120, -3.2, 3.2 );
144 status = m_thistsvc->regHist( "/VAL/PHY/ee_phi_em", m_ee_phi_em );
145
146 m_ee_nneu = new TH1F( "ee_nneu", "ee_nneu", 5, 0, 5 );
147 status = m_thistsvc->regHist( "/VAL/PHY/ee_nneu", m_ee_nneu );
148
149 m_ee_eemc_ep = new TH1F( "ee_eemc_ep", "ee_eemc_ep", 100, 1.0, 2.0 );
150 status = m_thistsvc->regHist( "/VAL/PHY/ee_eemc_ep", m_ee_eemc_ep );
151 m_ee_eemc_em = new TH1F( "ee_eemc_em", "ee_eemc_em", 100, 1.0, 2.0 );
152 status = m_thistsvc->regHist( "/VAL/PHY/ee_eemc_em", m_ee_eemc_em );
153 m_ee_x_ep = new TH1F( "ee_x_ep", "ee_x_ep", 100, -1.0, 1.0 );
154 status = m_thistsvc->regHist( "/VAL/PHY/ee_x_ep", m_ee_x_ep );
155 m_ee_y_ep = new TH1F( "ee_y_ep", "ee_y_ep", 100, -1.0, 1.0 );
156 status = m_thistsvc->regHist( "/VAL/PHY/ee_y_ep", m_ee_y_ep );
157 m_ee_z_ep = new TH1F( "ee_z_ep", "ee_z_ep", 100, -10.0, 10.0 );
158 status = m_thistsvc->regHist( "/VAL/PHY/ee_z_ep", m_ee_z_ep );
159 m_ee_x_em = new TH1F( "ee_x_em", "ee_x_em", 100, -1.0, 1.0 );
160 status = m_thistsvc->regHist( "/VAL/PHY/ee_x_em", m_ee_x_em );
161 m_ee_y_em = new TH1F( "ee_y_em", "ee_y_em", 100, -1.0, 1.0 );
162 status = m_thistsvc->regHist( "/VAL/PHY/ee_y_em", m_ee_y_em );
163 m_ee_z_em = new TH1F( "ee_z_em", "ee_z_em", 100, -10.0, 10.0 );
164 status = m_thistsvc->regHist( "/VAL/PHY/ee_z_em", m_ee_z_em );
165
166 m_ee_px_ep = new TH1F( "ee_px_ep", "ee_px_ep", 200, -2.0, 2.0 );
167 status = m_thistsvc->regHist( "/VAL/PHY/ee_px_ep", m_ee_px_ep );
168 m_ee_py_ep = new TH1F( "ee_py_ep", "ee_py_ep", 200, -2.0, 2.0 );
169 status = m_thistsvc->regHist( "/VAL/PHY/ee_py_ep", m_ee_py_ep );
170 m_ee_pz_ep = new TH1F( "ee_pz_ep", "ee_pz_ep", 200, -2.0, 2.0 );
171 status = m_thistsvc->regHist( "/VAL/PHY/ee_pz_ep", m_ee_pz_ep );
172 m_ee_p_ep = new TH1F( "ee_p_ep", "ee_p_ep", 100, 1.0, 2.0 );
173 status = m_thistsvc->regHist( "/VAL/PHY/ee_p_ep", m_ee_p_ep );
174 m_ee_px_em = new TH1F( "ee_px_em", "ee_px_em", 100, -2.0, 2.0 );
175 status = m_thistsvc->regHist( "/VAL/PHY/ee_px_em", m_ee_px_em );
176 m_ee_py_em = new TH1F( "ee_py_em", "ee_py_em", 100, -2.0, 2.0 );
177 status = m_thistsvc->regHist( "/VAL/PHY/ee_py_em", m_ee_py_em );
178 m_ee_pz_em = new TH1F( "ee_pz_em", "ee_pz_em", 100, -2.0, 2.0 );
179 status = m_thistsvc->regHist( "/VAL/PHY/ee_pz_em", m_ee_pz_em );
180 m_ee_p_em = new TH1F( "ee_p_em", "ee_p_em", 100, 1.0, 2.0 );
181 status = m_thistsvc->regHist( "/VAL/PHY/ee_p_em", m_ee_p_em );
182 m_ee_deltatof = new TH1F( "ee_deltatof", "ee_deltatof", 50, 0.0, 10.0 );
183 status = m_thistsvc->regHist( "/VAL/PHY/ee_deltatof", m_ee_deltatof );
184
185 m_ee_pidchidedx_ep = new TH1F( "ee_pidchidedx_ep", "ee_pidchidedx_ep", 160, -4, 4 );
186 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchidedx_ep", m_ee_pidchidedx_ep );
187 m_ee_pidchidedx_em = new TH1F( "ee_pidchidedx_em", "ee_pidchidedx_em", 160, -4, 4 );
188 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchidedx_em", m_ee_pidchidedx_em );
189 m_ee_pidchitof1_ep = new TH1F( "ee_pidchitof1_ep", "ee_pidchitof1_ep", 160, -4, 4 );
190 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchitof1_ep", m_ee_pidchitof1_ep );
191 m_ee_pidchitof1_em = new TH1F( "ee_pidchitof1_em", "ee_pidchitof1_em", 160, -4, 4 );
192 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchitof1_em", m_ee_pidchitof1_em );
193 m_ee_pidchitof2_ep = new TH1F( "ee_pidchitof2_ep", "ee_pidchitof2_ep", 160, -4, 4 );
194 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchitof2_ep", m_ee_pidchitof2_ep );
195 m_ee_pidchitof2_em = new TH1F( "ee_pidchitof2_em", "ee_pidchitof2_em", 160, -4, 4 );
196 status = m_thistsvc->regHist( "/VAL/PHY/ee_pidchitof2_em", m_ee_pidchitof2_em );
197 }
198
199 if ( m_tagDimu )
200 {
201
202 m_mumu_mass = new TH1F( "mumu_mass", "mumu_mass", 80, 2.8, 3.2 );
203 status = m_thistsvc->regHist( "/VAL/PHY/mumu_mass", m_mumu_mass );
204 m_mumu_acoll = new TH1F( "mumu_acoll", "mumu_acoll", 60, 0, 6 );
205 status = m_thistsvc->regHist( "/VAL/PHY/mumu_acoll", m_mumu_acoll );
206 m_mumu_eop_mup = new TH1F( "mumu_eop_mup", "mumu_eop_mup", 100, 0., 0.5 );
207 status = m_thistsvc->regHist( "/VAL/PHY/mumu_eop_mup", m_mumu_eop_mup );
208 m_mumu_eop_mum = new TH1F( "mumu_eop_mum", "mumu_eop_mum", 100, 0., 0.5 );
209 status = m_thistsvc->regHist( "/VAL/PHY/mumu_eop_mum", m_mumu_eop_mum );
210 m_mumu_costheta_mup = new TH1F( "mumu_costheta_mup", "mumu_costheta_mup", 100, -1, 1 );
211 status = m_thistsvc->regHist( "/VAL/PHY/mumu_costheta_mup", m_mumu_costheta_mup );
212 m_mumu_costheta_mum = new TH1F( "mumu_costheta_mum", "mumu_costheta_mum", 100, -1, 1 );
213 status = m_thistsvc->regHist( "/VAL/PHY/mumu_costheta_mum", m_mumu_costheta_mum );
214
215 m_mumu_phi_mup = new TH1F( "mumu_phi_mup", "mumu_phi_mup", 120, -3.2, 3.2 );
216 status = m_thistsvc->regHist( "/VAL/PHY/mumu_phi_mup", m_mumu_phi_mup );
217 m_mumu_phi_mum = new TH1F( "mumu_phi_mum", "mumu_phi_mum", 120, -3.2, 3.2 );
218 status = m_thistsvc->regHist( "/VAL/PHY/mumu_phi_mum", m_mumu_phi_mum );
219
220 m_mumu_nneu = new TH1F( "mumu_nneu", "mumu_nneu", 5, 0, 5 );
221 status = m_thistsvc->regHist( "/VAL/PHY/mumu_nneu", m_mumu_nneu );
222 m_mumu_nlay = new TH1F( "mumu_nlay", "mumu_nlay", 9, 0, 10 );
223 status = m_thistsvc->regHist( "/VAL/PHY/mumu_nlay", m_mumu_nlay );
224 m_mumu_nhit = new TH1F( "mumu_nhit", "mumu_nhit", 19, 0, 20 );
225 status = m_thistsvc->regHist( "/VAL/PHY/mumu_nhit", m_mumu_nhit );
226
227 m_mumu_eemc_mup = new TH1F( "mumu_eemc_mup", "mumu_eemc_mup", 100, 0.0, 1.0 );
228 status = m_thistsvc->regHist( "/VAL/PHY/mumu_eemc_mup", m_mumu_eemc_mup );
229 m_mumu_eemc_mum = new TH1F( "mumu_eemc_mum", "mumu_eemc_mum", 100, 0.0, 1.0 );
230 status = m_thistsvc->regHist( "/VAL/PHY/mumu_eemc_mum", m_mumu_eemc_mum );
231 m_mumu_x_mup = new TH1F( "mumu_x_mup", "mumu_x_mup", 100, -1.0, 1.0 );
232 status = m_thistsvc->regHist( "/VAL/PHY/mumu_x_mup", m_mumu_x_mup );
233 m_mumu_y_mup = new TH1F( "mumu_y_mup", "mumu_y_mup", 100, -1.0, 1.0 );
234 status = m_thistsvc->regHist( "/VAL/PHY/mumu_y_mup", m_mumu_y_mup );
235 m_mumu_z_mup = new TH1F( "mumu_z_mup", "mumu_z_mup", 100, -10.0, 10.0 );
236 status = m_thistsvc->regHist( "/VAL/PHY/mumu_z_mup", m_mumu_z_mup );
237 m_mumu_x_mum = new TH1F( "mumu_x_mum", "mumu_x_mum", 100, -1.0, 1.0 );
238 status = m_thistsvc->regHist( "/VAL/PHY/mumu_x_mum", m_mumu_x_mum );
239 m_mumu_y_mum = new TH1F( "mumu_y_mum", "mumu_y_mum", 100, -1.0, 1.0 );
240 status = m_thistsvc->regHist( "/VAL/PHY/mumu_y_mum", m_mumu_y_mum );
241 m_mumu_z_mum = new TH1F( "mumu_z_mum", "mumu_z_mum", 100, -10.0, 10.0 );
242 status = m_thistsvc->regHist( "/VAL/PHY/mumu_z_mum", m_mumu_z_mum );
243
244 m_mumu_px_mup = new TH1F( "mumu_px_mup", "mumu_px_mup", 200, -2.0, 2.0 );
245 status = m_thistsvc->regHist( "/VAL/PHY/mumu_px_mup", m_mumu_px_mup );
246 m_mumu_py_mup = new TH1F( "mumu_py_mup", "mumu_py_mup", 200, -2.0, 2.0 );
247 status = m_thistsvc->regHist( "/VAL/PHY/mumu_py_mup", m_mumu_py_mup );
248 m_mumu_pz_mup = new TH1F( "mumu_pz_mup", "mumu_pz_mup", 200, -2.0, 2.0 );
249 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pz_mup", m_mumu_pz_mup );
250 m_mumu_p_mup = new TH1F( "mumu_p_mup", "mumu_p_mup", 100, 1.0, 2.0 );
251 status = m_thistsvc->regHist( "/VAL/PHY/mumu_p_mup", m_mumu_p_mup );
252 m_mumu_px_mum = new TH1F( "mumu_px_mum", "mumu_px_mum", 100, -2.0, 2.0 );
253 status = m_thistsvc->regHist( "/VAL/PHY/mumu_px_mum", m_mumu_px_mum );
254 m_mumu_py_mum = new TH1F( "mumu_py_mum", "mumu_py_mum", 100, -2.0, 2.0 );
255 status = m_thistsvc->regHist( "/VAL/PHY/mumu_py_mum", m_mumu_py_mum );
256 m_mumu_pz_mum = new TH1F( "mumu_pz_mum", "mumu_pz_mum", 100, -2.0, 2.0 );
257 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pz_mum", m_mumu_pz_mum );
258 m_mumu_p_mum = new TH1F( "mumu_p_mum", "mumu_p_mum", 100, 1.0, 2.0 );
259 status = m_thistsvc->regHist( "/VAL/PHY/mumu_p_mum", m_mumu_p_mum );
260 m_mumu_deltatof = new TH1F( "mumu_deltatof", "mumu_deltatof", 50, 0.0, 10.0 );
261 status = m_thistsvc->regHist( "/VAL/PHY/mumu_deltatof", m_mumu_deltatof );
262
263 m_mumu_pidchidedx_mup =
264 new TH1F( "mumu_pidchidedx_mup", "mumu_pidchidedx_mup", 160, -4, 4 );
265 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchidedx_mup", m_mumu_pidchidedx_mup );
266 m_mumu_pidchidedx_mum =
267 new TH1F( "mumu_pidchidedx_mum", "mumu_pidchidedx_mum", 160, -4, 4 );
268 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchidedx_mum", m_mumu_pidchidedx_mum );
269 m_mumu_pidchitof1_mup =
270 new TH1F( "mumu_pidchitof1_mup", "mumu_pidchitof1_mup", 160, -4, 4 );
271 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchitof1_mup", m_mumu_pidchitof1_mup );
272 m_mumu_pidchitof1_mum =
273 new TH1F( "mumu_pidchitof1_mum", "mumu_pidchitof1_mum", 160, -4, 4 );
274 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchitof1_mum", m_mumu_pidchitof1_mum );
275 m_mumu_pidchitof2_mup =
276 new TH1F( "mumu_pidchitof2_mup", "mumu_pidchitof2_mup", 160, -4, 4 );
277 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchitof2_mup", m_mumu_pidchitof2_mup );
278 m_mumu_pidchitof2_mum =
279 new TH1F( "mumu_pidchitof2_mum", "mumu_pidchitof2_mum", 160, -4, 4 );
280 status = m_thistsvc->regHist( "/VAL/PHY/mumu_pidchitof2_mum", m_mumu_pidchitof2_mum );
281 }
282
283 if ( m_mcdump && ( m_tagBhabha || m_tagDimu ) )
284 {
285 NTuplePtr nt001( ntupleSvc(), "FILE1/mc" );
286 if ( nt001 ) m_tuple001 = nt001;
287 else
288 {
289 m_tuple001 =
290 ntupleSvc()->book( "FILE1/mc", CLID_ColumnWiseTuple, "Bhabha N-Tuple example" );
291 if ( m_tuple001 )
292 {
293 status = m_tuple001->addItem( "mc_ep_e", m_mc_ep_e );
294 status = m_tuple001->addItem( "mc_ep_px", m_mc_ep_px );
295 status = m_tuple001->addItem( "mc_ep_py", m_mc_ep_py );
296 status = m_tuple001->addItem( "mc_ep_pz", m_mc_ep_pz );
297 status = m_tuple001->addItem( "mc_ep_costheta", m_mc_ep_costheta );
298 status = m_tuple001->addItem( "mc_em_e", m_mc_em_e );
299 status = m_tuple001->addItem( "mc_em_px", m_mc_em_px );
300 status = m_tuple001->addItem( "mc_em_py", m_mc_em_py );
301 status = m_tuple001->addItem( "mc_em_pz", m_mc_em_pz );
302 status = m_tuple001->addItem( "mc_em_costheta", m_mc_em_costheta );
303 }
304 else
305
306 {
307 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple001 ) << endmsg;
308 return StatusCode::FAILURE;
309 }
310 }
311 }
312
313 NTuplePtr nt1( ntupleSvc(), "FILE1/signal" );
314 if ( nt1 ) m_tuple1 = nt1;
315 else
316 {
317 m_tuple1 = ntupleSvc()->book( "FILE1/signal", CLID_ColumnWiseTuple, "N-Tuple" );
318 if ( m_tuple1 )
319 {
320 status = m_tuple1->addItem( "irun", m_run );
321 status = m_tuple1->addItem( "ievent", m_event );
322 status = m_tuple1->addItem( "Nchrg", m_nchrg );
323 status = m_tuple1->addItem( "Nneu", m_nneu, 0, 40 );
324 status = m_tuple1->addItem( "NGch", m_ngch, 0, 40 );
325 status = m_tuple1->addItem( "NGam", m_nGam );
326
327 status = m_tuple1->addItem( "bhabhatag", m_bhabhatag );
328 status = m_tuple1->addItem( "dimutag", m_dimutag );
329 status = m_tuple1->addItem( "hadrontag", m_hadrontag );
330 status = m_tuple1->addItem( "sbhabhatag", m_sbhabhatag );
331 status = m_tuple1->addItem( "sdimutag", m_sdimutag );
332
333 status = m_tuple1->addItem( "acoll", m_acoll );
334 status = m_tuple1->addItem( "acopl", m_acopl );
335 status = m_tuple1->addItem( "deltatof", m_deltatof );
336 status = m_tuple1->addItem( "eop1", m_eop1 );
337 status = m_tuple1->addItem( "eop2", m_eop2 );
338 status = m_tuple1->addItem( "eoeb1", m_eoeb1 );
339 status = m_tuple1->addItem( "eoeb2", m_eoeb2 );
340 status = m_tuple1->addItem( "poeb1", m_poeb1 );
341 status = m_tuple1->addItem( "poeb2", m_poeb2 );
342 status = m_tuple1->addItem( "etoeb1", m_etoeb1 );
343 status = m_tuple1->addItem( "etoeb2", m_etoeb2 );
344 status = m_tuple1->addItem( "mucinfo1", m_mucinfo1 );
345 status = m_tuple1->addItem( "mucinfo2", m_mucinfo2 );
346
347 status = m_tuple1->addIndexedItem( "delang", m_nneu, m_delang );
348 status = m_tuple1->addIndexedItem( "delphi", m_nneu, m_delphi );
349 status = m_tuple1->addIndexedItem( "delthe", m_nneu, m_delthe );
350 status = m_tuple1->addIndexedItem( "nemchits", m_nneu, m_nemchits );
351 status = m_tuple1->addIndexedItem( "x", m_nneu, m_x );
352 status = m_tuple1->addIndexedItem( "y", m_nneu, m_y );
353 status = m_tuple1->addIndexedItem( "z", m_nneu, m_z );
354 status = m_tuple1->addIndexedItem( "theta", m_nneu, m_theta );
355 status = m_tuple1->addIndexedItem( "phi", m_nneu, m_phi );
356 status = m_tuple1->addIndexedItem( "dx", m_nneu, m_dx );
357 status = m_tuple1->addIndexedItem( "dy", m_nneu, m_dy );
358 status = m_tuple1->addIndexedItem( "dz", m_nneu, m_dz );
359 status = m_tuple1->addIndexedItem( "dtheta", m_nneu, m_dtheta );
360 status = m_tuple1->addIndexedItem( "dphi", m_nneu, m_dphi );
361 status = m_tuple1->addIndexedItem( "energy", m_nneu, m_energy );
362 status = m_tuple1->addIndexedItem( "dE", m_nneu, m_dE );
363 status = m_tuple1->addIndexedItem( "eSeed", m_nneu, m_eSeed );
364 status = m_tuple1->addIndexedItem( "e3x3", m_nneu, m_e3x3 );
365 status = m_tuple1->addIndexedItem( "e5x5", m_nneu, m_e5x5 );
366 status = m_tuple1->addIndexedItem( "secondMoment", m_nneu, m_secondMoment );
367 status = m_tuple1->addIndexedItem( "latMoment", m_nneu, m_latMoment );
368 status = m_tuple1->addIndexedItem( "a20Moment", m_nneu, m_a20Moment );
369 status = m_tuple1->addIndexedItem( "a42Moment", m_nneu, m_a42Moment );
370 status = m_tuple1->addIndexedItem( "getTime", m_nneu, m_getTime );
371 status = m_tuple1->addIndexedItem( "getEAll", m_nneu, m_getEAll );
372 // status = m_tuple1->addIndexedItem ("getELepton",m_nneu, m_getELepton);
373 // status = m_tuple1->addIndexedItem ("getETof2x1",m_nneu, m_getETof2x1);
374 // status = m_tuple1->addIndexedItem ("getETof2x3",m_nneu, m_getETof2x3);
375
376 // status = m_tuple1->addIndexedItem ("ThetaGap",m_nneu, m_ThetaGap);
377 // status = m_tuple1->addIndexedItem ("PhiGap",m_nneu, m_PhiGap);
378 status = m_tuple1->addIndexedItem( "charge", m_ngch, m_charge );
379 status = m_tuple1->addIndexedItem( "vx", m_ngch, m_vx0 );
380 status = m_tuple1->addIndexedItem( "vy", m_ngch, m_vy0 );
381 status = m_tuple1->addIndexedItem( "vz", m_ngch, m_vz0 );
382
383 status = m_tuple1->addIndexedItem( "px", m_ngch, m_px );
384 status = m_tuple1->addIndexedItem( "py", m_ngch, m_py );
385 status = m_tuple1->addIndexedItem( "pz", m_ngch, m_pz );
386 status = m_tuple1->addIndexedItem( "p", m_ngch, m_p );
387
388 status = m_tuple1->addIndexedItem( "kal_vx", m_ngch, m_kal_vx0 );
389 status = m_tuple1->addIndexedItem( "kal_vy", m_ngch, m_kal_vy0 );
390 status = m_tuple1->addIndexedItem( "kal_vz", m_ngch, m_kal_vz0 );
391
392 status = m_tuple1->addIndexedItem( "kal_px", m_ngch, m_kal_px );
393 status = m_tuple1->addIndexedItem( "kal_py", m_ngch, m_kal_py );
394 status = m_tuple1->addIndexedItem( "kal_pz", m_ngch, m_kal_pz );
395 status = m_tuple1->addIndexedItem( "kal_p", m_ngch, m_kal_p );
396
397 status = m_tuple1->addIndexedItem( "probPH", m_ngch, m_probPH );
398 status = m_tuple1->addIndexedItem( "normPH", m_ngch, m_normPH );
399 status = m_tuple1->addIndexedItem( "chie", m_ngch, m_chie );
400 status = m_tuple1->addIndexedItem( "chimu", m_ngch, m_chimu );
401 status = m_tuple1->addIndexedItem( "chipi", m_ngch, m_chipi );
402 status = m_tuple1->addIndexedItem( "chik", m_ngch, m_chik );
403 status = m_tuple1->addIndexedItem( "chip", m_ngch, m_chip );
404 status = m_tuple1->addIndexedItem( "ghit", m_ngch, m_ghit );
405 status = m_tuple1->addIndexedItem( "thit", m_ngch, m_thit );
406
407 status = m_tuple1->addIndexedItem( "e_emc", m_ngch, m_e_emc );
408 status = m_tuple1->addIndexedItem( "phi_emc", m_ngch, m_phi_emc );
409 status = m_tuple1->addIndexedItem( "theta_emc", m_ngch, m_theta_emc );
410
411 status = m_tuple1->addIndexedItem( "nhit_muc", m_ngch, m_nhit_muc );
412 status = m_tuple1->addIndexedItem( "nlay_muc", m_ngch, m_nlay_muc );
413 status = m_tuple1->addIndexedItem( "t_btof", m_ngch, m_t_btof );
414 status = m_tuple1->addIndexedItem( "t_etof", m_ngch, m_t_etof );
415 status = m_tuple1->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
416 status = m_tuple1->addIndexedItem( "tof_etof", m_ngch, m_tof_etof );
417 status = m_tuple1->addIndexedItem( "te_etof", m_ngch, m_te_etof );
418 status = m_tuple1->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
419 status = m_tuple1->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
420 status = m_tuple1->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
421 status = m_tuple1->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
422
423 status = m_tuple1->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
424 status = m_tuple1->addIndexedItem( "tof_btof1", m_ngch, m_tof_btof1 );
425 status = m_tuple1->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
426 status = m_tuple1->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
427 status = m_tuple1->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
428 status = m_tuple1->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
429 status = m_tuple1->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
430
431 status = m_tuple1->addIndexedItem( "qual_btof2", m_ngch, m_qual_btof2 );
432 status = m_tuple1->addIndexedItem( "tof_btof2", m_ngch, m_tof_btof2 );
433 status = m_tuple1->addIndexedItem( "te_btof2", m_ngch, m_te_btof2 );
434 status = m_tuple1->addIndexedItem( "tmu_btof2", m_ngch, m_tmu_btof2 );
435 status = m_tuple1->addIndexedItem( "tpi_btof2", m_ngch, m_tpi_btof2 );
436 status = m_tuple1->addIndexedItem( "tk_btof2", m_ngch, m_tk_btof2 );
437 status = m_tuple1->addIndexedItem( "tp_btof2", m_ngch, m_tp_btof2 );
438 status = m_tuple1->addIndexedItem( "pidcode", m_ngch, m_pidcode );
439 status = m_tuple1->addIndexedItem( "pidprob", m_ngch, m_pidprob );
440 status = m_tuple1->addIndexedItem( "pidchiDedx", m_ngch, m_pidchiDedx );
441 status = m_tuple1->addIndexedItem( "pidchiTof1", m_ngch, m_pidchiTof1 );
442 status = m_tuple1->addIndexedItem( "pidchiTof2", m_ngch, m_pidchiTof2 );
443
444 status = m_tuple1->addItem( "px_cms_ep", m_px_cms_ep ); // momentum of electron+
445 status = m_tuple1->addItem( "py_cms_ep", m_py_cms_ep ); // momentum of electron+
446 status = m_tuple1->addItem( "pz_cms_ep", m_pz_cms_ep ); // momentum of electron+
447 status = m_tuple1->addItem( "e_cms_ep", m_e_cms_ep ); // momentum of electron+
448 status = m_tuple1->addItem( "cos_ep", m_cos_ep ); // momentum of electron+
449 status = m_tuple1->addItem( "px_cms_em", m_px_cms_em ); // momentum of electron-
450 status = m_tuple1->addItem( "py_cms_em", m_py_cms_em ); // momentum of electron-
451 status = m_tuple1->addItem( "pz_cms_em", m_pz_cms_em ); // momentum of electron-
452 status = m_tuple1->addItem( "e_cms_em", m_e_cms_em ); // momentum of electron-
453 status = m_tuple1->addItem( "cos_em", m_cos_em ); // momentum of electron-
454 status = m_tuple1->addItem( "mass_ee", m_mass_ee ); //
455 status = m_tuple1->addItem( "emax", m_emax ); //
456 status = m_tuple1->addItem( "esum", m_esum ); //
457 status = m_tuple1->addItem( "npip", m_npip );
458 status = m_tuple1->addItem( "npim", m_npim );
459 status = m_tuple1->addItem( "nkp", m_nkp );
460 status = m_tuple1->addItem( "nkm", m_nkm );
461 status = m_tuple1->addItem( "np", m_np );
462 status = m_tuple1->addItem( "npb", m_npb );
463
464 status = m_tuple1->addItem( "nep", m_nep );
465 status = m_tuple1->addItem( "nem", m_nem );
466 status = m_tuple1->addItem( "nmup", m_nmup );
467 status = m_tuple1->addItem( "nmum", m_nmum );
468 }
469 else
470 {
471 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
472 return StatusCode::FAILURE;
473 }
474 }
475
476 //
477 //--------end of book--------
478 //
479
480 log << MSG::INFO << "successfully return from initialize()" << endmsg;
481 return StatusCode::SUCCESS;
482}
483
484// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
485// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
486StatusCode ValidPhyJPsill::execute() {
487 const double beamEnergy = m_ecms / 2.;
488 const HepLorentzVector p_cms( m_ecms * sin( m_beamangle * 0.5 ), 0.0, 0.0, m_ecms );
489 const Hep3Vector u_cms = -p_cms.boostVector();
490 MsgStream log( msgSvc(), name() );
491 log << MSG::INFO << "in execute()" << endmsg;
492
493 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
494 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
495 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
496
497 m_run = eventHeader->runNumber();
498 m_event = eventHeader->eventNumber();
499 m_nchrg = evtRecEvent->totalCharged();
500 m_nneu = evtRecEvent->totalNeutral();
501
502 log << MSG::INFO << "get event tag OK" << endmsg;
503
504 if ( m_mcdump && ( m_tagBhabha || m_tagDimu ) )
505 {
506 SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
507 if ( !mcParticleCol )
508 {
509 log << MSG::ERROR << "Can not find mcParticleCol" << endmsg;
510 return StatusCode::FAILURE;
511 }
512 // std::cout<<"here"<<std::endl;
513 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
514 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
515 {
516 if ( ( *iter_mc )->particleProperty() == -11 || ( *iter_mc )->particleProperty() == -13 )
517 {
518
519 m_mc_ep_px = ( *iter_mc )->initialFourMomentum().px() * 0.001;
520 m_mc_ep_py = ( *iter_mc )->initialFourMomentum().py() * 0.001;
521 m_mc_ep_pz = ( *iter_mc )->initialFourMomentum().pz() * 0.001;
522 m_mc_ep_e = ( *iter_mc )->initialFourMomentum().e() * 0.001;
523 m_mc_ep_costheta = cos( ( *iter_mc )->initialFourMomentum().theta() );
524 // std::cout<<"costheta= "<<cos((*iter_mc)->initialFourMomentum().theta())<<std::endl;
525 }
526 if ( ( *iter_mc )->particleProperty() == 11 || ( *iter_mc )->particleProperty() == 13 )
527 {
528
529 m_mc_em_px = ( *iter_mc )->initialFourMomentum().px() * 0.001;
530 m_mc_em_py = ( *iter_mc )->initialFourMomentum().py() * 0.001;
531 m_mc_em_pz = ( *iter_mc )->initialFourMomentum().pz() * 0.001;
532 m_mc_em_e = ( *iter_mc )->initialFourMomentum().e() * 0.001;
533 m_mc_em_costheta = cos( ( *iter_mc )->initialFourMomentum().theta() );
534 }
535 }
536 m_tuple001->write();
537 }
538
539 // cout <<"ncharg, nneu, tottks = "
540 // << evtRecEvent->totalCharged() << " , "
541 // << evtRecEvent->totalNeutral() << " , "
542 // << evtRecEvent->totalTracks() << endl ;
543 // std::cout<<"dbg_1"<<std::endl;
544 counter[0]++;
545
546 //
547 // Good charged track selection
548 //
549 Vint iGood;
550 iGood.clear();
551
552 int nCharge = 0;
553 int nneu = evtRecEvent->totalNeutral();
554
555 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
556 {
557 if ( i >= evtRecTrkCol->size() ) break;
558 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
559 if ( !( *itTrk )->isMdcTrackValid() ) continue;
560 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
561 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
562 double vz0 = mdcTrk->z();
563 double vr0 = mdcTrk->r();
564 if ( !( *itTrk )->isEmcShowerValid() ) continue;
565
566 if ( fabs( vz0 ) >= m_vz0cut ) continue;
567 if ( vr0 >= m_vr0cut ) continue;
568 double cost = cos( mdcTrk->theta() );
569 if ( fabs( cost ) >= m_coscut ) continue;
570
571 iGood.push_back( ( *itTrk )->trackId() );
572 nCharge += mdcTrk->charge();
573 }
574
575 //
576 // Finish Good Charged Track Selection
577 //
578 m_ngch = iGood.size();
579 if ( m_ngch == 0 ) n0prong++;
580 if ( m_ngch < 0 || m_ngch > 2 || ( abs( nCharge ) > 1 ) || m_nneu > 40 )
581 { return StatusCode::SUCCESS; }
582 // if((m_ngch != 2) || (nCharge != 0) ) { return StatusCode::SUCCESS; }
583 // std::cout << "ngood, totcharge = " << m_ngch << " , " << nCharge << std::endl;
584 // std::cout<<"dbg_2"<<std::endl;
585 counter[1]++;
586
587 if ( m_ngch == 2 && nCharge == 0 ) n2prong++;
588 if ( m_ngch == 4 && nCharge == 0 ) n4prong++;
589 if ( m_ngch > 4 ) ng4prong++;
590 //
591 // Particle ID
592 //
593 Vint ipip, ipim, ikp, ikm, ipp, ipm, iep, iem, imup, imum;
594 ipip.clear();
595 ipim.clear();
596 ikp.clear();
597 ikm.clear();
598 ipp.clear();
599 ipm.clear();
600 iep.clear();
601 iem.clear();
602 imup.clear();
603 imum.clear();
604 Vp4 p_pip, p_pim, p_kp, p_km, p_pp, p_pm;
605 p_pip.clear();
606 p_pim.clear();
607 p_kp.clear();
608 p_km.clear();
609 p_pp.clear();
610 p_pm.clear();
611
612 ParticleID* pid = ParticleID::instance();
613 for ( int i = 0; i < m_ngch; i++ )
614 {
615 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
616 // if(pid) delete pid;
617 pid->init();
618 pid->setMethod( pid->methodProbability() );
619 pid->setChiMinCut( 4 );
620 pid->setRecTrack( *itTrk );
621 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
622 // pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTofE() |
623 // pid->useTofQ()); // use PID sub-system
624 pid->identify( pid->onlyElectron() | pid->onlyMuon() | pid->onlyPion() | pid->onlyKaon() |
625 pid->onlyProton() ); // seperater Pion/Kaon/Proton
626 // pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
627 // pid->identify(pid->onlyPion());
628 // pid->identify(pid->onlyKaon());
629 pid->calculate();
630 if ( !( pid->IsPidInfoValid() ) ) continue;
631
632 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
633 RecMdcKalTrack* mdcKalTrk = nullptr;
634 if ( ( *itTrk )->isMdcKalTrackValid() ) mdcKalTrk = ( *itTrk )->mdcKalTrack();
635
636 double prob_pi = pid->probPion();
637 double prob_K = pid->probKaon();
638 double prob_p = pid->probProton();
639 double prob_e = pid->probElectron();
640 double prob_mu = pid->probMuon();
641 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
642 HepLorentzVector ptrk;
643 ptrk.setPx( mdcTrk->px() );
644 ptrk.setPy( mdcTrk->py() );
645 ptrk.setPz( mdcTrk->pz() );
646 double p3 = ptrk.mag();
647
648 if ( prob_pi > prob_K && prob_pi > prob_p )
649 {
650 m_pidcode[i] = 2;
651 m_pidprob[i] = pid->prob( 2 );
652 m_pidchiDedx[i] = pid->chiDedx( 2 );
653 m_pidchiTof1[i] = pid->chiTof1( 2 );
654 m_pidchiTof2[i] = pid->chiTof2( 2 );
655
656 ptrk.setE( sqrt( p3 * p3 + xmass[2] * xmass[2] ) );
657 if ( mdcTrk->charge() > 0 )
658 {
659 ipip.push_back( iGood[i] );
660 p_pip.push_back( ptrk );
661 }
662 if ( mdcTrk->charge() < 0 )
663 {
664 ipim.push_back( iGood[i] );
665 p_pim.push_back( ptrk );
666 }
667 }
668
669 if ( prob_K > prob_pi && prob_K > prob_p )
670 {
671 m_pidcode[i] = 3;
672 m_pidprob[i] = pid->prob( 3 );
673 m_pidchiDedx[i] = pid->chiDedx( 3 );
674 m_pidchiTof1[i] = pid->chiTof1( 3 );
675 m_pidchiTof2[i] = pid->chiTof2( 3 );
676 ptrk.setE( sqrt( p3 * p3 + xmass[3] * xmass[3] ) );
677 if ( mdcTrk->charge() > 0 )
678 {
679 ikp.push_back( iGood[i] );
680 p_kp.push_back( ptrk );
681 }
682 if ( mdcTrk->charge() < 0 )
683 {
684 ikm.push_back( iGood[i] );
685 p_km.push_back( ptrk );
686 }
687 }
688
689 if ( prob_p > prob_pi && prob_p > prob_K )
690 {
691 m_pidcode[i] = 4;
692 m_pidprob[i] = pid->prob( 4 );
693 m_pidchiDedx[i] = pid->chiDedx( 4 );
694 m_pidchiTof1[i] = pid->chiTof1( 4 );
695 m_pidchiTof2[i] = pid->chiTof2( 4 );
696 ptrk.setE( sqrt( p3 * p3 + xmass[4] * xmass[4] ) );
697 if ( mdcTrk->charge() > 0 )
698 {
699 ipp.push_back( iGood[i] );
700 p_pp.push_back( ptrk );
701 }
702 if ( mdcTrk->charge() < 0 )
703 {
704 ipm.push_back( iGood[i] );
705 p_pm.push_back( ptrk );
706 }
707 }
708 // if (prob_e > prob_mu) {
709 if ( m_tagBhabha )
710 {
711 m_pidcode[i] = 0;
712 m_pidprob[i] = pid->prob( 0 );
713 m_pidchiDedx[i] = pid->chiDedx( 0 );
714 m_pidchiTof1[i] = pid->chiTof1( 0 );
715 m_pidchiTof2[i] = pid->chiTof2( 0 );
716 if ( mdcTrk->charge() > 0 ) { iep.push_back( iGood[i] ); }
717 if ( mdcTrk->charge() < 0 ) { iem.push_back( iGood[i] ); }
718 }
719 // if (prob_mu > prob_e) {
720 if ( m_tagDimu )
721 {
722 m_pidcode[i] = 1;
723 m_pidprob[i] = pid->prob( 1 );
724 m_pidchiDedx[i] = pid->chiDedx( 1 );
725 m_pidchiTof1[i] = pid->chiTof1( 1 );
726 m_pidchiTof2[i] = pid->chiTof2( 1 );
727 if ( mdcTrk->charge() > 0 ) { imup.push_back( iGood[i] ); }
728 if ( mdcTrk->charge() < 0 ) { imum.push_back( iGood[i] ); }
729 }
730 }
731 m_npip = ipip.size();
732 m_npim = ipim.size();
733 m_nkp = ikp.size();
734 m_nkm = ikm.size();
735 m_np = ipp.size();
736 m_npb = ipm.size();
737 m_nep = iep.size();
738 m_nem = iem.size();
739 m_nmup = imup.size();
740 m_nmum = imum.size();
741 counter[2]++;
742 // std::cout<<"dbg_3"<<std::endl;
743
744 //
745 // Good neutral track selection
746 //
747 Vint iGam;
748 iGam.clear();
749 // ineu=0
750 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
751 {
752 if ( i >= evtRecTrkCol->size() ) break;
753 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
754 if ( !( *itTrk )->isEmcShowerValid() ) continue;
755 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
756 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
757
758 int n = emcTrk->numHits();
759 double x = emcTrk->x();
760 double y = emcTrk->y();
761 double z = emcTrk->z();
762 double dx = emcTrk->dx();
763 double dy = emcTrk->dy();
764 double dth = emcTrk->dtheta();
765 double dph = emcTrk->dphi();
766 double dz = emcTrk->dz();
767 double energy = emcTrk->energy();
768 double dE = emcTrk->dE();
769 double eSeed = emcTrk->eSeed();
770 double e3x3 = emcTrk->e3x3();
771 double e5x5 = emcTrk->e5x5();
772 double secondMoment = emcTrk->secondMoment();
773 double latMoment = emcTrk->latMoment();
774 double getTime = emcTrk->time();
775 double getEAll = emcTrk->getEAll();
776 double a20Moment = emcTrk->a20Moment();
777 double a42Moment = emcTrk->a42Moment();
778 // int phigap=emcTrk->PhiGap();
779 // int thetagap=emcTrk->ThetaGap();
780 // double getETof2x1 = emcTrk->getETof2x1();
781 // double getETof2x3 = emcTrk->getETof2x3();
782 // double getELepton = emcTrk->getELepton();
783
784 HepPoint3D EmcPos( x, y, z );
785 m_nemchits[i - m_nchrg] = n;
786 m_theta[i - m_nchrg] = EmcPos.theta();
787 m_phi[i - m_nchrg] = EmcPos.phi();
788 m_x[i - m_nchrg] = x;
789 m_y[i - m_nchrg] = y;
790 m_z[i - m_nchrg] = z;
791 m_dx[i - m_nchrg] = dx;
792 m_dy[i - m_nchrg] = dy;
793 m_dz[i - m_nchrg] = dz;
794 m_dtheta[i - m_nchrg] = dth;
795 m_dphi[i - m_nchrg] = dph;
796 m_energy[i - m_nchrg] = energy;
797 m_dE[i - m_nchrg] = dE;
798 m_eSeed[i - m_nchrg] = eSeed;
799 m_e3x3[i - m_nchrg] = e3x3;
800 m_e5x5[i - m_nchrg] = e5x5;
801 m_secondMoment[i - m_nchrg] = secondMoment;
802 m_latMoment[i - m_nchrg] = latMoment;
803 m_getTime[i - m_nchrg] = getTime;
804 m_getEAll[i - m_nchrg] = getEAll;
805 m_a20Moment[i - m_nchrg] = a20Moment;
806 m_a42Moment[i - m_nchrg] = a42Moment;
807
808 // m_getELepton[i]=getELepton;
809 // m_getETof2x1[i]=getETof2x1;
810 // m_getETof2x3[i]=getETof2x3;
811 // m_PhiGap[i]=phigap;
812 // m_ThetaGap[i]=thetagap;
813 double dthe = 200.;
814 double dphi = 200.;
815 double dang = 200.;
816
817 // find the nearest charged track
818 for ( int j = 0; j < m_ngch; j++ )
819 {
820
821 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + iGood[j];
822 if ( !( *jtTrk )->isMdcTrackValid() ) continue;
823 RecMdcTrack* jtmdcTrk = ( *jtTrk )->mdcTrack();
824 double jtcharge = jtmdcTrk->charge();
825 if ( !( *jtTrk )->isExtTrackValid() ) continue;
826 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
827 if ( extTrk->emcVolumeNumber() == -1 ) continue;
828 Hep3Vector extpos = extTrk->emcPosition();
829 // double ctht = extpos.cosTheta(emcpos);
830 double angd = extpos.angle( emcpos );
831 double thed = extpos.theta() - emcpos.theta();
832 double phid = extpos.deltaPhi( emcpos );
833 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
834 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
835
836 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
837 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
838 if ( angd < dang ) dang = angd;
839 }
840
841 //
842 // good photon cut will be set here
843 //
844
845 dthe = dthe * 180 / ( CLHEP::pi );
846 dphi = dphi * 180 / ( CLHEP::pi );
847 dang = dang * 180 / ( CLHEP::pi );
848 double eraw = emcTrk->energy();
849 double phi = emcTrk->phi();
850 double the = emcTrk->theta();
851
852 m_delphi[i - m_nchrg] = dphi;
853 m_delthe[i - m_nchrg] = dthe;
854 m_delang[i - m_nchrg] = dang;
855 if ( eraw < m_energyThreshold ) continue;
856 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
857 if ( dang < m_gammaTrkCut ) continue;
858 iGam.push_back( ( *itTrk )->trackId() );
859 }
860
861 int nGam = iGam.size();
862 m_nGam = nGam;
863 // std::cout << "num Good Photon " << m_nGam << " , "
864 // <<evtRecEvent->totalNeutral()<<std::endl;
865 // std::cout<<"dbg_4"<<std::endl;
866 counter[3]++;
867
868 double egam_ext = 0;
869 double ex_gam = 0;
870 double ey_gam = 0;
871 double ez_gam = 0;
872 double et_gam = 0;
873 double e_gam = 0;
874 for ( int i = 0; i < m_nGam; i++ )
875 {
876 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
877 if ( !( *itTrk )->isEmcShowerValid() ) continue;
878 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
879 double eraw = emcTrk->energy();
880 double phi = emcTrk->phi();
881 double the = emcTrk->theta();
882 HepLorentzVector ptrk;
883 ex_gam += eraw * sin( the ) * cos( phi );
884 ey_gam += eraw * sin( the ) * sin( phi );
885 ez_gam += eraw * cos( the );
886 et_gam += eraw * sin( the );
887 e_gam += eraw;
888 if ( eraw >= egam_ext ) { egam_ext = eraw; }
889 }
890
891 double px_had = 0;
892 double py_had = 0;
893 double pz_had = 0;
894 double pt_had = 0;
895 double p_had = 0;
896 double e_had = 0;
897 //
898 // check good charged track's infomation
899 //
900 int ii;
901 m_e_emc[0] = -0.1;
902 m_e_emc[1] = -0.1;
903 for ( int i = 0; i < m_ngch; i++ )
904 {
905
906 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
907
908 if ( !( *itTrk )->isMdcTrackValid() ) continue; // MDC information
909 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
910 // if(!(*itTrk)->isEmcShowerValid()) return StatusCode::SUCCESS;///dbg
911 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
912 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
913 ii = i;
914 // if ( m_ngch==2 &&mdcTrk->charge()>0) ii = 0 ;
915 // if ( m_ngch==2 &&mdcTrk->charge()<0) ii = 1 ;
916
917 m_charge[ii] = mdcTrk->charge();
918 m_vx0[ii] = mdcTrk->x();
919 m_vy0[ii] = mdcTrk->y();
920 m_vz0[ii] = mdcTrk->z();
921
922 m_px[ii] = mdcTrk->px();
923 m_py[ii] = mdcTrk->py();
924 m_pz[ii] = mdcTrk->pz();
925 m_p[ii] = mdcTrk->p();
926
927 if ( m_tagBhabha ) mdcKalTrk->setPidType( RecMdcKalTrack::electron );
928 if ( m_tagDimu ) mdcKalTrk->setPidType( RecMdcKalTrack::muon );
929
930 /// if(m_pidcode[ii]==3)mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
931 /// if(m_pidcode[ii]==4)mdcKalTrk->setPidType(RecMdcKalTrack::proton);
932 m_kal_vx0[ii] = mdcKalTrk->x();
933 m_kal_vy0[ii] = mdcKalTrk->y();
934 m_kal_vz0[ii] = mdcKalTrk->z();
935
936 m_kal_px[ii] = mdcKalTrk->px();
937 m_kal_py[ii] = mdcKalTrk->py();
938 m_kal_pz[ii] = mdcKalTrk->pz();
939 m_kal_p[ii] = mdcKalTrk->p();
940
941 px_had += mdcKalTrk->px();
942 py_had += mdcKalTrk->py();
943 pz_had += mdcKalTrk->pz();
944 pt_had += mdcKalTrk->pxy();
945 p_had += mdcKalTrk->p();
946 e_had += sqrt( mdcKalTrk->p() * mdcKalTrk->p() + mdcKalTrk->mass() * mdcKalTrk->mass() );
947
948 double ptrk = mdcKalTrk->p();
949
950 if ( ( *itTrk )->isMdcDedxValid() )
951 { // DEDX information
952
953 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
954 m_probPH[ii] = dedxTrk->probPH();
955 m_normPH[ii] = dedxTrk->normPH();
956
957 m_chie[ii] = dedxTrk->chiE();
958 m_chimu[ii] = dedxTrk->chiMu();
959 m_chipi[ii] = dedxTrk->chiPi();
960 m_chik[ii] = dedxTrk->chiK();
961 m_chip[ii] = dedxTrk->chiP();
962 m_ghit[ii] = dedxTrk->numGoodHits();
963 m_thit[ii] = dedxTrk->numTotalHits();
964 }
965
966 if ( ( *itTrk )->isEmcShowerValid() )
967 {
968
969 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
970 m_e_emc[ii] = emcTrk->energy();
971 m_phi_emc[ii] = emcTrk->phi();
972 m_theta_emc[ii] = emcTrk->theta();
973 }
974
975 if ( ( *itTrk )->isMucTrackValid() )
976 {
977
978 RecMucTrack* mucTrk = ( *itTrk )->mucTrack();
979 m_nhit_muc[ii] = mucTrk->numHits();
980 m_nlay_muc[ii] = mucTrk->numLayers();
981 }
982
983 if ( ( *itTrk )->isTofTrackValid() )
984 { // TOF information
985
986 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
987
988 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
989
990 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
991 {
992 TofHitStatus* status = new TofHitStatus;
993 status->setStatus( ( *iter_tof )->status() );
994
995 if ( !( status->is_barrel() ) )
996 { // endcap
997 if ( ( status->is_cluster() ) ) m_t_etof[ii] = ( *iter_tof )->tof();
998 if ( !( status->is_counter() ) ) continue; // ?
999 if ( status->layer() != 0 ) continue; // layer1
1000 double path = ( *iter_tof )->path(); // ?
1001 double tof = ( *iter_tof )->tof();
1002 double ph = ( *iter_tof )->ph();
1003 double rhit = ( *iter_tof )->zrhit();
1004 double qual = 0.0 + ( *iter_tof )->quality();
1005 double cntr = 0.0 + ( *iter_tof )->tofID();
1006 double texp[5];
1007 for ( int j = 0; j < 5; j++ )
1008 {
1009 double gb = ptrk / xmass[j];
1010 double beta = gb / sqrt( 1 + gb * gb );
1011 texp[j] = path / beta / velc;
1012 }
1013
1014 m_qual_etof[ii] = qual;
1015 m_tof_etof[ii] = tof;
1016 }
1017 else
1018 { // barrel
1019 if ( ( status->is_cluster() ) ) m_t_btof[ii] = ( *iter_tof )->tof();
1020 if ( !( status->is_counter() ) ) continue; // ?
1021 if ( status->layer() == 1 )
1022 { // layer1
1023 double path = ( *iter_tof )->path(); // ?
1024 double tof = ( *iter_tof )->tof();
1025 double ph = ( *iter_tof )->ph();
1026 double rhit = ( *iter_tof )->zrhit();
1027 double qual = 0.0 + ( *iter_tof )->quality();
1028 double cntr = 0.0 + ( *iter_tof )->tofID();
1029 double texp[5];
1030 for ( int j = 0; j < 5; j++ )
1031 {
1032 double gb = ptrk / xmass[j];
1033 double beta = gb / sqrt( 1 + gb * gb );
1034 texp[j] = path / beta / velc;
1035 }
1036
1037 m_qual_btof1[ii] = qual;
1038 m_tof_btof1[ii] = tof;
1039 }
1040
1041 if ( status->layer() == 2 )
1042 { // layer2
1043 double path = ( *iter_tof )->path(); // ?
1044 double tof = ( *iter_tof )->tof();
1045 double ph = ( *iter_tof )->ph();
1046 double rhit = ( *iter_tof )->zrhit();
1047 double qual = 0.0 + ( *iter_tof )->quality();
1048 double cntr = 0.0 + ( *iter_tof )->tofID();
1049 double texp[5];
1050 for ( int j = 0; j < 5; j++ )
1051 {
1052 double gb = ptrk / xmass[j];
1053 double beta = gb / sqrt( 1 + gb * gb );
1054 texp[j] = path / beta / velc;
1055 }
1056
1057 m_qual_btof2[ii] = qual;
1058 m_tof_btof2[ii] = tof;
1059 }
1060 }
1061 }
1062 }
1063 }
1064 counter[4]++;
1065 // std::cout << "dbg_energy "<< "run, event, E1, E2" << m_run << " " << m_event
1066 // << " " << m_e_emc[0] << " " << m_e_emc[1] << std::endl;
1067
1068 // std::cout<<"dbg_5"<<std::endl;
1069 // tag
1070 m_sbhabhatag = 0;
1071 m_sdimutag = 0;
1072 m_bhabhatag = 0;
1073 m_dimutag = 0;
1074 m_hadrontag = 0;
1075 if ( ( m_ngch == 1 ) || ( m_ngch == 2 && nCharge == 0 ) )
1076 {
1077 EvtRecTrackIterator itTrk1;
1078
1079 EvtRecTrackIterator itTrk2;
1080
1081 RecMdcKalTrack* mdcKalTrk1 = nullptr ;
1082
1083 RecMdcKalTrack* mdcKalTrk2 = nullptr;
1084
1085 HepLorentzVector p41e, p42e, p4le;
1086 Hep3Vector p31e, p32e, p3le;
1087 HepLorentzVector p41m, p42m, p4lm;
1088 Hep3Vector p31m, p32m, p3lm;
1089 HepLorentzVector p41h, p42h, p4lh;
1090 Hep3Vector p31h, p32h, p3lh;
1091 WTrackParameter w1_ini;
1092 WTrackParameter w2_ini;
1093 int iip = -1;
1094 int iim = -1;
1095 for ( int i = 0; i < m_ngch; i++ )
1096 {
1097 if ( m_charge[i] > 0 ) itTrk1 = evtRecTrkCol->begin() + iGood[i];
1098 if ( m_charge[i] < 0 ) itTrk2 = evtRecTrkCol->begin() + iGood[i];
1099 if ( m_charge[i] > 0 ) mdcKalTrk1 = ( *itTrk1 )->mdcKalTrack();
1100 if ( m_charge[i] < 0 ) mdcKalTrk2 = ( *itTrk2 )->mdcKalTrack();
1101 if ( m_charge[i] > 0 ) iip = i;
1102 if ( m_charge[i] < 0 ) iim = i;
1103 if ( m_tagBhabha )
1104 {
1105
1106 if ( m_charge[i] > 0 )
1107 w1_ini =
1108 WTrackParameter( xmass[0], mdcKalTrk1->getZHelixE(), mdcKalTrk1->getZErrorE() );
1109 if ( m_charge[i] < 0 )
1110 w2_ini =
1111 WTrackParameter( xmass[0], mdcKalTrk2->getZHelixE(), mdcKalTrk2->getZErrorE() );
1112 if ( m_charge[i] > 0 ) p41e = w1_ini.p();
1113 if ( m_charge[i] < 0 ) p42e = w2_ini.p();
1114 if ( m_charge[i] > 0 ) p41e.boost( u_cms );
1115 if ( m_charge[i] < 0 ) p42e.boost( u_cms );
1116 if ( m_charge[i] > 0 ) p31e = p41e.vect();
1117 if ( m_charge[i] < 0 ) p32e = p42e.vect();
1118 }
1119
1120 if ( m_tagDimu )
1121 {
1122
1123 if ( m_charge[i] > 0 )
1124 w1_ini = WTrackParameter( xmass[1], mdcKalTrk1->getZHelixMu(),
1125 mdcKalTrk1->getZErrorMu() );
1126 if ( m_charge[i] < 0 )
1127 w2_ini = WTrackParameter( xmass[1], mdcKalTrk2->getZHelixMu(),
1128 mdcKalTrk2->getZErrorMu() );
1129 if ( m_charge[i] > 0 ) p41m = w1_ini.p();
1130 if ( m_charge[i] < 0 ) p42m = w2_ini.p();
1131 if ( m_charge[i] > 0 ) p41m.boost( u_cms );
1132 if ( m_charge[i] < 0 ) p42m.boost( u_cms );
1133 if ( m_charge[i] > 0 ) p31m = p41m.vect();
1134 if ( m_charge[i] < 0 ) p32m = p42m.vect();
1135 }
1136
1137 if ( m_tagBhabha )
1138 {
1139 if ( m_charge[i] > 0 )
1140 {
1141 m_px_cms_ep = p41e.px();
1142 m_py_cms_ep = p41e.py();
1143 m_pz_cms_ep = p41e.pz();
1144 m_e_cms_ep = p41e.e();
1145 }
1146 if ( m_charge[i] < 0 )
1147 {
1148 m_px_cms_em = p42e.px();
1149 m_py_cms_em = p42e.py();
1150 m_pz_cms_em = p42e.pz();
1151 m_e_cms_em = p42e.e();
1152 }
1153 }
1154 if ( m_tagDimu )
1155 {
1156 if ( m_charge[i] > 0 )
1157 {
1158 m_px_cms_ep = p41m.px();
1159 m_py_cms_ep = p41m.py();
1160 m_pz_cms_ep = p41m.pz();
1161 m_e_cms_ep = p41m.e();
1162 }
1163 if ( m_charge[i] < 0 )
1164 {
1165 m_px_cms_em = p42m.px();
1166 m_py_cms_em = p42m.py();
1167 m_pz_cms_em = p42m.pz();
1168 m_e_cms_em = p42m.e();
1169 }
1170 }
1171 }
1172
1173 double e01 = ( iip != -1 ) ? m_e_emc[iip] : 0; // m_e_cms_ep;
1174 double e02 = ( iim != -1 ) ? m_e_emc[iim] : 0; // m_e_cms_em;
1175 int ilarge = ( e01 > e02 ) ? iip : iim;
1176 p4le = ( e01 > e02 ) ? p41e : p42e;
1177 p4lm = ( e01 > e02 ) ? p41m : p42m;
1178
1179 p3le = ( e01 > e02 ) ? p31e : p32e;
1180 p3lm = ( e01 > e02 ) ? p31m : p32m;
1181
1182 double acolle = 180. - p31e.angle( p32e ) * 180.0 / CLHEP::pi;
1183 double acople = 180. - ( p31e.perpPart() ).angle( p32e.perpPart() ) * 180.0 / CLHEP::pi;
1184 double poeb1e = p41e.rho() / beamEnergy;
1185 double poeb2e = p42e.rho() / beamEnergy;
1186 double poeble = p4le.rho() / beamEnergy;
1187 double acollm = 180. - p31m.angle( p32m ) * 180.0 / CLHEP::pi;
1188 double acoplm = 180. - ( p31m.perpPart() ).angle( p32m.perpPart() ) * 180.0 / CLHEP::pi;
1189 double poeb1m = p41m.rho() / beamEnergy;
1190 double poeb2m = p42m.rho() / beamEnergy;
1191 double poeblm = p4lm.rho() / beamEnergy;
1192
1193 double eemc1 = m_e_emc[iip];
1194 double eemc2 = m_e_emc[iim];
1195
1196 double ex1 = m_kal_vx0[iip];
1197 double ey1 = m_kal_vy0[iip];
1198 double ez1 = m_kal_vz0[iip];
1199 double epx1 = m_kal_px[iip];
1200 double epy1 = m_kal_py[iip];
1201 double epz1 = m_kal_pz[iip];
1202 double epp1 = m_kal_p[iip];
1203 double ex2 = m_kal_vx0[iim];
1204 double ey2 = m_kal_vy0[iim];
1205 double ez2 = m_kal_vz0[iim];
1206 double epx2 = m_kal_px[iim];
1207 double epy2 = m_kal_py[iim];
1208 double epz2 = m_kal_pz[iim];
1209 double epp2 = m_kal_p[iim];
1210
1211 double pidchidedx1 = m_pidchiDedx[iip];
1212 double pidchitof11 = m_pidchiTof1[iip];
1213 double pidchitof21 = m_pidchiTof2[iip];
1214 double pidchidedx2 = m_pidchiDedx[iim];
1215 double pidchitof12 = m_pidchiTof1[iim];
1216 double pidchitof22 = m_pidchiTof2[iim];
1217
1218 double eoeb1 = m_e_emc[iip] / beamEnergy;
1219 double eoeb2 = m_e_emc[iim] / beamEnergy;
1220
1221 double eope1 = 0;
1222 if ( p41e.rho() > 0 ) eope1 = m_e_emc[iip] / p41e.rho();
1223 double eope2 = 0;
1224 if ( p42e.rho() > 0 ) eope2 = m_e_emc[iim] / p42e.rho();
1225 double eopm1 = 0;
1226 if ( p41m.rho() > 0 ) eopm1 = m_e_emc[iip] / p41m.rho();
1227 double eopm2 = 0;
1228 if ( p42m.rho() > 0 ) eopm2 = m_e_emc[iim] / p42m.rho();
1229
1230 double exoeb1 =
1231 m_e_emc[iip] * sin( m_theta_emc[iip] ) * cos( m_phi_emc[iip] ) / beamEnergy;
1232 double eyoeb1 =
1233 m_e_emc[iip] * sin( m_theta_emc[iip] ) * sin( m_phi_emc[iip] ) / beamEnergy;
1234 double ezoeb1 = m_e_emc[iip] * cos( m_theta_emc[iip] ) / beamEnergy;
1235 double etoeb1 = m_e_emc[iip] * sin( m_theta_emc[iip] ) / beamEnergy;
1236
1237 double exoeb2 =
1238 m_e_emc[iim] * sin( m_theta_emc[iim] ) * cos( m_phi_emc[iim] ) / beamEnergy;
1239 double eyoeb2 =
1240 m_e_emc[iim] * sin( m_theta_emc[iim] ) * sin( m_phi_emc[iim] ) / beamEnergy;
1241 double ezoeb2 = m_e_emc[iim] * cos( m_theta_emc[iim] ) / beamEnergy;
1242 double etoeb2 = m_e_emc[iim] * sin( m_theta_emc[iim] ) / beamEnergy;
1243
1244 double eoebl = m_e_emc[ilarge] / beamEnergy;
1245
1246 double eopl = 0;
1247 if ( p4lh.rho() > 0 ) eopl = m_e_emc[ilarge] / p4lh.rho();
1248
1249 double exoebl =
1250 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * cos( m_phi_emc[ilarge] ) / beamEnergy;
1251 double eyoebl =
1252 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * sin( m_phi_emc[ilarge] ) / beamEnergy;
1253 double ezoebl = m_e_emc[ilarge] * cos( m_theta_emc[ilarge] ) / beamEnergy;
1254 double etoebl = m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) / beamEnergy;
1255
1256 int mucinfo1 = ( m_nhit_muc[iip] >= 2 && m_nlay_muc[iip] >= 2 ) ? 1 : 0;
1257 int mucinfo2 = ( m_nhit_muc[iim] >= 2 && m_nlay_muc[iim] >= 2 ) ? 1 : 0;
1258 int mucinfol = ( m_nhit_muc[ilarge] >= 2 && m_nlay_muc[ilarge] >= 2 ) ? 1 : 0;
1259 int pidel = ( e01 > e02 ) ? m_nep : m_nem;
1260 int pidmul = ( e01 > e02 ) ? m_nmup : m_nmum;
1261 double deltatof = 0;
1262 // if(m_e_emc[0]<0.4||m_e_emc[1]<0.4)return StatusCode::SUCCESS;///dbg
1263 // if(m_e_emc[0]<1&&m_e_emc[1]<1)return StatusCode::SUCCESS;///dbg
1264 // if(m_e_emc[0]+m_e_emc[1]<2.7)return StatusCode::SUCCESS;///dbg
1265 // double ddphi = (fabs(m_phi[0]-m_phi[1])-CLHEP::pi)*180./CLHEP::pi;///dbg
1266 // double ddthe = (fabs(m_theta[0]+m_theta[1])-CLHEP::pi)*180./CLHEP::pi;///dbg
1267
1268 // if(((ddphi>-25&&ddphi<-4)||(ddphi>2&&ddphi<20)) )return StatusCode::SUCCESS;///dbg
1269 // if(abs(ddthe)<3)return StatusCode::SUCCESS;///dbg
1270 // if(m_tof_btof2[iip]*m_tof_btof2[iim]!=0)
1271 // deltatof+=fabs(m_tof_btof2[iip]-m_tof_btof2[iim]);
1272 // if(m_tof_btof1[iip]*m_tof_btof1[iim]!=0)deltatof+=fabs(m_tof_btof1[iip]-m_tof_btof1[iim]);
1273 // if(m_tof_etof[iip]*m_tof_etof[iim]!=0)deltatof+=fabs(m_tof_etof[iip]-m_tof_etof[iim]);
1274
1275 // if(!m_endcap) {
1276 if ( m_t_btof[iip] * m_t_btof[iim] != 0 ) deltatof = fabs( m_t_btof[iip] - m_t_btof[iim] );
1277 // }
1278 // else {
1279 // if(m_t_etof[iip]*m_t_etof[iim]!=0)deltatof=fabs(m_t_etof[iip]-m_t_etof[iim]);
1280 // }
1281 if ( m_tagBhabha )
1282 {
1283 // if (acolle<m_acoll_e_cut) m_bhabhatag+=1;
1284 // if (acople<m_acopl_e_cut)m_bhabhatag+=10;
1285 // if
1286 // ((fabs(poeb1e-1)<m_poeb_e_cut)&&(fabs(poeb2e-1)<m_poeb_e_cut))m_bhabhatag+=100;
1287 // if (m_useTOF&&(deltatof<m_dtof_e_cut))m_bhabhatag+=1000;
1288 // if (m_useMUC&&(mucinfo1==0||mucinfo2==0))m_bhabhatag+=10000;
1289 // if (m_useEMC&&(eoeb1>m_eoeb_e_cut&&eoeb2>m_eoeb_e_cut))m_bhabhatag+=100000;
1290 // if (m_useEMC&&(eoeb1+eoeb2>m_etotal_e_cut))m_bhabhatag+=1000000;
1291 // if (m_usePID&&(m_nep==1&&m_nem==1))m_bhabhatag+=10000000;
1292
1293 // if (m_ngch==1||(m_ngch == 2 &&acolle<m_acoll_e_cut)) m_sbhabhatag+=1;
1294 // if (m_ngch==1||(m_ngch == 2 &&acople<m_acopl_e_cut))m_sbhabhatag+=10;
1295 // if ((fabs(poeble-1)<m_poeb_e_cut))m_sbhabhatag+=100;
1296 // if (m_ngch==1||(m_ngch == 2
1297 // &&m_useTOF&&(deltatof<m_dtof_e_cut)))m_sbhabhatag+=1000; if
1298 // (m_useMUC&&(mucinfol==0))m_sbhabhatag+=10000; if
1299 // (m_useEMC&&(eoebl>m_eoeb_e_cut))m_sbhabhatag+=100000; if (m_useEMC&&(m_ngch ==
1300 // 1||eoeb1+eoeb2>m_etotal_e_cut))m_sbhabhatag+=1000000; if
1301 // (m_usePID&&(pidel==1))m_sbhabhatag+=10000000;
1302
1303 if ( m_ngch == 2 )
1304 {
1305
1306 if ( acolle < m_acoll_e_cut ) m_bhabhatag += 1;
1307 if ( acople < m_acopl_e_cut ) m_bhabhatag += 10;
1308 if ( sqrt( ( eoeb1 - 1 ) * ( eoeb1 - 1 ) + ( eoeb2 - 1 ) * ( eoeb2 - 1 ) ) <
1309 m_tetotal_e_cut )
1310 m_bhabhatag += 100;
1311 if ( !m_useTOF || ( deltatof < m_dtof_e_cut ) ) m_bhabhatag += 1000;
1312 if ( poeb1e > m_tpoeb_e_cut || poeb2e > m_tpoeb_e_cut ||
1313 ( sqrt( ( poeb1e - 1 ) * ( poeb1e - 1 ) + ( poeb2e - 1 ) * ( poeb2e - 1 ) ) <
1314 m_tptotal_e_cut ) )
1315 m_bhabhatag += 10000;
1316 if ( !m_useMUC || ( mucinfo1 == 0 || mucinfo2 == 0 ) ) m_bhabhatag += 100000;
1317 if ( !m_usePID || ( m_nep == 1 && m_nem == 1 ) ) m_bhabhatag += 1000000;
1318 }
1319
1320 if ( m_ngch == 1 || ( m_ngch == 2 && acolle < m_acoll_e_cut ) ) m_sbhabhatag += 1;
1321 if ( m_ngch == 1 || ( m_ngch == 2 && acople < m_acopl_e_cut ) ) m_sbhabhatag += 10;
1322 if ( ( fabs( poeble - 1 ) < m_poeb_e_cut ) ) m_sbhabhatag += 100;
1323 if ( m_ngch == 1 || !m_useTOF || ( deltatof < m_dtof_e_cut ) ) m_sbhabhatag += 1000;
1324 if ( !m_useMUC || ( mucinfol == 0 ) ) m_sbhabhatag += 10000;
1325 if ( !m_useEMC || ( eoebl > m_eoeb_e_cut ) ) m_sbhabhatag += 100000;
1326 if ( !m_useEMC || ( m_ngch == 1 || eoeb1 + eoeb2 > m_etotal_e_cut ) )
1327 m_sbhabhatag += 1000000;
1328 if ( !m_usePID || ( pidel == 1 ) ) m_sbhabhatag += 10000000;
1329 }
1330
1331 if ( m_tagDimu )
1332 {
1333
1334 if ( m_ngch == 2 )
1335 {
1336
1337 if ( acollm < m_acoll_mu_cut ) m_dimutag += 1;
1338 if ( acoplm < m_acopl_mu_cut ) m_dimutag += 10;
1339 if ( ( sqrt( ( ( poeb1m - poeb2m ) / 0.35 ) * ( ( poeb1m - poeb2m ) / 0.35 ) +
1340 ( ( poeb1m + poeb2m - 1.68 ) / 0.125 ) *
1341 ( ( poeb1m + poeb2m - 1.68 ) / 0.125 ) ) > m_tptotal_mu_cut ) &&
1342 ( ( ( poeb1m >= m_tpoebh_mu_cut ) && ( poeb2m >= m_tpoebl_mu_cut ) ) ||
1343 ( ( poeb2m >= m_tpoebh_mu_cut ) && ( poeb1m >= m_tpoebl_mu_cut ) ) ) )
1344 m_dimutag += 100;
1345 if ( !m_useTOF || ( deltatof < m_dtof_mu_cut ) ) m_dimutag += 1000;
1346 if ( !m_useMUC || ( mucinfo1 == 1 || mucinfo2 == 1 ) ) m_dimutag += 10000;
1347 if ( etoeb1 < m_eoeb_mu_cut && eoeb2 < m_eoeb_mu_cut ) m_dimutag += 100000;
1348 if ( etoeb1 + etoeb2 < m_etotal_mu_cut ) m_dimutag += 1000000;
1349 if ( !m_usePID || ( m_nmup == 1 && m_nmum == 1 ) ) m_dimutag += 10000000;
1350 }
1351
1352 if ( m_ngch == 1 || ( m_ngch == 2 && acollm < m_acoll_mu_cut ) ) m_sdimutag += 1;
1353 if ( m_ngch == 1 || ( m_ngch == 2 && acoplm < m_acopl_mu_cut ) ) m_sdimutag += 10;
1354 if ( ( fabs( poeblm - 1 ) < m_poeb_mu_cut ) ) m_sdimutag += 100;
1355 if ( m_ngch == 1 || !m_useTOF || ( deltatof < m_dtof_mu_cut ) ) m_sdimutag += 1000;
1356 if ( !m_useMUC || ( mucinfo1 == 1 || mucinfo2 == 1 ) ) m_sdimutag += 10000;
1357 if ( !m_useEMC || ( etoebl < m_eoeb_mu_cut ) ) m_sdimutag += 100000;
1358 if ( !m_useEMC || ( m_ngch == 1 || etoeb1 + etoeb2 < m_etotal_mu_cut ) )
1359 m_sdimutag += 1000000;
1360 if ( !m_usePID || ( pidmul == 1 ) ) m_sdimutag += 10000000;
1361 }
1362
1363 if ( m_tagBhabha )
1364 {
1365 m_acoll = acolle;
1366 m_acopl = acople;
1367 m_poeb1 = poeb1e;
1368 m_poeb2 = poeb2e;
1369 m_eop1 = eope1;
1370 m_eop2 = eope2;
1371 m_cos_ep = p41e.cosTheta();
1372 m_cos_em = p42e.cosTheta();
1373 m_mass_ee = ( p41e + p42e ).m();
1374
1375 // if(m_sbhabhatag-int(m_sbhabhatag/10000000)*10000000==1111111||
1376 // m_sbhabhatag-int(m_sbhabhatag/10000000)*10000000==1101111){
1377 if ( m_bhabhatag == 1111111 )
1378 {
1379 nbhabha++;
1380 m_ee_mass->Fill( ( p41e + p42e ).m() );
1381 m_ee_acoll->Fill( acolle );
1382 m_ee_eop_ep->Fill( eope1 );
1383 m_ee_eop_em->Fill( eope2 );
1384 m_ee_costheta_ep->Fill( p41e.cosTheta() );
1385 m_ee_costheta_em->Fill( p42e.cosTheta() );
1386 m_ee_phi_ep->Fill( p41e.phi() );
1387 m_ee_phi_em->Fill( p42e.phi() );
1388 m_ee_nneu->Fill( nneu );
1389
1390 m_ee_eemc_ep->Fill( eemc1 );
1391 m_ee_eemc_em->Fill( eemc2 );
1392
1393 m_ee_x_ep->Fill( ex1 );
1394 m_ee_y_ep->Fill( ey1 );
1395 m_ee_z_ep->Fill( ez1 );
1396
1397 m_ee_x_em->Fill( ex2 );
1398 m_ee_y_em->Fill( ey2 );
1399 m_ee_z_em->Fill( ez2 );
1400
1401 m_ee_px_ep->Fill( epx1 );
1402 m_ee_py_ep->Fill( epy1 );
1403 m_ee_pz_ep->Fill( epz1 );
1404 m_ee_p_ep->Fill( epp1 );
1405
1406 m_ee_px_em->Fill( epx2 );
1407 m_ee_py_em->Fill( epy2 );
1408 m_ee_pz_em->Fill( epz2 );
1409 m_ee_p_em->Fill( epp2 );
1410
1411 m_ee_deltatof->Fill( deltatof );
1412
1413 m_ee_pidchidedx_ep->Fill( pidchidedx1 );
1414 m_ee_pidchidedx_em->Fill( pidchidedx2 );
1415 m_ee_pidchitof1_ep->Fill( pidchitof11 );
1416 m_ee_pidchitof1_em->Fill( pidchitof12 );
1417 m_ee_pidchitof2_ep->Fill( pidchitof21 );
1418 m_ee_pidchitof2_em->Fill( pidchitof22 );
1419 }
1420 }
1421 if ( m_tagDimu )
1422 {
1423 m_acoll = acollm;
1424 m_acopl = acoplm;
1425 m_poeb1 = poeb1m;
1426 m_poeb2 = poeb2m;
1427 m_eop1 = eopm1;
1428 m_eop2 = eopm2;
1429 m_cos_ep = p41m.cosTheta();
1430 m_cos_em = p42m.cosTheta();
1431 m_mass_ee = ( p41m + p42m ).m();
1432
1433 // if(m_sdimutag-int(m_sdimutag/10000000)*10000000==1111112||
1434 // m_sdimutag-int(m_sdimutag/10000000)*10000000==1101112){
1435 if ( m_dimutag == 11111111 )
1436 {
1437 ndimu++;
1438 m_mumu_mass->Fill( ( p41m + p42m ).m() );
1439 m_mumu_acoll->Fill( acollm );
1440 m_mumu_eop_mup->Fill( eopm1 );
1441 m_mumu_eop_mum->Fill( eopm2 );
1442 m_mumu_costheta_mup->Fill( p41m.cosTheta() );
1443 m_mumu_costheta_mum->Fill( p42m.cosTheta() );
1444 m_mumu_phi_mup->Fill( p41m.phi() );
1445 m_mumu_phi_mum->Fill( p42m.phi() );
1446 m_mumu_nneu->Fill( nneu );
1447 m_mumu_nhit->Fill( m_nhit_muc[ilarge] );
1448 m_mumu_nlay->Fill( m_nlay_muc[ilarge] );
1449
1450 m_mumu_eemc_mup->Fill( eemc1 );
1451 m_mumu_eemc_mum->Fill( eemc2 );
1452
1453 m_mumu_x_mup->Fill( ex1 );
1454 m_mumu_y_mup->Fill( ey1 );
1455 m_mumu_z_mup->Fill( ez1 );
1456
1457 m_mumu_x_mum->Fill( ex2 );
1458 m_mumu_y_mum->Fill( ey2 );
1459 m_mumu_z_mum->Fill( ez2 );
1460
1461 m_mumu_px_mup->Fill( epx1 );
1462 m_mumu_py_mup->Fill( epy1 );
1463 m_mumu_pz_mup->Fill( epz1 );
1464 m_mumu_p_mup->Fill( epp1 );
1465
1466 m_mumu_px_mum->Fill( epx2 );
1467 m_mumu_py_mum->Fill( epy2 );
1468 m_mumu_pz_mum->Fill( epz2 );
1469 m_mumu_p_mum->Fill( epp2 );
1470
1471 m_mumu_deltatof->Fill( deltatof );
1472
1473 m_mumu_pidchidedx_mup->Fill( pidchidedx1 );
1474 m_mumu_pidchidedx_mum->Fill( pidchidedx2 );
1475 m_mumu_pidchitof1_mup->Fill( pidchitof11 );
1476 m_mumu_pidchitof1_mum->Fill( pidchitof12 );
1477 m_mumu_pidchitof2_mup->Fill( pidchitof21 );
1478 m_mumu_pidchitof2_mum->Fill( pidchitof22 );
1479 }
1480 }
1481
1482 m_deltatof = deltatof;
1483
1484 m_eoeb1 = eoeb1;
1485 m_eoeb2 = eoeb2;
1486
1487 m_etoeb1 = etoeb1;
1488 m_etoeb2 = etoeb2;
1489 m_mucinfo1 = mucinfo1;
1490 m_mucinfo2 = mucinfo2;
1491 }
1492
1493 // std::cout<<"m_dimutag=="<<m_dimutag<<std::endl;
1494 // if( (m_tagBhabha&&m_bhabhatag==1111111))m_tuple1 -> write();
1495 if ( !m_studymode )
1496 {
1497 if ( ( m_tagDimu && m_dimutag == 11111111 ) || ( m_tagBhabha && m_bhabhatag == 1111111 ) )
1498 m_tuple1->write();
1499 }
1500 else m_tuple1->write();
1501
1502 return StatusCode::SUCCESS;
1503}
1504
1505// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1506// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1507StatusCode ValidPhyJPsill::finalize() {
1508
1509 MsgStream log( msgSvc(), name() );
1510 log << MSG::INFO << "in finalize()" << endmsg;
1511
1512 std::cout << "************ for Singal *******************" << std::endl;
1513 std::cout << "*******************************************" << std::endl;
1514 std::cout << "the total events is " << counter[0] << std::endl;
1515 std::cout << "Good charged tracks " << counter[1] << std::endl;
1516 std::cout << "particle ID " << counter[2] << std::endl;
1517 std::cout << "info. for good charged track " << counter[3] << std::endl;
1518 std::cout << "number of bhabha " << nbhabha << std::endl;
1519 std::cout << "number of dimu " << ndimu << std::endl;
1520
1521 std::cout << "*******************************************" << std::endl;
1522 std::cout << "[OVAL] efficiency of j2ee " << 1.0 * nbhabha / counter[0] << std::endl;
1523 std::cout << "[OVAL] efficiency of j2mumu " << 1.0 * ndimu / counter[0] << std::endl;
1524 return StatusCode::SUCCESS;
1525}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
u_cms
const Hep3Vector u_cms
Definition DQADtagAlg.cxx:56
mks0
const double mks0
Definition DQAJpsi2PPbarAlg.cxx:40
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
n
const Int_t n
Definition DataBase/tau_mode.c:57
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi
double mpi
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
xmass
const double xmass[5]
Definition Gam4pikp.cxx:35
velc
const double velc
Definition Gam4pikp.cxx:36
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:37
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_KeyIHVP INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition KK2f.h:50
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DstEmcShower::dy
double dy() const
Definition DstEmcShower.cxx:10
DstEmcShower::dz
double dz() const
Definition DstEmcShower.cxx:17
DstEmcShower::dx
double dx() const
Definition DstEmcShower.cxx:3
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::instance
static ParticleID * instance()
Definition ParticleID.cxx:21
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:26
ParticleID::chiDedx
double chiDedx(int n) const
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:115
ValidPhyJPsill::initialize
StatusCode initialize()
Definition ValidPhyJPsill.cxx:111
ValidPhyJPsill::finalize
StatusCode finalize()
Definition ValidPhyJPsill.cxx:1507
ValidPhyJPsill::execute
StatusCode execute()
Definition ValidPhyJPsill.cxx:486
ValidPhyJPsill::ValidPhyJPsill
ValidPhyJPsill(const std::string &name, ISvcLocator *pSvcLocator)
Definition ValidPhyJPsill.cxx:63