CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
DQAKsKpi.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/AlgFactory.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/PropertyMgr.h"
7#include "VertexFit/IVertexDbSvc.h"
8#include "GaudiKernel/Bootstrap.h"
9
10#include "EventModel/EventModel.h"
11#include "EventModel/Event.h"
12#include "EventModel/EventHeader.h"
13
14#include "EvtRecEvent/EvtRecEvent.h"
15#include "EvtRecEvent/EvtRecTrack.h"
16#include "DstEvent/TofHitStatus.h"
17
18#include "TMath.h"
19#include "CLHEP/Vector/ThreeVector.h"
20#include "CLHEP/Vector/LorentzVector.h"
21#include "CLHEP/Vector/TwoVector.h"
22
23#include "VertexFit/KinematicFit.h"
24#include "VertexFit/VertexFit.h"
25#include "VertexFit/SecondVertexFit.h"
26#include "VertexFit/IVertexDbSvc.h"
27#include "VertexFit/Helix.h"
28
29#include "ParticleID/ParticleID.h"
30
31#include "TMath.h"
32#include "TH1F.h"
33#include "GaudiKernel/INTupleSvc.h"
34#include "GaudiKernel/NTuple.h"
35#include "GaudiKernel/Bootstrap.h"
36#include "GaudiKernel/ISvcLocator.h"
37#include "GaudiKernel/ITHistSvc.h"
38//#include "GaudiKernel/IHistogramSvc.h"
39
40using CLHEP::Hep3Vector;
41using CLHEP::Hep2Vector;
42using CLHEP::HepLorentzVector;
43#include "CLHEP/Geometry/Point3D.h"
44#ifndef ENABLE_BACKWARDS_COMPATIBILITY
46#endif
47
48#include "VertexFit/KinematicFit.h"
49#include "VertexFit/VertexFit.h"
50#include "VertexFit/SecondVertexFit.h"
51#include "VertexFit/IVertexDbSvc.h"
52#include "ParticleID/ParticleID.h"
53#include <vector>
54
55#include "DQAKsKpiAlg/DQAKsKpi.h"
56
57const double mpi = 0.13957;
58const double mk = 0.493677;
59const double mks0 = 0.497648;
60const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
61const double velc = 299.792458; // tof path unit in mm
62typedef std::vector<int> Vint;
63typedef std::vector<HepLorentzVector> Vp4;
64
65//boost
66const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097);
67const Hep3Vector u_cms = -p_cms.boostVector();
68 //declare one counter
69static int counter[10]={0,0,0,0,0,0,0,0,0,0};
70
71/////////////////////////////////////////////////////////////////////////////
72
73DQAKsKpi::DQAKsKpi(const std::string& name, ISvcLocator* pSvcLocator) :
74 Algorithm(name, pSvcLocator) {
75
76 //Declare the properties
77 declareProperty("Vr0cut", m_vr0cut=5.0);
78 declareProperty("Vz0cut", m_vz0cut=20.0);
79 declareProperty("Vr1cut", m_vr1cut=1.0);
80 declareProperty("Vz1cut", m_vz1cut=5.0);
81 declareProperty("Vctcut", m_cthcut=0.93);
82 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
83 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
84}
85
86
87// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
89 MsgStream log(msgSvc(), name());
90
91 log << MSG::INFO << "in initialize()" << endmsg;
92 StatusCode status;
93
94 if(service("THistSvc", m_thsvc).isFailure()) {
95 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
96 return StatusCode::FAILURE;
97 }
98
99 // "DQAHist" is fixed
100 // "DQAKsKpi" is control sample name, just as ntuple case.
101 TH1F* hks_dl = new TH1F("ks_dl", "ks_dl", 300, -5.0, 25.0);
102 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_dl", hks_dl).isFailure()) {
103 log << MSG::ERROR << "Couldn't register ks_dl" << endreq;
104 }
105
106 TH1F* hks_m = new TH1F("ks_m", "ks_m", 200,0.4, 0.6);
107 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_m", hks_m).isFailure()) {
108 log << MSG::ERROR << "Couldn't register ks_m" << endreq;
109 }
110
111 TH1F* hkspi_m = new TH1F("kspi_m", "kspi_m", 200,0.6, 2.6);
112 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hkspi_m", hkspi_m).isFailure()) {
113 log << MSG::ERROR << "Couldn't register kspi_m" << endreq;
114 }
115
116 TH1F* hks_p = new TH1F("ks_p", "ks_p", 100,0.0, 1.5);
117 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_p", hks_p).isFailure()) {
118 log << MSG::ERROR << "Couldn't register ks_p" << endreq;
119 }
120
121 TH1F* hkpi_m = new TH1F("kpi_m", "kpi_m", 200,0.6, 2.6);
122 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hkpi_m", hkpi_m).isFailure()) {
123 log << MSG::ERROR << "Couldn't register kpi_m" << endreq;
124 }
125
126 // DQA
127 // The first directory specifier must be "DQAFILE"
128 // The second is the control sample name, the first letter is in upper format. eg. "Rhopi"
129 NTuplePtr nt(ntupleSvc(), "DQAFILE/KsKpi");
130 if ( nt ) m_tuple = nt;
131 else {
132 m_tuple = ntupleSvc()->book("DQAFILE/KsKpi", CLID_ColumnWiseTuple, "KsKpi ntuple");
133 if( m_tuple ) {
134 status = m_tuple->addItem("runNo", m_runNo);
135 status = m_tuple->addItem("event", m_event);
136// status = m_tuple->addItem("Nchrg", m_nchrg);
137// status = m_tuple->addItem("Nneu", m_nneu);
138
139 status = m_tuple->addItem("npip", m_npip);
140 status = m_tuple->addItem("npim", m_npim);
141 status = m_tuple->addItem("nkp", m_nkp);
142 status = m_tuple->addItem("nkm", m_nkm);
143 status = m_tuple->addItem("np", m_np);
144 status = m_tuple->addItem("npb", m_npb);
145
146 status = m_tuple->addItem("vfits_chi", m_vfits_chi);
147 status = m_tuple->addItem("vfits_vx", m_vfits_vx);
148 status = m_tuple->addItem("vfits_vy", m_vfits_vy);
149 status = m_tuple->addItem("vfits_vz", m_vfits_vz);
150 status = m_tuple->addItem("vfits_vr", m_vfits_vr);
151
152 status = m_tuple->addItem("vfitp_chi", m_vfitp_chi);
153 status = m_tuple->addItem("vfitp_vx", m_vfitp_vx);
154 status = m_tuple->addItem("vfitp_vy", m_vfitp_vy);
155 status = m_tuple->addItem("vfitp_vz", m_vfitp_vz);
156 status = m_tuple->addItem("vfitp_vr", m_vfitp_vr);
157
158 status = m_tuple->addItem("vfit2_chi", m_vfit2_chi);
159 status = m_tuple->addItem("vfit2_mks", m_vfit2_mks);
160 status = m_tuple->addItem("vfit2_ct", m_vfit2_ct);
161 status = m_tuple->addItem("vfit2_dl", m_vfit2_dl);
162 status = m_tuple->addItem("vfit2_dle", m_vfit2_dle);
163
164 status = m_tuple->addItem("chi2_fs4c", m_chi2_fs4c);
165 status = m_tuple->addItem("mks_fs4c", m_mks_fs4c);
166 status = m_tuple->addItem("mkspi_fs4c",m_mkspi_fs4c);
167 status = m_tuple->addItem("mksk_fs4c", m_mksk_fs4c);
168 status = m_tuple->addItem("mkpi_fs4c", m_mkpi_fs4c);
169
170 status = m_tuple->addItem("4c_chi2", m_4c_chi2);
171 status = m_tuple->addItem("4c_mks", m_4c_mks);
172 status = m_tuple->addItem("4c_mkspi", m_4c_mkspi);
173 status = m_tuple->addItem("4c_mksk", m_4c_mksk);
174 status = m_tuple->addItem("4c_mkpi", m_4c_mkpi);
175 status = m_tuple->addItem("4c_ks_px", m_4c_ks_px);
176 status = m_tuple->addItem("4c_ks_py", m_4c_ks_py);
177 status = m_tuple->addItem("4c_ks_pz", m_4c_ks_pz);
178 status = m_tuple->addItem("4c_ks_p", m_4c_ks_p);
179 status = m_tuple->addItem("4c_ks_cos", m_4c_ks_cos);
180
181 status = m_tuple->addItem("NGch", m_ngch, 0, 10);
182 status = m_tuple->addIndexedItem("pidcode", m_ngch, m_pidcode);
183 status = m_tuple->addIndexedItem("pidprob", m_ngch, m_pidprob);
184 status = m_tuple->addIndexedItem("pidchiDedx", m_ngch, m_pidchiDedx);
185 status = m_tuple->addIndexedItem("pidchiTof1", m_ngch, m_pidchiTof1);
186 status = m_tuple->addIndexedItem("pidchiTof2", m_ngch, m_pidchiTof2);
187
188 status = m_tuple->addIndexedItem("charge",m_ngch, m_charge);
189 status = m_tuple->addIndexedItem("vx0", m_ngch, m_vx0);
190 status = m_tuple->addIndexedItem("vy0", m_ngch, m_vy0);
191 status = m_tuple->addIndexedItem("vz0", m_ngch, m_vz0);
192 status = m_tuple->addIndexedItem("vr0", m_ngch, m_vr0);
193
194 status = m_tuple->addIndexedItem("vx", m_ngch, m_vx);
195 status = m_tuple->addIndexedItem("vy", m_ngch, m_vy);
196 status = m_tuple->addIndexedItem("vz", m_ngch, m_vz);
197 status = m_tuple->addIndexedItem("vr", m_ngch, m_vr);
198
199 status = m_tuple->addIndexedItem("px", m_ngch, m_px);
200 status = m_tuple->addIndexedItem("py", m_ngch, m_py);
201 status = m_tuple->addIndexedItem("pz", m_ngch, m_pz);
202 status = m_tuple->addIndexedItem("p", m_ngch, m_p);
203 status = m_tuple->addIndexedItem("cost", m_ngch, m_cost);
204
205 status = m_tuple->addIndexedItem("probPH", m_ngch, m_probPH);
206 status = m_tuple->addIndexedItem("normPH", m_ngch, m_normPH);
207 status = m_tuple->addIndexedItem("chie", m_ngch, m_chie);
208 status = m_tuple->addIndexedItem("chimu", m_ngch, m_chimu);
209 status = m_tuple->addIndexedItem("chipi", m_ngch, m_chipi);
210 status = m_tuple->addIndexedItem("chik", m_ngch, m_chik);
211 status = m_tuple->addIndexedItem("chip", m_ngch, m_chip);
212 status = m_tuple->addIndexedItem("ghit", m_ngch, m_ghit);
213 status = m_tuple->addIndexedItem("thit", m_ngch, m_thit);
214
215 status = m_tuple->addIndexedItem("e_emc", m_ngch, m_e_emc);
216
217 status = m_tuple->addIndexedItem("qual_etof", m_ngch, m_qual_etof);
218 status = m_tuple->addIndexedItem("tof_etof", m_ngch, m_tof_etof);
219 status = m_tuple->addIndexedItem("te_etof", m_ngch, m_te_etof);
220 status = m_tuple->addIndexedItem("tmu_etof", m_ngch, m_tmu_etof);
221 status = m_tuple->addIndexedItem("tpi_etof", m_ngch, m_tpi_etof);
222 status = m_tuple->addIndexedItem("tk_etof", m_ngch, m_tk_etof);
223 status = m_tuple->addIndexedItem("tp_etof", m_ngch, m_tp_etof);
224
225 status = m_tuple->addIndexedItem("qual_btof1", m_ngch, m_qual_btof1);
226 status = m_tuple->addIndexedItem("tof_btof1", m_ngch, m_tof_btof1);
227 status = m_tuple->addIndexedItem("te_btof1", m_ngch, m_te_btof1);
228 status = m_tuple->addIndexedItem("tmu_btof1", m_ngch, m_tmu_btof1);
229 status = m_tuple->addIndexedItem("tpi_btof1", m_ngch, m_tpi_btof1);
230 status = m_tuple->addIndexedItem("tk_btof1", m_ngch, m_tk_btof1);
231 status = m_tuple->addIndexedItem("tp_btof1", m_ngch, m_tp_btof1);
232
233 status = m_tuple->addIndexedItem("qual_btof2", m_ngch, m_qual_btof2);
234 status = m_tuple->addIndexedItem("tof_btof2", m_ngch, m_tof_btof2);
235 status = m_tuple->addIndexedItem("te_btof2", m_ngch, m_te_btof2);
236 status = m_tuple->addIndexedItem("tmu_btof2", m_ngch, m_tmu_btof2);
237 status = m_tuple->addIndexedItem("tpi_btof2", m_ngch, m_tpi_btof2);
238 status = m_tuple->addIndexedItem("tk_btof2", m_ngch, m_tk_btof2);
239 status = m_tuple->addIndexedItem("tp_btof2", m_ngch, m_tp_btof2);
240
241 } else {
242 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
243 }
244 }
245
246//--------end of book--------
247
248 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
249 return StatusCode::SUCCESS;
250
251}
252
253// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
254StatusCode DQAKsKpi::execute() {
255
256 MsgStream log(msgSvc(), name());
257 log << MSG::INFO << "in execute()" << endreq;
258
259 // DQA
260 // Add the line below at the beginning of execute()
261// setFilterPassed(false);
262
263 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
264 m_runNo = eventHeader->runNumber();
265 m_event=eventHeader->eventNumber();
266 log << MSG::DEBUG <<"run, evtnum = "
267 << m_runNo << " , "
268 << m_event <<endreq;
269
270 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
271// m_nchrg = evtRecEvent->totalCharged();
272// m_nneu = evtRecEvent->totalNeutral();
273 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
274 << evtRecEvent->totalCharged() << " , "
275 << evtRecEvent->totalNeutral() << " , "
276 << evtRecEvent->totalTracks() <<endreq;
277
278 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
279
280 if(evtRecEvent->totalNeutral()>100) {
281 return StatusCode::SUCCESS;
282 }
283
284 Vint iGood;
285 iGood.clear();
286
287 Hep3Vector xorigin(0,0,0);
288
289 IVertexDbSvc* vtxsvc;
290 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
291 if(vtxsvc->isVertexValid()){
292 double* dbv = vtxsvc->PrimaryVertex();
293 double* vv = vtxsvc->SigmaPrimaryVertex();
294 xorigin.setX(dbv[0]);
295 xorigin.setY(dbv[1]);
296 xorigin.setZ(dbv[2]);
297 log << MSG::INFO
298 <<"xorigin.x="<<xorigin.x()<<", "
299 <<"xorigin.y="<<xorigin.y()<<", "
300 <<"xorigin.z="<<xorigin.z()<<", "
301 <<endreq ;
302 }
303
304 int nCharge = 0;
305 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
306 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
307 if(!(*itTrk)->isMdcTrackValid()) continue;
308 if (!(*itTrk)->isMdcKalTrackValid()) continue;
309 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
310
311 double pch =mdcTrk->p();
312 double x0 =mdcTrk->x();
313 double y0 =mdcTrk->y();
314 double z0 =mdcTrk->z();
315 double phi0=mdcTrk->helix(1);
316 double xv=xorigin.x();
317 double yv=xorigin.y();
318 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
319 double vx0 = x0;
320 double vy0 = y0;
321 double vz0 = z0-xorigin.z();
322 double vr0 = Rxy;
323 double Vct=cos(mdcTrk->theta());
324
325 HepVector a = mdcTrk->helix();
326 HepSymMatrix Ea = mdcTrk->err();
327 HepPoint3D point0(0.,0.,0.); // the initial point for MDC recosntruction
328 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
329 VFHelix helixip(point0,a,Ea);
330 helixip.pivot(IP);
331 HepVector vecipa = helixip.a();
332 double Rvxy0=fabs(vecipa[0]); //the nearest distance to IP in xy plane
333 double Rvz0=vecipa[3]; //the nearest distance to IP in z direction
334 double Rvphi0=vecipa[1];
335// m_rvxy0=Rvxy0;
336// m_rvz0=Rvz0;
337// m_rvphi0=Rvphi0;
338
339 if(fabs(Rvz0) >= m_vz0cut) continue;
340 if(Rvxy0 >= m_vr0cut) continue;
341// if(fabs(Vct)>=m_cthcut) continue;
342// iGood.push_back((*itTrk)->trackId());
343 iGood.push_back(i);
344 nCharge += mdcTrk->charge();
345 }
346
347 //
348 // Finish Good Charged Track Selection
349 //
350 int m_ngch = iGood.size();
351 log << MSG::DEBUG << "ngood, totcharge = " << m_ngch << " , " << nCharge << endreq;
352 if((m_ngch != 4)||(nCharge!=0)){
353 return StatusCode::SUCCESS;
354 }
355
356 // Particle ID
357 //
358 Vint ipip, ipim, ikp, ikm, ipp, ipm;
359 ipip.clear();
360 ipim.clear();
361 ikp.clear();
362 ikm.clear();
363 ipp.clear();
364 ipm.clear();
365
366 Vp4 p_pip, p_pim, p_kp, p_km, p_pp, p_pm ;
367 p_pip.clear();
368 p_pim.clear();
369 p_kp.clear();
370 p_km.clear();
371 p_pp.clear();
372 p_pm.clear();
373
375 for(int i = 0; i < m_ngch; i++) {
376 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
377 // if(pid) delete pid;
378 pid->init();
379 pid->setMethod(pid->methodProbability());
380 pid->setChiMinCut(4);
381 pid->setRecTrack(*itTrk);
382 pid->usePidSys(pid->useDedx()); // just use dedx PID
383 // pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTofE() | pid->useTofQ()); // use PID sub-system
384 pid->identify(pid->onlyPionKaonProton()); // seperater Pion/Kaon/Proton
385 // pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
386 // pid->identify(pid->onlyPion());
387 // pid->identify(pid->onlyKaon());
388 pid->calculate();
389 if(!(pid->IsPidInfoValid())) continue;
390
391 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
392 RecMdcKalTrack* mdcKalTrk = 0 ;
393 if((*itTrk)->isMdcKalTrackValid()) mdcKalTrk = (*itTrk)->mdcKalTrack();
394
395 double prob_pi = pid->probPion();
396 double prob_K = pid->probKaon();
397 double prob_p = pid->probProton();
398 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
399 HepLorentzVector ptrk;
400 ptrk.setPx(mdcTrk->px()) ;
401 ptrk.setPy(mdcTrk->py()) ;
402 ptrk.setPz(mdcTrk->pz()) ;
403 double p3 = ptrk.mag() ;
404
405 if (prob_pi > prob_K && prob_pi > prob_p) {
406 m_pidcode[i]=2;
407 m_pidprob[i]=pid->prob(2);
408 m_pidchiDedx[i]=pid->chiDedx(2);
409 m_pidchiTof1[i]=pid->chiTof1(2);
410 m_pidchiTof2[i]=pid->chiTof2(2);
411 ptrk.setE(sqrt(p3*p3 + xmass[2]*xmass[2])) ;
412 if(mdcTrk->charge() > 0) {
413 ipip.push_back(iGood[i]);
414 p_pip.push_back(ptrk);
415 }
416 if (mdcTrk->charge() < 0) {
417 ipim.push_back(iGood[i]);
418 p_pim.push_back(ptrk);
419 }
420 }
421
422 if (prob_K > prob_pi && prob_K > prob_p) {
423 m_pidcode[i]=3;
424 m_pidprob[i]=pid->prob(3);
425 m_pidchiDedx[i]=pid->chiDedx(3);
426 m_pidchiTof1[i]=pid->chiTof1(3);
427 m_pidchiTof2[i]=pid->chiTof2(3);
428 ptrk.setE(sqrt(p3*p3 + xmass[3]*xmass[3])) ;
429 if(mdcTrk->charge() > 0) {
430 ikp.push_back(iGood[i]);
431 p_kp.push_back(ptrk);
432 }
433 if (mdcTrk->charge() < 0) {
434 ikm.push_back(iGood[i]);
435 p_km.push_back(ptrk);
436 }
437 }
438
439 if (prob_p > prob_pi && prob_p > prob_K) {
440 m_pidcode[i]=4;
441 m_pidprob[i]=pid->prob(4);
442 m_pidchiDedx[i]=pid->chiDedx(4);
443 m_pidchiTof1[i]=pid->chiTof1(4);
444 m_pidchiTof2[i]=pid->chiTof2(4);
445 ptrk.setE(sqrt(p3*p3 + xmass[4]*xmass[4])) ;
446 if(mdcTrk->charge() > 0) {
447 ipp.push_back(iGood[i]);
448 p_pp.push_back(ptrk);
449 }
450 if (mdcTrk->charge() < 0) {
451 ipm.push_back(iGood[i]);
452 p_pm.push_back(ptrk);
453 }
454 }
455 }
456
457 m_npip= ipip.size() ;
458 m_npim= ipim.size() ;
459 m_nkp = ikp.size() ;
460 m_nkm = ikm.size() ;
461 m_np = ipp.size() ;
462 m_npb = ipm.size() ;
463
464// cout<<"m_npip*m_npim: "<<m_npip*m_npim<<endl;
465 if( m_npip*m_npim != 2 ) {
466 return StatusCode::SUCCESS;
467 }
468// cout<<"m_nkp+m_nkm: "<<m_nkp+m_nkm<<endl;
469 if( m_nkp+m_nkm != 1 ) {
470 return StatusCode::SUCCESS;
471 }
472// cout<<"haha"<<endl;
473
474 // Vertex Fit
475 HepPoint3D vx(0., 0., 0.);
476 HepSymMatrix Evx(3, 0);
477 double bx = 1E+6;
478 double by = 1E+6;
479 double bz = 1E+6;
480 Evx[0][0] = bx*bx;
481 Evx[1][1] = by*by;
482 Evx[2][2] = bz*bz;
483
484 VertexParameter vxpar;
485 vxpar.setVx(vx);
486 vxpar.setEvx(Evx);
487
488 VertexFit *vtxfit_s = VertexFit::instance(); // S second vertex
489 VertexFit *vtxfit_p = VertexFit::instance(); // P primary vertex
491
492 RecMdcKalTrack *pi1KalTrk, *pi2KalTrk, *pi3KalTrk, *k1KalTrk;
493 RecMdcKalTrack *pipKalTrk, *pimKalTrk, *piKalTrk, *kKalTrk;
494 WTrackParameter wks,vwks;
495
496 double chi_temp = 999.0;
497 double mks_temp = 10.0 ;
498 bool okloop=false;
499 for(unsigned int i1 = 0; i1 < m_npip; i1++) { //pion plus
500 RecMdcKalTrack *pi1KalTrk = (*(evtRecTrkCol->begin()+ipip[i1]))-> mdcKalTrack();
502 WTrackParameter wpi1trk(xmass[2], pi1KalTrk->getZHelix(), pi1KalTrk->getZError());
503
504 for(unsigned int i2 = 0; i2 < m_npim; i2++) { //pion minus
505 RecMdcKalTrack *pi2KalTrk = (*(evtRecTrkCol->begin()+ipim[i2]))-> mdcKalTrack();
507 WTrackParameter wpi2trk(xmass[2], pi2KalTrk->getZHelix(), pi2KalTrk->getZError());
508
509 vtxfit_s->init();
510 vtxfit_s->AddTrack(0, wpi1trk);
511 vtxfit_s->AddTrack(1, wpi2trk);
512 vtxfit_s->AddVertex(0, vxpar, 0, 1);
513
514 if(!(vtxfit_s->Fit(0))) continue;
515 vtxfit_s->BuildVirtualParticle(0);
516 m_vfits_chi = vtxfit_s->chisq(0);
517 WTrackParameter wkshort = vtxfit_s->wVirtualTrack(0);
518 VertexParameter vparks = vtxfit_s->vpar(0);
519
520 m_vfits_vx = (vparks.Vx())[0];
521 m_vfits_vy = (vparks.Vx())[1];
522 m_vfits_vz = (vparks.Vx())[2];
523 m_vfits_vr = sqrt(m_vfits_vx*m_vfits_vx + m_vfits_vy*m_vfits_vy) ;
524
525 if ( m_npip == 2 ) {
526 int j = i1 ;
527 int jj = ( i1 == 1 ) ? 0 : 1;
528 pi3KalTrk = (*(evtRecTrkCol->begin()+ipip[jj]))->mdcKalTrack();
529 k1KalTrk = (*(evtRecTrkCol->begin()+ikm[0]))->mdcKalTrack();
530 }
531 if (m_npim == 2 ) {
532 int j = i2 ;
533 int jj = ( i2 == 1 ) ? 0 : 1;
534 pi3KalTrk = (*(evtRecTrkCol->begin()+ipim[jj]))->mdcKalTrack();
535 k1KalTrk = (*(evtRecTrkCol->begin()+ikp[0]))->mdcKalTrack();
536 }
537
539 WTrackParameter wpi3trk( xmass[2], pi3KalTrk->getZHelix(), pi3KalTrk->getZError());
541 WTrackParameter wk1trk( xmass[3], k1KalTrk->getZHelixK(), k1KalTrk->getZErrorK());
542
543 vtxfit_p->init();
544// vtxfit_p->AddTrack(0, wkshort);
545 vtxfit_p->AddTrack(0, wpi3trk);
546 vtxfit_p->AddTrack(1, wk1trk);
547 vtxfit_p->AddVertex(0, vxpar, 0, 1);
548 if(!(vtxfit_p->Fit(0))) continue;
549
550 m_vfitp_chi = vtxfit_p->chisq(0) ;
551
552 VertexParameter primaryVpar = vtxfit_p->vpar(0);
553 m_vfitp_vx = (primaryVpar.Vx())[0];
554 m_vfitp_vy = (primaryVpar.Vx())[1];
555 m_vfitp_vz = (primaryVpar.Vx())[2];
556 m_vfitp_vr = sqrt(m_vfitp_vx*m_vfitp_vx + m_vfitp_vy*m_vfitp_vy);
557
558 vtxfit2->init();
559 vtxfit2->setPrimaryVertex(vtxfit_p->vpar(0));
560 vtxfit2->AddTrack(0, wkshort);
561 vtxfit2->setVpar(vtxfit_s->vpar(0));
562 if(!vtxfit2->Fit()) continue;
563
564 if ( fabs(((vtxfit2->wpar()).p()).m()-mks0) < mks_temp ) {
565 mks_temp = fabs(((vtxfit2->wpar()).p()).m()-mks0) ;
566
567 okloop = true;
568
569 wks = vtxfit2->wpar();
570 m_vfit2_mks = (wks.p()).m();
571 m_vfit2_chi = vtxfit2->chisq();
572 m_vfit2_ct = vtxfit2->ctau();
573 m_vfit2_dl = vtxfit2->decayLength();
574 m_vfit2_dle = vtxfit2->decayLengthError();
575
576 pipKalTrk = pi1KalTrk ;
577 pimKalTrk = pi2KalTrk ;
578 piKalTrk = pi3KalTrk ;
579 kKalTrk = k1KalTrk ;
580
581 }
582 }
583 }
584
585 if (! okloop ) {
586 return StatusCode::SUCCESS;
587 }
588
593
594 WTrackParameter wpip(xmass[2], pipKalTrk->getZHelix(), pipKalTrk->getZError()); //pi+ from Ks
595 WTrackParameter wpim(xmass[2], pimKalTrk->getZHelix(), pimKalTrk->getZError()); //pi- from Ks
596
597 WTrackParameter wpi(xmass[2], piKalTrk->getZHelix(), piKalTrk->getZError());
598 WTrackParameter wk(xmass[3], kKalTrk->getZHelixK(), kKalTrk->getZErrorK());
599
600 //
601 // check good charged track's infomation
602 int ii ;
603 for(int j=0; j<m_ngch; j++){
604 m_charge[j] = 9999.0;
605 m_vx0[j] = 9999.0;
606 m_vy0[j] = 9999.0;
607 m_vz0[j] = 9999.0;
608 m_vr0[j] = 9999.0;
609
610 m_vx[j] = 9999.0;
611 m_vy[j] = 9999.0;
612 m_vz[j] = 9999.0;
613 m_vr[j] = 9999.0;
614
615 m_px[j] = 9999.0;
616 m_py[j] = 9999.0;
617 m_pz[j] = 9999.0;
618 m_p[j] = 9999.0;
619 m_cost[j] = 9999.0;
620
621 m_probPH[j] = 9999.0;
622 m_normPH[j] = 9999.0;
623 m_chie[j] = 9999.0;
624 m_chimu[j] = 9999.0;
625 m_chipi[j] = 9999.0;
626 m_chik[j] = 9999.0;
627 m_chip[j] = 9999.0;
628 m_ghit[j] = 9999.0;
629 m_thit[j] = 9999.0;
630
631 m_e_emc[j] = 9999.0;
632
633 m_qual_etof[j] = 9999.0;
634 m_tof_etof[j] = 9999.0;
635 m_te_etof[j] = 9999.0;
636 m_tmu_etof[j] = 9999.0;
637 m_tpi_etof[j] = 9999.0;
638 m_tk_etof[j] = 9999.0;
639 m_tp_etof[j] = 9999.0;
640
641 m_qual_btof1[j] = 9999.0;
642 m_tof_btof1[j] = 9999.0;
643 m_te_btof1[j] = 9999.0;
644 m_tmu_btof1[j] = 9999.0;
645 m_tpi_btof1[j] = 9999.0;
646 m_tk_btof1[j] = 9999.0;
647 m_tp_btof1[j] = 9999.0;
648
649 m_qual_btof2[j] = 9999.0;
650 m_tof_btof2[j] = 9999.0;
651 m_te_btof2[j] = 9999.0;
652 m_tmu_btof2[j] = 9999.0;
653 m_tpi_btof2[j] = 9999.0;
654 m_tk_btof2[j] = 9999.0;
655 m_tp_btof2[j] = 9999.0;
656
657 m_pidcode[j] = 9999.0;
658 m_pidprob[j] = 9999.0;
659 m_pidchiDedx[j] = 9999.0;
660 m_pidchiTof1[j] = 9999.0;
661 m_pidchiTof2[j] = 99999.0;
662 }
663
664 for(int i = 0; i < m_ngch; i++ ){
665
666 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
667 if(!(*itTrk)->isMdcTrackValid()) continue; // MDC information
668 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
669 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
670
671 if ( mdcKalTrk == pipKalTrk ) {
672 ii = 0 ;
674 }
675 if ( mdcKalTrk == pimKalTrk ) {
676 ii = 1 ;
678 }
679 if ( mdcKalTrk == piKalTrk ) {
680 ii = 2 ;
682 }
683 if ( mdcKalTrk == kKalTrk ) {
684 ii = 3 ;
686 }
687
688 m_charge[ii] = mdcTrk->charge();
689 double x0=mdcTrk->x();
690 double y0=mdcTrk->y();
691 double z0=mdcTrk->z();
692 double phi0=mdcTrk->helix(1);
693 double xv=xorigin.x();
694 double yv=xorigin.y();
695 double zv=xorigin.z();
696 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
697
698 m_vx0[ii] = x0-xv ;
699 m_vy0[ii] = y0-yv ;
700 m_vz0[ii] = z0-zv ;
701 m_vr0[ii] = rv ;
702
703 x0=mdcKalTrk->x();
704 y0=mdcKalTrk->y();
705 z0=mdcKalTrk->z();
706 rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
707
708 m_vx[ii] = x0-xv ;
709 m_vy[ii] = y0-yv ;
710 m_vz[ii] = z0-zv ;
711 m_vr[ii] = rv ;
712
713 m_px[ii] = mdcKalTrk->px();
714 m_py[ii] = mdcKalTrk->py();
715 m_pz[ii] = mdcKalTrk->pz();
716 m_p[ii] = mdcKalTrk->p();
717 m_cost[ii] = mdcKalTrk->pz()/mdcKalTrk->p();
718
719 double ptrk = mdcKalTrk->p() ;
720 if((*itTrk)->isMdcDedxValid()) { // DEDX information
721 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
722 m_probPH[ii]= dedxTrk->probPH();
723 m_normPH[ii]= dedxTrk->normPH();
724
725 m_chie[ii] = dedxTrk->chiE();
726 m_chimu[ii] = dedxTrk->chiMu();
727 m_chipi[ii] = dedxTrk->chiPi();
728 m_chik[ii] = dedxTrk->chiK();
729 m_chip[ii] = dedxTrk->chiP();
730 m_ghit[ii] = dedxTrk->numGoodHits();
731 m_thit[ii] = dedxTrk->numTotalHits();
732 }
733
734 if((*itTrk)->isEmcShowerValid()) {
735 RecEmcShower *emcTrk = (*itTrk)->emcShower();
736 m_e_emc[ii] = emcTrk->energy();
737 }
738
739 if((*itTrk)->isTofTrackValid()) { //TOF information
740 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
741
742 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
743
744 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
745 TofHitStatus *status = new TofHitStatus;
746 status->setStatus((*iter_tof)->status());
747
748 if(!(status->is_barrel())){//endcap
749 if( !(status->is_counter()) ) continue; // ?
750 if( status->layer()!=0 ) continue;//layer1
751 double path=(*iter_tof)->path(); // ?
752 double tof = (*iter_tof)->tof();
753 double ph = (*iter_tof)->ph();
754 double rhit = (*iter_tof)->zrhit();
755 double qual = 0.0 + (*iter_tof)->quality();
756 double cntr = 0.0 + (*iter_tof)->tofID();
757 double texp[5];
758 for(int j = 0; j < 5; j++) {
759 double gb = ptrk/xmass[j];
760 double beta = gb/sqrt(1+gb*gb);
761 texp[j] = path /beta/velc;
762 }
763
764 m_qual_etof[ii] = qual;
765 m_tof_etof[ii] = tof ;
766 }
767 else {//barrel
768 if( !(status->is_counter()) ) continue; // ?
769 if(status->layer()==1){ //layer1
770 double path=(*iter_tof)->path(); // ?
771 double tof = (*iter_tof)->tof();
772 double ph = (*iter_tof)->ph();
773 double rhit = (*iter_tof)->zrhit();
774 double qual = 0.0 + (*iter_tof)->quality();
775 double cntr = 0.0 + (*iter_tof)->tofID();
776 double texp[5];
777 for(int j = 0; j < 5; j++) {
778 double gb = ptrk/xmass[j];
779 double beta = gb/sqrt(1+gb*gb);
780 texp[j] = path /beta/velc;
781 }
782
783 m_qual_btof1[ii] = qual;
784 m_tof_btof1[ii] = tof ;
785 }
786
787 if(status->layer()==2){//layer2
788 double path=(*iter_tof)->path(); // ?
789 double tof = (*iter_tof)->tof();
790 double ph = (*iter_tof)->ph();
791 double rhit = (*iter_tof)->zrhit();
792 double qual = 0.0 + (*iter_tof)->quality();
793 double cntr = 0.0 + (*iter_tof)->tofID();
794 double texp[5];
795 for(int j = 0; j < 5; j++) {
796 double gb = ptrk/xmass[j];
797 double beta = gb/sqrt(1+gb*gb);
798 texp[j] = path /beta/velc;
799 }
800
801 m_qual_btof2[ii] = qual;
802 m_tof_btof2[ii] = tof ;
803 }
804 }
805 }
806 }
807 }
808
809 // Kinamatic Fit
811
812 double ecms = 3.097;
813 double chisq = 9999.;
814 m_4c_chi2 = 9999.;
815 m_4c_mks = 10.0;
816 m_4c_mkspi = 10.0;
817 m_4c_mksk = 10.0;
818 m_4c_mkpi = 10.0;
819 m_4c_ks_px = 10.0;
820 m_4c_ks_py = 10.0;
821 m_4c_ks_pz = 10.0;
822 m_4c_ks_p = 10.0;
823 m_4c_ks_cos= 10.0;
824
825 kmfit->init();
826 kmfit->AddTrack(0, wpi);
827 kmfit->AddTrack(1, wk);
828 kmfit->AddTrack(2, wks);
829 kmfit->AddFourMomentum(0, p_cms);
830 bool oksq = kmfit->Fit();
831 if(oksq) {
832 chisq = kmfit->chisq();
833
834 HepLorentzVector pk = kmfit->pfit(1);
835 HepLorentzVector pks = kmfit->pfit(2);
836 HepLorentzVector pkspi = kmfit->pfit(0) + kmfit->pfit(2);
837 HepLorentzVector pksk = kmfit->pfit(1) + kmfit->pfit(2);
838 HepLorentzVector pkpi = kmfit->pfit(0) + kmfit->pfit(1);
839
840 pks.boost(u_cms);
841 pkspi.boost(u_cms);
842 pksk.boost(u_cms);
843 pkpi.boost(u_cms);
844
845 m_4c_chi2 = chisq ;
846 m_4c_mks = pks.m();
847 m_4c_mkspi = pkspi.m();
848 m_4c_mksk = pksk.m();
849 m_4c_mkpi = pkpi.m();
850
851 m_4c_ks_px = pks.px() ;
852 m_4c_ks_py = pks.py() ;
853 m_4c_ks_pz = pks.pz() ;
854 m_4c_ks_p = (pks.vect()).mag() ;
855 m_4c_ks_cos = m_4c_ks_pz/m_4c_ks_p ;
856
857 }
858
859 chisq = 9999.;
860 m_chi2_fs4c = 9999.;
861 m_mks_fs4c = 10.0;
862 m_mkspi_fs4c = 10.0;
863 m_mksk_fs4c = 10.0;
864 m_mkpi_fs4c = 10.0;
865
866 kmfit->init();
867 kmfit->AddTrack(0, wpip);
868 kmfit->AddTrack(1, wpim);
869 kmfit->AddTrack(2, wpi);
870 kmfit->AddTrack(3, wk);
871 kmfit->AddFourMomentum(0, p_cms);
872 oksq = kmfit->Fit();
873 if(oksq) {
874 chisq = kmfit->chisq();
875
876 HepLorentzVector pks = kmfit->pfit(0) + kmfit->pfit(1);
877 HepLorentzVector pkspi = pks + kmfit->pfit(2);
878 HepLorentzVector pksk = pks + kmfit->pfit(3);
879 HepLorentzVector pkpi = kmfit->pfit(2) + kmfit->pfit(3);
880
881 pks.boost(u_cms);
882 pkspi.boost(u_cms);
883 pksk.boost(u_cms);
884 pkpi.boost(u_cms);
885
886 m_chi2_fs4c = chisq ;
887 m_mks_fs4c = pks.m();
888 m_mkspi_fs4c = pkspi.m();
889 m_mksk_fs4c = pksk.m();
890 m_mkpi_fs4c = pkpi.m();
891 }
892
893 // finale selection
894 if(chisq > 20) { return StatusCode::SUCCESS; } //4C chi2
895 if(m_vfit2_dl < 0.5) { return StatusCode::SUCCESS; } // Ks decay length
896 if(fabs(m_4c_mks-mks0) > 0.01) { return StatusCode::SUCCESS; } // Ks mass
897 if(m_4c_mkspi < 1.25) { return StatusCode::SUCCESS; } // Kspi mass
898
899 // DQA
900 TH1 *h(0);
901 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_dl", h).isSuccess()) {
902 h->Fill(m_vfit2_dl);
903 } else {
904 log << MSG::ERROR << "Couldn't retrieve hks_dl" << endreq;
905 }
906
907 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_m", h).isSuccess()) {
908 h->Fill(m_4c_mks);
909 } else {
910 log << MSG::ERROR << "Couldn't retrieve hks_m" << endreq;
911 }
912
913 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hkspi_m", h).isSuccess()) {
914 h->Fill(m_4c_mkspi);
915 } else {
916 log << MSG::ERROR << "Couldn't retrieve hkspi_m" << endreq;
917 }
918
919 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_p", h).isSuccess()) {
920 h->Fill(m_4c_ks_p);
921 } else {
922 log << MSG::ERROR << "Couldn't retrieve hks_p" << endreq;
923 }
924
925 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hkpi_m", h).isSuccess()) {
926 h->Fill(m_4c_mkpi);
927 } else {
928 log << MSG::ERROR << "Couldn't retrieve hkpi_m" << endreq;
929 }
930
931 ////////////////////////////////////////////////////////////
932 // DQA
933 // set tag and quality
934
935 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
936 if(m_npip==2 && m_npim==1){
937 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
938 (*(evtRecTrkCol->begin()+ipip[1]))->setPartId(2);
939 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
940 (*(evtRecTrkCol->begin()+ikm[0]))->setPartId(4);
941 }
942 if(m_npip==1 && m_npim==2){
943 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
944 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
945 (*(evtRecTrkCol->begin()+ipim[1]))->setPartId(2);
946 (*(evtRecTrkCol->begin()+ikp[0]))->setPartId(4);
947 }
948 // Quality: defined by whether dE/dx or TOF is used to identify particle
949 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
950 // 1 - only dE/dx (can be used for TOF calibration)
951 // 2 - only TOF (can be used for dE/dx calibration)
952 // 3 - Both dE/dx and TOF
953 if(m_npip==2 && m_npim==1){
954 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(1);
955 (*(evtRecTrkCol->begin()+ipip[1]))->setQuality(1);
956 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(1);
957 (*(evtRecTrkCol->begin()+ikm[0]))->setQuality(1);
958 }
959 if(m_npip==1 && m_npim==2){
960 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(1);
961 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(1);
962 (*(evtRecTrkCol->begin()+ipim[1]))->setQuality(1);
963 (*(evtRecTrkCol->begin()+ikp[0]))->setQuality(1);
964 }
965 // DQA
966 // Add the line below at the end of execute(), (before return)
967 //
968 setFilterPassed(true);
969 ////////////////////////////////////////////////////////////
970
971 m_tuple->write();
972
973 return StatusCode::SUCCESS;
974
975}
976
977// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
978StatusCode DQAKsKpi::finalize() {
979
980 MsgStream log(msgSvc(), name());
981 log << MSG::INFO << "in finalize()" << endmsg;
982 return StatusCode::SUCCESS;
983}
984
985
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
const double mks0
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
HepGeom::Point3D< double > HepPoint3D
Definition: DQAKsKpi.cxx:45
const double mks0
Definition: DQAKsKpi.cxx:59
const Hep3Vector u_cms
Definition: DQAKsKpi.cxx:67
std::vector< HepLorentzVector > Vp4
Definition: DQAKsKpi.cxx:63
const double xmass[5]
Definition: DQAKsKpi.cxx:60
const double velc
Definition: DQAKsKpi.cxx:61
const double mk
Definition: DQAKsKpi.cxx:58
const double mpi
Definition: DQAKsKpi.cxx:57
std::vector< int > Vint
Definition: DQAKsKpi.cxx:62
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
const double xmass[5]
Definition: Gam4pikp.cxx:50
const double velc
Definition: Gam4pikp.cxx:51
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
double sin(const BesAngle a)
double cos(const BesAngle a)
DQAKsKpi(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DQAKsKpi.cxx:73
StatusCode finalize()
Definition: DQAKsKpi.cxx:978
StatusCode execute()
Definition: DQAKsKpi.cxx:254
StatusCode initialize()
Definition: DQAKsKpi.cxx:88
const HepSymMatrix err() const
const HepVector helix() const
......
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
static KinematicFit * instance()
void AddFourMomentum(int number, HepLorentzVector p4)
double chiTof2(int n) const
static ParticleID * instance()
Definition: ParticleID.cxx:22
bool IsPidInfoValid() const
double chiTof1(int n) const
void calculate()
Definition: ParticleID.cxx:101
void init()
Definition: ParticleID.cxx:27
double chiDedx(int n) const
static SecondVertexFit * instance()
void setStatus(unsigned int status)
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
const HepPoint3D & pivot(void) const
returns pivot position.
const HepVector & a(void) const
returns helix parameters.
void init()
Definition: VertexFit.cxx:29
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition: VertexFit.cxx:89
static VertexFit * instance()
Definition: VertexFit.cxx:15
void BuildVirtualParticle(int number)
Definition: VertexFit.cxx:619
bool Fit()
Definition: VertexFit.cxx:301
const double ecms
Definition: inclkstar.cxx:42