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