BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
PPbar.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/EventHeader.h"
9#include "EventModel/EventModel.h"
10
11#include "DstEvent/TofHitStatus.h"
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14
15#include "VertexDbSvc/IVertexDbSvc.h"
16#include "VertexFit/KinematicFit.h"
17#include "VertexFit/SecondVertexFit.h"
18#include "VertexFit/VertexFit.h"
19
20#include "ParticleID/ParticleID.h"
21
22#include "CLHEP/Vector/LorentzVector.h"
23#include "CLHEP/Vector/ThreeVector.h"
24#include "CLHEP/Vector/TwoVector.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/IHistogramSvc.h"
27#include "GaudiKernel/INTupleSvc.h"
28#include "GaudiKernel/ISvcLocator.h"
29#include "GaudiKernel/NTuple.h"
30#include "TMath.h"
31using CLHEP::Hep2Vector;
32using CLHEP::Hep3Vector;
33using CLHEP::HepLorentzVector;
34#include "CLHEP/Geometry/Point3D.h"
35#ifndef ENABLE_BACKWARDS_COMPATIBILITY
36typedef HepGeom::Point3D<double> HepPoint3D;
37#endif
38#include "PPbar.h"
39
40#include "ParticleID/ParticleID.h"
41#include "VertexFit/KinematicFit.h"
42#include "VertexFit/VertexFit.h"
43
44#include <vector>
45// const double twopi = 6.2831853;
46// const double pi = 3.1415927;
47const double mpi0 = 0.134977;
48const double mks0 = 0.497648;
49const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
50// const double velc = 29.9792458; tof_path unit in cm.
51const double velc = 299.792458; // tof path unit in mm
52typedef std::vector<int> Vint;
53typedef std::vector<HepLorentzVector> Vp4;
54
55// boosted
56const HepLorentzVector p_cms( 0.034067, 0.0, 0.0, 3.097 );
57const Hep3Vector u_cms = -p_cms.boostVector();
58
59// declare one counter
60static int counter[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
61/////////////////////////////////////////////////////////////////////////////
62DECLARE_COMPONENT( PPbar )
63PPbar::PPbar( const std::string& name, ISvcLocator* pSvcLocator )
64 : Algorithm( name, pSvcLocator ) {
65
66 // Declare the properties
67 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
68 declareProperty( "Vz0cut", m_vz0cut = 10.0 );
69 declareProperty( "Coscut", m_coscut = 0.93 );
70 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
71 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
72 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
73 declareProperty( "Test4C", m_test4C = 1 );
74 declareProperty( "Test5C", m_test5C = 1 );
75 declareProperty( "CheckDedx", m_checkDedx = 1 );
76 declareProperty( "CheckTof", m_checkTof = 1 );
77
78 declareProperty( "tagPPbar", m_tagPPbar = false );
79}
80
81// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
82StatusCode PPbar::initialize() {
83 MsgStream log( msgSvc(), name() );
84
85 log << MSG::INFO << "in initialize()" << endmsg;
86
87 StatusCode status;
88
89 status = service( "THistSvc", m_thistsvc );
90 if ( status.isFailure() )
91 {
92 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
93 return status;
94 }
95
96 if ( m_tagPPbar )
97 {
98
99 m_pp_mass = new TH1F( "pp_mass", "pp_mass", 100, 3.0, 3.2 );
100 status = m_thistsvc->regHist( "/VAL/PHY/pp_mass", m_pp_mass );
101 m_pp_angle = new TH1F( "pp_angle", "pp_angle", 100, 170, 180 );
102 status = m_thistsvc->regHist( "/VAL/PHY/pp_angle", m_pp_angle );
103 m_pp_deltatof = new TH1F( "pp_deltatof", "pp_deltatof", 100, 0, 10 );
104 status = m_thistsvc->regHist( "/VAL/PHY/pp_deltatof", m_pp_deltatof );
105
106 m_pp_x_pp = new TH1F( "pp_x_pp", "pp_x_pp", 100, -0.5, 0.5 );
107 status = m_thistsvc->regHist( "/VAL/PHY/pp_x_pp", m_pp_x_pp );
108 m_pp_y_pp = new TH1F( "pp_y_pp", "pp_y_pp", 100, -0.5, 0.5 );
109 status = m_thistsvc->regHist( "/VAL/PHY/pp_y_pp", m_pp_y_pp );
110 m_pp_z_pp = new TH1F( "pp_z_pp", "pp_z_pp", 100, -10.0, 10. );
111 status = m_thistsvc->regHist( "/VAL/PHY/pp_z_pp", m_pp_z_pp );
112 m_pp_r_pp = new TH1F( "pp_r_pp", "pp_r_pp", 100, 0.0, 0.5 );
113 status = m_thistsvc->regHist( "/VAL/PHY/pp_r_pp", m_pp_r_pp );
114 m_pp_px_pp = new TH1F( "pp_px_pp", "pp_px_pp", 100, -1.5, 1.5 );
115 status = m_thistsvc->regHist( "/VAL/PHY/pp_px_pp", m_pp_px_pp );
116 m_pp_py_pp = new TH1F( "pp_py_pp", "pp_py_pp", 100, -1.5, 1.5 );
117 status = m_thistsvc->regHist( "/VAL/PHY/pp_py_pp", m_pp_py_pp );
118 m_pp_pz_pp = new TH1F( "pp_pz_pp", "pp_pz_pp", 100, -1.5, 1.5 );
119 status = m_thistsvc->regHist( "/VAL/PHY/pp_pz_pp", m_pp_pz_pp );
120 m_pp_p_pp = new TH1F( "pp_p_pp", "pp_p_pp", 100, 1.15, 1.35 );
121 status = m_thistsvc->regHist( "/VAL/PHY/pp_p_pp", m_pp_p_pp );
122 m_pp_cos_pp = new TH1F( "pp_cos_pp", "pp_cos_pp", 100, -1.0, 1.0 );
123 status = m_thistsvc->regHist( "/VAL/PHY/pp_cos_pp", m_pp_cos_pp );
124 m_pp_emc_pp = new TH1F( "pp_emc_pp", "pp_emc_pp", 100, 0.0, 2.0 );
125 status = m_thistsvc->regHist( "/VAL/PHY/pp_emc_pp", m_pp_emc_pp );
126 m_pp_pidchidedx_pp = new TH1F( "pp_pidchidedx_pp", "pp_pidchidedx_pp", 100, -10.0, 10.0 );
127 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchidedx_pp", m_pp_pidchidedx_pp );
128 m_pp_pidchitof1_pp = new TH1F( "pp_pidchitof1_pp", "pp_pidchitof1_pp", 100, -10.0, 10.0 );
129 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchitof1_pp", m_pp_pidchitof1_pp );
130 m_pp_pidchitof2_pp = new TH1F( "pp_pidchitof2_pp", "pp_pidchitof2_pp", 100, -10.0, 10.0 );
131 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchitof2_pp", m_pp_pidchitof2_pp );
132
133 m_pp_x_pb = new TH1F( "pp_x_pb", "pp_x_pb", 100, -0.5, 0.5 );
134 status = m_thistsvc->regHist( "/VAL/PHY/pp_x_pb", m_pp_x_pb );
135 m_pp_y_pb = new TH1F( "pp_y_pb", "pp_y_pb", 100, -0.5, 0.5 );
136 status = m_thistsvc->regHist( "/VAL/PHY/pp_y_pb", m_pp_y_pb );
137 m_pp_z_pb = new TH1F( "pp_z_pb", "pp_z_pb", 100, -10.0, 10. );
138 status = m_thistsvc->regHist( "/VAL/PHY/pp_z_pb", m_pp_z_pb );
139 m_pp_r_pb = new TH1F( "pp_r_pb", "pp_r_pb", 100, 0.0, 0.5 );
140 status = m_thistsvc->regHist( "/VAL/PHY/pp_r_pb", m_pp_r_pb );
141 m_pp_px_pb = new TH1F( "pp_px_pb", "pp_px_pb", 100, -1.5, 1.5 );
142 status = m_thistsvc->regHist( "/VAL/PHY/pp_px_pb", m_pp_px_pb );
143 m_pp_py_pb = new TH1F( "pp_py_pb", "pp_py_pb", 100, -1.5, 1.5 );
144 status = m_thistsvc->regHist( "/VAL/PHY/pp_py_pb", m_pp_py_pb );
145 m_pp_pz_pb = new TH1F( "pp_pz_pb", "pp_pz_pb", 100, -1.5, 1.5 );
146 status = m_thistsvc->regHist( "/VAL/PHY/pp_pz_pb", m_pp_pz_pb );
147 m_pp_p_pb = new TH1F( "pp_p_pb", "pp_p_pb", 100, 1.15, 1.35 );
148 status = m_thistsvc->regHist( "/VAL/PHY/pp_p_pb", m_pp_p_pb );
149 m_pp_cos_pb = new TH1F( "pp_cos_pb", "pp_cos_pb", 100, -1.0, 1.0 );
150 status = m_thistsvc->regHist( "/VAL/PHY/pp_cos_pb", m_pp_cos_pb );
151 m_pp_emc_pb = new TH1F( "pp_emc_pb", "pp_emc_pb", 100, 0.0, 2.0 );
152 status = m_thistsvc->regHist( "/VAL/PHY/pp_emc_pb", m_pp_emc_pb );
153 m_pp_pidchidedx_pb = new TH1F( "pp_pidchidedx_pb", "pp_pidchidedx_pb", 100, -10.0, 10.0 );
154 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchidedx_pb", m_pp_pidchidedx_pb );
155 m_pp_pidchitof1_pb = new TH1F( "pp_pidchitof1_pb", "pp_pidchitof1_pb", 100, -10.0, 10.0 );
156 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchitof1_pb", m_pp_pidchitof1_pb );
157 m_pp_pidchitof2_pb = new TH1F( "pp_pidchitof2_pb", "pp_pidchitof2_pb", 100, -10.0, 10.0 );
158 status = m_thistsvc->regHist( "/VAL/PHY/pp_pidchitof2_pb", m_pp_pidchitof2_pb );
159 }
160
161 NTuplePtr nt1( ntupleSvc(), "FILE1/testv" );
162 if ( nt1 ) m_tuple1 = nt1;
163 else
164 {
165 m_tuple1 = ntupleSvc()->book( "FILE1/testv", CLID_ColumnWiseTuple, "N-Tuple" );
166 if ( m_tuple1 )
167 {
168 status = m_tuple1->addItem( "irun", m_run );
169 status = m_tuple1->addItem( "ievent", m_event );
170 status = m_tuple1->addItem( "Nchrg", m_nchrg );
171 status = m_tuple1->addItem( "Nneu", m_nneu );
172 status = m_tuple1->addItem( "NGch", m_ngch, 0, 10 );
173
174 status = m_tuple1->addIndexedItem( "charge", m_ngch, m_charge );
175 status = m_tuple1->addIndexedItem( "vx", m_ngch, m_vx0 );
176 status = m_tuple1->addIndexedItem( "vy", m_ngch, m_vy0 );
177 status = m_tuple1->addIndexedItem( "vz", m_ngch, m_vz0 );
178 status = m_tuple1->addIndexedItem( "vr", m_ngch, m_vr0 );
179 status = m_tuple1->addIndexedItem( "px", m_ngch, m_px );
180 status = m_tuple1->addIndexedItem( "py", m_ngch, m_py );
181 status = m_tuple1->addIndexedItem( "pz", m_ngch, m_pz );
182 status = m_tuple1->addIndexedItem( "p", m_ngch, m_p );
183 status = m_tuple1->addIndexedItem( "cos", m_ngch, m_cos );
184
185 status = m_tuple1->addIndexedItem( "bst_px", m_ngch, m_bst_px );
186 status = m_tuple1->addIndexedItem( "bst_py", m_ngch, m_bst_py );
187 status = m_tuple1->addIndexedItem( "bst_pz", m_ngch, m_bst_pz );
188 status = m_tuple1->addIndexedItem( "bst_p", m_ngch, m_bst_p );
189 status = m_tuple1->addIndexedItem( "bst_cos", m_ngch, m_bst_cos );
190
191 status = m_tuple1->addIndexedItem( "kal_vx", m_ngch, m_kal_vx0 );
192 status = m_tuple1->addIndexedItem( "kal_vy", m_ngch, m_kal_vy0 );
193 status = m_tuple1->addIndexedItem( "kal_vz", m_ngch, m_kal_vz0 );
194 status = m_tuple1->addIndexedItem( "kal_vr", m_ngch, m_kal_vr0 );
195
196 status = m_tuple1->addIndexedItem( "kal_px", m_ngch, m_kal_px );
197 status = m_tuple1->addIndexedItem( "kal_py", m_ngch, m_kal_py );
198 status = m_tuple1->addIndexedItem( "kal_pz", m_ngch, m_kal_pz );
199 status = m_tuple1->addIndexedItem( "kal_p", m_ngch, m_kal_p );
200 status = m_tuple1->addIndexedItem( "kal_cos", m_ngch, m_kal_cos );
201
202 status = m_tuple1->addIndexedItem( "bst_kal_px", m_ngch, m_bst_kal_px );
203 status = m_tuple1->addIndexedItem( "bst_kal_py", m_ngch, m_bst_kal_py );
204 status = m_tuple1->addIndexedItem( "bst_kal_pz", m_ngch, m_bst_kal_pz );
205 status = m_tuple1->addIndexedItem( "bst_kal_p", m_ngch, m_bst_kal_p );
206 status = m_tuple1->addIndexedItem( "bst_kal_cos", m_ngch, m_bst_kal_cos );
207
208 status = m_tuple1->addIndexedItem( "vtx_px", m_ngch, m_vtx_px );
209 status = m_tuple1->addIndexedItem( "vtx_py", m_ngch, m_vtx_py );
210 status = m_tuple1->addIndexedItem( "vtx_pz", m_ngch, m_vtx_pz );
211 status = m_tuple1->addIndexedItem( "vtx_p", m_ngch, m_vtx_p );
212 status = m_tuple1->addIndexedItem( "vtx_cos", m_ngch, m_vtx_cos );
213
214 status = m_tuple1->addIndexedItem( "chie", m_ngch, m_chie );
215 status = m_tuple1->addIndexedItem( "chimu", m_ngch, m_chimu );
216 status = m_tuple1->addIndexedItem( "chipi", m_ngch, m_chipi );
217 status = m_tuple1->addIndexedItem( "chik", m_ngch, m_chik );
218 status = m_tuple1->addIndexedItem( "chip", m_ngch, m_chip );
219 status = m_tuple1->addIndexedItem( "ghit", m_ngch, m_ghit );
220 status = m_tuple1->addIndexedItem( "thit", m_ngch, m_thit );
221 status = m_tuple1->addIndexedItem( "probPH", m_ngch, m_probPH );
222 status = m_tuple1->addIndexedItem( "normPH", m_ngch, m_normPH );
223
224 status = m_tuple1->addIndexedItem( "e_emc", m_ngch, m_e_emc );
225
226 status = m_tuple1->addIndexedItem( "clus_tof", m_ngch, m_clus_tof );
227 status = m_tuple1->addIndexedItem( "clus_texp_e", m_ngch, m_clus_texp_e );
228 status = m_tuple1->addIndexedItem( "clus_texp_mu", m_ngch, m_clus_texp_mu );
229 status = m_tuple1->addIndexedItem( "clus_texp_pi", m_ngch, m_clus_texp_pi );
230 status = m_tuple1->addIndexedItem( "clus_texp_k", m_ngch, m_clus_texp_k );
231 status = m_tuple1->addIndexedItem( "clus_texp_p", m_ngch, m_clus_texp_p );
232 status = m_tuple1->addIndexedItem( "clus_toff_e", m_ngch, m_clus_toff_e );
233 status = m_tuple1->addIndexedItem( "clus_toff_mu", m_ngch, m_clus_toff_mu );
234 status = m_tuple1->addIndexedItem( "clus_toff_pi", m_ngch, m_clus_toff_pi );
235 status = m_tuple1->addIndexedItem( "clus_toff_k", m_ngch, m_clus_toff_k );
236 status = m_tuple1->addIndexedItem( "clus_toff_p", m_ngch, m_clus_toff_p );
237 status = m_tuple1->addIndexedItem( "clus_tsig_e", m_ngch, m_clus_tsig_e );
238 status = m_tuple1->addIndexedItem( "clus_tsig_mu", m_ngch, m_clus_tsig_mu );
239 status = m_tuple1->addIndexedItem( "clus_tsig_pi", m_ngch, m_clus_tsig_pi );
240 status = m_tuple1->addIndexedItem( "clus_tsig_k", m_ngch, m_clus_tsig_k );
241 status = m_tuple1->addIndexedItem( "clus_tsig_p", m_ngch, m_clus_tsig_p );
242 status = m_tuple1->addIndexedItem( "clus_path", m_ngch, m_clus_path );
243 status = m_tuple1->addIndexedItem( "clus_zrhit", m_ngch, m_clus_zrhit );
244 status = m_tuple1->addIndexedItem( "clus_ph", m_ngch, m_clus_ph );
245 status = m_tuple1->addIndexedItem( "clus_beta", m_ngch, m_clus_beta );
246 status = m_tuple1->addIndexedItem( "clus_qual", m_ngch, m_clus_qual );
247 status = m_tuple1->addIndexedItem( "clus_t0", m_ngch, m_clus_t0 );
248
249 status = m_tuple1->addIndexedItem( "cntr_etof", m_ngch, m_cntr_etof );
250 status = m_tuple1->addIndexedItem( "ptot_etof", m_ngch, m_ptot_etof );
251 status = m_tuple1->addIndexedItem( "ph_etof", m_ngch, m_ph_etof );
252 status = m_tuple1->addIndexedItem( "rhit_etof", m_ngch, m_rhit_etof );
253 status = m_tuple1->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
254 status = m_tuple1->addIndexedItem( "tof_etof", m_ngch, m_tof_etof );
255 status = m_tuple1->addIndexedItem( "te_etof", m_ngch, m_te_etof );
256 status = m_tuple1->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
257 status = m_tuple1->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
258 status = m_tuple1->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
259 status = m_tuple1->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
260
261 status = m_tuple1->addIndexedItem( "cntr_btof1", m_ngch, m_cntr_btof1 );
262 status = m_tuple1->addIndexedItem( "ptot_btof1", m_ngch, m_ptot_btof1 );
263 status = m_tuple1->addIndexedItem( "ph_btof1", m_ngch, m_ph_btof1 );
264 status = m_tuple1->addIndexedItem( "zhit_btof1", m_ngch, m_zhit_btof1 );
265 status = m_tuple1->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
266 status = m_tuple1->addIndexedItem( "tof_btof1", m_ngch, m_tof_btof1 );
267 status = m_tuple1->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
268 status = m_tuple1->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
269 status = m_tuple1->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
270 status = m_tuple1->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
271 status = m_tuple1->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
272
273 status = m_tuple1->addIndexedItem( "cntr_btof2", m_ngch, m_cntr_btof2 );
274 status = m_tuple1->addIndexedItem( "ptot_btof2", m_ngch, m_ptot_btof2 );
275 status = m_tuple1->addIndexedItem( "ph_btof2", m_ngch, m_ph_btof2 );
276 status = m_tuple1->addIndexedItem( "zhit_btof2", m_ngch, m_zhit_btof2 );
277 status = m_tuple1->addIndexedItem( "qual_btof2", m_ngch, m_qual_btof2 );
278 status = m_tuple1->addIndexedItem( "tof_btof2", m_ngch, m_tof_btof2 );
279 status = m_tuple1->addIndexedItem( "te_btof2", m_ngch, m_te_btof2 );
280 status = m_tuple1->addIndexedItem( "tmu_btof2", m_ngch, m_tmu_btof2 );
281 status = m_tuple1->addIndexedItem( "tpi_btof2", m_ngch, m_tpi_btof2 );
282 status = m_tuple1->addIndexedItem( "tk_btof2", m_ngch, m_tk_btof2 );
283 status = m_tuple1->addIndexedItem( "tp_btof2", m_ngch, m_tp_btof2 );
284
285 status = m_tuple1->addIndexedItem( "m_ptrk_pid", m_ngch, m_ptrk_pid );
286 status = m_tuple1->addIndexedItem( "m_cost_pid", m_ngch, m_cost_pid );
287 status = m_tuple1->addIndexedItem( "m_dedx_pid", m_ngch, m_dedx_pid );
288 status = m_tuple1->addIndexedItem( "m_tof1_pid", m_ngch, m_tof1_pid );
289 status = m_tuple1->addIndexedItem( "m_tof2_pid", m_ngch, m_tof2_pid );
290 status = m_tuple1->addIndexedItem( "m_prob_pi", m_ngch, m_prob_pi );
291 status = m_tuple1->addIndexedItem( "m_prob_k", m_ngch, m_prob_k );
292 status = m_tuple1->addIndexedItem( "m_prob_p", m_ngch, m_prob_p );
293
294 status = m_tuple1->addItem( "np", m_np );
295 status = m_tuple1->addItem( "npb", m_npb );
296
297 status = m_tuple1->addItem( "m2p", m_m2p );
298 status = m_tuple1->addItem( "angle", m_angle );
299 status = m_tuple1->addItem( "deltatof", m_deltatof );
300
301 status = m_tuple1->addItem( "kal_m2p", m_kal_m2p );
302 status = m_tuple1->addItem( "kal_angle", m_kal_angle );
303
304 status = m_tuple1->addItem( "vtx_m2p", m_vtx_m2p );
305 status = m_tuple1->addItem( "vtx_angle", m_vtx_angle );
306
307 //
308 // status = m_tuple1->addItem ( "Jpx", m_Jpx );
309 // status = m_tuple1->addItem ( "Jpy", m_Jpy );
310 // status = m_tuple1->addItem ( "Jpz", m_Jpz );
311 // status = m_tuple1->addItem ( "Jpp", m_Jpp );
312 // status = m_tuple1->addItem ( "m2p", m_m2p );
313 //
314 // status = m_tuple1->addItem ( "p1pp", m_p1pp );
315 // status = m_tuple1->addItem ( "p1px", m_p1px );
316 // status = m_tuple1->addItem ( "p1py", m_p1py );
317 // status = m_tuple1->addItem ( "p1pz", m_p1pz );
318 // status = m_tuple1->addItem ( "p1cos", m_p1cos );
319 //
320 // status = m_tuple1->addItem ( "p2pp", m_p2pp );
321 // status = m_tuple1->addItem ( "p2px", m_p2px );
322 // status = m_tuple1->addItem ( "p2py", m_p2py );
323 // status = m_tuple1->addItem ( "p2pz", m_p2pz );
324 // status = m_tuple1->addItem ( "p2cos",m_p2cos);
325 //
326 // status = m_tuple1->addItem ( "angle", m_angle );
327 //
328 }
329 else
330 {
331 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
332 return StatusCode::FAILURE;
333 }
334 }
335
336 //
337 //--------end of book--------
338 //
339
340 log << MSG::INFO << "successfully return from initialize()" << endmsg;
341 return StatusCode::SUCCESS;
342}
343
344// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
345StatusCode PPbar::execute() {
346
347 MsgStream log( msgSvc(), name() );
348 log << MSG::INFO << "in execute()" << endmsg;
349
350 counter[0]++;
351
352 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
353 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
354 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
355
356 m_run = eventHeader->runNumber();
357 m_event = eventHeader->eventNumber();
358 m_nchrg = evtRecEvent->totalCharged();
359 m_nneu = evtRecEvent->totalNeutral();
360
361 log << MSG::INFO << "get event tag OK" << endmsg;
362 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
363 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
364
365 //
366 // check x0, y0, z0, r0
367 // suggest cut: |z0|<10 && r0<1
368 //
369 Hep3Vector xorigin( 0, 0, 0 );
370 IVertexDbSvc* vtxsvc;
371 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
372 if ( vtxsvc->isVertexValid() )
373 {
374 double* dbv = vtxsvc->PrimaryVertex();
375 double* vv = vtxsvc->SigmaPrimaryVertex();
376 xorigin.setX( dbv[0] );
377 xorigin.setY( dbv[1] );
378 xorigin.setZ( dbv[2] );
379 }
380
381 Vint iGood;
382 iGood.clear();
383 int nCharge = 0;
384 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
385 {
386 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
387 if ( !( *itTrk )->isMdcTrackValid() ) continue;
388 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
389 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
390 double x0 = mdcTrk->x();
391 double y0 = mdcTrk->y();
392 double z0 = mdcTrk->z();
393 double phi0 = mdcTrk->helix( 1 );
394 double xv = xorigin.x();
395 double yv = xorigin.y();
396 double zv = xorigin.z();
397 double rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
398 double cost = cos( mdcTrk->theta() );
399 if ( fabs( z0 - zv ) >= m_vz0cut ) continue;
400 if ( fabs( rv ) >= m_vr0cut ) continue;
401 if ( fabs( cost ) >= m_coscut ) continue;
402
403 iGood.push_back( ( *itTrk )->trackId() );
404 nCharge += mdcTrk->charge();
405 }
406
407 //
408 // Finish Good Charged Track Selection
409 //
410 int nGood = iGood.size();
411 m_ngch = iGood.size();
412 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
413 //
414 // Charge track number cut
415 //
416
417 if ( ( nGood != 2 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
418 counter[1]++;
419
420 //
421 // Assign 4-momentum to each charged track
422 //
423
424 Vint ipp, ipm;
425 ipp.clear();
426 ipm.clear();
427 Vp4 p_pp, p_pm;
428 p_pp.clear();
429 p_pm.clear();
430 Vp4 p_kal_pp, p_kal_pm;
431 p_kal_pp.clear();
432 p_kal_pm.clear();
433 RecMdcTrack * ppTrk, *pmTrk;
434 RecMdcKalTrack *ppKalTrk = 0, *pmKalTrk = 0;
435 int ii = -1;
436
437 ParticleID* pid = ParticleID::instance();
438 for ( int i = 0; i < nGood; i++ )
439 {
440
441 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
442 if ( !( *itTrk )->isMdcTrackValid() ) continue;
443 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
444 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
445 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
446
447 // if(pid) delete pid;
448 pid->init();
449 pid->setMethod( pid->methodProbability() );
450 pid->setChiMinCut( 4 );
451 pid->setRecTrack( *itTrk );
452 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
453 pid->identify( pid->onlyPionKaonProton() ); // seperater Pion/Kaon/Proton
454 // pid->identify(pid->onlyProton());
455 // pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
456 // pid->identify(pid->onlyPion());
457 // pid->identify(pid->onlyKaon());
458 pid->calculate();
459 if ( !( pid->IsPidInfoValid() ) ) continue;
460
461 double prob_pi = pid->probPion();
462 double prob_k = pid->probKaon();
463 double prob_p = pid->probProton();
464
465 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
466 // if(pid->probPion() < 0.001) continue;
467 // if (prob_p > prob_pi && prob_p > prob_k) {
468 if ( prob_p > prob_pi && prob_p > prob_k )
469 {
470
471 HepLorentzVector ptrk, pkaltrk;
472
473 ptrk.setPx( mdcTrk->px() );
474 ptrk.setPy( mdcTrk->py() );
475 ptrk.setPz( mdcTrk->pz() );
476 double p3 = ptrk.mag();
477 ptrk.setE( sqrt( p3 * p3 + xmass[4] * xmass[4] ) );
478
479 mdcKalTrk->setPidType( RecMdcKalTrack::proton );
480 pkaltrk.setPx( mdcKalTrk->px() );
481 pkaltrk.setPy( mdcKalTrk->py() );
482 pkaltrk.setPz( mdcKalTrk->pz() );
483 p3 = pkaltrk.mag();
484 pkaltrk.setE( sqrt( p3 * p3 + xmass[4] * xmass[4] ) );
485
486 if ( mdcTrk->charge() > 0 )
487 {
488 ii = 0;
489
490 ipp.push_back( iGood[i] );
491 p_pp.push_back( ptrk );
492 p_kal_pp.push_back( pkaltrk );
493 ppTrk = mdcTrk;
494 ppKalTrk = mdcKalTrk;
495 }
496 if ( mdcTrk->charge() < 0 )
497 {
498 ii = 1;
499
500 ipm.push_back( iGood[i] );
501 p_pm.push_back( ptrk );
502 p_kal_pm.push_back( pkaltrk );
503 pmTrk = mdcTrk;
504 pmKalTrk = mdcKalTrk;
505 }
506 m_ptrk_pid[ii] = mdcTrk->p();
507 m_cost_pid[ii] = cos( mdcTrk->theta() );
508 m_dedx_pid[ii] = pid->chiDedx( 4 );
509 m_tof1_pid[ii] = pid->chiTof1( 4 );
510 m_tof2_pid[ii] = pid->chiTof2( 4 );
511 m_prob_pi[ii] = pid->probPion();
512 m_prob_k[ii] = pid->probKaon();
513 m_prob_p[ii] = pid->probProton();
514 }
515 }
516 m_np = ipp.size();
517 m_npb = ipm.size();
518
519 if ( m_np * m_npb != 1 ) return StatusCode::SUCCESS;
520
521 counter[2]++;
522
523 // boosted mdcTrk
524 p_pp[0].boost( u_cms );
525 p_pm[0].boost( u_cms );
526 for ( int i = 0; i <= 1; i++ )
527 {
528 HepLorentzVector p;
529 if ( i == 0 ) p = p_pp[0];
530 if ( i == 1 ) p = p_pm[0];
531
532 m_bst_px[i] = p.px();
533 m_bst_py[i] = p.py();
534 m_bst_pz[i] = p.pz();
535 m_bst_p[i] = p.rho();
536 m_bst_cos[i] = p.cosTheta();
537 }
538
539 m_m2p = ( p_pp[0] + p_pm[0] ).m();
540 m_angle = ( p_pp[0].vect() ).angle( ( p_pm[0] ).vect() ) * 180.0 / ( CLHEP::pi );
541
542 // boosted mdcKalTrk
543 p_kal_pp[0].boost( u_cms );
544 p_kal_pm[0].boost( u_cms );
545
546 for ( int i = 0; i <= 1; i++ )
547 {
548 HepLorentzVector p;
549 if ( i == 0 ) p = p_kal_pp[0];
550 if ( i == 1 ) p = p_kal_pm[0];
551
552 m_bst_kal_px[i] = p.px();
553 m_bst_kal_py[i] = p.py();
554 m_bst_kal_pz[i] = p.pz();
555 m_bst_kal_p[i] = p.rho();
556 m_bst_kal_cos[i] = p.cosTheta();
557 }
558 m_kal_m2p = ( p_kal_pp[0] + p_kal_pm[0] ).m();
559 m_kal_angle = ( p_kal_pp[0].vect() ).angle( ( p_kal_pm[0] ).vect() ) * 180.0 / ( CLHEP::pi );
560
561 //
562 //
563 // full information for charged tracks
564 //
565 //
566
567 ////////////////////////////////////////////
568 // Get all information for good charged tracks
569 ////////////////////////////////////////////
570
571 for ( int i = 0; i < m_ngch; i++ )
572 {
573 ////////////////////////////////////////////
574 // MDC information
575 ////////////////////////////////////////////
576 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
577 if ( !( *itTrk )->isMdcTrackValid() ) continue;
578 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
579 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
580 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
581
582 if ( mdcTrk->charge() > 0 ) { ii = 0; }
583 else { ii = 1; }
584
585 m_charge[ii] = mdcTrk->charge();
586
587 double x0 = mdcTrk->x();
588 double y0 = mdcTrk->y();
589 double z0 = mdcTrk->z();
590 double phi0 = mdcTrk->helix( 1 );
591 double xv = xorigin.x();
592 double yv = xorigin.y();
593 double zv = xorigin.z();
594 double rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
595
596 m_vx0[ii] = x0 - xv;
597 m_vy0[ii] = y0 - yv;
598 m_vz0[ii] = z0 - zv;
599 m_vr0[ii] = rv;
600
601 m_px[ii] = mdcTrk->px();
602 m_py[ii] = mdcTrk->py();
603 m_pz[ii] = mdcTrk->pz();
604 m_p[ii] = mdcTrk->p();
605 m_cos[ii] = mdcTrk->pz() / mdcTrk->p();
606
607 mdcKalTrk->setPidType( RecMdcKalTrack::proton );
608 x0 = mdcKalTrk->x();
609 y0 = mdcKalTrk->y();
610 z0 = mdcKalTrk->z();
611 rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
612
613 m_kal_vx0[ii] = x0 - xv;
614 m_kal_vy0[ii] = y0 - yv;
615 m_kal_vz0[ii] = z0 - zv;
616 m_kal_vr0[ii] = rv;
617
618 m_kal_px[ii] = mdcKalTrk->px();
619 m_kal_py[ii] = mdcKalTrk->py();
620 m_kal_pz[ii] = mdcKalTrk->pz();
621 m_kal_p[ii] = mdcKalTrk->p();
622 m_kal_cos[ii] = mdcKalTrk->pz() / mdcKalTrk->p();
623
624 double ptrk = mdcKalTrk->p();
625
626 ////////////////////////////////////////////
627 // DEDX information
628 ////////////////////////////////////////////
629 if ( ( *itTrk )->isMdcDedxValid() )
630 { // Dedx information
631 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
632 m_chie[ii] = dedxTrk->chiE();
633 m_chimu[ii] = dedxTrk->chiMu();
634 m_chipi[ii] = dedxTrk->chiPi();
635 m_chik[ii] = dedxTrk->chiK();
636 m_chip[ii] = dedxTrk->chiP();
637 m_ghit[ii] = dedxTrk->numGoodHits();
638 m_thit[ii] = dedxTrk->numTotalHits();
639 m_probPH[ii] = dedxTrk->probPH();
640 m_normPH[ii] = dedxTrk->normPH();
641 }
642
643 ////////////////////////////////////////////
644 // TOF information
645 ////////////////////////////////////////////
646 log << MSG::INFO << "TOF info " << endmsg;
647 m_clus_tof[ii] = 100.0;
648 m_clus_texp_e[ii] = 100.0;
649 m_clus_texp_mu[ii] = 100.0;
650 m_clus_texp_pi[ii] = 100.0;
651 m_clus_texp_k[ii] = 100.0;
652 m_clus_texp_p[ii] = 100.0;
653
654 m_clus_toff_e[ii] = 100.0;
655 m_clus_toff_mu[ii] = 100.0;
656 m_clus_toff_pi[ii] = 100.0;
657 m_clus_toff_k[ii] = 100.0;
658 m_clus_toff_p[ii] = 100.0;
659
660 m_clus_tsig_e[ii] = 100.0;
661 m_clus_tsig_mu[ii] = 100.0;
662 m_clus_tsig_pi[ii] = 100.0;
663 m_clus_tsig_k[ii] = 100.0;
664 m_clus_tsig_p[ii] = 100.0;
665
666 m_clus_path[ii] = 100.0;
667 m_clus_zrhit[ii] = 100.0;
668 m_clus_ph[ii] = 100.0;
669 m_clus_beta[ii] = 100.0;
670 m_clus_qual[ii] = 100.0;
671 m_clus_t0[ii] = 100.0;
672
673 if ( ( *itTrk )->isEmcShowerValid() )
674 {
675 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
676 m_e_emc[ii] = emcTrk->energy();
677 }
678
679 if ( ( *itTrk )->isTofTrackValid() )
680 {
681
682 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
683
684 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
685 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
686 {
687 TofHitStatus* status = new TofHitStatus;
688 status->setStatus( ( *iter_tof )->status() );
689
690 if ( !( status->is_barrel() ) )
691 { // endcap
692 if ( !( status->is_counter() ) ) continue; // ?
693 if ( status->layer() != 0 ) continue; // layer1
694 double path = ( *iter_tof )->path(); // ?
695 double tof = ( *iter_tof )->tof();
696 double ph = ( *iter_tof )->ph();
697 double rhit = ( *iter_tof )->zrhit();
698 double qual = 0.0 + ( *iter_tof )->quality();
699 double cntr = 0.0 + ( *iter_tof )->tofID();
700 double texp[5];
701 for ( int j = 0; j < 5; j++ )
702 {
703 double gb = ptrk / xmass[j];
704 double beta = gb / sqrt( 1 + gb * gb );
705 texp[j] = path / beta / velc;
706 }
707
708 m_cntr_etof[ii] = cntr;
709 m_ptot_etof[ii] = ptrk;
710 m_ph_etof[ii] = ph;
711 m_rhit_etof[ii] = rhit;
712 m_qual_etof[ii] = qual;
713 m_tof_etof[ii] = tof;
714 m_te_etof[ii] = tof - texp[0];
715 m_tmu_etof[ii] = tof - texp[1];
716 m_tpi_etof[ii] = tof - texp[2];
717 m_tk_etof[ii] = tof - texp[3];
718 m_tp_etof[ii] = tof - texp[4];
719 }
720 else
721 { // barrel
722 if ( ( status->is_cluster() ) )
723 {
724 //////////////////////////////////////////////////////////
725 /// please use is_cluster rather is_counter
726 /// please use texp(i) rahter than calculations by yourself
727 // if( !(status->is_counter()) ) continue; // ?
728 // if(status->layer()==1){ //barrel TOF layer1
729 m_clus_tof[ii] = ( *iter_tof )->tof();
730 m_clus_texp_e[ii] = ( *iter_tof )->texp( 0 );
731 m_clus_texp_mu[ii] = ( *iter_tof )->texp( 1 );
732 m_clus_texp_pi[ii] = ( *iter_tof )->texp( 2 );
733 m_clus_texp_k[ii] = ( *iter_tof )->texp( 3 );
734 m_clus_texp_p[ii] = ( *iter_tof )->texp( 4 );
735
736 m_clus_toff_e[ii] = ( *iter_tof )->toffset( 0 );
737 m_clus_toff_mu[ii] = ( *iter_tof )->toffset( 1 );
738 m_clus_toff_pi[ii] = ( *iter_tof )->toffset( 2 );
739 m_clus_toff_k[ii] = ( *iter_tof )->toffset( 3 );
740 m_clus_toff_p[ii] = ( *iter_tof )->toffset( 4 );
741
742 m_clus_tsig_e[ii] = ( *iter_tof )->sigma( 0 );
743 m_clus_tsig_mu[ii] = ( *iter_tof )->sigma( 1 );
744 m_clus_tsig_pi[ii] = ( *iter_tof )->sigma( 2 );
745 m_clus_tsig_k[ii] = ( *iter_tof )->sigma( 3 );
746 m_clus_tsig_p[ii] = ( *iter_tof )->sigma( 4 );
747
748 m_clus_path[ii] = ( *iter_tof )->path();
749 m_clus_zrhit[ii] = ( *iter_tof )->zrhit();
750 m_clus_ph[ii] = ( *iter_tof )->ph();
751 m_clus_beta[ii] = ( *iter_tof )->beta();
752 m_clus_qual[ii] = ( *iter_tof )->quality();
753 m_clus_t0[ii] = ( *iter_tof )->t0();
754 }
755 if ( !( status->is_counter() ) ) continue; // ?
756 if ( status->layer() == 1 )
757 { // layer1
758 double path = ( *iter_tof )->path(); // ?
759 double tof = ( *iter_tof )->tof();
760 double ph = ( *iter_tof )->ph();
761 double rhit = ( *iter_tof )->zrhit();
762 double qual = 0.0 + ( *iter_tof )->quality();
763 double cntr = 0.0 + ( *iter_tof )->tofID();
764 double texp[5];
765 for ( int j = 0; j < 5; j++ )
766 {
767 double gb = ptrk / xmass[j];
768 double beta = gb / sqrt( 1 + gb * gb );
769 texp[j] = path / beta / velc;
770 }
771
772 m_cntr_btof1[ii] = cntr;
773 m_ptot_btof1[ii] = ptrk;
774 m_ph_btof1[ii] = ph;
775 m_zhit_btof1[ii] = rhit;
776 m_qual_btof1[ii] = qual;
777 m_tof_btof1[ii] = tof;
778 m_te_btof1[ii] = tof - texp[0];
779 m_tmu_btof1[ii] = tof - texp[1];
780 m_tpi_btof1[ii] = tof - texp[2];
781 m_tk_btof1[ii] = tof - texp[3];
782 m_tp_btof1[ii] = tof - texp[4];
783 }
784
785 if ( status->layer() == 2 )
786 { // layer2
787 double path = ( *iter_tof )->path(); // ?
788 double tof = ( *iter_tof )->tof();
789 double ph = ( *iter_tof )->ph();
790 double rhit = ( *iter_tof )->zrhit();
791 double qual = 0.0 + ( *iter_tof )->quality();
792 double cntr = 0.0 + ( *iter_tof )->tofID();
793 double texp[5];
794 for ( int j = 0; j < 5; j++ )
795 {
796 double gb = ptrk / xmass[j];
797 double beta = gb / sqrt( 1 + gb * gb );
798 texp[j] = path / beta / velc;
799 }
800
801 m_cntr_btof2[ii] = cntr;
802 m_ptot_btof2[ii] = ptrk;
803 m_ph_btof2[ii] = ph;
804 m_zhit_btof2[ii] = rhit;
805 m_qual_btof2[ii] = qual;
806 m_tof_btof2[ii] = tof;
807 m_te_btof2[ii] = tof - texp[0];
808 m_tmu_btof2[ii] = tof - texp[1];
809 m_tpi_btof2[ii] = tof - texp[2];
810 m_tk_btof2[ii] = tof - texp[3];
811 m_tp_btof2[ii] = tof - texp[4];
812 }
813 }
814 delete status;
815 }
816 }
817 double deltatof = 0.0;
818 if ( m_tof_btof1[0] * m_tof_btof1[1] != 0.0 )
819 deltatof += fabs( m_tof_btof1[0] - m_tof_btof1[1] );
820 if ( m_tof_btof2[0] * m_tof_btof2[1] != 0.0 )
821 deltatof += fabs( m_tof_btof2[0] - m_tof_btof2[1] );
822 if ( m_tof_etof[0] * m_tof_etof[1] != 0.0 )
823 deltatof += fabs( m_tof_etof[0] - m_tof_etof[1] );
824 m_deltatof = deltatof;
825 }
826
827 counter[3]++;
828
829 if ( m_tagPPbar )
830 {
831 m_pp_mass->Fill( m_kal_m2p );
832 m_pp_angle->Fill( m_kal_angle );
833 m_pp_deltatof->Fill( m_deltatof );
834
835 m_pp_x_pp->Fill( m_vx0[0] );
836 m_pp_y_pp->Fill( m_vy0[0] );
837 m_pp_z_pp->Fill( m_vz0[0] );
838 m_pp_r_pp->Fill( m_vr0[0] );
839 m_pp_px_pp->Fill( m_bst_kal_px[0] );
840 m_pp_py_pp->Fill( m_bst_kal_py[0] );
841 m_pp_pz_pp->Fill( m_bst_kal_pz[0] );
842 m_pp_p_pp->Fill( m_bst_kal_p[0] );
843 m_pp_cos_pp->Fill( m_bst_kal_cos[0] );
844 m_pp_emc_pp->Fill( m_e_emc[0] );
845 m_pp_pidchidedx_pp->Fill( m_dedx_pid[0] );
846 m_pp_pidchitof1_pp->Fill( m_tof1_pid[0] );
847 m_pp_pidchitof2_pp->Fill( m_tof2_pid[0] );
848
849 m_pp_x_pb->Fill( m_vx0[1] );
850 m_pp_y_pb->Fill( m_vy0[1] );
851 m_pp_z_pb->Fill( m_vz0[1] );
852 m_pp_r_pb->Fill( m_vr0[1] );
853 m_pp_px_pb->Fill( m_bst_kal_px[1] );
854 m_pp_py_pb->Fill( m_bst_kal_py[1] );
855 m_pp_pz_pb->Fill( m_bst_kal_pz[1] );
856 m_pp_p_pb->Fill( m_bst_kal_p[1] );
857 m_pp_cos_pb->Fill( m_bst_kal_cos[1] );
858 m_pp_emc_pb->Fill( m_e_emc[1] );
859 m_pp_pidchidedx_pb->Fill( m_dedx_pid[1] );
860 m_pp_pidchitof1_pb->Fill( m_tof1_pid[1] );
861 m_pp_pidchitof2_pb->Fill( m_tof2_pid[1] );
862 }
863 //
864 // Test vertex fit
865 //
866
867 HepPoint3D vx( 0., 0., 0. );
868 HepSymMatrix Evx( 3, 0 );
869 double bx = 1E+6;
870 double by = 1E+6;
871 double bz = 1E+6;
872 Evx[0][0] = bx * bx;
873 Evx[1][1] = by * by;
874 Evx[2][2] = bz * bz;
875
876 VertexParameter vxpar;
877 vxpar.setVx( vx );
878 vxpar.setEvx( Evx );
879
880 VertexFit* vtxfit = VertexFit::instance();
881 vtxfit->init();
882
883 WTrackParameter wp( xmass[4], ppKalTrk->getZHelixP(), ppKalTrk->getZErrorP() );
884 WTrackParameter wpb( xmass[4], pmKalTrk->getZHelixP(), pmKalTrk->getZErrorP() );
885
886 vtxfit->AddTrack( 0, wp );
887 vtxfit->AddTrack( 1, wpb );
888 vtxfit->AddVertex( 0, vxpar, 0, 1 );
889 if ( !vtxfit->Fit( 0 ) ) return StatusCode::SUCCESS;
890 vtxfit->Swim( 0 );
891
892 counter[4]++;
893
894 wp = vtxfit->wtrk( 0 );
895 wpb = vtxfit->wtrk( 1 );
896
897 //
898 // boosted
899 //
900 HepLorentzVector p = wp.p();
901 HepLorentzVector pb = wpb.p();
902
903 p.boost( u_cms );
904 pb.boost( u_cms );
905
906 HepLorentzVector ptot = p + pb;
907
908 for ( int i = 0; i < 2; i++ )
909 {
910 HepLorentzVector p_vtx;
911 if ( m_charge[i] > 0 )
912 {
913 p_vtx = p;
914 ii = 0;
915 }
916 else
917 {
918 p_vtx = pb;
919 ii = 1;
920 }
921
922 m_vtx_px[ii] = p_vtx.px();
923 m_vtx_py[ii] = p_vtx.py();
924 m_vtx_pz[ii] = p_vtx.pz();
925 m_vtx_p[ii] = p_vtx.rho();
926 m_vtx_cos[ii] = p_vtx.cosTheta();
927 }
928
929 m_vtx_m2p = ( p + pb ).m();
930 m_vtx_angle = ( p.vect() ).angle( pb.vect() ) * 180.0 / ( CLHEP::pi );
931
932 counter[5]++;
933 m_tuple1->write();
934
935 return StatusCode::SUCCESS;
936}
937
938// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
939StatusCode PPbar::finalize() {
940
941 MsgStream log( msgSvc(), name() );
942 log << MSG::INFO << "in finalize()" << endmsg;
943
944 std::cout << "************* PPbar.cxx ****************" << std::endl;
945 std::cout << "the totale events is " << counter[0] << std::endl;
946 std::cout << "select good charged track " << counter[1] << std::endl;
947 std::cout << "particle ID " << counter[2] << std::endl;
948 std::cout << "inform. for charged track " << counter[3] << std::endl;
949 std::cout << "vertex fit " << counter[4] << std::endl;
950 std::cout << "boosted information " << counter[5] << std::endl;
951 std::cout << "****************************************" << std::endl;
952
953 return StatusCode::SUCCESS;
954}
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)
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi0
double mpi0
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
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
PPbar
Definition PPbar.h:10
PPbar::PPbar
PPbar(const std::string &name, ISvcLocator *pSvcLocator)
Definition PPbar.cxx:63
PPbar::finalize
StatusCode finalize()
Definition PPbar.cxx:939
PPbar::execute
StatusCode execute()
Definition PPbar.cxx:345
PPbar::initialize
StatusCode initialize()
Definition PPbar.cxx:82
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
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:91
VertexFit::init
void init()
Definition VertexFit.cxx:27
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288