BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
KsKpi.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 "GaudiKernel/Bootstrap.h"
8#include "GaudiKernel/IHistogramSvc.h"
9#include "GaudiKernel/INTupleSvc.h"
10#include "GaudiKernel/NTuple.h"
11
12#include "EventModel/Event.h"
13#include "EventModel/EventHeader.h"
14#include "EventModel/EventModel.h"
15
16#include "DstEvent/TofHitStatus.h"
17#include "EvtRecEvent/EvtRecEvent.h"
18#include "EvtRecEvent/EvtRecTrack.h"
19
20#include "McTruth/McParticle.h"
21
22#include "VertexDbSvc/IVertexDbSvc.h"
23#include "VertexFit/KinematicFit.h"
24#include "VertexFit/SecondVertexFit.h"
25#include "VertexFit/VertexFit.h"
26
27#include "ParticleID/ParticleID.h"
28
29#include "TMath.h"
30
31#include "CLHEP/Vector/LorentzVector.h"
32#include "CLHEP/Vector/ThreeVector.h"
33#include "CLHEP/Vector/TwoVector.h"
34#include "KsKpi.h"
35
36using CLHEP::Hep2Vector;
37using CLHEP::Hep3Vector;
38using CLHEP::HepLorentzVector;
39#include "CLHEP/Geometry/Point3D.h"
40#ifndef ENABLE_BACKWARDS_COMPATIBILITY
41typedef HepGeom::Point3D<double> HepPoint3D;
42#endif
43
44#include <vector>
45const double mpi0 = 0.134977;
46const double mks0 = 0.497648;
47const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
48// const double velc = 29.9792458; tof_path unit in cm.
49const double velc = 299.792458; // tof path unit in mm
50typedef std::vector<int> Vint;
51typedef std::vector<HepLorentzVector> Vp4;
52
53// boost
54const HepLorentzVector p_cms( 0.034067, 0.0, 0.0, 3.097 );
55const Hep3Vector u_cms = -p_cms.boostVector();
56
57// declare one counter
58static int counter[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
59/////////////////////////////////////////////////////////////////////////////
60DECLARE_COMPONENT( KsKpi )
61KsKpi::KsKpi( const std::string& name, ISvcLocator* pSvcLocator )
62 : Algorithm( name, pSvcLocator ) {
63
64 // Declare the properties
65 declareProperty( "Vr0cut", m_vr0cut = 5.0 );
66 declareProperty( "Vz0cut", m_vz0cut = 20.0 );
67 declareProperty( "Coscut", m_coscut = 0.93 );
68
69 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
70 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
71 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
72 declareProperty( "Test4C", m_test4C = 1 );
73 declareProperty( "Test5C", m_test5C = 1 );
74 declareProperty( "CheckDedx", m_checkDedx = 1 );
75 declareProperty( "CheckTof", m_checkTof = 1 );
76
77 declareProperty( "tagKsKpi", m_tagKsKpi = false );
78}
79
80// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
81// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
82StatusCode KsKpi::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_tagKsKpi )
97 {
98
99 m_kskpi_vx_pi1 = new TH1F( "kskpi_vx_pi1", "kskpi_vx_pi1", 100, -5.0, 5.0 );
100 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vx_pi1", m_kskpi_vx_pi1 );
101 m_kskpi_vy_pi1 = new TH1F( "kskpi_vy_pi1", "kskpi_vy_pi1", 100, -5.0, 5.0 );
102 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vy_pi1", m_kskpi_vy_pi1 );
103 m_kskpi_vz_pi1 = new TH1F( "kskpi_vz_pi1", "kskpi_vz_pi1", 100, -20.0, 20.0 );
104 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vz_pi1", m_kskpi_vz_pi1 );
105 m_kskpi_vr_pi1 = new TH1F( "kskpi_vr_pi1", "kskpi_vr_pi1", 100, 0.0, 5.0 );
106 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vr_pi1", m_kskpi_vr_pi1 );
107 m_kskpi_px_pi1 = new TH1F( "kskpi_px_pi1", "kskpi_px_pi1", 100, -1.5, 1.5 );
108 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_px_pi1", m_kskpi_px_pi1 );
109 m_kskpi_py_pi1 = new TH1F( "kskpi_py_pi1", "kskpi_py_pi1", 100, -1.5, 1.5 );
110 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_py_pi1", m_kskpi_py_pi1 );
111 m_kskpi_pz_pi1 = new TH1F( "kskpi_pz_pi1", "kskpi_pz_pi1", 100, -1.5, 1.5 );
112 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pz_pi1", m_kskpi_pz_pi1 );
113 m_kskpi_pp_pi1 = new TH1F( "kskpi_pp_pi1", "kskpi_pp_pi1", 100, 0.0, 1.5 );
114 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pp_pi1", m_kskpi_pp_pi1 );
115 m_kskpi_cos_pi1 = new TH1F( "kskpi_cos_pi1", "kskpi_cos_pi1", 100, -1.0, 1.0 );
116 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_cos_pi1", m_kskpi_cos_pi1 );
117 m_kskpi_emc_pi1 = new TH1F( "kskpi_emc_pi1", "kskpi_emc_pi1", 100, 0.0, 1.5 );
118 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_emc_pi1", m_kskpi_emc_pi1 );
119
120 m_kskpi_vx_pi2 = new TH1F( "kskpi_vx_pi2", "kskpi_vx_pi2", 100, -5.0, 5.0 );
121 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vx_pi2", m_kskpi_vx_pi2 );
122 m_kskpi_vy_pi2 = new TH1F( "kskpi_vy_pi2", "kskpi_vy_pi2", 100, -5.0, 5.0 );
123 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vy_pi2", m_kskpi_vy_pi2 );
124 m_kskpi_vz_pi2 = new TH1F( "kskpi_vz_pi2", "kskpi_vz_pi2", 100, -20.0, 20.0 );
125 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vz_pi2", m_kskpi_vz_pi2 );
126 m_kskpi_vr_pi2 = new TH1F( "kskpi_vr_pi2", "kskpi_vr_pi2", 100, 0.0, 5.0 );
127 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vr_pi2", m_kskpi_vr_pi2 );
128 m_kskpi_px_pi2 = new TH1F( "kskpi_px_pi2", "kskpi_px_pi2", 100, -1.5, 1.5 );
129 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_px_pi2", m_kskpi_px_pi2 );
130 m_kskpi_py_pi2 = new TH1F( "kskpi_py_pi2", "kskpi_py_pi2", 100, -1.5, 1.5 );
131 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_py_pi2", m_kskpi_py_pi2 );
132 m_kskpi_pz_pi2 = new TH1F( "kskpi_pz_pi2", "kskpi_pz_pi2", 100, -1.5, 1.5 );
133 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pz_pi2", m_kskpi_pz_pi2 );
134 m_kskpi_pp_pi2 = new TH1F( "kskpi_pp_pi2", "kskpi_pp_pi2", 100, 0.0, 1.5 );
135 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pp_pi2", m_kskpi_pp_pi2 );
136 m_kskpi_cos_pi2 = new TH1F( "kskpi_cos_pi2", "kskpi_cos_pi2", 100, -1.0, 1.0 );
137 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_cos_pi2", m_kskpi_cos_pi2 );
138 m_kskpi_emc_pi2 = new TH1F( "kskpi_emc_pi2", "kskpi_emc_pi2", 100, 0.0, 1.5 );
139 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_emc_pi2", m_kskpi_emc_pi2 );
140
141 m_kskpi_vx_pi = new TH1F( "kskpi_vx_pi", "kskpi_vx_pi", 100, -1.0, 1.0 );
142 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vx_pi", m_kskpi_vx_pi );
143 m_kskpi_vy_pi = new TH1F( "kskpi_vy_pi", "kskpi_vy_pi", 100, -1.0, 1.0 );
144 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vy_pi", m_kskpi_vy_pi );
145 m_kskpi_vz_pi = new TH1F( "kskpi_vz_pi", "kskpi_vz_pi", 100, -10.0, 10.0 );
146 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vz_pi", m_kskpi_vz_pi );
147 m_kskpi_vr_pi = new TH1F( "kskpi_vr_pi", "kskpi_vr_pi", 100, 0.0, 1.0 );
148 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vr_pi", m_kskpi_vr_pi );
149 m_kskpi_px_pi = new TH1F( "kskpi_px_pi", "kskpi_px_pi", 100, -1.5, 1.5 );
150 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_px_pi", m_kskpi_px_pi );
151 m_kskpi_py_pi = new TH1F( "kskpi_py_pi", "kskpi_py_pi", 100, -1.5, 1.5 );
152 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_py_pi", m_kskpi_py_pi );
153 m_kskpi_pz_pi = new TH1F( "kskpi_pz_pi", "kskpi_pz_pi", 100, -1.5, 1.5 );
154 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pz_pi", m_kskpi_pz_pi );
155 m_kskpi_pp_pi = new TH1F( "kskpi_pp_pi", "kskpi_pp_pi", 100, 0.0, 1.5 );
156 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pp_pi", m_kskpi_pp_pi );
157 m_kskpi_cos_pi = new TH1F( "kskpi_cos_pi", "kskpi_cos_pi", 100, -1.0, 1.0 );
158 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_cos_pi", m_kskpi_cos_pi );
159 m_kskpi_emc_pi = new TH1F( "kskpi_emc_pi", "kskpi_emc_pi", 100, 0.0, 1.5 );
160 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_emc_pi", m_kskpi_emc_pi );
161
162 m_kskpi_vx_k = new TH1F( "kskpi_vx_k", "kskpi_vx_k", 100, -1.0, 1.0 );
163 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vx_k", m_kskpi_vx_k );
164 m_kskpi_vy_k = new TH1F( "kskpi_vy_k", "kskpi_vy_k", 100, -1.0, 1.0 );
165 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vy_k", m_kskpi_vy_k );
166 m_kskpi_vz_k = new TH1F( "kskpi_vz_k", "kskpi_vz_k", 100, -10.0, 10.0 );
167 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vz_k", m_kskpi_vz_k );
168 m_kskpi_vr_k = new TH1F( "kskpi_vr_k", "kskpi_vr_k", 100, 0.0, 1.0 );
169 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vr_k", m_kskpi_vr_k );
170 m_kskpi_px_k = new TH1F( "kskpi_px_k", "kskpi_px_k", 100, -1.5, 1.5 );
171 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_px_k", m_kskpi_px_k );
172 m_kskpi_py_k = new TH1F( "kskpi_py_k", "kskpi_py_k", 100, -1.5, 1.5 );
173 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_py_k", m_kskpi_py_k );
174 m_kskpi_pz_k = new TH1F( "kskpi_pz_k", "kskpi_pz_k", 100, -1.5, 1.5 );
175 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pz_k", m_kskpi_pz_k );
176 m_kskpi_pp_k = new TH1F( "kskpi_pp_k", "kskpi_pp_k", 100, 0.0, 1.5 );
177 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pp_k", m_kskpi_pp_k );
178 m_kskpi_cos_k = new TH1F( "kskpi_cos_k", "kskpi_cos_k", 100, -1.0, 1.0 );
179 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_cos_k", m_kskpi_cos_k );
180 m_kskpi_emc_k = new TH1F( "kskpi_emc_k", "kskpi_emc_k", 100, 0.0, 1.5 );
181 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_emc_k", m_kskpi_emc_k );
182
183 m_kskpi_pidchidedx_1 =
184 new TH1F( "kskpi_pidchidedx_1", "kskpi_pidchidedx_1", 100, -10.0, 10.0 );
185 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchidedx_1", m_kskpi_pidchidedx_1 );
186 m_kskpi_pidchitof1_1 =
187 new TH1F( "kskpi_pidchitof1_1", "kskpi_pidchitof1_1", 100, -10.0, 10.0 );
188 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof1_1", m_kskpi_pidchitof1_1 );
189 m_kskpi_pidchitof2_1 =
190 new TH1F( "kskpi_pidchitof2_1", "kskpi_pidchitof2_1", 100, -10.0, 10.0 );
191 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof2_1", m_kskpi_pidchitof2_1 );
192
193 m_kskpi_pidchidedx_2 =
194 new TH1F( "kskpi_pidchidedx_2", "kskpi_pidchidedx_2", 100, -10.0, 10.0 );
195 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchidedx_2", m_kskpi_pidchidedx_2 );
196 m_kskpi_pidchitof1_2 =
197 new TH1F( "kskpi_pidchitof1_2", "kskpi_pidchitof1_2", 100, -10.0, 10.0 );
198 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof1_2", m_kskpi_pidchitof1_2 );
199 m_kskpi_pidchitof2_2 =
200 new TH1F( "kskpi_pidchitof2_2", "kskpi_pidchitof2_2", 100, -10.0, 10.0 );
201 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof2_2", m_kskpi_pidchitof2_2 );
202
203 m_kskpi_pidchidedx_3 =
204 new TH1F( "kskpi_pidchidedx_3", "kskpi_pidchidedx_3", 100, -10.0, 10.0 );
205 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchidedx_3", m_kskpi_pidchidedx_3 );
206 m_kskpi_pidchitof1_3 =
207 new TH1F( "kskpi_pidchitof1_3", "kskpi_pidchitof1_3", 100, -10.0, 10.0 );
208 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof1_3", m_kskpi_pidchitof1_3 );
209 m_kskpi_pidchitof2_3 =
210 new TH1F( "kskpi_pidchitof2_3", "kskpi_pidchitof2_3", 100, -10.0, 10.0 );
211 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof2_3", m_kskpi_pidchitof2_3 );
212
213 m_kskpi_pidchidedx_4 =
214 new TH1F( "kskpi_pidchidedx_4", "kskpi_pidchidedx_4", 100, -10.0, 10.0 );
215 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchidedx_4", m_kskpi_pidchidedx_4 );
216 m_kskpi_pidchitof1_4 =
217 new TH1F( "kskpi_pidchitof1_4", "kskpi_pidchitof1_4", 100, -10.0, 10.0 );
218 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof1_4", m_kskpi_pidchitof1_4 );
219 m_kskpi_pidchitof2_4 =
220 new TH1F( "kskpi_pidchitof2_4", "kskpi_pidchitof2_4", 100, -10.0, 10.0 );
221 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_pidchitof2_4", m_kskpi_pidchitof2_4 );
222
223 m_kskpi_vfits_chi = new TH1F( "kskpi_vfits_chi", "kskpi_vfits_chi", 100, 0.0, 20.0 );
224 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfits_chi", m_kskpi_vfits_chi );
225 m_kskpi_vfits_vx = new TH1F( "kskpi_vfits_vx", "kskpi_vfits_vx", 100, -20.0, 20.0 );
226 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfits_vx", m_kskpi_vfits_vx );
227 m_kskpi_vfits_vy = new TH1F( "kskpi_vfits_vy", "kskpi_vfits_vy", 100, -20.0, 20.0 );
228 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfits_vy", m_kskpi_vfits_vy );
229 m_kskpi_vfits_vz = new TH1F( "kskpi_vfits_vz", "kskpi_vfits_vz", 100, -20.0, 20.0 );
230 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfits_vz", m_kskpi_vfits_vz );
231 m_kskpi_vfits_vr = new TH1F( "kskpi_vfits_vr", "kskpi_vfits_vr", 100, 0.0, 20.0 );
232 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfits_vr", m_kskpi_vfits_vr );
233
234 m_kskpi_vfitp_chi = new TH1F( "kskpi_vfitp_chi", "kskpi_vfitp_chi", 100, 0.0, 50.0 );
235 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfitp_chi", m_kskpi_vfitp_chi );
236 m_kskpi_vfitp_vx = new TH1F( "kskpi_vfitp_vx", "kskpi_vfitp_vx", 100, -1.0, 1.0 );
237 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfitp_vx", m_kskpi_vfitp_vx );
238 m_kskpi_vfitp_vy = new TH1F( "kskpi_vfitp_vy", "kskpi_vfitp_vy", 100, -1.0, 1.0 );
239 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfitp_vy", m_kskpi_vfitp_vy );
240 m_kskpi_vfitp_vz = new TH1F( "kskpi_vfitp_vz", "kskpi_vfitp_vz", 100, -5.0, 5.0 );
241 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfitp_vz", m_kskpi_vfitp_vz );
242 m_kskpi_vfitp_vr = new TH1F( "kskpi_vfitp_vr", "kskpi_vfitp_vr", 100, 0.0, 1.0 );
243 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfitp_vr", m_kskpi_vfitp_vr );
244
245 m_kskpi_vfit2_chi = new TH1F( "kskpi_vfit2_chi", "kskpi_vfit2_chi", 100, 0.0, 20.0 );
246 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfit2_chi", m_kskpi_vfit2_chi );
247 m_kskpi_vfit2_mks = new TH1F( "kskpi_vfit2_mks", "kskpi_vfit2_mks", 100, 0.4, 0.6 );
248 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfit2_mks", m_kskpi_vfit2_mks );
249 m_kskpi_vfit2_ct = new TH1F( "kskpi_vfit2_ct", "kskpi_vfit2_ct", 100, -3.0, 13.0 );
250 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfit2_ct", m_kskpi_vfit2_ct );
251 m_kskpi_vfit2_dl = new TH1F( "kskpi_vfit2_dl", "kskpi_vfit2_dl", 100, -5.0, 25.0 );
252 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfit2_dl", m_kskpi_vfit2_dl );
253 m_kskpi_vfit2_dle = new TH1F( "kskpi_vfit2_dle", "kskpi_vfit2_dle", 100, 0.0, 1.0 );
254 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_vfit2_dle", m_kskpi_vfit2_dle );
255
256 m_kskpi_4c_chi = new TH1F( "kskpi_4c_chi", "kskpi_4c_chi", 100, 0.0, 50 );
257 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_chi", m_kskpi_4c_chi );
258 m_kskpi_4c_mks = new TH1F( "kskpi_4c_mks", "kskpi_4c_mks", 100, 0.4, 0.6 );
259 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_mks", m_kskpi_4c_mks );
260 m_kskpi_4c_mksk = new TH1F( "kskpi_4c_mksk", "kskpi_4c_mksk", 100, 1.0, 3.0 );
261 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_mksk", m_kskpi_4c_mksk );
262 m_kskpi_4c_mkspi = new TH1F( "kskpi_4c_mkspi", "kskpi_4c_mkspi", 100, 0.6, 2.6 );
263 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_mkspi", m_kskpi_4c_mkspi );
264 m_kskpi_4c_mkpi = new TH1F( "kskpi_4c_mkpi", "kskpi_4c_mkpi", 100, 1.0, 3.0 );
265 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_mkpi", m_kskpi_4c_mkpi );
266 m_kskpi_4c_ks_px = new TH1F( "kskpi_4c_ks_px", "kskpi_4c_ks_px", 100, -1.5, 1.5 );
267 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_ks_px", m_kskpi_4c_ks_px );
268 m_kskpi_4c_ks_py = new TH1F( "kskpi_4c_ks_py", "kskpi_4c_ks_py", 100, -1.5, 1.5 );
269 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_ks_py", m_kskpi_4c_ks_py );
270 m_kskpi_4c_ks_pz = new TH1F( "kskpi_4c_ks_pz", "kskpi_4c_ks_pz", 100, -1.5, 1.5 );
271 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_ks_pz", m_kskpi_4c_ks_pz );
272 m_kskpi_4c_ks_p = new TH1F( "kskpi_4c_ks_p", "kskpi_4c_ks_p", 100, 0.0, 1.5 );
273 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_ks_p", m_kskpi_4c_ks_p );
274 m_kskpi_4c_ks_cos = new TH1F( "kskpi_4c_ks_cos", "kskpi_4c_ks_cos", 100, -1.0, 1.0 );
275 status = m_thistsvc->regHist( "/VAL/PHY/kskpi_4c_ks_cos", m_kskpi_4c_ks_cos );
276 }
277
278 NTuplePtr nt1( ntupleSvc(), "FILE1/signal" );
279 if ( nt1 ) m_tuple1 = nt1;
280 else
281 {
282 m_tuple1 = ntupleSvc()->book( "FILE1/signal", CLID_ColumnWiseTuple, "N-Tuple" );
283 if ( m_tuple1 )
284 {
285 status = m_tuple1->addItem( "irun", m_run );
286 status = m_tuple1->addItem( "ievent", m_event );
287 status = m_tuple1->addItem( "Nchrg", m_nchrg );
288 status = m_tuple1->addItem( "Nneu", m_nneu );
289 status = m_tuple1->addItem( "NGch", m_ngch, 0, 10 );
290
291 status = m_tuple1->addIndexedItem( "pidcode", m_ngch, m_pidcode );
292 status = m_tuple1->addIndexedItem( "pidprob", m_ngch, m_pidprob );
293 status = m_tuple1->addIndexedItem( "pidchiDedx", m_ngch, m_pidchiDedx );
294 status = m_tuple1->addIndexedItem( "pidchiTof1", m_ngch, m_pidchiTof1 );
295 status = m_tuple1->addIndexedItem( "pidchiTof2", m_ngch, m_pidchiTof2 );
296
297 status = m_tuple1->addItem( "npip", m_npip );
298 status = m_tuple1->addItem( "npim", m_npim );
299 status = m_tuple1->addItem( "nkp", m_nkp );
300 status = m_tuple1->addItem( "nkm", m_nkm );
301 status = m_tuple1->addItem( "np", m_np );
302 status = m_tuple1->addItem( "npb", m_npb );
303
304 status = m_tuple1->addItem( "vfits_chi", m_vfits_chi );
305 status = m_tuple1->addItem( "vfits_vx", m_vfits_vx );
306 status = m_tuple1->addItem( "vfits_vy", m_vfits_vy );
307 status = m_tuple1->addItem( "vfits_vz", m_vfits_vz );
308 status = m_tuple1->addItem( "vfits_vr", m_vfits_vr );
309
310 status = m_tuple1->addItem( "vfitp_chi", m_vfitp_chi );
311 status = m_tuple1->addItem( "vfitp_vx", m_vfitp_vx );
312 status = m_tuple1->addItem( "vfitp_vy", m_vfitp_vy );
313 status = m_tuple1->addItem( "vfitp_vz", m_vfitp_vz );
314 status = m_tuple1->addItem( "vfitp_vr", m_vfitp_vr );
315
316 status = m_tuple1->addItem( "vfit2_chi", m_vfit2_chi );
317 status = m_tuple1->addItem( "vfit2_mks", m_vfit2_mks );
318 status = m_tuple1->addItem( "vfit2_ct", m_vfit2_ct );
319 status = m_tuple1->addItem( "vfit2_dl", m_vfit2_dl );
320 status = m_tuple1->addItem( "vfit2_dle", m_vfit2_dle );
321
322 status = m_tuple1->addIndexedItem( "charge", m_ngch, m_charge );
323 status = m_tuple1->addIndexedItem( "vx0", m_ngch, m_vx0 );
324 status = m_tuple1->addIndexedItem( "vy0", m_ngch, m_vy0 );
325 status = m_tuple1->addIndexedItem( "vz0", m_ngch, m_vz0 );
326 status = m_tuple1->addIndexedItem( "vr0", m_ngch, m_vr0 );
327
328 status = m_tuple1->addIndexedItem( "vx", m_ngch, m_vx );
329 status = m_tuple1->addIndexedItem( "vy", m_ngch, m_vy );
330 status = m_tuple1->addIndexedItem( "vz", m_ngch, m_vz );
331 status = m_tuple1->addIndexedItem( "vr", m_ngch, m_vr );
332
333 status = m_tuple1->addIndexedItem( "px", m_ngch, m_px );
334 status = m_tuple1->addIndexedItem( "py", m_ngch, m_py );
335 status = m_tuple1->addIndexedItem( "pz", m_ngch, m_pz );
336 status = m_tuple1->addIndexedItem( "p", m_ngch, m_p );
337 status = m_tuple1->addIndexedItem( "cost", m_ngch, m_cost );
338
339 status = m_tuple1->addIndexedItem( "probPH", m_ngch, m_probPH );
340 status = m_tuple1->addIndexedItem( "normPH", m_ngch, m_normPH );
341 status = m_tuple1->addIndexedItem( "chie", m_ngch, m_chie );
342 status = m_tuple1->addIndexedItem( "chimu", m_ngch, m_chimu );
343 status = m_tuple1->addIndexedItem( "chipi", m_ngch, m_chipi );
344 status = m_tuple1->addIndexedItem( "chik", m_ngch, m_chik );
345 status = m_tuple1->addIndexedItem( "chip", m_ngch, m_chip );
346 status = m_tuple1->addIndexedItem( "ghit", m_ngch, m_ghit );
347 status = m_tuple1->addIndexedItem( "thit", m_ngch, m_thit );
348
349 status = m_tuple1->addIndexedItem( "e_emc", m_ngch, m_e_emc );
350
351 status = m_tuple1->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
352 status = m_tuple1->addIndexedItem( "tof_etof", m_ngch, m_tof_etof );
353 status = m_tuple1->addIndexedItem( "te_etof", m_ngch, m_te_etof );
354 status = m_tuple1->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
355 status = m_tuple1->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
356 status = m_tuple1->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
357 status = m_tuple1->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
358
359 status = m_tuple1->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
360 status = m_tuple1->addIndexedItem( "tof_btof1", m_ngch, m_tof_btof1 );
361 status = m_tuple1->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
362 status = m_tuple1->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
363 status = m_tuple1->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
364 status = m_tuple1->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
365 status = m_tuple1->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
366
367 status = m_tuple1->addIndexedItem( "qual_btof2", m_ngch, m_qual_btof2 );
368 status = m_tuple1->addIndexedItem( "tof_btof2", m_ngch, m_tof_btof2 );
369 status = m_tuple1->addIndexedItem( "te_btof2", m_ngch, m_te_btof2 );
370 status = m_tuple1->addIndexedItem( "tmu_btof2", m_ngch, m_tmu_btof2 );
371 status = m_tuple1->addIndexedItem( "tpi_btof2", m_ngch, m_tpi_btof2 );
372 status = m_tuple1->addIndexedItem( "tk_btof2", m_ngch, m_tk_btof2 );
373 status = m_tuple1->addIndexedItem( "tp_btof2", m_ngch, m_tp_btof2 );
374
375 status = m_tuple1->addItem( "chi2_fs4c", m_chi2_fs4c );
376 status = m_tuple1->addItem( "mks_fs4c", m_mks_fs4c );
377 status = m_tuple1->addItem( "mkspi_fs4c", m_mkspi_fs4c );
378 status = m_tuple1->addItem( "mksk_fs4c", m_mksk_fs4c );
379 status = m_tuple1->addItem( "mkpi_fs4c", m_mkpi_fs4c );
380
381 status = m_tuple1->addItem( "4c_chi2", m_4c_chi2 );
382 status = m_tuple1->addItem( "4c_mks", m_4c_mks );
383 status = m_tuple1->addItem( "4c_mkspi", m_4c_mkspi );
384 status = m_tuple1->addItem( "4c_mksk", m_4c_mksk );
385 status = m_tuple1->addItem( "4c_mkpi", m_4c_mkpi );
386 status = m_tuple1->addItem( "4c_ks_px", m_4c_ks_px );
387 status = m_tuple1->addItem( "4c_ks_py", m_4c_ks_py );
388 status = m_tuple1->addItem( "4c_ks_pz", m_4c_ks_pz );
389 status = m_tuple1->addItem( "4c_ks_p", m_4c_ks_p );
390 status = m_tuple1->addItem( "4c_ks_cos", m_4c_ks_cos );
391 }
392 else
393 {
394 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
395 return StatusCode::FAILURE;
396 }
397 }
398
399 //
400 //--------end of book--------
401 //
402
403 log << MSG::INFO << "successfully return from initialize()" << endmsg;
404 return StatusCode::SUCCESS;
405}
406
407// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
408// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
409StatusCode KsKpi::execute() {
410
411 MsgStream log( msgSvc(), name() );
412 log << MSG::INFO << "in execute()" << endmsg;
413
414 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
415 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
416 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
417
418 m_run = eventHeader->runNumber();
419 m_event = eventHeader->eventNumber();
420 m_nchrg = evtRecEvent->totalCharged();
421 m_nneu = evtRecEvent->totalNeutral();
422
423 log << MSG::INFO << "get event tag OK" << endmsg;
424 // cout <<"ncharg, nneu, tottks = "
425 // << evtRecEvent->totalCharged() << " , "
426 // << evtRecEvent->totalNeutral() << " , "
427 // << evtRecEvent->totalTracks() << endl ;
428
429 counter[0]++;
430
431 //
432 // Good charged track selection
433 //
434 Vint iGood;
435 iGood.clear();
436
437 int nCharge = 0;
438 Hep3Vector xorigin( 0, 0, 0 );
439
440 // if (m_reader.isRunNumberValid(runNo)) {
441 IVertexDbSvc* vtxsvc;
442 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
443 if ( vtxsvc->isVertexValid() )
444 {
445 double* dbv = vtxsvc->PrimaryVertex();
446 double* vv = vtxsvc->SigmaPrimaryVertex();
447 // HepVector dbv = m_reader.PrimaryVertex(runNo);
448 // HepVector vv = m_reader.SigmaPrimaryVertex(runNo);
449 xorigin.setX( dbv[0] );
450 xorigin.setY( dbv[1] );
451 xorigin.setZ( dbv[2] );
452 }
453
454 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
455 {
456 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
457 if ( !( *itTrk )->isMdcTrackValid() ) continue;
458 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
459 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
460 double x0 = mdcTrk->x();
461 double y0 = mdcTrk->y();
462 double z0 = mdcTrk->z();
463 double phi0 = mdcTrk->helix( 1 );
464 double xv = xorigin.x();
465 double yv = xorigin.y();
466 double zv = xorigin.z();
467 double rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
468 double cost = cos( mdcTrk->theta() );
469 if ( fabs( z0 - zv ) >= m_vz0cut ) continue;
470 if ( fabs( rv ) >= m_vr0cut ) continue;
471 if ( fabs( cost ) >= m_coscut ) continue;
472
473 iGood.push_back( ( *itTrk )->trackId() );
474 nCharge += mdcTrk->charge();
475 }
476
477 //
478 // Finish Good Charged Track Selection
479 //
480 m_ngch = iGood.size();
481 if ( ( m_ngch != 4 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
482 // std::cout << "ngood, totcharge = " << m_ngch << " , " << nCharge << std::endl;
483
484 counter[1]++;
485
486 //
487 // Particle ID
488 //
489 Vint ipip, ipim, ikp, ikm, ipp, ipm;
490 ipip.clear();
491 ipim.clear();
492 ikp.clear();
493 ikm.clear();
494 ipp.clear();
495 ipm.clear();
496
497 Vp4 p_pip, p_pim, p_kp, p_km, p_pp, p_pm;
498 p_pip.clear();
499 p_pim.clear();
500 p_kp.clear();
501 p_km.clear();
502 p_pp.clear();
503 p_pm.clear();
504
505 ParticleID* pid = ParticleID::instance();
506 for ( int i = 0; i < m_ngch; i++ )
507 {
508 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
509 // if(pid) delete pid;
510 pid->init();
511 pid->setMethod( pid->methodProbability() );
512 pid->setChiMinCut( 4 );
513 pid->setRecTrack( *itTrk );
514 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
515 // pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTofE() |
516 // pid->useTofQ()); // use PID sub-system
517 pid->identify( pid->onlyPionKaonProton() ); // seperater Pion/Kaon/Proton
518 // pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
519 // pid->identify(pid->onlyPion());
520 // pid->identify(pid->onlyKaon());
521 pid->calculate();
522 if ( !( pid->IsPidInfoValid() ) ) continue;
523
524 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
525 RecMdcKalTrack* mdcKalTrk = 0;
526 if ( ( *itTrk )->isMdcKalTrackValid() ) mdcKalTrk = ( *itTrk )->mdcKalTrack();
527
528 double prob_pi = pid->probPion();
529 double prob_K = pid->probKaon();
530 double prob_p = pid->probProton();
531 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
532 HepLorentzVector ptrk;
533 ptrk.setPx( mdcTrk->px() );
534 ptrk.setPy( mdcTrk->py() );
535 ptrk.setPz( mdcTrk->pz() );
536 double p3 = ptrk.mag();
537
538 if ( prob_pi > prob_K && prob_pi > prob_p )
539 {
540 m_pidcode[i] = 2;
541 m_pidprob[i] = pid->prob( 2 );
542 m_pidchiDedx[i] = pid->chiDedx( 2 );
543 m_pidchiTof1[i] = pid->chiTof1( 2 );
544 m_pidchiTof2[i] = pid->chiTof2( 2 );
545 ptrk.setE( sqrt( p3 * p3 + xmass[2] * xmass[2] ) );
546 if ( mdcTrk->charge() > 0 )
547 {
548 ipip.push_back( iGood[i] );
549 p_pip.push_back( ptrk );
550 }
551 if ( mdcTrk->charge() < 0 )
552 {
553 ipim.push_back( iGood[i] );
554 p_pim.push_back( ptrk );
555 }
556 }
557
558 if ( prob_K > prob_pi && prob_K > prob_p )
559 {
560 m_pidcode[i] = 3;
561 m_pidprob[i] = pid->prob( 3 );
562 m_pidchiDedx[i] = pid->chiDedx( 3 );
563 m_pidchiTof1[i] = pid->chiTof1( 3 );
564 m_pidchiTof2[i] = pid->chiTof2( 3 );
565 ptrk.setE( sqrt( p3 * p3 + xmass[3] * xmass[3] ) );
566 if ( mdcTrk->charge() > 0 )
567 {
568 ikp.push_back( iGood[i] );
569 p_kp.push_back( ptrk );
570 }
571 if ( mdcTrk->charge() < 0 )
572 {
573 ikm.push_back( iGood[i] );
574 p_km.push_back( ptrk );
575 }
576 }
577
578 if ( prob_p > prob_pi && prob_p > prob_K )
579 {
580 m_pidcode[i] = 4;
581 m_pidprob[i] = pid->prob( 4 );
582 m_pidchiDedx[i] = pid->chiDedx( 4 );
583 m_pidchiTof1[i] = pid->chiTof1( 4 );
584 m_pidchiTof2[i] = pid->chiTof2( 4 );
585 ptrk.setE( sqrt( p3 * p3 + xmass[4] * xmass[4] ) );
586 if ( mdcTrk->charge() > 0 )
587 {
588 ipp.push_back( iGood[i] );
589 p_pp.push_back( ptrk );
590 }
591 if ( mdcTrk->charge() < 0 )
592 {
593 ipm.push_back( iGood[i] );
594 p_pm.push_back( ptrk );
595 }
596 }
597 }
598 m_npip = ipip.size();
599 m_npim = ipim.size();
600 m_nkp = ikp.size();
601 m_nkm = ikm.size();
602 m_np = ipp.size();
603 m_npb = ipm.size();
604
605 if ( m_npip * m_npim != 2 ) { return StatusCode::SUCCESS; }
606 if ( m_nkp + m_nkm != 1 ) { return StatusCode::SUCCESS; }
607
608 counter[2]++;
609
610 //
611 // Test vertex fit
612 //
613
614 HepPoint3D vx( 0., 0., 0. );
615 HepSymMatrix Evx( 3, 0 );
616 double bx = 1E+6;
617 double by = 1E+6;
618 double bz = 1E+6;
619 Evx[0][0] = bx * bx;
620 Evx[1][1] = by * by;
621 Evx[2][2] = bz * bz;
622
623 VertexParameter vxpar;
624 vxpar.setVx( vx );
625 vxpar.setEvx( Evx );
626
627 // HepPoint3D bvx(0., 0., 0.);
628 // HepSymMatrix bEvx(3, 0);
629 // bEvx[0][0] = 0.038*0.038;;
630 // bEvx[1][1] = 0.00057*0.00057;
631 // bEvx[2][2] = 1.5*1.5;
632 // VertexParameter bs;
633 // bs.setVx(bvx);
634 // bs.setEvx(bEvx);
635
636 VertexFit* vtxfit_s = VertexFit::instance(); // S second vertex
637 VertexFit* vtxfit_p = VertexFit::instance(); // P primary vertex
638 SecondVertexFit* vtxfit2 = SecondVertexFit::instance();
639
640 RecMdcKalTrack *pi1KalTrk = nullptr, *pi2KalTrk = nullptr, *pi3KalTrk = nullptr, *k1KalTrk = nullptr;
641 RecMdcKalTrack *pipKalTrk = nullptr, *pimKalTrk = nullptr, *piKalTrk = nullptr, *kKalTrk = nullptr;
642 WTrackParameter wks, vwks;
643
644 double chi_temp = 999.0;
645 double mks_temp = 10.0;
646 bool okloop = false;
647 for ( unsigned int i1 = 0; i1 < m_npip; i1++ )
648 {
649 RecMdcKalTrack* pi1KalTrk = ( *( evtRecTrkCol->begin() + ipip[i1] ) )->mdcKalTrack();
650 pi1KalTrk->setPidType( RecMdcKalTrack::pion );
651 WTrackParameter wpi1trk( xmass[2], pi1KalTrk->getZHelix(), pi1KalTrk->getZError() );
652
653 for ( unsigned int i2 = 0; i2 < m_npim; i2++ )
654 {
655 RecMdcKalTrack* pi2KalTrk = ( *( evtRecTrkCol->begin() + ipim[i2] ) )->mdcKalTrack();
656 pi2KalTrk->setPidType( RecMdcKalTrack::pion );
657 WTrackParameter wpi2trk( xmass[2], pi2KalTrk->getZHelix(), pi2KalTrk->getZError() );
658
659 vtxfit_s->init();
660 vtxfit_s->AddTrack( 0, wpi1trk );
661 vtxfit_s->AddTrack( 1, wpi2trk );
662 vtxfit_s->AddVertex( 0, vxpar, 0, 1 );
663
664 if ( !( vtxfit_s->Fit( 0 ) ) ) continue;
665 vtxfit_s->BuildVirtualParticle( 0 );
666 m_vfits_chi = vtxfit_s->chisq( 0 );
667 WTrackParameter wkshort = vtxfit_s->wVirtualTrack( 0 );
668 VertexParameter vparks = vtxfit_s->vpar( 0 );
669
670 m_vfits_vx = ( vparks.Vx() )[0];
671 m_vfits_vy = ( vparks.Vx() )[1];
672 m_vfits_vz = ( vparks.Vx() )[2];
673 m_vfits_vr = sqrt( m_vfits_vx * m_vfits_vx + m_vfits_vy * m_vfits_vy );
674
675 if ( m_npip == 2 )
676 {
677 int j = i1;
678 int jj = ( i1 == 1 ) ? 0 : 1;
679 pi3KalTrk = ( *( evtRecTrkCol->begin() + ipip[jj] ) )->mdcKalTrack();
680 k1KalTrk = ( *( evtRecTrkCol->begin() + ikm[0] ) )->mdcKalTrack();
681 }
682 if ( m_npim == 2 )
683 {
684 int j = i2;
685 int jj = ( i2 == 1 ) ? 0 : 1;
686 pi3KalTrk = ( *( evtRecTrkCol->begin() + ipim[jj] ) )->mdcKalTrack();
687 k1KalTrk = ( *( evtRecTrkCol->begin() + ikp[0] ) )->mdcKalTrack();
688 }
689
690 pi3KalTrk->setPidType( RecMdcKalTrack::pion );
691 WTrackParameter wpi3trk( xmass[2], pi3KalTrk->getZHelix(), pi3KalTrk->getZError() );
692 k1KalTrk->setPidType( RecMdcKalTrack::kaon );
693 WTrackParameter wk1trk( xmass[3], k1KalTrk->getZHelixK(), k1KalTrk->getZErrorK() );
694
695 vtxfit_p->init();
696 // vtxfit_p->AddTrack(0, wkshort);
697 vtxfit_p->AddTrack( 0, wpi3trk );
698 vtxfit_p->AddTrack( 1, wk1trk );
699 vtxfit_p->AddVertex( 0, vxpar, 0, 1 );
700 if ( !( vtxfit_p->Fit( 0 ) ) ) continue;
701
702 m_vfitp_chi = vtxfit_p->chisq( 0 );
703
704 VertexParameter primaryVpar = vtxfit_p->vpar( 0 );
705 m_vfitp_vx = ( primaryVpar.Vx() )[0];
706 m_vfitp_vy = ( primaryVpar.Vx() )[1];
707 m_vfitp_vz = ( primaryVpar.Vx() )[2];
708 m_vfitp_vr = sqrt( m_vfitp_vx * m_vfitp_vx + m_vfitp_vy * m_vfitp_vy );
709
710 vtxfit2->init();
711 vtxfit2->setPrimaryVertex( vtxfit_p->vpar( 0 ) );
712 vtxfit2->AddTrack( 0, wkshort );
713 vtxfit2->setVpar( vtxfit_s->vpar( 0 ) );
714 if ( !vtxfit2->Fit() ) continue;
715
716 if ( fabs( ( ( vtxfit2->wpar() ).p() ).m() - mks0 ) < mks_temp )
717 {
718 mks_temp = fabs( ( ( vtxfit2->wpar() ).p() ).m() - mks0 );
719
720 okloop = true;
721
722 wks = vtxfit2->wpar();
723 m_vfit2_mks = ( wks.p() ).m();
724 m_vfit2_chi = vtxfit2->chisq();
725 m_vfit2_ct = vtxfit2->ctau();
726 m_vfit2_dl = vtxfit2->decayLength();
727 m_vfit2_dle = vtxfit2->decayLengthError();
728
729 pipKalTrk = pi1KalTrk;
730 pimKalTrk = pi2KalTrk;
731 piKalTrk = pi3KalTrk;
732 kKalTrk = k1KalTrk;
733 }
734 }
735 }
736
737 if ( !okloop ) { return StatusCode::SUCCESS; }
738
739 pipKalTrk->setPidType( RecMdcKalTrack::pion );
740 pimKalTrk->setPidType( RecMdcKalTrack::pion );
741 piKalTrk->setPidType( RecMdcKalTrack::pion );
742 kKalTrk->setPidType( RecMdcKalTrack::kaon );
743
744 WTrackParameter wpip( xmass[2], pipKalTrk->getZHelix(),
745 pipKalTrk->getZError() ); // pi+ from Ks
746 WTrackParameter wpim( xmass[2], pimKalTrk->getZHelix(),
747 pimKalTrk->getZError() ); // pi- from Ks
748 WTrackParameter wpi( xmass[2], piKalTrk->getZHelix(), piKalTrk->getZError() );
749 WTrackParameter wk( xmass[3], kKalTrk->getZHelixK(), kKalTrk->getZErrorK() );
750
751 counter[3]++;
752
753 //
754 // check good charged track's infomation
755 //
756 int ii = -1;
757 for ( int i = 0; i < m_ngch; i++ )
758 {
759
760 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
761 if ( !( *itTrk )->isMdcTrackValid() ) continue; // MDC information
762 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
763 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
764
765 if ( mdcKalTrk == pipKalTrk )
766 {
767 ii = 0;
768 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
769 }
770 if ( mdcKalTrk == pimKalTrk )
771 {
772 ii = 1;
773 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
774 }
775 if ( mdcKalTrk == piKalTrk )
776 {
777 ii = 2;
778 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
779 }
780 if ( mdcKalTrk == kKalTrk )
781 {
782 ii = 3;
783 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
784 }
785
786 m_charge[ii] = mdcTrk->charge();
787 double x0 = mdcTrk->x();
788 double y0 = mdcTrk->y();
789 double z0 = mdcTrk->z();
790 double phi0 = mdcTrk->helix( 1 );
791 double xv = xorigin.x();
792 double yv = xorigin.y();
793 double zv = xorigin.z();
794 double rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
795
796 m_vx0[ii] = x0 - xv;
797 m_vy0[ii] = y0 - yv;
798 m_vz0[ii] = z0 - zv;
799 m_vr0[ii] = rv;
800
801 x0 = mdcKalTrk->x();
802 y0 = mdcKalTrk->y();
803 z0 = mdcKalTrk->z();
804 rv = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
805
806 m_vx[ii] = x0 - xv;
807 m_vy[ii] = y0 - yv;
808 m_vz[ii] = z0 - zv;
809 m_vr[ii] = rv;
810
811 m_px[ii] = mdcKalTrk->px();
812 m_py[ii] = mdcKalTrk->py();
813 m_pz[ii] = mdcKalTrk->pz();
814 m_p[ii] = mdcKalTrk->p();
815 m_cost[ii] = mdcKalTrk->pz() / mdcKalTrk->p();
816
817 double ptrk = mdcKalTrk->p();
818 if ( ( *itTrk )->isMdcDedxValid() )
819 { // DEDX information
820 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
821 m_probPH[ii] = dedxTrk->probPH();
822 m_normPH[ii] = dedxTrk->normPH();
823
824 m_chie[ii] = dedxTrk->chiE();
825 m_chimu[ii] = dedxTrk->chiMu();
826 m_chipi[ii] = dedxTrk->chiPi();
827 m_chik[ii] = dedxTrk->chiK();
828 m_chip[ii] = dedxTrk->chiP();
829 m_ghit[ii] = dedxTrk->numGoodHits();
830 m_thit[ii] = dedxTrk->numTotalHits();
831 }
832
833 if ( ( *itTrk )->isEmcShowerValid() )
834 {
835 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
836 m_e_emc[ii] = emcTrk->energy();
837 }
838
839 if ( ( *itTrk )->isTofTrackValid() )
840 { // TOF information
841 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
842
843 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
844
845 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
846 {
847 TofHitStatus* status = new TofHitStatus;
848 status->setStatus( ( *iter_tof )->status() );
849
850 if ( !( status->is_barrel() ) )
851 { // endcap
852 if ( !( status->is_counter() ) ) continue; // ?
853 if ( status->layer() != 0 ) continue; // layer1
854 double path = ( *iter_tof )->path(); // ?
855 double tof = ( *iter_tof )->tof();
856 double ph = ( *iter_tof )->ph();
857 double rhit = ( *iter_tof )->zrhit();
858 double qual = 0.0 + ( *iter_tof )->quality();
859 double cntr = 0.0 + ( *iter_tof )->tofID();
860 double texp[5];
861 for ( int j = 0; j < 5; j++ )
862 {
863 double gb = ptrk / xmass[j];
864 double beta = gb / sqrt( 1 + gb * gb );
865 texp[j] = path / beta / velc;
866 }
867
868 m_qual_etof[ii] = qual;
869 m_tof_etof[ii] = tof;
870 }
871 else
872 { // barrel
873 if ( !( status->is_counter() ) ) continue; // ?
874 if ( status->layer() == 1 )
875 { // layer1
876 double path = ( *iter_tof )->path(); // ?
877 double tof = ( *iter_tof )->tof();
878 double ph = ( *iter_tof )->ph();
879 double rhit = ( *iter_tof )->zrhit();
880 double qual = 0.0 + ( *iter_tof )->quality();
881 double cntr = 0.0 + ( *iter_tof )->tofID();
882 double texp[5];
883 for ( int j = 0; j < 5; j++ )
884 {
885 double gb = ptrk / xmass[j];
886 double beta = gb / sqrt( 1 + gb * gb );
887 texp[j] = path / beta / velc;
888 }
889
890 m_qual_btof1[ii] = qual;
891 m_tof_btof1[ii] = tof;
892 }
893
894 if ( status->layer() == 2 )
895 { // layer2
896 double path = ( *iter_tof )->path(); // ?
897 double tof = ( *iter_tof )->tof();
898 double ph = ( *iter_tof )->ph();
899 double rhit = ( *iter_tof )->zrhit();
900 double qual = 0.0 + ( *iter_tof )->quality();
901 double cntr = 0.0 + ( *iter_tof )->tofID();
902 double texp[5];
903 for ( int j = 0; j < 5; j++ )
904 {
905 double gb = ptrk / xmass[j];
906 double beta = gb / sqrt( 1 + gb * gb );
907 texp[j] = path / beta / velc;
908 }
909
910 m_qual_btof2[ii] = qual;
911 m_tof_btof2[ii] = tof;
912 }
913 }
914 }
915 }
916 }
917 counter[4]++;
918
919 //////////
920 //
921 // Apply Kinematic 4C fit
922 //
923
924 KinematicFit* kmfit = KinematicFit::instance();
925
926 if ( m_test4C == 1 )
927 {
928 double ecms = 3.097;
929 double chisq = 9999.;
930 m_4c_chi2 = 9999.;
931 m_4c_mks = 10.0;
932 m_4c_mkspi = 10.0;
933 m_4c_mksk = 10.0;
934 m_4c_mkpi = 10.0;
935 m_4c_ks_px = 10.0;
936 m_4c_ks_py = 10.0;
937 m_4c_ks_pz = 10.0;
938 m_4c_ks_p = 10.0;
939 m_4c_ks_cos = 10.0;
940
941 kmfit->init();
942 kmfit->AddTrack( 0, wpi );
943 kmfit->AddTrack( 1, wk );
944 kmfit->AddTrack( 2, wks );
945 kmfit->AddFourMomentum( 0, p_cms );
946 bool oksq = kmfit->Fit();
947 if ( oksq )
948 {
949 chisq = kmfit->chisq();
950
951 HepLorentzVector pks = kmfit->pfit( 2 );
952 HepLorentzVector pkspi = kmfit->pfit( 0 ) + kmfit->pfit( 2 );
953 HepLorentzVector pksk = kmfit->pfit( 1 ) + kmfit->pfit( 2 );
954 HepLorentzVector pkpi = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
955
956 pks.boost( u_cms );
957 pkspi.boost( u_cms );
958 pksk.boost( u_cms );
959 pkpi.boost( u_cms );
960
961 m_4c_chi2 = chisq;
962 m_4c_mks = pks.m();
963 m_4c_mkspi = pkspi.m();
964 m_4c_mksk = pksk.m();
965 m_4c_mkpi = pkpi.m();
966
967 m_4c_ks_px = pks.px();
968 m_4c_ks_py = pks.py();
969 m_4c_ks_pz = pks.pz();
970 m_4c_ks_p = ( pks.vect() ).mag();
971 m_4c_ks_cos = m_4c_ks_pz / m_4c_ks_p;
972 }
973
974 chisq = 9999.;
975 m_chi2_fs4c = 9999.;
976 m_mks_fs4c = 10.0;
977 m_mkspi_fs4c = 10.0;
978 m_mksk_fs4c = 10.0;
979 m_mkpi_fs4c = 10.0;
980
981 kmfit->init();
982 kmfit->AddTrack( 0, wpip );
983 kmfit->AddTrack( 1, wpim );
984 kmfit->AddTrack( 2, wpi );
985 kmfit->AddTrack( 3, wk );
986 kmfit->AddFourMomentum( 0, p_cms );
987 oksq = kmfit->Fit();
988 if ( oksq )
989 {
990 chisq = kmfit->chisq();
991
992 HepLorentzVector pks = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
993 HepLorentzVector pkspi = pks + kmfit->pfit( 2 );
994 HepLorentzVector pksk = pks + kmfit->pfit( 3 );
995 HepLorentzVector pkpi = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
996
997 pks.boost( u_cms );
998 pkspi.boost( u_cms );
999 pksk.boost( u_cms );
1000 pkpi.boost( u_cms );
1001
1002 m_chi2_fs4c = chisq;
1003 m_mks_fs4c = pks.m();
1004 m_mkspi_fs4c = pkspi.m();
1005 m_mksk_fs4c = pksk.m();
1006 m_mkpi_fs4c = pkpi.m();
1007 }
1008 }
1009
1010 counter[5]++;
1011
1012 if ( m_tagKsKpi )
1013 {
1014 m_kskpi_vx_pi1->Fill( m_vx[0] );
1015 m_kskpi_vy_pi1->Fill( m_vy[0] );
1016 m_kskpi_vz_pi1->Fill( m_vz[0] );
1017 m_kskpi_vr_pi1->Fill( m_vr[0] );
1018 m_kskpi_px_pi1->Fill( m_px[0] );
1019 m_kskpi_py_pi1->Fill( m_py[0] );
1020 m_kskpi_pz_pi1->Fill( m_pz[0] );
1021 m_kskpi_pp_pi1->Fill( m_p[0] );
1022 m_kskpi_cos_pi1->Fill( m_cost[0] );
1023 m_kskpi_emc_pi1->Fill( m_e_emc[0] );
1024
1025 m_kskpi_vx_pi2->Fill( m_vx[1] );
1026 m_kskpi_vy_pi2->Fill( m_vy[1] );
1027 m_kskpi_vz_pi2->Fill( m_vz[1] );
1028 m_kskpi_vr_pi2->Fill( m_vr[1] );
1029 m_kskpi_px_pi2->Fill( m_px[1] );
1030 m_kskpi_py_pi2->Fill( m_py[1] );
1031 m_kskpi_pz_pi2->Fill( m_pz[1] );
1032 m_kskpi_pp_pi2->Fill( m_p[1] );
1033 m_kskpi_cos_pi2->Fill( m_cost[1] );
1034 m_kskpi_emc_pi2->Fill( m_e_emc[1] );
1035
1036 m_kskpi_vx_pi->Fill( m_vx[2] );
1037 m_kskpi_vy_pi->Fill( m_vy[2] );
1038 m_kskpi_vz_pi->Fill( m_vz[2] );
1039 m_kskpi_vr_pi->Fill( m_vr[2] );
1040 m_kskpi_px_pi->Fill( m_px[2] );
1041 m_kskpi_py_pi->Fill( m_py[2] );
1042 m_kskpi_pz_pi->Fill( m_pz[2] );
1043 m_kskpi_pp_pi->Fill( m_p[2] );
1044 m_kskpi_cos_pi->Fill( m_cost[2] );
1045 m_kskpi_emc_pi->Fill( m_e_emc[2] );
1046
1047 m_kskpi_vx_k->Fill( m_vx[3] );
1048 m_kskpi_vy_k->Fill( m_vy[3] );
1049 m_kskpi_vz_k->Fill( m_vz[3] );
1050 m_kskpi_vr_k->Fill( m_vr[3] );
1051 m_kskpi_px_k->Fill( m_px[3] );
1052 m_kskpi_py_k->Fill( m_py[3] );
1053 m_kskpi_pz_k->Fill( m_pz[3] );
1054 m_kskpi_pp_k->Fill( m_p[3] );
1055 m_kskpi_cos_k->Fill( m_cost[3] );
1056 m_kskpi_emc_k->Fill( m_e_emc[3] );
1057
1058 m_kskpi_pidchidedx_1->Fill( m_pidchiDedx[0] );
1059 m_kskpi_pidchitof1_1->Fill( m_pidchiTof1[0] );
1060 m_kskpi_pidchitof2_1->Fill( m_pidchiTof2[0] );
1061 m_kskpi_pidchidedx_2->Fill( m_pidchiDedx[1] );
1062 m_kskpi_pidchitof1_2->Fill( m_pidchiTof1[1] );
1063 m_kskpi_pidchitof2_2->Fill( m_pidchiTof2[1] );
1064 m_kskpi_pidchidedx_3->Fill( m_pidchiDedx[2] );
1065 m_kskpi_pidchitof1_3->Fill( m_pidchiTof1[2] );
1066 m_kskpi_pidchitof2_3->Fill( m_pidchiTof2[2] );
1067 m_kskpi_pidchidedx_4->Fill( m_pidchiDedx[3] );
1068 m_kskpi_pidchitof1_4->Fill( m_pidchiTof1[3] );
1069 m_kskpi_pidchitof2_4->Fill( m_pidchiTof2[3] );
1070
1071 m_kskpi_vfits_chi->Fill( m_vfits_chi );
1072 m_kskpi_vfits_vx->Fill( m_vfits_vx );
1073 m_kskpi_vfits_vy->Fill( m_vfits_vy );
1074 m_kskpi_vfits_vz->Fill( m_vfits_vz );
1075 m_kskpi_vfits_vr->Fill( m_vfits_vr );
1076
1077 m_kskpi_vfitp_chi->Fill( m_vfitp_chi );
1078 m_kskpi_vfitp_vx->Fill( m_vfitp_vx );
1079 m_kskpi_vfitp_vy->Fill( m_vfitp_vy );
1080 m_kskpi_vfitp_vz->Fill( m_vfitp_vz );
1081 m_kskpi_vfitp_vr->Fill( m_vfitp_vr );
1082
1083 m_kskpi_vfit2_chi->Fill( m_vfit2_chi );
1084 m_kskpi_vfit2_mks->Fill( m_vfit2_mks );
1085 m_kskpi_vfit2_ct->Fill( m_vfit2_ct );
1086 m_kskpi_vfit2_dl->Fill( m_vfit2_dl );
1087 m_kskpi_vfit2_dle->Fill( m_vfit2_dle );
1088
1089 m_kskpi_4c_ks_px->Fill( m_4c_ks_px );
1090 m_kskpi_4c_ks_py->Fill( m_4c_ks_py );
1091 m_kskpi_4c_ks_pz->Fill( m_4c_ks_pz );
1092 m_kskpi_4c_ks_p->Fill( m_4c_ks_p );
1093 m_kskpi_4c_ks_cos->Fill( m_4c_ks_cos );
1094
1095 m_kskpi_4c_chi->Fill( m_4c_chi2 );
1096 m_kskpi_4c_mks->Fill( m_4c_mks );
1097 m_kskpi_4c_mksk->Fill( m_4c_mksk );
1098 m_kskpi_4c_mkspi->Fill( m_4c_mkspi );
1099 m_kskpi_4c_mkpi->Fill( m_4c_mkpi );
1100 }
1101
1102 m_tuple1->write();
1103
1104 return StatusCode::SUCCESS;
1105}
1106
1107// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1108// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1109StatusCode KsKpi::finalize() {
1110
1111 MsgStream log( msgSvc(), name() );
1112 log << MSG::INFO << "in finalize()" << endmsg;
1113
1114 std::cout << "************ for KsKpi *******************" << std::endl;
1115 std::cout << "*******************************************" << std::endl;
1116 std::cout << "the total events is " << counter[0] << std::endl;
1117 std::cout << "Good charged tracks " << counter[1] << std::endl;
1118 std::cout << "particle ID " << counter[2] << std::endl;
1119 std::cout << "resort tracks considering Ks " << counter[3] << std::endl;
1120 std::cout << "info. for good charged track " << counter[4] << std::endl;
1121 std::cout << "kinematic fit 4C " << counter[5] << std::endl;
1122 std::cout << "kinematic fit 5C " << counter[6] << std::endl;
1123 std::cout << "*******************************************" << std::endl;
1124
1125 return StatusCode::SUCCESS;
1126}
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
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:18
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
KsKpi
Definition KsKpi.h:9
KsKpi::initialize
StatusCode initialize()
Definition KsKpi.cxx:82
KsKpi::KsKpi
KsKpi(const std::string &name, ISvcLocator *pSvcLocator)
Definition KsKpi.cxx:61
KsKpi::execute
StatusCode execute()
Definition KsKpi.cxx:409
KsKpi::finalize
StatusCode finalize()
Definition KsKpi.cxx:1109
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
SecondVertexFit::setPrimaryVertex
void setPrimaryVertex(const VertexParameter vpar)
Definition Analysis/VertexFit/include/VertexFit/SecondVertexFit.h:38
SecondVertexFit::instance
static SecondVertexFit * instance()
Definition SecondVertexFit.cxx:45
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::wVirtualTrack
WTrackParameter wVirtualTrack(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:104
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::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition VertexFit.cxx:600
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288
ecms
const double ecms
Definition inclkstar.cxx:26