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