BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQASelDimu.cxx
Go to the documentation of this file.
1#include "CLHEP/Vector/LorentzVector.h"
2#include "CLHEP/Vector/ThreeVector.h"
3#include "CLHEP/Vector/TwoVector.h"
4#include "GaudiKernel/Bootstrap.h"
5#include "GaudiKernel/INTupleSvc.h"
6#include "GaudiKernel/ISvcLocator.h"
7#include "GaudiKernel/MsgStream.h"
8#include "GaudiKernel/NTuple.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "TH1F.h"
11#include "TMath.h"
12#include <vector>
13
14#include "DstEvent/TofHitStatus.h"
15#include "EmcRecEventModel/RecEmcEventModel.h"
16#include "EmcRecGeoSvc/EmcRecBarrelGeo.h"
17#include "EmcRecGeoSvc/IEmcRecGeoSvc.h"
18#include "EvTimeEvent/RecEsTime.h"
19#include "EventModel/Event.h"
20#include "EventModel/EventHeader.h"
21#include "EventModel/EventModel.h"
22#include "EvtRecEvent/EvtRecEvent.h"
23#include "EvtRecEvent/EvtRecTrack.h"
24#include "McTruth/McParticle.h"
25#include "ParticleID/ParticleID.h"
26#include "VertexDbSvc/IVertexDbSvc.h"
27#include "VertexFit/Helix.h"
28#include "VertexFit/KinematicFit.h"
29#include "VertexFit/SecondVertexFit.h"
30#include "VertexFit/VertexFit.h"
31
32using CLHEP::Hep2Vector;
33using CLHEP::Hep3Vector;
34using CLHEP::HepLorentzVector;
35#include "CLHEP/Geometry/Point3D.h"
36#ifndef ENABLE_BACKWARDS_COMPATIBILITY
37typedef HepGeom::Point3D<double> HepPoint3D;
38#endif
39
40#include "DQASelDimu.h"
41
42#ifndef ENABLE_BACKWARDS_COMPATIBILITY
43typedef HepGeom::Point3D<double> HepPoint3D;
44#endif
45using CLHEP::HepLorentzVector;
46const double mpsi2s = 3.68609;
47const double mpi = 0.13957;
48const double mk = 0.493677;
49const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
50const double velc = 299.792458; // tof path unit in mm
51typedef std::vector<int> Vint;
52typedef std::vector<HepLorentzVector> Vp4;
53// declare one counter
54static int counter[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
55static int ndimu = 0;
56/////////////////////////////////////////////////////////////////////////////
57DECLARE_COMPONENT( DQASelDimu )
58DQASelDimu::DQASelDimu( const std::string& name, ISvcLocator* pSvcLocator )
59 : Algorithm( name, pSvcLocator ) {
60
61 // Declare the properties
62 declareProperty( "writentuple", m_writentuple = false );
63 declareProperty( "useVertexDB", m_useVertexDB = false );
64 declareProperty( "ecms", m_ecms = 3.097 );
65 declareProperty( "beamangle", m_beamangle = 0.022 );
66 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
67 declareProperty( "Vz0cut", m_vz0cut = 8.0 );
68 declareProperty( "Coscut", m_coscut = 0.93 );
69
70 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
71 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
72 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
73 declareProperty( "GammaTrkCut", m_gammaTrkCut = 20.0 );
74 declareProperty( "GammaTLCut", m_gammatlCut = 0 );
75 declareProperty( "GammaTHCut", m_gammathCut = 60 );
76
77 declareProperty( "acoll_mu_cut", m_acoll_mu_cut = 6. );
78 declareProperty( "acopl_mu_cut", m_acopl_mu_cut = 6. );
79 declareProperty( "poeb_mu_cut", m_poeb_mu_cut = 0.15 );
80 declareProperty( "dtof_mu_cut", m_dtof_mu_cut = 4. );
81 declareProperty( "eoeb_mu_cut", m_eoeb_mu_cut = 0.35 );
82 declareProperty( "etotal_mu_cut", m_etotal_mu_cut = 0.6 );
83 declareProperty( "tpoebh_mu_cut", m_tpoebh_mu_cut = 0.9 );
84 declareProperty( "tpoebl_mu_cut", m_tpoebl_mu_cut = 0.7 );
85 declareProperty( "tptotal_mu_cut", m_tptotal_mu_cut = 1.0 );
86
87 // normally, MDC+EMC, otherwise EMC only
88 declareProperty( "m_useEMConly", m_useEMConly = false );
89 declareProperty( "m_usePID", m_usePID = false ); // sub-system is under study
90 declareProperty( "m_useMDC", m_useMDC = true );
91 declareProperty( "m_useDEDX", m_useDEDX = false ); // not used
92 declareProperty( "m_useTOF", m_useTOF = false ); // sub-system is under study
93 declareProperty( "m_useEMC", m_useEMC = true );
94 declareProperty( "m_useMUC", m_useMUC = false ); // efficiency
95}
96
97// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
98StatusCode DQASelDimu::initialize() {
99 MsgStream log( msgSvc(), name() );
100
101 log << MSG::INFO << "in initialize()" << endmsg;
102 StatusCode status;
103 status = service( "THistSvc", m_thistsvc );
104 if ( status.isFailure() )
105 {
106 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
107 return status;
108 }
109
110 m_mumu_mass = new TH1F( "mumu_mass", "mumu_mass", 80, m_ecms - 0.3, m_ecms + 0.5 );
111 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_mass", m_mumu_mass );
112 m_mumu_acoll = new TH1F( "mumu_acoll", "mumu_acoll", 60, 0, 6 );
113 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_acoll", m_mumu_acoll );
114 m_mumu_eop_mup = new TH1F( "mumu_eop_mup", "mumu_eop_mup", 100, 0., 0.5 );
115 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_eop_mup", m_mumu_eop_mup );
116 m_mumu_eop_mum = new TH1F( "mumu_eop_mum", "mumu_eop_mum", 100, 0., 0.5 );
117 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_eop_mum", m_mumu_eop_mum );
118 m_mumu_costheta_mup = new TH1F( "mumu_costheta_mup", "mumu_costheta_mup", 100, -1, 1 );
119 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_costheta_mup", m_mumu_costheta_mup );
120 m_mumu_costheta_mum = new TH1F( "mumu_costheta_mum", "mumu_costheta_mum", 100, -1, 1 );
121 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_costheta_mum", m_mumu_costheta_mum );
122
123 m_mumu_phi_mup = new TH1F( "mumu_phi_mup", "mumu_phi_mup", 120, -3.2, 3.2 );
124 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_phi_mup", m_mumu_phi_mup );
125 m_mumu_phi_mum = new TH1F( "mumu_phi_mum", "mumu_phi_mum", 120, -3.2, 3.2 );
126 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_phi_mum", m_mumu_phi_mum );
127
128 m_mumu_nneu = new TH1F( "mumu_nneu", "mumu_nneu", 5, 0, 5 );
129 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_nneu", m_mumu_nneu );
130 m_mumu_nlay = new TH1F( "mumu_nlay", "mumu_nlay", 9, 0, 10 );
131 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_nlay", m_mumu_nlay );
132 m_mumu_nhit = new TH1F( "mumu_nhit", "mumu_nhit", 19, 0, 20 );
133 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_nhit", m_mumu_nhit );
134
135 m_mumu_eemc_mup = new TH1F( "mumu_eemc_mup", "mumu_eemc_mup", 100, 0.0, 1.0 );
136 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_eemc_mup", m_mumu_eemc_mup );
137 m_mumu_eemc_mum = new TH1F( "mumu_eemc_mum", "mumu_eemc_mum", 100, 0.0, 1.0 );
138 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_eemc_mum", m_mumu_eemc_mum );
139 m_mumu_x_mup = new TH1F( "mumu_x_mup", "mumu_x_mup", 100, -1.0, 1.0 );
140 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_x_mup", m_mumu_x_mup );
141 m_mumu_y_mup = new TH1F( "mumu_y_mup", "mumu_y_mup", 100, -1.0, 1.0 );
142 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_y_mup", m_mumu_y_mup );
143 m_mumu_z_mup = new TH1F( "mumu_z_mup", "mumu_z_mup", 100, -10.0, 10.0 );
144 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_z_mup", m_mumu_z_mup );
145 m_mumu_x_mum = new TH1F( "mumu_x_mum", "mumu_x_mum", 100, -1.0, 1.0 );
146 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_x_mum", m_mumu_x_mum );
147 m_mumu_y_mum = new TH1F( "mumu_y_mum", "mumu_y_mum", 100, -1.0, 1.0 );
148 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_y_mum", m_mumu_y_mum );
149 m_mumu_z_mum = new TH1F( "mumu_z_mum", "mumu_z_mum", 100, -10.0, 10.0 );
150 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_z_mum", m_mumu_z_mum );
151
152 m_mumu_px_mup = new TH1F( "mumu_px_mup", "mumu_px_mup", 200, -2.0, 2.0 );
153 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_px_mup", m_mumu_px_mup );
154 m_mumu_py_mup = new TH1F( "mumu_py_mup", "mumu_py_mup", 200, -2.0, 2.0 );
155 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_py_mup", m_mumu_py_mup );
156 m_mumu_pz_mup = new TH1F( "mumu_pz_mup", "mumu_pz_mup", 200, -2.0, 2.0 );
157 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pz_mup", m_mumu_pz_mup );
158 m_mumu_p_mup = new TH1F( "mumu_p_mup", "mumu_p_mup", 100, 1.0, 2.0 );
159 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_p_mup", m_mumu_p_mup );
160 m_mumu_px_mum = new TH1F( "mumu_px_mum", "mumu_px_mum", 100, -2.0, 2.0 );
161 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_px_mum", m_mumu_px_mum );
162 m_mumu_py_mum = new TH1F( "mumu_py_mum", "mumu_py_mum", 100, -2.0, 2.0 );
163 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_py_mum", m_mumu_py_mum );
164 m_mumu_pz_mum = new TH1F( "mumu_pz_mum", "mumu_pz_mum", 100, -2.0, 2.0 );
165 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pz_mum", m_mumu_pz_mum );
166 m_mumu_p_mum = new TH1F( "mumu_p_mum", "mumu_p_mum", 100, 1.0, 2.0 );
167 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_p_mum", m_mumu_p_mum );
168 m_mumu_deltatof = new TH1F( "mumu_deltatof", "mumu_deltatof", 50, 0.0, 10.0 );
169 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_deltatof", m_mumu_deltatof );
170
171 m_mumu_pidchidedx_mup = new TH1F( "mumu_pidchidedx_mup", "mumu_pidchidedx_mup", 160, -4, 4 );
172 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchidedx_mup", m_mumu_pidchidedx_mup );
173 m_mumu_pidchidedx_mum = new TH1F( "mumu_pidchidedx_mum", "mumu_pidchidedx_mum", 160, -4, 4 );
174 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchidedx_mum", m_mumu_pidchidedx_mum );
175 m_mumu_pidchitof1_mup = new TH1F( "mumu_pidchitof1_mup", "mumu_pidchitof1_mup", 160, -4, 4 );
176 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchitof1_mup", m_mumu_pidchitof1_mup );
177 m_mumu_pidchitof1_mum = new TH1F( "mumu_pidchitof1_mum", "mumu_pidchitof1_mum", 160, -4, 4 );
178 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchitof1_mum", m_mumu_pidchitof1_mum );
179 m_mumu_pidchitof2_mup = new TH1F( "mumu_pidchitof2_mup", "mumu_pidchitof2_mup", 160, -4, 4 );
180 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchitof2_mup", m_mumu_pidchitof2_mup );
181 m_mumu_pidchitof2_mum = new TH1F( "mumu_pidchitof2_mum", "mumu_pidchitof2_mum", 160, -4, 4 );
182 status = m_thistsvc->regHist( "/DQAHist/Dimu/mumu_pidchitof2_mum", m_mumu_pidchitof2_mum );
183
184 NTuplePtr nt1( ntupleSvc(), "DQAFILE/Dimu" );
185 if ( nt1 ) m_tuple1 = nt1;
186 else
187 {
188 m_tuple1 = ntupleSvc()->book( "DQAFILE/Dimu", CLID_ColumnWiseTuple, "N-Tuple" );
189 if ( m_tuple1 )
190 {
191 status = m_tuple1->addItem( "run", m_run );
192 status = m_tuple1->addItem( "rec", m_rec );
193 status = m_tuple1->addItem( "Nchrg", m_ncharg );
194 status = m_tuple1->addItem( "Nneu", m_nneu, 0, 40 );
195 status = m_tuple1->addItem( "NGch", m_ngch, 0, 40 );
196 status = m_tuple1->addItem( "NGam", m_nGam );
197
198 status = m_tuple1->addItem( "dimutag", m_dimutag );
199
200 status = m_tuple1->addItem( "acoll", m_acoll );
201 status = m_tuple1->addItem( "acopl", m_acopl );
202 status = m_tuple1->addItem( "deltatof", m_deltatof );
203 status = m_tuple1->addItem( "eop1", m_eop1 );
204 status = m_tuple1->addItem( "eop2", m_eop2 );
205 status = m_tuple1->addItem( "eoeb1", m_eoeb1 );
206 status = m_tuple1->addItem( "eoeb2", m_eoeb2 );
207 status = m_tuple1->addItem( "poeb1", m_poeb1 );
208 status = m_tuple1->addItem( "poeb2", m_poeb2 );
209 status = m_tuple1->addItem( "etoeb1", m_etoeb1 );
210 status = m_tuple1->addItem( "etoeb2", m_etoeb2 );
211 status = m_tuple1->addItem( "mucinfo1", m_mucinfo1 );
212 status = m_tuple1->addItem( "mucinfo2", m_mucinfo2 );
213
214 status = m_tuple1->addIndexedItem( "delang", m_nneu, m_delang );
215 status = m_tuple1->addIndexedItem( "delphi", m_nneu, m_delphi );
216 status = m_tuple1->addIndexedItem( "delthe", m_nneu, m_delthe );
217 status = m_tuple1->addIndexedItem( "npart", m_nneu, m_npart );
218 status = m_tuple1->addIndexedItem( "nemchits", m_nneu, m_nemchits );
219 status = m_tuple1->addIndexedItem( "module", m_nneu, m_module );
220 status = m_tuple1->addIndexedItem( "x", m_nneu, m_x );
221 status = m_tuple1->addIndexedItem( "y", m_nneu, m_y );
222 status = m_tuple1->addIndexedItem( "z", m_nneu, m_z );
223 status = m_tuple1->addIndexedItem( "px", m_nneu, m_px );
224 status = m_tuple1->addIndexedItem( "py", m_nneu, m_py );
225 status = m_tuple1->addIndexedItem( "pz", m_nneu, m_pz );
226 status = m_tuple1->addIndexedItem( "theta", m_nneu, m_theta );
227 status = m_tuple1->addIndexedItem( "phi", m_nneu, m_phi );
228 status = m_tuple1->addIndexedItem( "dx", m_nneu, m_dx );
229 status = m_tuple1->addIndexedItem( "dy", m_nneu, m_dy );
230 status = m_tuple1->addIndexedItem( "dz", m_nneu, m_dz );
231 status = m_tuple1->addIndexedItem( "dtheta", m_nneu, m_dtheta );
232 status = m_tuple1->addIndexedItem( "dphi", m_nneu, m_dphi );
233 status = m_tuple1->addIndexedItem( "energy", m_nneu, m_energy );
234 status = m_tuple1->addIndexedItem( "dE", m_nneu, m_dE );
235 status = m_tuple1->addIndexedItem( "eSeed", m_nneu, m_eSeed );
236 status = m_tuple1->addIndexedItem( "nSeed", m_nneu, m_nSeed );
237 status = m_tuple1->addIndexedItem( "e3x3", m_nneu, m_e3x3 );
238 status = m_tuple1->addIndexedItem( "e5x5", m_nneu, m_e5x5 );
239 status = m_tuple1->addIndexedItem( "secondMoment", m_nneu, m_secondMoment );
240 status = m_tuple1->addIndexedItem( "latMoment", m_nneu, m_latMoment );
241 status = m_tuple1->addIndexedItem( "a20Moment", m_nneu, m_a20Moment );
242 status = m_tuple1->addIndexedItem( "a42Moment", m_nneu, m_a42Moment );
243 status = m_tuple1->addIndexedItem( "getTime", m_nneu, m_getTime );
244 status = m_tuple1->addIndexedItem( "getEAll", m_nneu, m_getEAll );
245
246 status = m_tuple1->addIndexedItem( "charge", m_ngch, m_charge );
247 status = m_tuple1->addIndexedItem( "vx", m_ngch, m_vx0 );
248 status = m_tuple1->addIndexedItem( "vy", m_ngch, m_vy0 );
249 status = m_tuple1->addIndexedItem( "vz", m_ngch, m_vz0 );
250
251 status = m_tuple1->addIndexedItem( "px", m_ngch, m_px );
252 status = m_tuple1->addIndexedItem( "py", m_ngch, m_py );
253 status = m_tuple1->addIndexedItem( "pz", m_ngch, m_pz );
254 status = m_tuple1->addIndexedItem( "p", m_ngch, m_p );
255
256 status = m_tuple1->addIndexedItem( "kal_vx", m_ngch, m_kal_vx0 );
257 status = m_tuple1->addIndexedItem( "kal_vy", m_ngch, m_kal_vy0 );
258 status = m_tuple1->addIndexedItem( "kal_vz", m_ngch, m_kal_vz0 );
259
260 status = m_tuple1->addIndexedItem( "kal_px", m_ngch, m_kal_px );
261 status = m_tuple1->addIndexedItem( "kal_py", m_ngch, m_kal_py );
262 status = m_tuple1->addIndexedItem( "kal_pz", m_ngch, m_kal_pz );
263 status = m_tuple1->addIndexedItem( "kal_p", m_ngch, m_kal_p );
264
265 status = m_tuple1->addIndexedItem( "probPH", m_ngch, m_probPH );
266 status = m_tuple1->addIndexedItem( "normPH", m_ngch, m_normPH );
267 status = m_tuple1->addIndexedItem( "chie", m_ngch, m_chie );
268 status = m_tuple1->addIndexedItem( "chimu", m_ngch, m_chimu );
269 status = m_tuple1->addIndexedItem( "chipi", m_ngch, m_chipi );
270 status = m_tuple1->addIndexedItem( "chik", m_ngch, m_chik );
271 status = m_tuple1->addIndexedItem( "chip", m_ngch, m_chip );
272 status = m_tuple1->addIndexedItem( "ghit", m_ngch, m_ghit );
273 status = m_tuple1->addIndexedItem( "thit", m_ngch, m_thit );
274
275 status = m_tuple1->addIndexedItem( "e_emc", m_ngch, m_e_emc );
276 status = m_tuple1->addIndexedItem( "phi_emc", m_ngch, m_phi_emc );
277 status = m_tuple1->addIndexedItem( "theta_emc", m_ngch, m_theta_emc );
278
279 status = m_tuple1->addIndexedItem( "nhit_muc", m_ngch, m_nhit_muc );
280 status = m_tuple1->addIndexedItem( "nlay_muc", m_ngch, m_nlay_muc );
281 status = m_tuple1->addIndexedItem( "t_btof", m_ngch, m_t_btof );
282 status = m_tuple1->addIndexedItem( "t_etof", m_ngch, m_t_etof );
283 status = m_tuple1->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
284 status = m_tuple1->addIndexedItem( "tof_etof", m_ngch, m_tof_etof );
285 status = m_tuple1->addIndexedItem( "te_etof", m_ngch, m_te_etof );
286 status = m_tuple1->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
287 status = m_tuple1->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
288 status = m_tuple1->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
289 status = m_tuple1->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
290
291 status = m_tuple1->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
292 status = m_tuple1->addIndexedItem( "tof_btof1", m_ngch, m_tof_btof1 );
293 status = m_tuple1->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
294 status = m_tuple1->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
295 status = m_tuple1->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
296 status = m_tuple1->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
297 status = m_tuple1->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
298
299 status = m_tuple1->addIndexedItem( "qual_btof2", m_ngch, m_qual_btof2 );
300 status = m_tuple1->addIndexedItem( "tof_btof2", m_ngch, m_tof_btof2 );
301 status = m_tuple1->addIndexedItem( "te_btof2", m_ngch, m_te_btof2 );
302 status = m_tuple1->addIndexedItem( "tmu_btof2", m_ngch, m_tmu_btof2 );
303 status = m_tuple1->addIndexedItem( "tpi_btof2", m_ngch, m_tpi_btof2 );
304 status = m_tuple1->addIndexedItem( "tk_btof2", m_ngch, m_tk_btof2 );
305 status = m_tuple1->addIndexedItem( "tp_btof2", m_ngch, m_tp_btof2 );
306 status = m_tuple1->addIndexedItem( "pidcode", m_ngch, m_pidcode );
307 status = m_tuple1->addIndexedItem( "pidprob", m_ngch, m_pidprob );
308 status = m_tuple1->addIndexedItem( "pidchiDedx", m_ngch, m_pidchiDedx );
309 status = m_tuple1->addIndexedItem( "pidchiTof1", m_ngch, m_pidchiTof1 );
310 status = m_tuple1->addIndexedItem( "pidchiTof2", m_ngch, m_pidchiTof2 );
311
312 status = m_tuple1->addItem( "px_cms_ep", m_px_cms_ep ); // momentum of electron+
313 status = m_tuple1->addItem( "py_cms_ep", m_py_cms_ep ); // momentum of electron+
314 status = m_tuple1->addItem( "pz_cms_ep", m_pz_cms_ep ); // momentum of electron+
315 status = m_tuple1->addItem( "e_cms_ep", m_e_cms_ep ); // momentum of electron+
316 status = m_tuple1->addItem( "cos_ep", m_cos_ep ); // momentum of electron+
317 status = m_tuple1->addItem( "px_cms_em", m_px_cms_em ); // momentum of electron-
318 status = m_tuple1->addItem( "py_cms_em", m_py_cms_em ); // momentum of electron-
319 status = m_tuple1->addItem( "pz_cms_em", m_pz_cms_em ); // momentum of electron-
320 status = m_tuple1->addItem( "e_cms_em", m_e_cms_em ); // momentum of electron-
321 status = m_tuple1->addItem( "cos_em", m_cos_em ); // momentum of electron-
322 status = m_tuple1->addItem( "mass_ee", m_mass_ee ); //
323 status = m_tuple1->addItem( "emax", m_emax ); //
324 status = m_tuple1->addItem( "esum", m_esum ); //
325 status = m_tuple1->addItem( "npip", m_npip );
326 status = m_tuple1->addItem( "npim", m_npim );
327 status = m_tuple1->addItem( "nkp", m_nkp );
328 status = m_tuple1->addItem( "nkm", m_nkm );
329 status = m_tuple1->addItem( "np", m_np );
330 status = m_tuple1->addItem( "npb", m_npb );
331
332 status = m_tuple1->addItem( "nep", m_nep );
333 status = m_tuple1->addItem( "nem", m_nem );
334 status = m_tuple1->addItem( "nmup", m_nmup );
335 status = m_tuple1->addItem( "nmum", m_nmum );
336 }
337 else
338 {
339 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
340 return StatusCode::FAILURE;
341 }
342 }
343
344 //
345 //--------end of book--------
346 //
347
348 log << MSG::INFO << "successfully return from initialize()" << endmsg;
349 return StatusCode::SUCCESS;
350}
351
352// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
353StatusCode DQASelDimu::execute() {
354 setFilterPassed( false );
355 const double beamEnergy = m_ecms / 2.;
356 const HepLorentzVector p_cms( m_ecms * sin( m_beamangle * 0.5 ), 0.0, 0.0, m_ecms );
357 const Hep3Vector u_cms = -p_cms.boostVector();
358 MsgStream log( msgSvc(), name() );
359 log << MSG::INFO << "in execute()" << endmsg;
360
361 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
362 if ( !eventHeader )
363 {
364 log << MSG::FATAL << "Could not find Event Header" << endmsg;
365 return StatusCode::SUCCESS;
366 }
367
368 m_run = eventHeader->runNumber();
369 m_rec = eventHeader->eventNumber();
370
371 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
372 if ( !evtRecEvent )
373 {
374 log << MSG::FATAL << "Could not find EvtRecEvent" << endmsg;
375 return StatusCode::SUCCESS;
376 }
377 log << MSG::INFO << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
378 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
379 // if(evtRecEvent->totalNeutral()>30)return sc;
380 m_ncharg = evtRecEvent->totalCharged();
381
382 m_nneu = evtRecEvent->totalNeutral();
383
384 HepPoint3D vx( 0., 0., 0. );
385 HepSymMatrix Evx( 3, 0 );
386 if ( m_useVertexDB )
387 {
388 IVertexDbSvc* vtxsvc;
389 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
390 if ( vtxsvc->isVertexValid() )
391 {
392 double* dbv = vtxsvc->PrimaryVertex();
393 double* vv = vtxsvc->SigmaPrimaryVertex();
394 // if (m_reader.isRunNumberValid( m_run)) {
395 // HepVector dbv = m_reader.PrimaryVertex( m_run);
396 // HepVector vv = m_reader.SigmaPrimaryVertex( m_run);
397 vx.setX( dbv[0] );
398 vx.setY( dbv[1] );
399 vx.setZ( dbv[2] );
400 Evx[0][0] = vv[0] * vv[0];
401 Evx[0][1] = vv[0] * vv[1];
402 Evx[1][1] = vv[1] * vv[1];
403 Evx[1][2] = vv[1] * vv[2];
404 Evx[2][2] = vv[2] * vv[2];
405 }
406 }
407
408 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
409 if ( !evtRecTrkCol )
410 {
411 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endmsg;
412 return StatusCode::SUCCESS;
413 }
414 Vint iGood;
415 iGood.clear();
416
417 int nCharge = 0;
418
419 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
420 {
421 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
422 if ( !( *itTrk )->isMdcTrackValid() ) continue;
423 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
424
425 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
426 double pch = mdcTrk->p();
427 double x0 = mdcTrk->x();
428 double y0 = mdcTrk->y();
429 double z0 = mdcTrk->z();
430 double phi0 = mdcTrk->helix( 1 );
431 double xv = vx.x();
432 double yv = vx.y();
433 double zv = vx.z();
434 double Rxy = ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 );
435 double m_vx0 = x0;
436 double m_vy0 = y0;
437 double m_vz0 = z0;
438 double m_vr0 = Rxy;
439 if ( fabs( z0 ) >= m_vz0cut ) continue;
440 if ( fabs( Rxy ) >= m_vr0cut ) continue;
441
442 if ( fabs( m_vz0 ) >= m_vz0cut ) continue;
443 if ( m_vr0 >= m_vr0cut ) continue;
444
445 // double cost = cos(mdcTrk->theta());
446 // if(fabs(cost) >= m_coscut ) continue;
447 // iGood.push_back((*itTrk)->trackId());
448 iGood.push_back( i );
449 nCharge += mdcTrk->charge();
450 }
451
452 //
453 // Finish Good Charged Track Selection
454 //
455 int nGood = iGood.size();
456 m_ngch = nGood;
457 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
458
459 if ( ( nGood != 2 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
460
461 counter[1]++;
462
463 //
464 // Particle ID
465 //
466 Vint ipip, ipim, iep, iem, imup, imum;
467 ipip.clear();
468 ipim.clear();
469 iep.clear();
470 iem.clear();
471 imup.clear();
472 imum.clear();
473
474 ParticleID* pid = ParticleID::instance();
475 for ( int i = 0; i < m_ngch; i++ )
476 {
477 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
478 // if(pid) delete pid;
479 pid->init();
480 pid->setMethod( pid->methodProbability() );
481 pid->setChiMinCut( 4 );
482 pid->setRecTrack( *itTrk );
483 pid->usePidSys( pid->useDedx() | pid->useTof1() |
484 pid->useTof2() ); //|pid->useEmc()|pid->useMuc());
485 //// use PID sub-system
486 pid->identify( pid->onlyElectron() | pid->onlyMuon() |
487 pid->onlyPion() ); // seperater Pion/Kaon/Proton
488 pid->calculate();
489 if ( !( pid->IsPidInfoValid() ) ) continue;
490 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
491 /// RecMdcKalTrack* mdcKalTrk = 0 ;
492 /// if((*itTrk)->isMdcKalTrackValid()) mdcKalTrk = (*itTrk)->mdcKalTrack();
493 double prob_pi = pid->probPion();
494 double prob_K = pid->probKaon();
495 double prob_p = pid->probProton();
496 double prob_e = pid->probElectron();
497 double prob_mu = pid->probMuon();
498 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
499 HepLorentzVector ptrk;
500 ptrk.setPx( mdcTrk->px() );
501 ptrk.setPy( mdcTrk->py() );
502 ptrk.setPz( mdcTrk->pz() );
503 double p3 = ptrk.mag();
504
505 m_pidcode[i] = 1;
506 m_pidprob[i] = pid->prob( 1 );
507 m_pidchiDedx[i] = pid->chiDedx( 1 );
508 m_pidchiTof1[i] = pid->chiTof1( 1 );
509 m_pidchiTof2[i] = pid->chiTof2( 1 );
510 if ( mdcTrk->charge() > 0 ) { imup.push_back( iGood[i] ); }
511 if ( mdcTrk->charge() < 0 ) { imum.push_back( iGood[i] ); }
512 }
513
514 m_nep = iep.size();
515 m_nem = iem.size();
516 m_nmup = imup.size();
517 m_nmum = imum.size();
518
519 counter[2]++;
520
521 //
522 // Good neutral track selection
523 //
524 Vint iGam;
525 iGam.clear();
526 int iphoton = 0;
527 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
528 {
529 if ( i >= evtRecTrkCol->size() ) break;
530 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
531 if ( !( *itTrk )->isEmcShowerValid() ) continue;
532 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
533 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
534
535 RecEmcID showerId = emcTrk->getShowerId();
536 unsigned int npart = EmcID::barrel_ec( showerId );
537 int n = emcTrk->numHits();
538 int module = emcTrk->module();
539 double x = emcTrk->x();
540 double y = emcTrk->y();
541 double z = emcTrk->z();
542 double dx = emcTrk->dx();
543 double dy = emcTrk->dy();
544 double dth = emcTrk->dtheta();
545 double dph = emcTrk->dphi();
546 double dz = emcTrk->dz();
547 double energy = emcTrk->energy();
548 double dE = emcTrk->dE();
549 double eSeed = emcTrk->eSeed();
550 double e3x3 = emcTrk->e3x3();
551 double e5x5 = emcTrk->e5x5();
552 double secondMoment = emcTrk->secondMoment();
553 double latMoment = emcTrk->latMoment();
554 double getTime = emcTrk->time();
555 double getEAll = emcTrk->getEAll();
556 double a20Moment = emcTrk->a20Moment();
557 double a42Moment = emcTrk->a42Moment();
558 // int phigap=emcTrk->PhiGap();
559 // int thetagap=emcTrk->ThetaGap();
560 // double getETof2x1 = emcTrk->getETof2x1();
561 // double getETof2x3 = emcTrk->getETof2x3();
562 // double getELepton = emcTrk->getELepton();
563 double nseed = 0; //(emcTrk->getCluster() )->getSeedSize() ;
564 HepPoint3D EmcPos( x, y, z );
565 m_nemchits[iphoton] = n;
566 m_npart[iphoton] = npart;
567 m_module[iphoton] = module;
568 m_theta[iphoton] = EmcPos.theta();
569 m_phi[iphoton] = EmcPos.phi();
570 m_x[iphoton] = x;
571 m_y[iphoton] = y;
572 m_z[iphoton] = z;
573 m_dx[iphoton] = dx;
574 m_dy[iphoton] = dy;
575 m_dz[iphoton] = dz;
576 m_dtheta[iphoton] = dth;
577 m_dphi[iphoton] = dph;
578 m_energy[iphoton] = energy;
579 m_dE[iphoton] = dE;
580 m_eSeed[iphoton] = eSeed;
581 m_nSeed[iphoton] = nseed;
582 m_e3x3[iphoton] = e3x3;
583 m_e5x5[iphoton] = e5x5;
584 m_secondMoment[iphoton] = secondMoment;
585 m_latMoment[iphoton] = latMoment;
586 m_getTime[iphoton] = getTime;
587 m_getEAll[iphoton] = getEAll;
588 m_a20Moment[iphoton] = a20Moment;
589 m_a42Moment[iphoton] = a42Moment;
590
591 // m_getELepton[iphoton]=getELepton;
592 // m_getETof2x1[iphoton]=getETof2x1;
593 // m_getETof2x3[iphoton]=getETof2x3;
594 // m_PhiGap[iphoton]=phigap;
595 // m_ThetaGap[iphoton]=thetagap;
596 double dthe = 200.;
597 double dphi = 200.;
598 double dang = 200.;
599
600 // find the nearest charged track
601 for ( int j = 0; j < nGood; j++ )
602 {
603
604 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + iGood[j];
605 if ( !( *jtTrk )->isMdcTrackValid() ) continue;
606 RecMdcTrack* jtmdcTrk = ( *jtTrk )->mdcTrack();
607 double jtcharge = jtmdcTrk->charge();
608 if ( !( *jtTrk )->isExtTrackValid() ) continue;
609 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
610 if ( extTrk->emcVolumeNumber() == -1 ) continue;
611 Hep3Vector extpos = extTrk->emcPosition();
612 // double ctht = extpos.cosTheta(emcpos);
613 double angd = extpos.angle( emcpos );
614 double thed = extpos.theta() - emcpos.theta();
615 double phid = extpos.deltaPhi( emcpos );
616 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
617 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
618
619 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
620 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
621 if ( angd < dang ) dang = angd;
622 }
623
624 //
625 // good photon cut will be set here
626 //
627
628 dthe = dthe * 180 / ( CLHEP::pi );
629 dphi = dphi * 180 / ( CLHEP::pi );
630 dang = dang * 180 / ( CLHEP::pi );
631 double eraw = emcTrk->energy();
632 double phi = emcTrk->phi();
633 double the = emcTrk->theta();
634
635 m_delphi[iphoton] = dphi;
636 m_delthe[iphoton] = dthe;
637 m_delang[iphoton] = dang;
638 if ( energy < m_energyThreshold ) continue;
639 if ( getTime > m_gammathCut || getTime < m_gammatlCut ) continue;
640 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
641 if ( dang < m_gammaTrkCut ) continue;
642 iphoton++;
643 iGam.push_back( i );
644 if ( iphoton >= 40 ) return StatusCode::SUCCESS;
645 }
646
647 int nGam = iGam.size();
648 m_nGam = nGam;
649 // std::cout << "num Good Photon " << m_nGam << " , "
650 // <<evtRecEvent->totalNeutral()<<std::endl;
651 // std::cout<<"dbg_4"<<std::endl;
652 counter[3]++;
653
654 double egam_ext = 0;
655 double ex_gam = 0;
656 double ey_gam = 0;
657 double ez_gam = 0;
658 double et_gam = 0;
659 double e_gam = 0;
660 for ( int i = 0; i < m_nGam; i++ )
661 {
662 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
663 if ( !( *itTrk )->isEmcShowerValid() ) continue;
664 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
665 double eraw = emcTrk->energy();
666 double phi = emcTrk->phi();
667 double the = emcTrk->theta();
668 HepLorentzVector ptrk;
669 ex_gam += eraw * sin( the ) * cos( phi );
670 ey_gam += eraw * sin( the ) * sin( phi );
671 ez_gam += eraw * cos( the );
672 et_gam += eraw * sin( the );
673 e_gam += eraw;
674 if ( eraw >= egam_ext ) { egam_ext = eraw; }
675 }
676
677 double px_had = 0;
678 double py_had = 0;
679 double pz_had = 0;
680 double pt_had = 0;
681 double t_pxy2 = 0;
682 double p_had = 0;
683 double e_had = 0;
684 //
685 // check good charged track's infomation
686 //
687 int ii;
688 m_e_emc[0] = -0.1;
689 m_e_emc[1] = -0.1;
690 for ( int i = 0; i < m_ngch; i++ )
691 {
692
693 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
694
695 if ( !( *itTrk )->isMdcTrackValid() ) continue; // MDC information
696 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
697 // if(!(*itTrk)->isEmcShowerValid()) return StatusCode::SUCCESS;///dbg
698 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
699 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
700 ii = i;
701
702 m_charge[ii] = mdcTrk->charge();
703 m_vx0[ii] = mdcTrk->x();
704 m_vy0[ii] = mdcTrk->y();
705 m_vz0[ii] = mdcTrk->z();
706
707 m_px[ii] = mdcTrk->px();
708 m_py[ii] = mdcTrk->py();
709 m_pz[ii] = mdcTrk->pz();
710 m_p[ii] = mdcTrk->p();
711
712 mdcKalTrk->setPidType( RecMdcKalTrack::muon );
713
714 /// if(m_pidcode[ii]==3)mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
715 /// if(m_pidcode[ii]==4)mdcKalTrk->setPidType(RecMdcKalTrack::proton);
716 m_kal_vx0[ii] = mdcKalTrk->x();
717 m_kal_vy0[ii] = mdcKalTrk->y();
718 m_kal_vz0[ii] = mdcKalTrk->z();
719
720 m_kal_px[ii] = mdcKalTrk->px();
721 m_kal_py[ii] = mdcKalTrk->py();
722 m_kal_pz[ii] = mdcKalTrk->pz();
723 // pxy() and p() are not filled in the reconstruction algorithm
724 t_pxy2 = m_kal_px[i] * m_kal_px[i] + m_kal_py[i] * m_kal_py[i];
725 m_kal_p[i] = sqrt( t_pxy2 + m_kal_pz[i] * m_kal_pz[i] );
726 double ptrk = m_kal_p[i];
727 px_had += m_kal_px[i];
728 py_had += m_kal_py[i];
729 pz_had += m_kal_pz[i];
730 pt_had += sqrt( t_pxy2 );
731 p_had += m_kal_p[i];
732 e_had += sqrt( m_kal_p[i] * m_kal_p[i] + mdcKalTrk->mass() * mdcKalTrk->mass() );
733
734 if ( ( *itTrk )->isMdcDedxValid() )
735 { // DEDX information
736
737 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
738 m_probPH[ii] = dedxTrk->probPH();
739 m_normPH[ii] = dedxTrk->normPH();
740
741 m_chie[ii] = dedxTrk->chiE();
742 m_chimu[ii] = dedxTrk->chiMu();
743 m_chipi[ii] = dedxTrk->chiPi();
744 m_chik[ii] = dedxTrk->chiK();
745 m_chip[ii] = dedxTrk->chiP();
746 m_ghit[ii] = dedxTrk->numGoodHits();
747 m_thit[ii] = dedxTrk->numTotalHits();
748 }
749
750 if ( ( *itTrk )->isEmcShowerValid() )
751 {
752
753 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
754 m_e_emc[ii] = emcTrk->energy();
755 m_phi_emc[ii] = emcTrk->phi();
756 m_theta_emc[ii] = emcTrk->theta();
757 }
758
759 if ( ( *itTrk )->isMucTrackValid() )
760 {
761
762 RecMucTrack* mucTrk = ( *itTrk )->mucTrack();
763 m_nhit_muc[ii] = mucTrk->numHits();
764 m_nlay_muc[ii] = mucTrk->numLayers();
765 }
766 else
767 {
768 m_nhit_muc[ii] = 0;
769 m_nlay_muc[ii] = 0;
770 }
771 m_t_btof[ii] = 0;
772 m_t_etof[ii] = 0;
773
774 if ( ( *itTrk )->isTofTrackValid() )
775 { // TOF information
776
777 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
778
779 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
780
781 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
782 {
783 TofHitStatus* status = new TofHitStatus;
784 status->setStatus( ( *iter_tof )->status() );
785
786 if ( !( status->is_barrel() ) )
787 { // endcap
788 if ( ( status->is_cluster() ) ) m_t_etof[ii] = ( *iter_tof )->tof();
789 if ( !( status->is_counter() ) )
790 {
791 if ( status ) delete status;
792 continue;
793 } // ?
794 if ( status->layer() != 0 )
795 {
796 if ( status ) delete status;
797 continue;
798 } // layer1
799 double path = ( *iter_tof )->path(); // ?
800 double tof = ( *iter_tof )->tof();
801 double ph = ( *iter_tof )->ph();
802 double rhit = ( *iter_tof )->zrhit();
803 double qual = 0.0 + ( *iter_tof )->quality();
804 double cntr = 0.0 + ( *iter_tof )->tofID();
805 double texp[5];
806 for ( int j = 0; j < 5; j++ )
807 {
808 double gb = ptrk / xmass[j];
809 double beta = gb / sqrt( 1 + gb * gb );
810 texp[j] = path / beta / velc;
811 }
812
813 m_qual_etof[ii] = qual;
814 m_tof_etof[ii] = tof;
815 }
816 else
817 { // barrel
818 if ( ( status->is_cluster() ) ) m_t_btof[ii] = ( *iter_tof )->tof();
819 if ( !( status->is_counter() ) )
820 {
821 if ( status ) delete status;
822 continue;
823 } // ?
824 if ( status->layer() == 1 )
825 { // layer1
826 double path = ( *iter_tof )->path(); // ?
827 double tof = ( *iter_tof )->tof();
828 double ph = ( *iter_tof )->ph();
829 double rhit = ( *iter_tof )->zrhit();
830 double qual = 0.0 + ( *iter_tof )->quality();
831 double cntr = 0.0 + ( *iter_tof )->tofID();
832 double texp[5];
833 for ( int j = 0; j < 5; j++ )
834 {
835 double gb = ptrk / xmass[j];
836 double beta = gb / sqrt( 1 + gb * gb );
837 texp[j] = path / beta / velc;
838 }
839
840 m_qual_btof1[ii] = qual;
841 m_tof_btof1[ii] = tof;
842 }
843
844 if ( status->layer() == 2 )
845 { // layer2
846 double path = ( *iter_tof )->path(); // ?
847 double tof = ( *iter_tof )->tof();
848 double ph = ( *iter_tof )->ph();
849 double rhit = ( *iter_tof )->zrhit();
850 double qual = 0.0 + ( *iter_tof )->quality();
851 double cntr = 0.0 + ( *iter_tof )->tofID();
852 double texp[5];
853 for ( int j = 0; j < 5; j++ )
854 {
855 double gb = ptrk / xmass[j];
856 double beta = gb / sqrt( 1 + gb * gb );
857 texp[j] = path / beta / velc;
858 }
859
860 m_qual_btof2[ii] = qual;
861 m_tof_btof2[ii] = tof;
862 }
863 }
864 if ( status ) delete status;
865 }
866 }
867 }
868 counter[4]++;
869
870 // std::cout<<"dbg_5"<<std::endl;
871 // tag
872
873 m_dimutag = 0;
874 if ( m_ngch != 2 || nCharge != 0 ) return StatusCode::SUCCESS;
875
876 EvtRecTrackIterator itTrk1;
877
878 EvtRecTrackIterator itTrk2;
879
880 RecMdcKalTrack* mdcKalTrk1{ nullptr };
881 RecMdcKalTrack* mdcKalTrk2{ nullptr };
882
883 HepLorentzVector p41e, p42e, p4le;
884 Hep3Vector p31e, p32e, p3le;
885 HepLorentzVector p41m, p42m, p4lm;
886 Hep3Vector p31m, p32m, p3lm;
887 HepLorentzVector p41h, p42h, p4lh;
888 Hep3Vector p31h, p32h, p3lh;
889 WTrackParameter w1_ini;
890 WTrackParameter w2_ini;
891 int iip = -1;
892 int iim = -1;
893 for ( int i = 0; i < m_ngch; i++ )
894 {
895 if ( m_charge[i] > 0 ) itTrk1 = evtRecTrkCol->begin() + iGood[i];
896 if ( m_charge[i] < 0 ) itTrk2 = evtRecTrkCol->begin() + iGood[i];
897 if ( m_charge[i] > 0 ) mdcKalTrk1 = ( *itTrk1 )->mdcKalTrack();
898 if ( m_charge[i] < 0 ) mdcKalTrk2 = ( *itTrk2 )->mdcKalTrack();
899 if ( m_charge[i] > 0 ) iip = i;
900 if ( m_charge[i] < 0 ) iim = i;
901
902 if ( m_charge[i] > 0 )
903 w1_ini =
904 WTrackParameter( xmass[1], mdcKalTrk1->getZHelixMu(), mdcKalTrk1->getZErrorMu() );
905 if ( m_charge[i] < 0 )
906 w2_ini =
907 WTrackParameter( xmass[1], mdcKalTrk2->getZHelixMu(), mdcKalTrk2->getZErrorMu() );
908
909 if ( m_charge[i] > 0 ) p41m = w1_ini.p();
910 if ( m_charge[i] < 0 ) p42m = w2_ini.p();
911 if ( m_charge[i] > 0 ) p41m.boost( u_cms );
912 if ( m_charge[i] < 0 ) p42m.boost( u_cms );
913 if ( m_charge[i] > 0 ) p31m = p41m.vect();
914 if ( m_charge[i] < 0 ) p32m = p42m.vect();
915
916 if ( m_charge[i] > 0 )
917 {
918 m_px_cms_ep = p41m.px();
919 m_py_cms_ep = p41m.py();
920 m_pz_cms_ep = p41m.pz();
921 m_e_cms_ep = p41m.e();
922 }
923 if ( m_charge[i] < 0 )
924 {
925 m_px_cms_em = p42m.px();
926 m_py_cms_em = p42m.py();
927 m_pz_cms_em = p42m.pz();
928 m_e_cms_em = p42m.e();
929 }
930 }
931
932 double e01 = ( iip != -1 ) ? m_e_emc[iip] : 0; // m_e_cms_ep;
933 double e02 = ( iim != -1 ) ? m_e_emc[iim] : 0; // m_e_cms_em;
934 int ilarge = ( e01 > e02 ) ? iip : iim;
935
936 p4lm = ( e01 > e02 ) ? p41m : p42m;
937
938 p3lm = ( e01 > e02 ) ? p31m : p32m;
939
940 double acollm = 180. - p31m.angle( p32m ) * 180.0 / CLHEP::pi;
941 double acoplm = 180. - ( p31m.perpPart() ).angle( p32m.perpPart() ) * 180.0 / CLHEP::pi;
942 double poeb1m = p41m.rho() / beamEnergy;
943 double poeb2m = p42m.rho() / beamEnergy;
944 double poeblm = p4lm.rho() / beamEnergy;
945
946 double eemc1 = m_e_emc[iip];
947 double eemc2 = m_e_emc[iim];
948
949 double ex1 = m_kal_vx0[iip];
950 double ey1 = m_kal_vy0[iip];
951 double ez1 = m_kal_vz0[iip];
952 double epx1 = m_kal_px[iip];
953 double epy1 = m_kal_py[iip];
954 double epz1 = m_kal_pz[iip];
955 double epp1 = m_kal_p[iip];
956 double ex2 = m_kal_vx0[iim];
957 double ey2 = m_kal_vy0[iim];
958 double ez2 = m_kal_vz0[iim];
959 double epx2 = m_kal_px[iim];
960 double epy2 = m_kal_py[iim];
961 double epz2 = m_kal_pz[iim];
962 double epp2 = m_kal_p[iim];
963
964 double pidchidedx1 = m_pidchiDedx[iip];
965 double pidchitof11 = m_pidchiTof1[iip];
966 double pidchitof21 = m_pidchiTof2[iip];
967 double pidchidedx2 = m_pidchiDedx[iim];
968 double pidchitof12 = m_pidchiTof1[iim];
969 double pidchitof22 = m_pidchiTof2[iim];
970
971 double eoeb1 = m_e_emc[iip] / beamEnergy;
972 double eoeb2 = m_e_emc[iim] / beamEnergy;
973
974 double eopm1 = 0;
975 if ( p41m.rho() > 0 ) eopm1 = m_e_emc[iip] / p41m.rho();
976 double eopm2 = 0;
977 if ( p42m.rho() > 0 ) eopm2 = m_e_emc[iim] / p42m.rho();
978
979 double exoeb1 = m_e_emc[iip] * sin( m_theta_emc[iip] ) * cos( m_phi_emc[iip] ) / beamEnergy;
980 double eyoeb1 = m_e_emc[iip] * sin( m_theta_emc[iip] ) * sin( m_phi_emc[iip] ) / beamEnergy;
981 double ezoeb1 = m_e_emc[iip] * cos( m_theta_emc[iip] ) / beamEnergy;
982 double etoeb1 = m_e_emc[iip] * sin( m_theta_emc[iip] ) / beamEnergy;
983
984 double exoeb2 = m_e_emc[iim] * sin( m_theta_emc[iim] ) * cos( m_phi_emc[iim] ) / beamEnergy;
985 double eyoeb2 = m_e_emc[iim] * sin( m_theta_emc[iim] ) * sin( m_phi_emc[iim] ) / beamEnergy;
986 double ezoeb2 = m_e_emc[iim] * cos( m_theta_emc[iim] ) / beamEnergy;
987 double etoeb2 = m_e_emc[iim] * sin( m_theta_emc[iim] ) / beamEnergy;
988
989 double eoebl = m_e_emc[ilarge] / beamEnergy;
990
991 double eopl = 0;
992 if ( p4lm.rho() > 0 ) eopl = m_e_emc[ilarge] / p4lh.rho();
993
994 double exoebl =
995 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * cos( m_phi_emc[ilarge] ) / beamEnergy;
996 double eyoebl =
997 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * sin( m_phi_emc[ilarge] ) / beamEnergy;
998 double ezoebl = m_e_emc[ilarge] * cos( m_theta_emc[ilarge] ) / beamEnergy;
999 double etoebl = m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) / beamEnergy;
1000
1001 int mucinfo1 = ( m_nhit_muc[iip] >= 2 && m_nlay_muc[iip] >= 2 ) ? 1 : 0;
1002 int mucinfo2 = ( m_nhit_muc[iim] >= 2 && m_nlay_muc[iim] >= 2 ) ? 1 : 0;
1003 int mucinfol = ( m_nhit_muc[ilarge] >= 2 && m_nlay_muc[ilarge] >= 2 ) ? 1 : 0;
1004 int pidel = ( e01 > e02 ) ? m_nep : m_nem;
1005 int pidmul = ( e01 > e02 ) ? m_nmup : m_nmum;
1006 double deltatof = 0;
1007
1008 if ( m_t_btof[iip] * m_t_btof[iim] != 0 ) deltatof += fabs( m_t_btof[iip] - m_t_btof[iim] );
1009 if ( m_t_etof[iip] * m_t_etof[iim] != 0 ) deltatof += fabs( m_t_etof[iip] - m_t_etof[iim] );
1010
1011 if ( acollm < m_acoll_mu_cut ) m_dimutag += 1;
1012 if ( acoplm < m_acopl_mu_cut ) m_dimutag += 10;
1013 if ( fabs( m_ecms - mpsi2s ) < 0.001 )
1014 {
1015 if ( ( sqrt( ( ( poeb1m - poeb2m ) / 0.35 ) * ( ( poeb1m - poeb2m ) / 0.35 ) +
1016 ( ( poeb1m + poeb2m - 1.68 ) / 0.125 ) *
1017 ( ( poeb1m + poeb2m - 1.68 ) / 0.125 ) ) > m_tptotal_mu_cut ) &&
1018 ( ( ( poeb1m >= m_tpoebh_mu_cut ) && ( poeb2m >= m_tpoebl_mu_cut ) ) ||
1019 ( ( poeb2m >= m_tpoebh_mu_cut ) && ( poeb1m >= m_tpoebl_mu_cut ) ) ) )
1020 m_dimutag += 100;
1021 }
1022 else
1023 {
1024 if ( ( fabs( poeb1m - 1 ) < m_poeb_mu_cut ) && ( fabs( poeb2m - 1 ) < m_poeb_mu_cut ) )
1025 m_dimutag += 100;
1026 }
1027 if ( !m_useTOF || ( deltatof < m_dtof_mu_cut ) ) m_dimutag += 1000;
1028 if ( !m_useMUC || ( mucinfo1 == 1 && mucinfo2 == 1 ) ) m_dimutag += 10000;
1029 if ( etoeb1 < m_eoeb_mu_cut && eoeb2 < m_eoeb_mu_cut ) m_dimutag += 100000;
1030 if ( etoeb1 + etoeb2 < m_etotal_mu_cut ) m_dimutag += 1000000;
1031 if ( !m_usePID || ( m_nmup == 1 && m_nmum == 1 ) ) m_dimutag += 10000000;
1032
1033 m_acoll = acollm;
1034 m_acopl = acoplm;
1035 m_poeb1 = poeb1m;
1036 m_poeb2 = poeb2m;
1037 m_eop1 = eopm1;
1038 m_eop2 = eopm2;
1039 m_cos_ep = p41m.cosTheta();
1040 m_cos_em = p42m.cosTheta();
1041 m_mass_ee = ( p41m + p42m ).m();
1042
1043 m_deltatof = deltatof;
1044
1045 m_eoeb1 = eoeb1;
1046 m_eoeb2 = eoeb2;
1047
1048 m_etoeb1 = etoeb1;
1049 m_etoeb2 = etoeb2;
1050 m_mucinfo1 = mucinfo1;
1051 m_mucinfo2 = mucinfo2;
1052
1053 if ( m_dimutag == 11111111 )
1054 {
1055 ndimu++;
1056 ////////////////////////////////////////////////////////////
1057 // DQA
1058 // set tag and quality
1059
1060 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
1061 ( *itTrk1 )->tagMuon();
1062 ( *itTrk2 )->tagMuon();
1063
1064 // Quality: defined by whether dE/dx or TOF is used to identify particle
1065 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
1066 // 1 - only dE/dx (can be used for TOF calibration)
1067 // 2 - only TOF (can be used for dE/dx calibration)
1068 // 3 - Both dE/dx and TOF
1069
1070 ( *itTrk1 )->setQuality( 0 );
1071 ( *itTrk2 )->setQuality( 0 );
1072 // DQA
1073 // Add the line below at the end of execute(), (before return)
1074 //
1075 setFilterPassed( true );
1076 ////////////////////////////////////////////////////////////
1077
1078 m_mumu_mass->Fill( ( p41m + p42m ).m() );
1079 m_mumu_acoll->Fill( acollm );
1080 m_mumu_eop_mup->Fill( eopm1 );
1081 m_mumu_eop_mum->Fill( eopm2 );
1082 m_mumu_costheta_mup->Fill( p41m.cosTheta() );
1083 m_mumu_costheta_mum->Fill( p42m.cosTheta() );
1084 m_mumu_phi_mup->Fill( p41m.phi() );
1085 m_mumu_phi_mum->Fill( p42m.phi() );
1086 m_mumu_nneu->Fill( m_nGam );
1087 m_mumu_nhit->Fill( m_nhit_muc[ilarge] );
1088 m_mumu_nlay->Fill( m_nlay_muc[ilarge] );
1089
1090 m_mumu_eemc_mup->Fill( eemc1 );
1091 m_mumu_eemc_mum->Fill( eemc2 );
1092
1093 m_mumu_x_mup->Fill( ex1 );
1094 m_mumu_y_mup->Fill( ey1 );
1095 m_mumu_z_mup->Fill( ez1 );
1096
1097 m_mumu_x_mum->Fill( ex2 );
1098 m_mumu_y_mum->Fill( ey2 );
1099 m_mumu_z_mum->Fill( ez2 );
1100
1101 m_mumu_px_mup->Fill( epx1 );
1102 m_mumu_py_mup->Fill( epy1 );
1103 m_mumu_pz_mup->Fill( epz1 );
1104 m_mumu_p_mup->Fill( epp1 );
1105
1106 m_mumu_px_mum->Fill( epx2 );
1107 m_mumu_py_mum->Fill( epy2 );
1108 m_mumu_pz_mum->Fill( epz2 );
1109 m_mumu_p_mum->Fill( epp2 );
1110
1111 m_mumu_deltatof->Fill( deltatof );
1112
1113 m_mumu_pidchidedx_mup->Fill( pidchidedx1 );
1114 m_mumu_pidchidedx_mum->Fill( pidchidedx2 );
1115 m_mumu_pidchitof1_mup->Fill( pidchitof11 );
1116 m_mumu_pidchitof1_mum->Fill( pidchitof12 );
1117 m_mumu_pidchitof2_mup->Fill( pidchitof21 );
1118 m_mumu_pidchitof2_mum->Fill( pidchitof22 );
1119
1120 if ( m_writentuple ) m_tuple1->write().ignore();
1121 }
1122
1123 return StatusCode::SUCCESS;
1124}
1125
1126// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1127StatusCode DQASelDimu::finalize() {
1128
1129 MsgStream log( msgSvc(), name() );
1130 log << MSG::INFO << "in finalize()" << endmsg;
1131 return StatusCode::SUCCESS;
1132}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
u_cms
const Hep3Vector u_cms
Definition DQADtagAlg.cxx:56
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
mpsi2s
const double mpsi2s
Definition DQASelBhabha.cxx:47
n
const Int_t n
Definition DataBase/tau_mode.c:57
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi
double mpi
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
xmass
const double xmass[5]
Definition Gam4pikp.cxx:35
velc
const double velc
Definition Gam4pikp.cxx:36
mk
const double mk
Definition Gam4pikp.cxx:33
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:37
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
energy
************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
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DQASelDimu::DQASelDimu
DQASelDimu(const std::string &name, ISvcLocator *pSvcLocator)
Definition DQASelDimu.cxx:58
DQASelDimu::initialize
StatusCode initialize()
Definition DQASelDimu.cxx:98
DQASelDimu::finalize
StatusCode finalize()
Definition DQASelDimu.cxx:1127
DQASelDimu::execute
StatusCode execute()
Definition DQASelDimu.cxx:353
DstEmcShower::dy
double dy() const
Definition DstEmcShower.cxx:10
DstEmcShower::dz
double dz() const
Definition DstEmcShower.cxx:17
DstEmcShower::dx
double dx() const
Definition DstEmcShower.cxx:3
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::instance
static ParticleID * instance()
Definition ParticleID.cxx:21
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition ParticleID.cxx:97
ParticleID::init
void init()
Definition ParticleID.cxx:26
ParticleID::chiDedx
double chiDedx(int n) const
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:115