BOSS 7.0.1
BESIII Offline Software System
Loading...
Searching...
No Matches
PrimaryVertex.cxx
Go to the documentation of this file.
1#include "CLHEP/Vector/LorentzVector.h"
2#include "CLHEP/Geometry/Point3D.h"
3#ifndef ENABLE_BACKWARDS_COMPATIBILITY
4 typedef HepGeom::Point3D<double> HepPoint3D;
5#endif
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/SmartDataPtr.h"
8#include "GaudiKernel/IDataManagerSvc.h"
9#include "EventModel/EventModel.h"
10#include "EventModel/Event.h"
11#include "EventModel/EventHeader.h"
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14#include "EvtRecEvent/EvtRecPrimaryVertex.h"
15#include "VertexFit/KinematicFit.h"
16#include "VertexFit/VertexFit.h"
17#include "VertexFit/HTrackParameter.h"
18#include "VertexFit/KalmanVertexFit.h"
19#include "VertexFit/BField.h"
20#include "VertexFit/FastVertexFit.h"
21#include "ParticleID/ParticleID.h"
22#include "PrimaryVertexAlg/PrimaryVertex.h"
23
24#include <assert.h>
25#include "TMath.h"
26#include "TH2D.h"
27#include "TH1D.h"
28#include "TF1.h"
29#include <map>
30#include <iostream>
31#include <fstream>
32
33using CLHEP::HepLorentzVector;
34using namespace std;
35
36const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272}; //GeV
37const double ecms = 3.097;
38const double mpi0 = 0.134976;
39const double momega = 0.78265;
40
41typedef std::vector<int> Vint;
42typedef std::vector<HepLorentzVector> Vp4;
43
44//*******************************************************************************
45PrimaryVertex::PrimaryVertex(const std::string& name, ISvcLocator* pSvcLocator) :
46 Algorithm (name, pSvcLocator)
47{
48 // Declare the properties
49 declareProperty("Output", m_output = 0);
50 declareProperty("TrackNumberCut", m_trackNumberCut = 1);
51 declareProperty("CosThetaCut", m_cosThetaCut=0.93);
52 declareProperty("Vz0Cut", m_vz0Cut = 40.0);
53 // fit Method
54 declareProperty("FitMethod", m_fitMethod = 0);
55 // for global method
56 declareProperty("Chi2CutOfGlobal", m_globalChisqCut = 200.0);
57
58 // for Kalman method
59 declareProperty("ChisqCut", m_chisqCut = 200.0);
60 declareProperty("TrackIteration", m_trackIteration = 5);
61 declareProperty("VertexIteration", m_vertexIteration = 5);
62 declareProperty("Chi2CutforTrkIter", m_chi2CutforTrkIter = 0.1);
63 declareProperty("FreedomCut", m_freedomCut = 1);
64 declareProperty("Chi2CutforSmooth", m_chi2CutforSmooth = 200.0);
65 //declareProperty("MinDistance", m_minDistance = 100.0);
66 //declareProperty("MinPointX", m_minPointX = 0.1);
67 //declareProperty("MinPointY", m_minPointY = 0.1);
68}
69
70//*******************************************************************************
71StatusCode PrimaryVertex::initialize()
72{
73 MsgStream log(msgSvc(), name());
74 log << MSG::INFO << "in initialize()" << endmsg;
75
76 for(int i = 0; i < 15; i++){
77 m_sel_number[i] = 0;
78 }
79
80 StatusCode status;
81
82 if (m_output == 1) {
83 NTuplePtr nt1(ntupleSvc(), "FILE999/glbvtx");
84 if(nt1) m_tuple1 = nt1;
85 else {
86 m_tuple1 = ntupleSvc()->book ("FILE999/glbvtx", CLID_ColumnWiseTuple, "global vertex");
87 if(m_tuple1) {
88 status = m_tuple1->addItem("gvx", m_gvx);
89 status = m_tuple1->addItem("gvy", m_gvy);
90 status = m_tuple1->addItem("gvz", m_gvz);
91 status = m_tuple1->addItem("chig", m_chig);
92 status = m_tuple1->addItem("ndofg", m_ndofg);
93 status = m_tuple1->addItem("probg", m_probg);
94 } else {
95 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple1) << endmsg;
96 return StatusCode::FAILURE;
97 }
98 }
99
100 NTuplePtr nt2(ntupleSvc(), "FILE999/chisq"); //chisq of Kalman filter
101 if(nt2) m_tuple2 = nt2;
102 else {
103 m_tuple2 = ntupleSvc()->book ("FILE999/chisq", CLID_ColumnWiseTuple, "chi-square of ");
104 if(m_tuple2) {
105 status = m_tuple2->addItem("chis", m_chis);
106 status = m_tuple2->addItem("probs", m_probs);
107 status = m_tuple2->addItem("chif", m_chif);
108 status = m_tuple2->addItem("probf", m_probf);
109 } else {
110 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple2) << endmsg;
111 return StatusCode::FAILURE;
112 }
113 }
114
115 NTuplePtr nt3(ntupleSvc(), "FILE999/Pull");
116 if(nt3) m_tuple3 = nt3;
117 else {
118 m_tuple3 = ntupleSvc()->book ("FILE999/Pull", CLID_ColumnWiseTuple, "Pull");
119 if(m_tuple3) {
120 status = m_tuple3->addItem("pull_drho", m_pull_drho);
121 status = m_tuple3->addItem("pull_phi", m_pull_phi);
122 status = m_tuple3->addItem("pull_kapha", m_pull_kapha);
123 status = m_tuple3->addItem("pull_dz", m_pull_dz);
124 status = m_tuple3->addItem("pull_lamb", m_pull_lamb);
125 status = m_tuple3->addItem("pull_momentum", m_pull_momentum);
126 } else {
127 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple3) << endmsg;
128 return StatusCode::FAILURE;
129 }
130 }
131
132 NTuplePtr nt4(ntupleSvc(), "FILE999/kalvtx");
133 if(nt4) m_tuple4 = nt4;
134 else {
135 m_tuple4 = ntupleSvc()->book ("FILE999/kalvtx", CLID_ColumnWiseTuple, "kalman vertex");
136 if(m_tuple4) {
137 status = m_tuple4->addItem("kvx", m_kvx);
138 status = m_tuple4->addItem("kvy", m_kvy);
139 status = m_tuple4->addItem("kvz", m_kvz);
140 status = m_tuple4->addItem("chik", m_chik);
141 status = m_tuple4->addItem("ndofk", m_ndofk);
142 status = m_tuple4->addItem("probk", m_probk);
143 } else {
144 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple4) << endmsg;
145 return StatusCode::FAILURE;
146 }
147 }
148 } //end if output
149
150 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
151 return StatusCode::SUCCESS;
152}
153//*******************************************************************************
154
155StatusCode PrimaryVertex::RegisterEvtRecPrimaryVertex(
156 EvtRecPrimaryVertex* aNewEvtRecPrimaryVertex, MsgStream& log) {
157 DataObject *aEvtRecEvent;
158 eventSvc()->findObject("/Event/EvtRec",aEvtRecEvent);
159 if (aEvtRecEvent==NULL) {
160 aEvtRecEvent = new EvtRecEvent();
161 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec",aEvtRecEvent);
162 if(sc!=StatusCode::SUCCESS) {
163 log << MSG::FATAL << "Could not register EvtRecEvent" <<endreq;
164 return StatusCode::FAILURE;
165 }
166 }
167 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
168 DataObject *aEvtRecPrimaryVertex;
169 eventSvc()->findObject("/Event/EvtRec/EvtRecPrimaryVertex",aEvtRecPrimaryVertex);
170 if(aEvtRecPrimaryVertex != NULL) {
171 //IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
172 dataManSvc->clearSubTree("/Event/EvtRec/EvtRecPrimaryVertex");
173 eventSvc()->unregisterObject("/Event/EvtRec/EvtRecPrimaryVertex");
174 }
175 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecPrimaryVertex",
176 aNewEvtRecPrimaryVertex);
177 if (sc != StatusCode::SUCCESS) {
178 log << MSG::FATAL << "Could not register EvtRecPrimaryVertex in TDS!" << endreq;
179 return StatusCode::FAILURE;
180 }
181
182 return StatusCode::SUCCESS;
183}
184
185///////////////////////////////////////////////////
186// Select good charged tracks in MDC ( no PID )
187//////////////////////////////////////////////////
188void PrimaryVertex::SelectGoodChargedTracks(
189 SmartDataPtr<EvtRecEvent>& recEvent,
190 SmartDataPtr<EvtRecTrackCol>& recTrackCol,
191 Vint& icp, Vint& icm, Vint& iGood) {
192 assert(recEvent->totalCharged() <= recTrackCol->size());
193 for (unsigned int i = 0; i < recEvent->totalCharged(); i++) {
194 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
195 if(!(*itTrk)->isMdcTrackValid()) continue;
196 if(!(*itTrk)->isMdcKalTrackValid()) continue;
197 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
198 if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
199 if (fabs(mdcTrk->z()) >= m_vz0Cut) continue;
200 iGood.push_back((*itTrk)->trackId());
201 if (mdcTrk->charge() > 0) {
202 icp.push_back((*itTrk)->trackId());
203 } else if (mdcTrk->charge() < 0) {
204 icm.push_back((*itTrk)->trackId());
205 }
206 }
207}
208
209void InitVertexParameter(VertexParameter& vxpar) {
210 HepPoint3D vx(0., 0., 0.);
211 HepSymMatrix Evx(3, 0);
212 double bx = 1E+6;
213 double by = 1E+6;
214 double bz = 1E+6;
215 Evx[0][0] = bx*bx;
216 Evx[1][1] = by*by;
217 Evx[2][2] = bz*bz;
218 vxpar.setVx(vx);
219 vxpar.setEvx(Evx);
220}
221
222//***************************************************************************
223StatusCode PrimaryVertex::execute()
224{
225 MsgStream log(msgSvc(), name());
226 log << MSG::INFO << "in execute()" << endmsg;
227
228 cout.precision(20);
229 ///////////////////////////////////////////
230 // Read REC information
231 //////////////////////////////////////////
232 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
233 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
234 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
235 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
236 << " " << eventHeader->eventNumber() << endreq;
237 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
238 << recEvent->totalCharged() << " , " << recEvent->totalNeutral() << " , "
239 << recEvent->totalTracks() <<endreq;
240 m_sel_number[0]++;
241 /*
242 if (eventHeader->eventNumber() % 1000 == 0) {
243 cout << "Event number = " << eventHeader->eventNumber() << endl;
244 }*/
245
246 EvtRecPrimaryVertex* aNewEvtRecPrimaryVertex = new EvtRecPrimaryVertex;
247
248 StatusCode sc = RegisterEvtRecPrimaryVertex(aNewEvtRecPrimaryVertex, log);
249 if (sc != StatusCode::SUCCESS) {
250 return sc;
251 }
252
253 // Cut 1 : for anything sample, at least 3 charged tracks
254 if (recEvent->totalCharged() < m_trackNumberCut) {
255 return StatusCode::SUCCESS;
256 }
257 m_sel_number[1]++;
258
259 Vint icp, icm, iGood;
260 SelectGoodChargedTracks(recEvent, recTrackCol, icp, icm, iGood);
261
262 // Cut 2 : for anything sample, at least 3 good charged tracks
263 if (iGood.size() < m_trackNumberCut) {
264 return StatusCode::SUCCESS;
265 }
266 m_sel_number[2]++;
267
268 // Vertex Fit
269 VertexParameter vxpar;
270 InitVertexParameter(vxpar);
271
272 // For Global Vertex Fitting
273 if (m_fitMethod == 0) {
274 VertexFit* vtxfit = VertexFit::instance();
275 vtxfit->init();
276 Vint tlis, trkidlis;
277 for (int i = 0; i < iGood.size(); i++) {
278 int trk_id = iGood[i];
279 EvtRecTrackIterator itTrk = recTrackCol->begin() + trk_id;
280 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
281 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
282 WTrackParameter wtrk(xmass[2], kalTrk->helix(), kalTrk->err());
283 vtxfit->AddTrack(i ,wtrk);
284 tlis.push_back(i);
285 trkidlis.push_back(trk_id);
286 }
287 vtxfit->AddVertex(0, vxpar, tlis);
288 if(!vtxfit->Fit(0)) return StatusCode::SUCCESS; //FIXME
289 if(vtxfit->chisq(0) > m_globalChisqCut) return StatusCode::SUCCESS; //FIXME
290 if (m_output == 1) {
291 m_chig = vtxfit->chisq(0);
292 m_ndofg = 2 * iGood.size() - 3;
293 m_probg = TMath::Prob(m_chig, m_ndofg);
294 HepVector glvtx = vtxfit->Vx(0);
295 m_gvx = glvtx[0];
296 m_gvy = glvtx[1];
297 m_gvz = glvtx[2];
298 m_tuple1->write();
299 }
300/*
301 cout << "======== After global vertex fitting =============" << endl;
302 cout << "Event number = " << eventHeader->eventNumber() << endl;
303 cout << "NTracks: " << iGood.size() << endl;
304 cout << "Chisq: " << vtxfit->chisq(0) << endl;
305 cout << "Ndof: " << 2 * iGood.size() - 3 << endl;
306 cout << "FitMethod: "<< " 0 " << endl;
307 cout << "Vertex position: " << vtxfit->Vx(0)<< endl;
308 cout << "Vertex resolution: " << vtxfit->Evx(0) << endl;
309 cout << "track id list : " << endl;
310 for (int j = 0; j < trkidlis.size(); j++) {
311 cout << trkidlis[j] << " ";
312 }
313 cout << " " << endl; */
314
315 aNewEvtRecPrimaryVertex->setNTracks(iGood.size());
316 aNewEvtRecPrimaryVertex->setTrackIdList(trkidlis);
317 aNewEvtRecPrimaryVertex->setChi2(vtxfit->chisq(0));
318 aNewEvtRecPrimaryVertex->setNdof(2 * iGood.size() - 3);
319 aNewEvtRecPrimaryVertex->setFitMethod(0);
320 aNewEvtRecPrimaryVertex->setVertex(vtxfit->Vx(0));
321 aNewEvtRecPrimaryVertex->setErrorVertex(vtxfit->Evx(0));
322 aNewEvtRecPrimaryVertex->setIsValid(true);
323
324 } else if (m_fitMethod == 1) {
325 // For Kalman Vertex Fitting
326 KalmanVertexFit *kalvtx = KalmanVertexFit::instance();
327 kalvtx->init();
328 kalvtx->initVertex(vxpar);
329 kalvtx->setChisqCut(m_chisqCut);
330 kalvtx->setTrackIteration(m_trackIteration);
331 kalvtx->setVertexIteration(m_vertexIteration);
332 kalvtx->setTrackIterationCut(m_chi2CutforTrkIter);
333 for(int i = 0; i < iGood.size(); i++) {
334 int trk_id = iGood[i];
335 EvtRecTrackIterator itTrk = recTrackCol->begin()+trk_id;
336 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
337 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
338 HTrackParameter htrk(kalTrk->helix(), kalTrk->err(), trk_id, 2);
339 kalvtx->addTrack(htrk);
340 }
341 kalvtx->filter();
342
343 //int freedomCut = -3 + 2*m_trackNumberCut;
344 if(kalvtx->ndof() >= m_freedomCut) { //Cut : add 2 when add every track
345 //for(int i = 0; i < kalvtx->numTrack(); i++) {
346 for(int i = 0; i < iGood.size(); i++) {
347 if (m_output == 1) {
348 m_chif = kalvtx->chiF(i);
349 m_chis = kalvtx->chiS(i);
350 m_probs = TMath::Prob(m_chis, 2);
351 m_probf = TMath::Prob(m_chif, 2);
352 m_tuple2->write();
353 }
354 if(kalvtx->chiS(i) > m_chi2CutforSmooth) kalvtx->remove(i);
355 }
356 }
357 if(kalvtx->ndof() >= m_freedomCut) { //FIXME to Rhopi is 0 , others are 1
358 for(int i = 0; i < iGood.size(); i++) {
359 kalvtx->smooth(i);
360
361 HepVector Pull(5, 0);
362 Pull = kalvtx->pull(i);
363 if (m_output == 1) {
364 m_pull_drho = Pull[0];
365 m_pull_phi = Pull[1];
366 m_pull_kapha = Pull[2];
367 m_pull_dz = Pull[3];
368 m_pull_lamb = Pull[4];
369 m_pull_momentum = kalvtx->pullmomentum(i);
370 m_tuple3->write();
371 }
372 }
373 if (m_output == 1) {
374 m_chik = kalvtx->chisq();
375 m_ndofk = kalvtx->ndof();
376 m_probk = TMath::Prob(m_chik, m_ndofk);
377 HepVector xp(3, 0);
378 xp = kalvtx->x();
379 m_kvx = xp[0];
380 m_kvy = xp[1];
381 m_kvz = xp[2];
382 m_tuple4->write();
383 }
384
385 m_sel_number[3]++;
386 /*
387 cout << "======== After Kalman vertex fitting =============" << endl;
388 cout << "NTracks: " << kalvtx->numTrack() << endl;
389 cout << "Chisq: " << kalvtx->chisq() << endl;
390 cout << "Ndof: " << kalvtx->ndof() << endl;
391 cout << "FitMethod: "<< " 1 " << endl;
392 cout << "Vertex position: " << kalvtx->x() << endl;
393 cout << "Vertex resolution: " << kalvtx->Ex() << endl;
394 for (int j = 0; j < (kalvtx->trackIDList()).size(); j++) {
395 cout << kalvtx->trackIDList()[j] << " ";
396 }
397 cout << " " << endl;*/
398
399 aNewEvtRecPrimaryVertex->setNTracks(kalvtx->numTrack());
400 aNewEvtRecPrimaryVertex->setTrackIdList(kalvtx->trackIDList());
401 aNewEvtRecPrimaryVertex->setChi2(kalvtx->chisq());
402 aNewEvtRecPrimaryVertex->setNdof(kalvtx->ndof());
403 aNewEvtRecPrimaryVertex->setFitMethod(1);
404 aNewEvtRecPrimaryVertex->setVertex(kalvtx->x());
405 aNewEvtRecPrimaryVertex->setErrorVertex(kalvtx->Ex());
406 aNewEvtRecPrimaryVertex->setIsValid(true);
407
408 }
409 }
410 return StatusCode::SUCCESS;
411}
412
413//***************************************************************************
414
415StatusCode PrimaryVertex::finalize()
416{
417 MsgStream log(msgSvc(), name());
418 log << MSG::INFO << "in finalize()" << endmsg;
419
420 log << MSG::ALWAYS << "==================================" << endreq;
421 log << MSG::ALWAYS << "survived event :";
422 for(int i = 0; i < 10; i++){
423 log << MSG::ALWAYS << m_sel_number[i] << " ";
424 }
425 log << MSG::ALWAYS << endreq;
426 log << MSG::ALWAYS << "==================================" << endreq;
427 return StatusCode::SUCCESS;
428}
429
xmass
const double xmass[5]
Definition: Gam4pikp.cxx:50
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: PrimaryVertex.cxx:4
Vp4
std::vector< HepLorentzVector > Vp4
Definition: PrimaryVertex.cxx:42
xmass
const double xmass[5]
Definition: PrimaryVertex.cxx:36
mpi0
const double mpi0
Definition: PrimaryVertex.cxx:38
momega
const double momega
Definition: PrimaryVertex.cxx:39
ecms
const double ecms
Definition: PrimaryVertex.cxx:37
Vint
std::vector< int > Vint
Definition: PrimaryVertex.cxx:41
InitVertexParameter
void InitVertexParameter(VertexParameter &vxpar)
Definition: PrimaryVertex.cxx:209
KalmanVertexFit::pullmomentum
double pullmomentum(const int k)
Definition: KalmanVertexFit.cxx:480
KalmanVertexFit::Ex
HepSymMatrix Ex() const
Definition: KalmanVertexFit.cxx:76
KalmanVertexFit::addTrack
void addTrack(const HTrackParameter)
Definition: KalmanVertexFit.cxx:81
KalmanVertexFit::smooth
void smooth(const int k)
Definition: KalmanVertexFit.cxx:281
KalmanVertexFit::filter
int filter(const int k)
Definition: KalmanVertexFit.cxx:135
KalmanVertexFit::trackIDList
std::vector< int > trackIDList() const
Definition: KalmanVertexFit.cxx:126
KalmanVertexFit::chiS
double chiS(const int k) const
Definition: KalmanVertexFit.cxx:263
KalmanVertexFit::initVertex
void initVertex(const VertexParameter vtx)
Definition: KalmanVertexFit.cxx:59
KalmanVertexFit::numTrack
int numTrack() const
Definition: KalmanVertexFit.cxx:118
KalmanVertexFit::pull
HepVector pull(const int k)
Definition: KalmanVertexFit.cxx:468
KalmanVertexFit::instance
static KalmanVertexFit * instance()
Definition: KalmanVertexFit.cxx:19
KalmanVertexFit::remove
void remove(const int k)
Definition: KalmanVertexFit.cxx:355
PrimaryVertex::finalize
StatusCode finalize()
Definition: PrimaryVertex.cxx:415
PrimaryVertex::initialize
StatusCode initialize()
Definition: PrimaryVertex.cxx:71
PrimaryVertex::PrimaryVertex
PrimaryVertex(const std::string &name, ISvcLocator *pSvcLocator)
Definition: PrimaryVertex.cxx:45
PrimaryVertex::execute
StatusCode execute()
Definition: PrimaryVertex.cxx:223
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::Fit
bool Fit()
Definition: VertexFit.cxx:301