BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
LambdaReconstruction.cxx
Go to the documentation of this file.
1//
2// Lambda -> p+ pi- Reconstruction
3// Anti-Lambda -> p- pi+
4//
5
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/AlgFactory.h"
8#include "GaudiKernel/ISvcLocator.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "GaudiKernel/IDataProviderSvc.h"
11#include "GaudiKernel/PropertyMgr.h"
12
13#include "EventModel/Event.h"
14#include "EvtRecEvent/EvtRecEvent.h"
15#include "EvtRecEvent/EvtRecTrack.h"
16#include "EvtRecEvent/EvtRecVeeVertex.h"
17#include "EventModel/EventHeader.h"
18#include "EventModel/Event.h"
19#include "EventModel/EventModel.h"
20#include "McTruth/McParticle.h"
21#include "ParticleID/ParticleID.h"
22#include "VertexFit/VertexFit.h"
23#include "VertexFitRefine/VertexFitRefine.h"
24#include "VertexFit/SecondVertexFit.h"
25#include "VeeVertexAlg/LambdaReconstruction.h"
26#include <vector>
27
28typedef std::vector<int> Vint;
29const double mpi = 0.13957;
30const double mp = 0.938272;
31//*******************************************************************************************
32DECLARE_COMPONENT(LambdaReconstruction)
33LambdaReconstruction::LambdaReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
34Algorithm(name, pSvcLocator) {
35 //Declare the properties
36
37 declareProperty("Vz0Cut", m_vzCut = 20.0);
38 declareProperty("CosThetaCut", m_cosThetaCut=0.93);
39 declareProperty("ChisqCut", m_chisqCut = 100.0);
40 declareProperty("UseVFRefine", m_useVFrefine = false);
41 /*
42 declareProperty("PidUseDedx", m_useDedx = true);
43 declareProperty("PidUseTof1", m_useTof1 = true);
44 declareProperty("PidUseTof2", m_useTof2 = true);
45 declareProperty("PidUseTofE", m_useTofE = false);
46 declareProperty("PidUseTofQ", m_useTofQ = false);
47 declareProperty("PidUseEmc", m_useEmc = false);
48 declareProperty("PidProbCut", m_PidProbCut= 0.001);
49 declareProperty("MassRangeLower", m_massRangeLower = 0.45);
50 declareProperty("MassRangeHigher", m_massRangeHigher = 0.55);
51 */
52}
53
54//******************************************************************************************
55StatusCode LambdaReconstruction::initialize() {
56 MsgStream log(msgSvc(), name());
57 log << MSG::INFO << "in initialize()" <<endreq;
58
59 return StatusCode::SUCCESS;
60}
61
62//********************************************************************************************
63StatusCode LambdaReconstruction::finalize() {
64 MsgStream log(msgSvc(), name());
65 log << MSG::INFO << "in finalize()" << endreq;
66
67 return StatusCode::SUCCESS;
68}
69
70StatusCode LambdaReconstruction::registerEvtRecVeeVertexCol(
71 EvtRecVeeVertexCol* aNewEvtRecVeeVertexCol, MsgStream& log) {
72 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecVeeVertexCol",
73 aNewEvtRecVeeVertexCol);
74 if (sc != StatusCode::SUCCESS) {
75 log << MSG::FATAL << "Could not register EvtRecVeeVertexCol in TDS!" << endreq;
76 }
77 return sc;
78}
79
80//*********************************************************************************************
81StatusCode LambdaReconstruction::execute() {
82 MsgStream log(msgSvc(), name());
83 log << MSG::INFO << "in execute()" << endreq;
84
85 StatusCode sc;
86
87 //////////////////
88 // Read REC data
89 /////////////////
90 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
91 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
92 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
93 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
94 << " " << eventHeader->eventNumber() << endreq;
95 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
96 << recEvent->totalCharged() << " , "
97 << recEvent->totalNeutral() << " , "
98 << recEvent->totalTracks() <<endreq;
99 int evtNo = eventHeader->eventNumber();
100 //if (eventHeader->eventNumber() % 1000 == 0) {
101 // cout << "Event number = " << evtNo << endl;
102 //}
103
104 SmartDataPtr<EvtRecVeeVertexCol> veeVertexCol(eventSvc(),
105 EventModel::EvtRec::EvtRecVeeVertexCol);
106 if (!veeVertexCol) {
107 veeVertexCol = new EvtRecVeeVertexCol;
108 sc = registerEvtRecVeeVertexCol(veeVertexCol, log);
109 if (sc != StatusCode::SUCCESS) {
110 return sc;
111 }
112 }
113
114 //////////////////////////////////////////////
115 // No PID : identify positive and negative
116 ///////////////////////////////////////////////
117 Vint icp, icm, iGood;
118 for (unsigned int i = 0; i < recEvent->totalCharged(); i++) {
119 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
120 if(!(*itTrk)->isMdcTrackValid()) continue;
121 if(!(*itTrk)->isMdcKalTrackValid()) continue;
122 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
123 if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
124 if (fabs(mdcTrk->z()) >= m_vzCut) continue;
125 iGood.push_back(i);
126 if (mdcTrk->charge() > 0) icp.push_back(i);
127 if (mdcTrk->charge() < 0) icm.push_back(i);
128 }
129 // Cut on good charged tracks
130 //if (iGood.size() < 2) return StatusCode::SUCCESS;
131 if (icp.size() < 1 || icm.size() < 1) return StatusCode::SUCCESS;
132
133 ///////////////////////////////////
134 // Vertex Fit
135 //////////////////////////////////
136 HepPoint3D vx(0., 0., 0.);
137 HepSymMatrix Evx(3, 0);
138 double bx = 1E+6;
139 double by = 1E+6;
140 double bz = 1E+6;
141 Evx[0][0] = bx*bx;
142 Evx[1][1] = by*by;
143 Evx[2][2] = bz*bz;
144
145 //VertexFit *vtxfit0 = VertexFit::instance();
146
147 // For Lambda reconstruction
148 for(unsigned int i1 = 0; i1 < icp.size(); i1++) {
149 int ip1 = icp[i1];
150 RecMdcKalTrack *ppKalTrk = (*(recTrackCol->begin()+ip1))->mdcKalTrack();
151 RecMdcKalTrack::setPidType(RecMdcKalTrack::proton);
152 WTrackParameter wpptrk(mp, ppKalTrk->helix(), ppKalTrk->err());
153
154 for(unsigned int i2 = 0; i2 < icm.size(); i2++) {
155 int ip2 = icm[i2];
156 RecMdcKalTrack *pimKalTrk = (*(recTrackCol->begin()+ip2))->mdcKalTrack();
157 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
158 WTrackParameter wpimtrk(mpi, pimKalTrk->helix(), pimKalTrk->err());
159
160 if(m_useVFrefine){
161 VertexParameter vxpar;
162 vxpar.setVx(vx);
163 vxpar.setEvx(Evx);
164
165 VertexFitRefine *vtxfit0 = VertexFitRefine::instance();
166 vtxfit0->init();
167 //vtxfit0->AddTrack(0, wpptrk);
168 //vtxfit0->AddTrack(1, wpimtrk);
169 vtxfit0->AddTrack(0, ppKalTrk, RecMdcKalTrack::proton);
170 vtxfit0->AddTrack(1, pimKalTrk, RecMdcKalTrack::pion);
171 vtxfit0->AddVertex(0, vxpar, 0, 1);
172
173 bool fitok = vtxfit0->Fit();
174 if(!fitok) continue;
175
176 //if(!(vtxfit0->Fit(0))) continue;
177 //vtxfit0->BuildVirtualParticle(0);
178 if(vtxfit0->chisq(0) > m_chisqCut) continue;
179 WTrackParameter wLamb = vtxfit0->wVirtualTrack(0); //FIXME
180 std::pair<int, int> pair;
181 pair.first = 5;
182 pair.second = 3;
183
184 EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
185 EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
186
187 EvtRecVeeVertex* LambdaVertex = new EvtRecVeeVertex;
188 LambdaVertex->setVertexId(3122);
189 LambdaVertex->setVertexType(1);
190 LambdaVertex->setChi2(vtxfit0->chisq(0));
191 LambdaVertex->setNdof(1);
192 LambdaVertex->setMass(wLamb.p().m());
193 LambdaVertex->setW(wLamb.w());
194 LambdaVertex->setEw(wLamb.Ew());
195 LambdaVertex->setPair(pair);
196 LambdaVertex->setNCharge(0);
197 LambdaVertex->setNTracks(2);
198 LambdaVertex->addDaughter(track0, 0);
199 LambdaVertex->addDaughter(track1, 1);
200 veeVertexCol->push_back(LambdaVertex);
201
202 /*
203 cout << "============= After Vertex fitting ===========" << endl;
204 cout << "Event number = " << evtNo << endl;
205 cout << "chi2 = " << vtxfit0->chisq(0) << endl;
206 cout << "mass = " << wLamb.p().m() << endl;
207 cout << "w = " << wLamb.w() << endl;
208 cout << "Ew = " << wLamb.Ew() << endl;
209 cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
210 cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
211
212 }else{
213 VertexParameter vxpar;
214 vxpar.setVx(vx);
215 vxpar.setEvx(Evx);
216
217 VertexFit *vtxfit0 = VertexFit::instance();
218 vtxfit0->init();
219 vtxfit0->AddTrack(0, wpptrk);
220 vtxfit0->AddTrack(1, wpimtrk);
221 vtxfit0->AddVertex(0, vxpar, 0, 1);
222 if(!(vtxfit0->Fit(0))) continue;
223 vtxfit0->BuildVirtualParticle(0);
224 if(vtxfit0->chisq(0) > m_chisqCut) continue;
225 WTrackParameter wLamb = vtxfit0->wVirtualTrack(0); //FIXME
226 std::pair<int, int> pair;
227 pair.first = 5;
228 pair.second = 3;
229
230 EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
231 EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
232
233 EvtRecVeeVertex* LambdaVertex = new EvtRecVeeVertex;
234 LambdaVertex->setVertexId(3122);
235 LambdaVertex->setVertexType(1);
236 LambdaVertex->setChi2(vtxfit0->chisq(0));
237 LambdaVertex->setNdof(1);
238 LambdaVertex->setMass(wLamb.p().m());
239 LambdaVertex->setW(wLamb.w());
240 LambdaVertex->setEw(wLamb.Ew());
241 LambdaVertex->setPair(pair);
242 LambdaVertex->setNCharge(0);
243 LambdaVertex->setNTracks(2);
244 LambdaVertex->addDaughter(track0, 0);
245 LambdaVertex->addDaughter(track1, 1);
246 veeVertexCol->push_back(LambdaVertex);
247
248 /*
249 cout << "============= After Vertex fitting ===========" << endl;
250 cout << "Event number = " << evtNo << endl;
251 cout << "chi2 = " << vtxfit0->chisq(0) << endl;
252 cout << "mass = " << wLamb.p().m() << endl;
253 cout << "w = " << wLamb.w() << endl;
254 cout << "Ew = " << wLamb.Ew() << endl;
255 cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
256 cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
257
258 }
259
260 }
261 }
262
263 // For Anti-Lambda reconstruction
264 for(unsigned int i1 = 0; i1 < icp.size(); i1++) {
265 int ip1 = icp[i1];
266 RecMdcKalTrack *pipKalTrk = (*(recTrackCol->begin()+ip1))->mdcKalTrack();
267 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
268 WTrackParameter wpiptrk(mpi, pipKalTrk->helix(), pipKalTrk->err());
269
270 for(unsigned int i2 = 0; i2 < icm.size(); i2++) {
271 int ip2 = icm[i2];
272 RecMdcKalTrack *pmKalTrk = (*(recTrackCol->begin()+ip2))->mdcKalTrack();
273 RecMdcKalTrack::setPidType(RecMdcKalTrack::proton);
274 WTrackParameter wpmtrk(mp, pmKalTrk->helix(), pmKalTrk->err());
275
276 if(m_useVFrefine){
277 VertexParameter vxpar;
278 vxpar.setVx(vx);
279 vxpar.setEvx(Evx);
280
281 VertexFitRefine *vtxfit0 = VertexFitRefine::instance();
282 vtxfit0->init();
283 //vtxfit0->AddTrack(0, wpiptrk);
284 //vtxfit0->AddTrack(1, wpmtrk);
285 vtxfit0->AddTrack(0, pipKalTrk, RecMdcKalTrack::pion);
286 vtxfit0->AddTrack(1, pmKalTrk, RecMdcKalTrack::proton);
287 vtxfit0->AddVertex(0, vxpar, 0, 1);
288
289 bool fitok = vtxfit0->Fit();
290 if(!fitok) continue;
291
292 //if(!(vtxfit0->Fit(0))) continue;
293 //vtxfit0->BuildVirtualParticle(0);
294 if(vtxfit0->chisq(0) > m_chisqCut) continue;
295 WTrackParameter wALamb = vtxfit0->wVirtualTrack(0); //FIXME
296 std::pair<int, int> pair;
297 pair.first = 3;
298 pair.second = 5;
299
300 EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
301 EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
302
303 EvtRecVeeVertex* ALambdaVertex = new EvtRecVeeVertex;
304 ALambdaVertex->setVertexId(-3122);
305 ALambdaVertex->setVertexType(1);
306 ALambdaVertex->setChi2(vtxfit0->chisq(0));
307 ALambdaVertex->setNdof(1);
308 ALambdaVertex->setMass(wALamb.p().m());
309 ALambdaVertex->setW(wALamb.w());
310 ALambdaVertex->setEw(wALamb.Ew());
311 ALambdaVertex->setPair(pair);
312 ALambdaVertex->setNCharge(0);
313 ALambdaVertex->setNTracks(2);
314 ALambdaVertex->addDaughter(track0, 0);
315 ALambdaVertex->addDaughter(track1, 1);
316 veeVertexCol->push_back(ALambdaVertex);
317
318 /*
319 cout << "============= After Vertex fitting ===========" << endl;
320 cout << "Event number = " << evtNo << endl;
321 cout << "chi2 = " << vtxfit0->chisq(0) << endl;
322 cout << "mass = " << wALamb.p().m() << endl;
323 cout << "w = " << wALamb.w() << endl;
324 cout << "Ew = " << wALamb.Ew() << endl;
325 cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
326 cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
327
328 }else{
329 VertexParameter vxpar;
330 vxpar.setVx(vx);
331 vxpar.setEvx(Evx);
332
333 VertexFit *vtxfit0 = VertexFit::instance();
334 vtxfit0->init();
335 vtxfit0->AddTrack(0, wpiptrk);
336 vtxfit0->AddTrack(1, wpmtrk);
337 vtxfit0->AddVertex(0, vxpar, 0, 1);
338 if(!(vtxfit0->Fit(0))) continue;
339 vtxfit0->BuildVirtualParticle(0);
340 if(vtxfit0->chisq(0) > m_chisqCut) continue;
341 WTrackParameter wALamb = vtxfit0->wVirtualTrack(0); //FIXME
342 std::pair<int, int> pair;
343 pair.first = 3;
344 pair.second = 5;
345
346 EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
347 EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
348
349 EvtRecVeeVertex* ALambdaVertex = new EvtRecVeeVertex;
350 ALambdaVertex->setVertexId(-3122);
351 ALambdaVertex->setVertexType(1);
352 ALambdaVertex->setChi2(vtxfit0->chisq(0));
353 ALambdaVertex->setNdof(1);
354 ALambdaVertex->setMass(wALamb.p().m());
355 ALambdaVertex->setW(wALamb.w());
356 ALambdaVertex->setEw(wALamb.Ew());
357 ALambdaVertex->setPair(pair);
358 ALambdaVertex->setNCharge(0);
359 ALambdaVertex->setNTracks(2);
360 ALambdaVertex->addDaughter(track0, 0);
361 ALambdaVertex->addDaughter(track1, 1);
362 veeVertexCol->push_back(ALambdaVertex);
363
364 /*
365 cout << "============= After Vertex fitting ===========" << endl;
366 cout << "Event number = " << evtNo << endl;
367 cout << "chi2 = " << vtxfit0->chisq(0) << endl;
368 cout << "mass = " << wALamb.p().m() << endl;
369 cout << "w = " << wALamb.w() << endl;
370 cout << "Ew = " << wALamb.Ew() << endl;
371 cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
372 cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
373
374 }
375
376 }
377 }
378
379
380
381 return StatusCode::SUCCESS;
382}
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
EvtRecVeeVertexCol
ObjectVector< EvtRecVeeVertex > EvtRecVeeVertexCol
Definition: EvtRecVeeVertex.h:107
mpi
const double mpi
Definition: Gam4pikp.cxx:47
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
mp
const double mp
Definition: LambdaReconstruction.cxx:30
mpi
const double mpi
Definition: LambdaReconstruction.cxx:29
Vint
std::vector< int > Vint
Definition: LambdaReconstruction.cxx:28
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
DstMdcKalTrack::helix
const HepVector & helix() const
Definition: DstMdcKalTrack.h:97
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition: DstMdcKalTrack.h:59
DstMdcKalTrack::err
const HepSymMatrix & err() const
Definition: DstMdcKalTrack.h:98
DstMdcTrack::theta
const double theta() const
Definition: DstMdcTrack.h:59
DstMdcTrack::charge
const int charge() const
Definition: DstMdcTrack.h:53
DstMdcTrack::z
const double z() const
Definition: DstMdcTrack.h:63
EvtRecVeeVertex::setVertexType
void setVertexType(int vtxType)
Definition: EvtRecVeeVertex.h:52
EvtRecVeeVertex::setChi2
void setChi2(double chi2)
Definition: EvtRecVeeVertex.h:53
EvtRecVeeVertex::addDaughter
void addDaughter(const SmartRef< EvtRecTrack > &track, int i)
Definition: EvtRecVeeVertex.h:97
EvtRecVeeVertex::setEw
void setEw(const HepSymMatrix &Ew)
Definition: EvtRecVeeVertex.h:57
EvtRecVeeVertex::setMass
void setMass(double mass)
Definition: EvtRecVeeVertex.h:55
EvtRecVeeVertex::setNdof
void setNdof(int ndof)
Definition: EvtRecVeeVertex.h:54
EvtRecVeeVertex::setNTracks
void setNTracks(int nTracks)
Definition: EvtRecVeeVertex.h:60
EvtRecVeeVertex::setVertexId
void setVertexId(int vtxId)
Definition: EvtRecVeeVertex.h:51
EvtRecVeeVertex::setPair
void setPair(const std::pair< int, int > &pair)
Definition: EvtRecVeeVertex.h:58
EvtRecVeeVertex::setW
void setW(const HepVector &w)
Definition: EvtRecVeeVertex.h:56
EvtRecVeeVertex::setNCharge
void setNCharge(int nCharge)
Definition: EvtRecVeeVertex.h:59
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
VertexFitRefine::wVirtualTrack
WTrackParameter wVirtualTrack(int n) const
Definition: VertexFitRefine.cc:177
VertexFitRefine::chisq
double chisq() const
Definition: VertexFitRefine.cc:117
VertexFitRefine::AddTrack
void AddTrack(const int index, RecMdcKalTrack *p, const RecMdcKalTrack::PidType pid)
Definition: VertexFitRefine.cxx:76
VertexFitRefine::instance
static VertexFitRefine * instance()
Definition: VertexFitRefine.cxx:34
VertexFitRefine::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition: VertexFitRefine.cc:17
VertexFit::chisq
double chisq() const
Definition: VertexFit.h:66
VertexFit::wVirtualTrack
WTrackParameter wVirtualTrack(int n) const
Definition: VertexFit.h:92
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
VertexParameter::setEvx
void setEvx(const HepSymMatrix &eVx)
Definition: VertexParameter.h:29
VertexParameter::setVx
void setVx(const HepPoint3D &vx)
Definition: VertexParameter.h:54
WTrackParameter::p
HepLorentzVector p() const
Definition: WTrackParameter.h:69
WTrackParameter::w
HepVector w() const
Definition: WTrackParameter.h:67
WTrackParameter::Ew
HepSymMatrix Ew() const
Definition: WTrackParameter.h:68
mp
const double mp
Definition: incllambda.cxx:45
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
EventModel::EvtRec::EvtRecVeeVertexCol
_EXTERN_ std::string EvtRecVeeVertexCol
Definition: EventModel.h:121
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:117