BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
JsiLL.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 "GaudiKernel/Bootstrap.h"
10#include "GaudiKernel/ISvcLocator.h"
11
12#include "GaudiKernel/INTupleSvc.h"
13#include "GaudiKernel/NTuple.h"
14#include "GaudiKernel/ITHistSvc.h"
15//#include "GaudiKernel/IHistogramSvc.h"
16
17#include "EventModel/EventModel.h"
18#include "EventModel/Event.h"
19
20#include "EvtRecEvent/EvtRecEvent.h"
21#include "EvtRecEvent/EvtRecTrack.h"
22#include "DstEvent/TofHitStatus.h"
23#include "EventModel/EventHeader.h"
24
25#include "TH1F.h"
26// #include "TMath.h"
27#include "CLHEP/Vector/ThreeVector.h"
28#include "CLHEP/Vector/LorentzVector.h"
29#include "CLHEP/Vector/TwoVector.h"
30
31using CLHEP::Hep3Vector;
32using CLHEP::Hep2Vector;
33using CLHEP::HepLorentzVector;
34#include "CLHEP/Geometry/Point3D.h"
35
36#include "VertexFit/IVertexDbSvc.h"
37#include "VertexFit/KinematicFit.h"
38#include "VertexFit/VertexFit.h"
39#include "VertexFit/SecondVertexFit.h"
40#include "ParticleID/ParticleID.h"
41#include "TH1F.h"
42//
43#include "JsiLL/JsiLL.h"
44
45#ifndef ENABLE_BACKWARDS_COMPATIBILITY
46 typedef HepGeom::Point3D<double> HepPoint3D;
47#endif
48using CLHEP::HepLorentzVector;
49
50const double mpi = 0.13957;
51const double mk = 0.493677;
52const double mproton = 0.938272;
53const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
54const double velc = 299.792458; // tof path unit in mm
55typedef std::vector<int> Vint;
56typedef std::vector<HepLorentzVector> Vp4;
57
58/////////////////////////////////////////////////////////////////////////////
59
60JsiLL::JsiLL(const std::string& name, ISvcLocator* pSvcLocator) :
61 Algorithm(name, pSvcLocator) {
62
63 //Declare the properties
64 declareProperty("Vr0cut", m_vr0cut=5.0);
65 declareProperty("Vz0cut", m_vz0cut=20.0);
66 declareProperty("Vr1cut", m_vr1cut=1.0);
67 declareProperty("Vz1cut", m_vz1cut=5.0);
68 declareProperty("Vctcut", m_cthcut=0.93);
69 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
70 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
71}
72
73
74// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
75StatusCode JsiLL::initialize(){
76 MsgStream log(msgSvc(), name());
77
78 log << MSG::INFO << "in initialize()" << endmsg;
79 StatusCode status;
80
81 // DQA
82 // The first directory specifier must be "DQAFILE"
83 // The second is the control sample name, the first letter is in upper format. eg. "Rhopi"
84 NTuplePtr nt(ntupleSvc(), "DQAFILE/JsiLL");
85 if ( nt ) m_tuple = nt;
86 else {
87 m_tuple = ntupleSvc()->book("DQAFILE/JsiLL", CLID_ColumnWiseTuple, "JsiLL ntuple");
88 if( m_tuple ) {
89 status = m_tuple->addItem("runNo", m_runNo);
90 status = m_tuple->addItem("event", m_event);
91 status = m_tuple->addItem("chisq", m_chisq);
92 status = m_tuple->addItem("mLambda", m_mLambda);
93 status = m_tuple->addItem("mLambdabar", m_mLambdabar);
94 status = m_tuple->addItem("pLambda", m_pLambda);
95 status = m_tuple->addItem("pLambdabar", m_pLambdabar);
96
97 } else {
98 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
99 }
100 }
101
102 if(service("THistSvc", m_thsvc).isFailure()) {
103 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
104 return StatusCode::FAILURE;
105 }
106
107 // "DQAHist" is fixed
108 // "Rhopi" is control sample name, just as ntuple case.
109 TH1F* hrxy = new TH1F("Rxy", "Rxy distribution", 110, -1., 10.);
110 if(m_thsvc->regHist("/DQAHist/JsiLL/hrxy", hrxy).isFailure()) {
111 log << MSG::ERROR << "Couldn't register Rxy" << endreq;
112 }
113 TH1F* hz = new TH1F("z", "z distribution", 200, -100., 100.);
114 if(m_thsvc->regHist("/DQAHist/JsiLL/hz", hz).isFailure()) {
115 log << MSG::ERROR << "Couldn't register z" << endreq;
116 }
117
118 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
119 return StatusCode::SUCCESS;
120
121
122}
123
124// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
125StatusCode JsiLL::execute() {
126
127 MsgStream log(msgSvc(), name());
128 log << MSG::INFO << "in execute()" << endreq;
129
130 // DQA
131 // Add the line below at the beginning of execute()
132 //
133 setFilterPassed(false);
134
135 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
136 m_runNo=eventHeader->runNumber();
137 m_event=eventHeader->eventNumber();
138 log << MSG::DEBUG <<"run, evtnum = "
139 << m_runNo << " , "
140 << m_event <<endreq;
141
142
143 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
144 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
145 << evtRecEvent->totalCharged() << " , "
146 << evtRecEvent->totalNeutral() << " , "
147 << evtRecEvent->totalTracks() <<endreq;
148
149 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
150
151 if(evtRecEvent->totalNeutral()>100) {
152 return StatusCode::SUCCESS;
153 }
154
155 Vint iGood, iplus, iminus;
156 iGood.clear();
157 iplus.clear();
158 iminus.clear();
159 Vp4 ppip, ppim;
160 ppip.clear();
161 ppim.clear();
162
163 Hep3Vector xorigin(0,0,0);
164
165 IVertexDbSvc* vtxsvc;
166 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
167 if(vtxsvc->isVertexValid()){
168 double* dbv = vtxsvc->PrimaryVertex();
169 double* vv = vtxsvc->SigmaPrimaryVertex();
170 xorigin.setX(dbv[0]);
171 xorigin.setY(dbv[1]);
172 xorigin.setZ(dbv[2]);
173 }
174
175 int nCharge = 0;
176 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
177 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
178 if(!(*itTrk)->isMdcTrackValid()) continue;
179 if (!(*itTrk)->isMdcKalTrackValid()) continue;
180 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
181
182 double pch =mdcTrk->p();
183 double x0 =mdcTrk->x();
184 double y0 =mdcTrk->y();
185 double z0 =mdcTrk->z();
186 double phi0=mdcTrk->helix(1);
187 double xv=xorigin.x();
188 double yv=xorigin.y();
189 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
190
191 if(fabs(z0) >= m_vz0cut) continue;
192 if(Rxy >= m_vr0cut) continue;
193// if(fabs(m_Vct)>=m_cthcut) continue;
194
195 // DQA
196 TH1 *h(0);
197 if (m_thsvc->getHist("/DQAHist/JsiLL/hrxy", h).isSuccess()) {
198 h->Fill(Rxy);
199 } else {
200 log << MSG::ERROR << "Couldn't retrieve hrxy" << endreq;
201 }
202 if (m_thsvc->getHist("/DQAHist/JsiLL/hz", h).isSuccess()) {
203 h->Fill(z0);
204 } else {
205 log << MSG::ERROR << "Couldn't retrieve hz" << endreq;
206 }
207
208// iGood.push_back((*itTrk)->trackId());
209 iGood.push_back(i);
210 nCharge += mdcTrk->charge();
211 if (mdcTrk->charge() > 0) {
212 iplus.push_back(i);
213 } else {
214 iminus.push_back(i);
215 }
216 }
217
218 //
219 // Finish Good Charged Track Selection
220 //
221 int nGood = iGood.size();
222 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
223 if((nGood != 4)||(nCharge!=0)){
224 return StatusCode::SUCCESS;
225 }
226
227
228// for(int i = 0; i < nGood; i++) {
229// EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
230// if(!(*itTrk)->isMdcKalTrackValid())
231// return StatusCode::SUCCESS;
232// }
233
234 int pidp0 = 5, pidp1 = 3, pidm0 = 5, pidm1 = 3; //PID
235
236 RecMdcKalTrack *itTrkp = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
237 RecMdcKalTrack *itTrkpip = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
238 RecMdcKalTrack *itTrkpb = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
239 RecMdcKalTrack *itTrkpim = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
240
241// p() is not filled during reconstruction
242 double tmppp, tmpppb, tmpppip, tmpppim;
243 tmppp = sqrt(itTrkp->px()*itTrkp->px() + itTrkp->py()*itTrkp->py() + itTrkp->pz()*itTrkp->pz());
244 tmpppb = sqrt(itTrkpb->px()*itTrkpb->px() + itTrkpb->py()*itTrkpb->py() + itTrkpb->pz()*itTrkpb->pz());
245 tmpppip = sqrt(itTrkpip->px()*itTrkpip->px() + itTrkpip->py()*itTrkpip->py() + itTrkpip->pz()*itTrkpip->pz());
246 tmpppim = sqrt(itTrkpim->px()*itTrkpim->px() + itTrkpim->py()*itTrkpim->py() + itTrkpim->pz()*itTrkpim->pz());
247
248
249 if ( tmppp < tmpppip ) {
250 itTrkp = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
251 itTrkpip = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
252 pidp0 = 3;
253 pidp1 = 5;
254 double tmp;
255 tmp = tmppp;
256 tmppp = tmpppip;
257 tmpppip = tmp;
258 }
259 if ( tmpppb < tmpppim ) {
260 itTrkpb = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
261 itTrkpim = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
262 pidm0 = 3;
263 pidm1 = 5;
264 double tmp;
265 tmp = tmpppb;
266 tmpppb = tmpppim;
267 tmpppim = tmp;
268 }
269
270 if ( tmppp < 0.7 || tmppp >1.1) return StatusCode::SUCCESS;
271 if ( tmpppb < 0.7 || tmpppb >1.1) return StatusCode::SUCCESS;
272 if ( tmpppip > 0.35) return StatusCode::SUCCESS;
273 if ( tmpppim > 0.35) return StatusCode::SUCCESS;
274
275
276
277 itTrkp->setPidType(RecMdcKalTrack::proton);
278 itTrkpb->setPidType(RecMdcKalTrack::proton);
279 itTrkpip->setPidType(RecMdcKalTrack::pion);
280 itTrkpim->setPidType(RecMdcKalTrack::pion);
281
282
283 m_chisq = 9999.0;
284 // Vertex Fit
285 HepPoint3D vx(0., 0., 0.);
286 HepSymMatrix Evx(3, 0);
287 double bx = 1E+6;
288 double by = 1E+6;
289 double bz = 1E+6;
290 Evx[0][0] = bx*bx;
291 Evx[1][1] = by*by;
292 Evx[2][2] = bz*bz;
293
294 VertexParameter vxpar;
295 vxpar.setVx(vx);
296 vxpar.setEvx(Evx);
297
298 VertexFit* vtxfita0 = VertexFit::instance();
299 SecondVertexFit *vtxfita = SecondVertexFit::instance();
300 VertexFit* vtxfitb0 = VertexFit::instance();
301 SecondVertexFit *vtxfitb = SecondVertexFit::instance();
302 VertexFit* vtxfit = VertexFit::instance();
303
304 WTrackParameter wpip = WTrackParameter(mpi, itTrkpip->getZHelix(), itTrkpip->getZError());
305 WTrackParameter wpim = WTrackParameter(mpi, itTrkpim->getZHelix(), itTrkpim->getZError());
306 WTrackParameter wp = WTrackParameter(mproton, itTrkp->getZHelixP(), itTrkp->getZErrorP());
307 WTrackParameter wpb = WTrackParameter(mproton, itTrkpb->getZHelixP(), itTrkpb->getZErrorP());
308
309 vtxfita0->init();
310 vtxfita0->AddTrack(0, wp);
311 vtxfita0->AddTrack(1, wpim);
312 vtxfita0->AddVertex(0, vxpar, 0, 1);
313 if(!vtxfita0->Fit(0)) return StatusCode::SUCCESS;
314 vtxfita0->Swim(0);
315 vtxfita0->BuildVirtualParticle(0);
316 vtxfita->init();
317 vtxfita->AddTrack(0, vtxfita0->wVirtualTrack(0));
318 vtxfita->setVpar(vtxfita0->vpar(0));
319 if(!vtxfita->Fit()) return StatusCode::SUCCESS;
320
321 WTrackParameter wLambda = vtxfita->wpar();
322
323 vtxfitb0->init();
324 vtxfitb0->AddTrack(0, wpb);
325 vtxfitb0->AddTrack(1, wpip);
326 vtxfitb0->AddVertex(0, vxpar, 0, 1);
327 if(!vtxfitb0->Fit(0)) return StatusCode::SUCCESS;
328 vtxfitb0->Swim(0);
329 vtxfitb0->BuildVirtualParticle(0);
330
331 vtxfitb->init();
332 vtxfitb->AddTrack(0, vtxfitb0->wVirtualTrack(0));
333 vtxfitb->setVpar(vtxfitb0->vpar(0));
334 if(!vtxfitb->Fit()) return StatusCode::SUCCESS;
335
336 WTrackParameter wLambdabar = vtxfitb->wpar();
337
338 vtxfit->init();
339 vtxfit->AddTrack(0, wLambda);
340 vtxfit->AddTrack(1, wLambdabar);
341 vtxfit->AddVertex(0, vxpar,0, 1);
342 if(!vtxfit->Fit(0)) return StatusCode::SUCCESS;
343 vtxfit->Swim(0);
344 WTrackParameter wwLambda = vtxfit->wtrk(0);
345 WTrackParameter wwLambdabar = vtxfit->wtrk(1);
346
347
348 // Kinamatic Fit
349
350 KinematicFit* kmfit = KinematicFit::instance();
351
352
353 HepLorentzVector ecms(0.034065,0.0,0.0,3.0969);
354 const Hep3Vector u_cms = -ecms.boostVector();
355
356 kmfit->init();
357 kmfit->AddTrack(0, wwLambda);
358 kmfit->AddTrack(1, wwLambdabar);
359 kmfit->AddFourMomentum(0, ecms);
360
361 if(!kmfit->Fit()) return StatusCode::SUCCESS;
362 m_chisq = kmfit->chisq();
363 if(m_chisq > 50) return StatusCode::SUCCESS;
364 HepLorentzVector kmf_pLambda = kmfit->pfit(0);
365 HepLorentzVector kmf_pLambdabar = kmfit->pfit(1);
366
367 kmf_pLambda.boost(u_cms);
368 kmf_pLambdabar.boost(u_cms);
369 m_mLambda = kmf_pLambda.m();
370 m_mLambdabar = kmf_pLambdabar.m();
371 m_pLambda = kmf_pLambda.rho();
372 m_pLambdabar = kmf_pLambdabar.rho();
373
374 if(fabs(m_mLambda-1.1157)>0.003||fabs(m_mLambdabar-1.1157)>0.003) return StatusCode::SUCCESS;
375 // finale selection
376
377 m_tuple->write();
378// return StatusCode::SUCCESS;
379
380 ////////////////////////////////////////////////////////////
381 // DQA
382 // set tag and quality
383
384 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
385 (*(evtRecTrkCol->begin()+iplus[0]))->setPartId(pidp0);
386 (*(evtRecTrkCol->begin()+iplus[1]))->setPartId(pidp1);
387 (*(evtRecTrkCol->begin()+iminus[0]))->setPartId(pidm0);
388 (*(evtRecTrkCol->begin()+iminus[1]))->setPartId(pidm1);
389 // Quality: defined by whether dE/dx or TOF is used to identify particle
390 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
391 // 1 - only dE/dx (can be used for TOF calibration)
392 // 2 - only TOF (can be used for dE/dx calibration)
393 // 3 - Both dE/dx and TOF
394 (*(evtRecTrkCol->begin()+iplus[0]))->setQuality(0);
395 (*(evtRecTrkCol->begin()+iplus[1]))->setQuality(0);
396 (*(evtRecTrkCol->begin()+iminus[0]))->setQuality(0);
397 (*(evtRecTrkCol->begin()+iminus[1]))->setQuality(0);
398
399 // DQA
400 // Add the line below at the end of execute(), (before return)
401 //
402 setFilterPassed(true);
403 ////////////////////////////////////////////////////////////
404
405 return StatusCode::SUCCESS;
406
407}
408
409// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
410StatusCode JsiLL::finalize() {
411
412 MsgStream log(msgSvc(), name());
413 log << MSG::INFO << "in finalize()" << endmsg;
414 return StatusCode::SUCCESS;
415}
416
417
u_cms
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
Vp4
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
mpi
const double mpi
Definition: Gam4pikp.cxx:47
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
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: JsiLL.cxx:46
mproton
const double mproton
Definition: JsiLL.cxx:52
Vp4
std::vector< HepLorentzVector > Vp4
Definition: JsiLL.cxx:56
xmass
const double xmass[5]
Definition: JsiLL.cxx:53
velc
const double velc
Definition: JsiLL.cxx:54
mk
const double mk
Definition: JsiLL.cxx:51
mpi
const double mpi
Definition: JsiLL.cxx:50
Vint
std::vector< int > Vint
Definition: JsiLL.cxx:55
mproton
const double mproton
Definition: PipiJpsi.cxx:50
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
JsiLL::finalize
StatusCode finalize()
Definition: JsiLL.cxx:410
JsiLL::execute
StatusCode execute()
Definition: JsiLL.cxx:125
JsiLL::JsiLL
JsiLL(const std::string &name, ISvcLocator *pSvcLocator)
Definition: JsiLL.cxx:60
JsiLL::initialize
StatusCode initialize()
Definition: JsiLL.cxx:75
KinematicFit::instance
static KinematicFit * instance()
Definition: KinematicFit.cxx:20
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition: KinematicFit.cxx:354
SecondVertexFit::instance
static SecondVertexFit * instance()
Definition: SecondVertexFit.cxx:45
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
VertexFit::init
void init()
Definition: VertexFit.cxx:29
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition: VertexFit.cxx:89
VertexFit::instance
static VertexFit * instance()
Definition: VertexFit.cxx:15
VertexFit::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition: VertexFit.cxx:619
VertexFit::Fit
bool Fit()
Definition: VertexFit.cxx:301
ecms
const double ecms
Definition: inclkstar.cxx:42