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