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