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