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