00001
00002
00003
00004 #include "CLHEP/Vector/LorentzVector.h"
00005 #include "CLHEP/Geometry/Point3D.h"
00006 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00007 typedef HepGeom::Point3D<double> HepPoint3D;
00008 #endif
00009 #include "GaudiKernel/MsgStream.h"
00010 #include "GaudiKernel/SmartDataPtr.h"
00011 #include "GaudiKernel/IDataManagerSvc.h"
00012 #include "EventModel/EventModel.h"
00013 #include "EventModel/Event.h"
00014 #include "EventModel/EventHeader.h"
00015 #include "EvtRecEvent/EvtRecEvent.h"
00016 #include "EvtRecEvent/EvtRecTrack.h"
00017 #include "VertexFit/VertexFit.h"
00018 #include "VertexFit/HTrackParameter.h"
00019 #include "VertexFit/KalmanVertexFit.h"
00020 #include "ParticleID/ParticleID.h"
00021 #include "BeamParamsAlg/BeamParams.h"
00022 #include "GaudiKernel/ITHistSvc.h"
00023 #include <assert.h>
00024 #include "TMath.h"
00025 #include "TH2D.h"
00026 #include "TH1D.h"
00027 #include "TF1.h"
00028 #include "TCanvas.h"
00029 #include "TPad.h"
00030 #include <map>
00031 #include <iostream>
00032 #include <fstream>
00033
00034 using CLHEP::HepLorentzVector;
00035 using namespace std;
00036
00037 const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
00038 const double ecms = 3.097;
00039
00040 typedef std::vector<int> Vint;
00041 typedef std::vector<HepLorentzVector> Vp4;
00042
00043
00044 BeamParams::BeamParams(const std::string& name, ISvcLocator* pSvcLocator) :
00045 Algorithm (name, pSvcLocator)
00046 {
00047
00048
00049
00050 declareProperty("ParVer", m_parVer = 1);
00051 declareProperty("TrackNumberCut", m_trackNumberCut = 1);
00052 declareProperty("Vz0Cut", m_vz0Cut = 50.0);
00053 declareProperty("CosThetaCut", m_cosThetaCut=0.93);
00054
00055
00056 declareProperty("OptionPID", m_pid = 0);
00057 declareProperty("PidUseDedx", m_useDedx = true);
00058 declareProperty("PidUseTof1", m_useTof1 = true);
00059 declareProperty("PidUseTof2", m_useTof2 = true);
00060 declareProperty("PidUseTofE", m_useTofE = false);
00061 declareProperty("PidUseTofQ", m_useTofQ = false);
00062 declareProperty("PidUseEmc", m_useEmc = false);
00063 declareProperty("PidUseMuc", m_useMuc = false);
00064
00065
00066 declareProperty("OptionVertexFind", m_vertexFind = 0);
00067 declareProperty("MinDistance", m_minDistance = 100.0);
00068 declareProperty("MinPointX", m_minPointX = 0.1);
00069 declareProperty("MinPointY", m_minPointY = 0.1);
00070 declareProperty("MeanPointX", m_meanPointX = 0.1);
00071 declareProperty("MeanPointY", m_meanPointY = -0.25);
00072
00073
00074 declareProperty("ChisqCut", m_chisqCut = 200.0);
00075 declareProperty("TrackIteration", m_trackIteration = 5);
00076 declareProperty("VertexIteration", m_vertexIteration = 5);
00077 declareProperty("Chi2CutforTrkIter", m_chi2CutforTrkIter = 0.1);
00078 declareProperty("Chi2CutforSmooth", m_chi2CutforSmooth = 200.0);
00079
00080
00081 declareProperty("HadronFile", m_hadronFile = 0);
00082 declareProperty("DstFileName", m_fileNameDst = std::string("dst.txt"));
00083 declareProperty("HadronFileName", m_fileNameHadron = std::string("hadron.txt"));
00084 declareProperty("FigsName", m_figsName = std::string("figs.ps"));
00085 }
00086
00087
00088 StatusCode BeamParams::initialize()
00089 {
00090 MsgStream log(msgSvc(), name());
00091 log << MSG::INFO << "in initialize()" << endmsg;
00092
00093 StatusCode sc=service("THistSvc", m_thsvc);
00094 if(sc.isFailure()) {
00095 log << MSG::FATAL << "Couldn't get THistSvc" << endreq;
00096 exit(1);
00097 }
00098
00099 for(int i = 0; i < 15; i++){
00100 m_sel_number[i] = 0;
00101 }
00102
00104
00106 m_vertex_x = new TH1D("x_of_vertex", "x of vertex", 200, -5, 5);
00107 m_vertex_y = new TH1D("y_of_vertex", "y of vertex", 200, -5, 5);
00108 m_vertex_z = new TH1D("z_of_vertex", "z of vertex", 200, -10, 10);
00109 m_vertex_x_kal = new TH1D("x_of_vertex_in_kal", "x of vertex in kal", 200, -5, 5);
00110 m_vertex_y_kal = new TH1D("y_of_vertex_in_kal", "y of vertex in kal", 200, -5, 5);
00111 m_vertex_z_kal = new TH1D("z_of_vertex_in_kal", "z of vertex in kal", 200, -10, 10);
00112 if(m_thsvc->regHist("/DQAHist/zhsVER/x_of_vertex",m_vertex_x).isFailure()){
00113 log << MSG::FATAL << "Couldn't register x of vertex " << endreq;
00114 exit(1);
00115 }
00116 if(m_thsvc->regHist("/DQAHist/zhsVER/y_of_vertex",m_vertex_y).isFailure()){
00117 log << MSG::FATAL << "Couldn't register y of vertex" << endreq;
00118 exit(1);
00119 }
00120 if(m_thsvc->regHist("/DQAHist/zhsVER/z_of_vertex",m_vertex_z).isFailure()){
00121 log << MSG::FATAL << "Couldn't register z of vertex" << endreq;
00122 exit(1);
00123 }
00124 if(m_thsvc->regHist("/DQAHist/zhsVER/x_of_vertex_in_kal",m_vertex_x_kal).isFailure()){
00125 log << MSG::FATAL << "Couldn't register x of vertex in kal" << endreq;
00126 exit(1);
00127 }
00128 if(m_thsvc->regHist("/DQAHist/zhsVER/y_of_vertex_in_kal",m_vertex_y_kal).isFailure()){
00129 log << MSG::FATAL << "Couldn't register y of vertex in kal" << endreq;
00130 exit(1);
00131 }
00132 if(m_thsvc->regHist("/DQAHist/zhsVER/z_of_vertex_in_kal",m_vertex_z_kal).isFailure()){
00133 log << MSG::FATAL << "Couldn't register z of vertex in kal" << endreq;
00134 exit(1);
00135 }
00136
00137
00138 StatusCode status;
00139
00140 NTuplePtr nt1(ntupleSvc(), "FILE1/minid");
00141 if(nt1) m_tuple1 = nt1;
00142 else {
00143 m_tuple1 = ntupleSvc()->book ("FILE1/minid", CLID_ColumnWiseTuple, "minimal distance");
00144 if(m_tuple1) {
00145 status = m_tuple1->addItem("xc", m_xc);
00146 status = m_tuple1->addItem("yc", m_yc);
00147 status = m_tuple1->addItem("zc", m_zc);
00148 status = m_tuple1->addItem("mind", m_mind);
00149 } else {
00150 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00151 return StatusCode::FAILURE;
00152 }
00153 }
00154
00155 NTuplePtr nt2(ntupleSvc(), "FILE1/chisq");
00156 if(nt2) m_tuple2 = nt2;
00157 else {
00158 m_tuple2 = ntupleSvc()->book ("FILE1/chisq", CLID_ColumnWiseTuple, "chi-square of ");
00159 if(m_tuple2) {
00160 status = m_tuple2->addItem("chis", m_chis);
00161 status = m_tuple2->addItem("probs", m_probs);
00162 status = m_tuple2->addItem("chif", m_chif);
00163 status = m_tuple2->addItem("probf", m_probf);
00164 } else {
00165 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00166 return StatusCode::FAILURE;
00167 }
00168 }
00169
00170 NTuplePtr nt3(ntupleSvc(), "FILE1/kalvtx");
00171 if(nt3) m_tuple3 = nt3;
00172 else {
00173 m_tuple3 = ntupleSvc()->book ("FILE1/kalvtx", CLID_ColumnWiseTuple, "kalman vertex");
00174 if(m_tuple3) {
00175 status = m_tuple3->addItem("kvx", m_kvx);
00176 status = m_tuple3->addItem("kvy", m_kvy);
00177 status = m_tuple3->addItem("kvz", m_kvz);
00178 status = m_tuple3->addItem("chik", m_chik);
00179 status = m_tuple3->addItem("ndofk", m_ndofk);
00180 status = m_tuple3->addItem("probk", m_probk);
00181 } else {
00182 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00183 return StatusCode::FAILURE;
00184 }
00185 }
00186
00187 NTuplePtr nt4(ntupleSvc(), "FILE1/glbvtx");
00188 if(nt4) m_tuple4 = nt4;
00189 else {
00190 m_tuple4 = ntupleSvc()->book ("FILE1/glbvtx", CLID_ColumnWiseTuple, "global vertex");
00191 if(m_tuple4) {
00192 status = m_tuple4->addItem("gvx", m_gvx);
00193 status = m_tuple4->addItem("gvy", m_gvy);
00194 status = m_tuple4->addItem("gvz", m_gvz);
00195 status = m_tuple4->addItem("chig", m_chig);
00196 status = m_tuple4->addItem("ndofg", m_ndofg);
00197 status = m_tuple4->addItem("probg", m_probg);
00198 } else {
00199 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00200 return StatusCode::FAILURE;
00201 }
00202 }
00203
00204 NTuplePtr nt5(ntupleSvc(), "FILE1/Pull");
00205 if(nt5) m_tuple5 = nt5;
00206 else {
00207 m_tuple5 = ntupleSvc()->book ("FILE1/Pull", CLID_ColumnWiseTuple, "Pull");
00208 if(m_tuple5) {
00209 status = m_tuple5->addItem("pull_drho", m_pull_drho);
00210 status = m_tuple5->addItem("pull_phi", m_pull_phi);
00211 status = m_tuple5->addItem("pull_kapha", m_pull_kapha);
00212 status = m_tuple5->addItem("pull_dz", m_pull_dz);
00213 status = m_tuple5->addItem("pull_lamb", m_pull_lamb);
00214 status = m_tuple5->addItem("pull_momentum", m_pull_momentum);
00215 } else {
00216 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple5) << endmsg;
00217 return StatusCode::FAILURE;
00218 }
00219 }
00220
00221 NTuplePtr nt6(ntupleSvc(), "FILE1/MdcTrack");
00222 if(nt6) m_tuple6 = nt6;
00223 else {
00224 m_tuple6 = ntupleSvc()->book ("FILE1/MdcTrack", CLID_ColumnWiseTuple, "MdcTrack");
00225 if(m_tuple6) {
00226 status = m_tuple6->addItem("MdcTrkX", m_mdcTrk_x);
00227 status = m_tuple6->addItem("MdcTrkY", m_mdcTrk_y);
00228 status = m_tuple6->addItem("MdcTrkZ", m_mdcTrk_z);
00229 status = m_tuple6->addItem("MdcTrkR", m_mdcTrk_r);
00230 status = m_tuple6->addItem("MdcTrkDrho", m_mdcTrk_dr);
00231 status = m_tuple6->addItem("Rxy", m_rxy);
00232 status = m_tuple6->addItem("MdcKalTrkZ", m_mdcKalTrk_z);
00233 } else {
00234 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple6) << endmsg;
00235 return StatusCode::FAILURE;
00236 }
00237 }
00238
00239 NTuplePtr nt7(ntupleSvc(), "FILE1/PullG");
00240 if(nt7) m_tuple7 = nt7;
00241 else {
00242 m_tuple7 = ntupleSvc()->book ("FILE1/PullG", CLID_ColumnWiseTuple, "Pull");
00243 if(m_tuple7) {
00244 status = m_tuple7->addItem("gpull_drho", m_gpull_drho);
00245 status = m_tuple7->addItem("gpull_phi", m_gpull_phi);
00246 status = m_tuple7->addItem("gpull_kapha", m_gpull_kapha);
00247 status = m_tuple7->addItem("gpull_dz", m_gpull_dz);
00248 status = m_tuple7->addItem("gpull_lamb", m_gpull_lamb);
00249 } else {
00250 log << MSG::FATAL << "Cannot book N-tuple:" << long(m_tuple7) << endmsg;
00251 return StatusCode::FAILURE;
00252 }
00253 }
00254
00255 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00256 return StatusCode::SUCCESS;
00257 }
00258
00259 #define DEB cout << __FILE__ << ":" << __LINE__ << endl;
00260
00261
00262 StatusCode BeamParams::execute()
00263 {
00264 MsgStream log(msgSvc(), name());
00265 log << MSG::INFO << "in execute()" << endmsg;
00266
00268
00270 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00271 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00272 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00273 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00274 << " " << eventHeader->eventNumber() << endreq;
00275 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00276 << recEvent->totalCharged() << " , "
00277 << recEvent->totalNeutral() << " , "
00278 << recEvent->totalTracks() <<endreq;
00279
00280 m_runNo = eventHeader->runNumber();
00281
00282 if (eventHeader->eventNumber() % 5000 == 0) {
00283 cout << "Event number = " << eventHeader->eventNumber()<<" run number = "<< m_runNo << endl;
00284 }
00285 m_sel_number[0]++;
00286
00287
00288 if (recEvent->totalCharged() < m_trackNumberCut) return StatusCode::SUCCESS;
00289 m_sel_number[1]++;
00290
00292
00294 Vint icp, icm, iGood, jGood;
00295 icp.clear();
00296 icm.clear();
00297 iGood.clear();
00298 jGood.clear();
00299
00300 map<int, int> ParticleType;
00301 ParticleID* pid = ParticleID::instance();
00302
00303 for(unsigned int i = 0; i < recEvent->totalCharged(); i++) {
00304 jGood.push_back(0);
00305 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
00306 if(!(*itTrk)->isMdcTrackValid()) continue;
00307 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00308 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00309 if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
00310 if (fabs(mdcTrk->z()) >= m_vz0Cut) continue;
00311
00312 iGood.push_back(i);
00313
00314 if (m_pid == 0) {
00315 if (mdcTrk->charge() > 0) {
00316
00317 icp.push_back(i);
00318 } else if (mdcTrk->charge() < 0) {
00319
00320 icm.push_back(i);
00321 }
00322 } else if (m_pid == 1) {
00323 ParticleType[i] = 2;
00324
00325 pid->init();
00326 pid->setChiMinCut(4);
00327 pid->setRecTrack(*itTrk);
00328 pid->setMethod(pid->methodProbability());
00329
00330
00331 if(m_useDedx)pid->usePidSys(pid->useDedx());
00332 if(m_useTof1)pid->usePidSys(pid->useTof1());
00333 if(m_useTof2)pid->usePidSys(pid->useTof2());
00334 if(m_useTofE)pid->usePidSys(pid->useTofE());
00335 if(m_useTofQ)pid->usePidSys(pid->useTofQ());
00336 if(m_useEmc)pid->usePidSys(pid->useEmc());
00337 if(m_useMuc)pid->usePidSys(pid->useMuc());
00338 pid->identify(pid->onlyElectron());
00339 pid->identify(pid->onlyMuon());
00340 pid->identify(pid->onlyPion());
00341 pid->identify(pid->onlyKaon());
00342 pid->identify(pid->onlyProton());
00343 pid->calculate();
00344 if(!pid->IsPidInfoValid()) continue;
00345 double prob_value[5];
00346 prob_value[0] = pid->probElectron();
00347 prob_value[1] = pid->probMuon();
00348 prob_value[2] = pid->probPion();
00349 prob_value[3] = pid->probKaon();
00350 prob_value[4] = pid->probProton();
00351
00352 int max = 0;
00353 for (int i = 1; i < 5; i++) {
00354 if (prob_value[i] > prob_value[max]) {
00355 max = i;
00356 }
00357 }
00358
00359
00360
00361 if(max == 0) ParticleType[i] = 0;
00362 if(max == 1) ParticleType[i] = 1;
00363 if(max == 2) ParticleType[i] = 2;
00364 if(max == 3) ParticleType[i] = 3;
00365 if(max == 4) ParticleType[i] =4;
00366 }
00367
00368
00369 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
00370 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00371 m_mdcTrk_x = mdcTrk->x();
00372 m_mdcTrk_y = mdcTrk->y();
00373 m_mdcTrk_z = mdcTrk->helix(3);
00374 m_mdcTrk_r = mdcTrk->r();
00375 m_mdcTrk_dr = mdcTrk->helix(0);
00376 m_mdcKalTrk_z = kalTrk->helix()[3];
00377 m_tuple6->write();
00378 }
00379
00380
00381 if (iGood.size() < m_trackNumberCut) return StatusCode::SUCCESS;
00382 m_sel_number[2]++;
00383
00385
00387 if (m_vertexFind == 1) {
00388 int nGood = iGood.size();
00389 for(int i1 = 0; i1 < nGood; i1++) {
00390 int id_trk1 = iGood[i1];
00391 EvtRecTrackIterator trk1 = (recTrackCol->begin()+id_trk1);
00392 RecMdcKalTrack *kalTrk1 = (*trk1)->mdcKalTrack();
00393
00394 HTrackParameter htrk1, htrk2;
00395 if (m_pid == 0) {
00396 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00397 htrk1 = HTrackParameter(kalTrk1->helix(), kalTrk1->err(), id_trk1, 2);
00398 } else if (m_pid == 1) {
00399 RecMdcKalTrack::setPidType((RecMdcKalTrack::PidType)ParticleType[id_trk1]);
00400 htrk1 = HTrackParameter(kalTrk1->helix(),kalTrk1->err(), id_trk1, ParticleType[id_trk1]);
00401 }
00402 for(int i2 = i1+1; i2 < nGood; i2++) {
00403 int id_trk2 = iGood[i2];
00404 EvtRecTrackIterator trk2 = (recTrackCol->begin()+id_trk2);
00405 RecMdcKalTrack *kalTrk2 = (*trk2)->mdcKalTrack();
00406
00407 if (m_pid == 0) {
00408 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00409 htrk2 = HTrackParameter(kalTrk2->helix(), kalTrk2->err(), id_trk2, 2);
00410 } else if (m_pid == 1) {
00411 RecMdcKalTrack::setPidType((RecMdcKalTrack::PidType)ParticleType[id_trk2]);
00412 htrk2 = HTrackParameter(kalTrk2->helix(),kalTrk2->err(),id_trk2,ParticleType[id_trk2]);
00413 }
00414 HepPoint3D cross(0., 0., 0.);
00415 double dist = htrk1.minDistanceTwoHelix(htrk2, cross);
00416 m_mind = dist;
00417 m_xc = cross.x();
00418 m_yc = cross.y();
00419 m_zc = cross.z();
00420 m_tuple1->write();
00421
00422 if(sqrt(cross.x()*cross.x()+cross.y()*cross.y())>2)continue;
00423
00424 if(fabs(dist) > m_minDistance) continue;
00425
00426
00427
00428
00429
00430 double cross_cut=0.;
00431 if(cross_cut < 3.5 * 3.5) {
00432 jGood[id_trk1] = 1;
00433 jGood[id_trk2] = 1;
00434 }
00435 }
00436 }
00437
00438 iGood.clear();
00439 for(int i =0; i < jGood.size(); i++) {
00440 if(jGood[i] == 1) iGood.push_back(i);
00441 }
00442
00443 }
00444
00445 if (iGood.size() >= 2) m_sel_number[3]++;
00446 if (iGood.size() >= 3) m_sel_number[4]++;
00447
00449
00451 HepPoint3D vx(0., 0., 0.);
00452 HepSymMatrix Evx(3, 0);
00453 double bx = 1E+6;
00454 double by = 1E+6;
00455 double bz = 1E+6;
00456 Evx[0][0] = bx*bx;
00457 Evx[1][1] = by*by;
00458 Evx[2][2] = bz*bz;
00459 VertexParameter vxpar;
00460 vxpar.setVx(vx);
00461 vxpar.setEvx(Evx);
00462
00463 VertexFit* vtxfit = VertexFit::instance();
00464
00465 vtxfit->init();
00466 Vint tlis;
00467 tlis.clear();
00468 for(int i = 0; i < iGood.size(); i++) {
00469 int trk_id = iGood[i];
00470 EvtRecTrackIterator itTrk = recTrackCol->begin()+trk_id;
00471 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
00472 if (m_pid == 0) {
00473 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00474 WTrackParameter wtrk(xmass[2], kalTrk->helix(), kalTrk->err());
00475 vtxfit->AddTrack(i ,wtrk);
00476 } else if (m_pid == 1) {
00477 RecMdcKalTrack::setPidType((RecMdcKalTrack::PidType)ParticleType[trk_id]);
00478 WTrackParameter wtrk(xmass[ParticleType[trk_id]], kalTrk->helix(), kalTrk->err());
00479 vtxfit->AddTrack(i ,wtrk);
00480 }
00481 tlis.push_back(i);
00482 }
00483
00484 vtxfit->AddVertex(0, vxpar, tlis);
00485 if(!vtxfit->Fit(0)) return StatusCode::SUCCESS;
00486 m_chig = vtxfit->chisq(0);
00487 m_ndofg = 2 * iGood.size() - 3;
00488 m_probg = TMath::Prob(m_chig, m_ndofg);
00489 HepVector glvtx = vtxfit->Vx(0);
00490 m_gvx = glvtx[0];
00491 m_gvy = glvtx[1];
00492 m_gvz = glvtx[2];
00493 m_tuple4->write();
00494
00495 if(!(m_vertex_x->Fill(glvtx[0]))){
00496 log << MSG::FATAL << "Couldn't Fill x of vertex" << endreq;
00497 exit(1);
00498 }
00499 if(!(m_vertex_y->Fill(glvtx[1]))){
00500 log << MSG::FATAL << "Couldn't Fill y of vertex" << endreq;
00501 exit(1);
00502 }
00503 if(!(m_vertex_z->Fill(glvtx[2]))){
00504 log << MSG::FATAL << "Couldn't Fill z of vertex" << endreq;
00505 exit(1);
00506 }
00507
00508 for (int j = 0; j < iGood.size(); j++) {
00509 HepVector pull(5, 0);
00510 bool success = vtxfit->pull(0, j, pull);
00511 if (success) {
00512 m_gpull_drho = pull[0];
00513 m_gpull_phi = pull[1];
00514 m_gpull_kapha = pull[2];
00515 m_gpull_dz = pull[3];
00516 m_gpull_lamb = pull[4];
00517 m_tuple7->write();
00518 }
00519 }
00520
00521
00522 KalmanVertexFit *kalvtx = KalmanVertexFit::instance();
00523 kalvtx->init();
00524 kalvtx->initVertex(vxpar);
00525 kalvtx->setChisqCut(m_chisqCut);
00526 kalvtx->setTrackIteration(m_trackIteration);
00527 kalvtx->setVertexIteration(m_vertexIteration);
00528 kalvtx->setTrackIterationCut(m_chi2CutforTrkIter);
00529 for(int i = 0; i < iGood.size(); i++) {
00530 int trk_id = iGood[i];
00531 EvtRecTrackIterator itTrk = recTrackCol->begin()+trk_id;
00532 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
00533 if (m_pid == 0) {
00534 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00535 HTrackParameter htrk(kalTrk->helix(), kalTrk->err(), trk_id, 2);
00536 kalvtx->addTrack(htrk);
00537 } else if (m_pid == 1) {
00538 RecMdcKalTrack::setPidType((RecMdcKalTrack::PidType)ParticleType[trk_id]);
00539 HTrackParameter htrk(kalTrk->helix(), kalTrk->err(), trk_id, ParticleType[trk_id]);
00540 kalvtx->addTrack(htrk);
00541 }
00542 }
00543 kalvtx->filter();
00544
00545 int freedomCut = 1;
00546 if(kalvtx->ndof() >= freedomCut) {
00547
00548 for(int i = 0; i < iGood.size(); i++) {
00549 m_chif = kalvtx->chiF(i);
00550 m_chis = kalvtx->chiS(i);
00551 m_probs = TMath::Prob(m_chis, 2);
00552 m_probf = TMath::Prob(m_chif, 2);
00553 m_tuple2->write();
00554
00555 if(kalvtx->chiS(i) > m_chi2CutforSmooth) kalvtx->remove(i);
00556 }
00557 }
00558 if(kalvtx->ndof() >= freedomCut) {
00559 for(int i = 0; i < iGood.size(); i++) {
00560 kalvtx->smooth(i);
00561
00562 HepVector Pull(5, 0);
00563 Pull = kalvtx->pull(i);
00564 m_pull_drho = Pull[0];
00565 m_pull_phi = Pull[1];
00566 m_pull_kapha = Pull[2];
00567 m_pull_dz = Pull[3];
00568 m_pull_lamb = Pull[4];
00569 m_pull_momentum = kalvtx->pullmomentum(i);
00570 m_tuple5->write();
00571 }
00572
00573 m_chik = kalvtx->chisq();
00574 m_ndofk = kalvtx->ndof();
00575 m_probk = TMath::Prob(m_chik, m_ndofk);
00576 HepVector xp(3, 0);
00577 xp = kalvtx->x();
00578 m_kvx = xp[0];
00579 m_kvy = xp[1];
00580 m_kvz = xp[2];
00581 m_tuple3->write();
00582
00583 if(!(m_vertex_x_kal->Fill(xp[0]))){
00584 log << MSG::FATAL << "Couldn't Fill kal x of vertex" << endreq;
00585 exit(1);
00586 }
00587 if(!(m_vertex_y_kal->Fill(xp[1]))){
00588 log << MSG::FATAL << "Couldn't Fill kal y of vertex" << endreq;
00589 exit(1);
00590 }
00591 if(!(m_vertex_z_kal->Fill(xp[2]))){
00592 log << MSG::FATAL << "Couldn't Fill kal z of vertex" << endreq;
00593 exit(1);
00594 }
00595
00596 m_sel_number[3]++;
00597 }
00598 return StatusCode::SUCCESS;
00599 }
00600
00601
00602
00603 StatusCode BeamParams::finalize()
00604 {
00605 MsgStream log(msgSvc(), name());
00606 log << MSG::INFO << "in finalize()" << endmsg;
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00661
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00689 log << MSG::ALWAYS << "==================================" << endreq;
00690 log << MSG::ALWAYS << "survived event :";
00691 for(int i = 0; i < 10; i++){
00692 log << MSG::ALWAYS << m_sel_number[i] << " ";
00693 }
00694 log << MSG::ALWAYS << endreq;
00695 log << MSG::ALWAYS << "==================================" << endreq;
00696 return StatusCode::SUCCESS;
00697 }
00698