00001 #include <vector>
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/AlgFactory.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "GaudiKernel/SmartDataPtr.h"
00007 #include "GaudiKernel/IDataProviderSvc.h"
00008 #include "GaudiKernel/PropertyMgr.h"
00009 #include "GaudiKernel/Bootstrap.h"
00010 #include "GaudiKernel/ISvcLocator.h"
00011
00012 #include "GaudiKernel/INTupleSvc.h"
00013 #include "GaudiKernel/NTuple.h"
00014 #include "GaudiKernel/ITHistSvc.h"
00015
00016
00017 #include "EventModel/EventModel.h"
00018 #include "EventModel/Event.h"
00019
00020 #include "EvtRecEvent/EvtRecEvent.h"
00021 #include "EvtRecEvent/EvtRecTrack.h"
00022 #include "DstEvent/TofHitStatus.h"
00023 #include "EventModel/EventHeader.h"
00024
00025 #include "TH1F.h"
00026
00027 #include "CLHEP/Vector/ThreeVector.h"
00028 #include "CLHEP/Vector/LorentzVector.h"
00029 #include "CLHEP/Vector/TwoVector.h"
00030
00031 using CLHEP::Hep3Vector;
00032 using CLHEP::Hep2Vector;
00033 using CLHEP::HepLorentzVector;
00034 #include "CLHEP/Geometry/Point3D.h"
00035
00036 #include "VertexFit/IVertexDbSvc.h"
00037 #include "VertexFit/KinematicFit.h"
00038 #include "VertexFit/VertexFit.h"
00039 #include "VertexFit/SecondVertexFit.h"
00040 #include "ParticleID/ParticleID.h"
00041 #include "TH1F.h"
00042
00043 #include "JsiLL/JsiLL.h"
00044
00045 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00046 typedef HepGeom::Point3D<double> HepPoint3D;
00047 #endif
00048 using CLHEP::HepLorentzVector;
00049
00050 const double mpi = 0.13957;
00051 const double mk = 0.493677;
00052 const double mproton = 0.938272;
00053 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00054 const double velc = 299.792458;
00055 typedef std::vector<int> Vint;
00056 typedef std::vector<HepLorentzVector> Vp4;
00057
00059
00060 JsiLL::JsiLL(const std::string& name, ISvcLocator* pSvcLocator) :
00061 Algorithm(name, pSvcLocator) {
00062
00063
00064 declareProperty("Vr0cut", m_vr0cut=5.0);
00065 declareProperty("Vz0cut", m_vz0cut=20.0);
00066 declareProperty("Vr1cut", m_vr1cut=1.0);
00067 declareProperty("Vz1cut", m_vz1cut=5.0);
00068 declareProperty("Vctcut", m_cthcut=0.93);
00069 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00070 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
00071 }
00072
00073
00074
00075 StatusCode JsiLL::initialize(){
00076 MsgStream log(msgSvc(), name());
00077
00078 log << MSG::INFO << "in initialize()" << endmsg;
00079 StatusCode status;
00080
00081
00082
00083
00084 NTuplePtr nt(ntupleSvc(), "DQAFILE/JsiLL");
00085 if ( nt ) m_tuple = nt;
00086 else {
00087 m_tuple = ntupleSvc()->book("DQAFILE/JsiLL", CLID_ColumnWiseTuple, "JsiLL ntuple");
00088 if( m_tuple ) {
00089 status = m_tuple->addItem("runNo", m_runNo);
00090 status = m_tuple->addItem("event", m_event);
00091 status = m_tuple->addItem("chisq", m_chisq);
00092 status = m_tuple->addItem("mLambda", m_mLambda);
00093 status = m_tuple->addItem("mLambdabar", m_mLambdabar);
00094 status = m_tuple->addItem("pLambda", m_pLambda);
00095 status = m_tuple->addItem("pLambdabar", m_pLambdabar);
00096
00097 } else {
00098 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
00099 }
00100 }
00101
00102 if(service("THistSvc", m_thsvc).isFailure()) {
00103 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00104 return StatusCode::FAILURE;
00105 }
00106
00107
00108
00109 TH1F* hrxy = new TH1F("Rxy", "Rxy distribution", 110, -1., 10.);
00110 if(m_thsvc->regHist("/DQAHist/JsiLL/hrxy", hrxy).isFailure()) {
00111 log << MSG::ERROR << "Couldn't register Rxy" << endreq;
00112 }
00113 TH1F* hz = new TH1F("z", "z distribution", 200, -100., 100.);
00114 if(m_thsvc->regHist("/DQAHist/JsiLL/hz", hz).isFailure()) {
00115 log << MSG::ERROR << "Couldn't register z" << endreq;
00116 }
00117
00118 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00119 return StatusCode::SUCCESS;
00120
00121
00122 }
00123
00124
00125 StatusCode JsiLL::execute() {
00126
00127 MsgStream log(msgSvc(), name());
00128 log << MSG::INFO << "in execute()" << endreq;
00129
00130
00131
00132
00133 setFilterPassed(false);
00134
00135 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00136 m_runNo=eventHeader->runNumber();
00137 m_event=eventHeader->eventNumber();
00138 log << MSG::DEBUG <<"run, evtnum = "
00139 << m_runNo << " , "
00140 << m_event <<endreq;
00141
00142
00143 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00144 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00145 << evtRecEvent->totalCharged() << " , "
00146 << evtRecEvent->totalNeutral() << " , "
00147 << evtRecEvent->totalTracks() <<endreq;
00148
00149 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00150
00151 if(evtRecEvent->totalNeutral()>100) {
00152 return StatusCode::SUCCESS;
00153 }
00154
00155 Vint iGood, iplus, iminus;
00156 iGood.clear();
00157 iplus.clear();
00158 iminus.clear();
00159 Vp4 ppip, ppim;
00160 ppip.clear();
00161 ppim.clear();
00162
00163 Hep3Vector xorigin(0,0,0);
00164
00165 IVertexDbSvc* vtxsvc;
00166 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00167 if(vtxsvc->isVertexValid()){
00168 double* dbv = vtxsvc->PrimaryVertex();
00169 double* vv = vtxsvc->SigmaPrimaryVertex();
00170 xorigin.setX(dbv[0]);
00171 xorigin.setY(dbv[1]);
00172 xorigin.setZ(dbv[2]);
00173 }
00174
00175 int nCharge = 0;
00176 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00177 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00178 if(!(*itTrk)->isMdcTrackValid()) continue;
00179 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00180 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00181
00182 double pch =mdcTrk->p();
00183 double x0 =mdcTrk->x();
00184 double y0 =mdcTrk->y();
00185 double z0 =mdcTrk->z();
00186 double phi0=mdcTrk->helix(1);
00187 double xv=xorigin.x();
00188 double yv=xorigin.y();
00189 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
00190
00191 if(fabs(z0) >= m_vz0cut) continue;
00192 if(Rxy >= m_vr0cut) continue;
00193
00194
00195
00196 TH1 *h(0);
00197 if (m_thsvc->getHist("/DQAHist/JsiLL/hrxy", h).isSuccess()) {
00198 h->Fill(Rxy);
00199 } else {
00200 log << MSG::ERROR << "Couldn't retrieve hrxy" << endreq;
00201 }
00202 if (m_thsvc->getHist("/DQAHist/JsiLL/hz", h).isSuccess()) {
00203 h->Fill(z0);
00204 } else {
00205 log << MSG::ERROR << "Couldn't retrieve hz" << endreq;
00206 }
00207
00208
00209 iGood.push_back(i);
00210 nCharge += mdcTrk->charge();
00211 if (mdcTrk->charge() > 0) {
00212 iplus.push_back(i);
00213 } else {
00214 iminus.push_back(i);
00215 }
00216 }
00217
00218
00219
00220
00221 int nGood = iGood.size();
00222 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00223 if((nGood != 4)||(nCharge!=0)){
00224 return StatusCode::SUCCESS;
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234 int pidp0 = 5, pidp1 = 3, pidm0 = 5, pidm1 = 3;
00235
00236 RecMdcKalTrack *itTrkp = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
00237 RecMdcKalTrack *itTrkpip = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
00238 RecMdcKalTrack *itTrkpb = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
00239 RecMdcKalTrack *itTrkpim = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
00240
00241
00242
00243 if (itTrkp->p() < itTrkpip->p()){
00244 itTrkp = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
00245 itTrkpip = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
00246 pidp0 = 3;
00247 pidp1 = 5;
00248 }
00249 if (itTrkpb->p() < itTrkpim->p()){
00250 itTrkpb = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
00251 itTrkpim = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
00252 pidm0 = 3;
00253 pidm1 = 5;
00254 }
00255 if (itTrkp->p() < 0.7 || itTrkp->p() >1.1) return StatusCode::SUCCESS;
00256 if (itTrkpb->p() < 0.7 || itTrkpb->p() >1.1) return StatusCode::SUCCESS;
00257 if (itTrkpip->p() > 0.35) return StatusCode::SUCCESS;
00258 if (itTrkpim->p() > 0.35) return StatusCode::SUCCESS;
00259
00260
00261
00262 itTrkp->setPidType(RecMdcKalTrack::proton);
00263 itTrkpb->setPidType(RecMdcKalTrack::proton);
00264 itTrkpip->setPidType(RecMdcKalTrack::pion);
00265 itTrkpim->setPidType(RecMdcKalTrack::pion);
00266
00267
00268 m_chisq = 9999.0;
00269
00270 HepPoint3D vx(0., 0., 0.);
00271 HepSymMatrix Evx(3, 0);
00272 double bx = 1E+6;
00273 double by = 1E+6;
00274 double bz = 1E+6;
00275 Evx[0][0] = bx*bx;
00276 Evx[1][1] = by*by;
00277 Evx[2][2] = bz*bz;
00278
00279 VertexParameter vxpar;
00280 vxpar.setVx(vx);
00281 vxpar.setEvx(Evx);
00282
00283 VertexFit* vtxfita0 = VertexFit::instance();
00284 SecondVertexFit *vtxfita = SecondVertexFit::instance();
00285 VertexFit* vtxfitb0 = VertexFit::instance();
00286 SecondVertexFit *vtxfitb = SecondVertexFit::instance();
00287 VertexFit* vtxfit = VertexFit::instance();
00288
00289 WTrackParameter wpip = WTrackParameter(mpi, itTrkpip->getZHelix(), itTrkpip->getZError());
00290 WTrackParameter wpim = WTrackParameter(mpi, itTrkpim->getZHelix(), itTrkpim->getZError());
00291 WTrackParameter wp = WTrackParameter(mproton, itTrkp->getZHelixP(), itTrkp->getZErrorP());
00292 WTrackParameter wpb = WTrackParameter(mproton, itTrkpb->getZHelixP(), itTrkpb->getZErrorP());
00293
00294 vtxfita0->init();
00295 vtxfita0->AddTrack(0, wp);
00296 vtxfita0->AddTrack(1, wpim);
00297 vtxfita0->AddVertex(0, vxpar, 0, 1);
00298 if(!vtxfita0->Fit(0)) return StatusCode::SUCCESS;
00299 vtxfita0->Swim(0);
00300 vtxfita0->BuildVirtualParticle(0);
00301 vtxfita->init();
00302 vtxfita->AddTrack(0, vtxfita0->wVirtualTrack(0));
00303 vtxfita->setVpar(vtxfita0->vpar(0));
00304 if(!vtxfita->Fit()) return StatusCode::SUCCESS;
00305
00306 WTrackParameter wLambda = vtxfita->wpar();
00307
00308 vtxfitb0->init();
00309 vtxfitb0->AddTrack(0, wpb);
00310 vtxfitb0->AddTrack(1, wpip);
00311 vtxfitb0->AddVertex(0, vxpar, 0, 1);
00312 if(!vtxfitb0->Fit(0)) return StatusCode::SUCCESS;
00313 vtxfitb0->Swim(0);
00314 vtxfitb0->BuildVirtualParticle(0);
00315
00316 vtxfitb->init();
00317 vtxfitb->AddTrack(0, vtxfitb0->wVirtualTrack(0));
00318 vtxfitb->setVpar(vtxfitb0->vpar(0));
00319 if(!vtxfitb->Fit()) return StatusCode::SUCCESS;
00320
00321 WTrackParameter wLambdabar = vtxfitb->wpar();
00322
00323 vtxfit->init();
00324 vtxfit->AddTrack(0, wLambda);
00325 vtxfit->AddTrack(1, wLambdabar);
00326 vtxfit->AddVertex(0, vxpar,0, 1);
00327 if(!vtxfit->Fit(0)) return StatusCode::SUCCESS;
00328 vtxfit->Swim(0);
00329 WTrackParameter wwLambda = vtxfit->wtrk(0);
00330 WTrackParameter wwLambdabar = vtxfit->wtrk(1);
00331
00332
00333
00334
00335 KinematicFit* kmfit = KinematicFit::instance();
00336
00337
00338 HepLorentzVector ecms(0.034065,0.0,0.0,3.0969);
00339 const Hep3Vector u_cms = -ecms.boostVector();
00340
00341 kmfit->init();
00342 kmfit->AddTrack(0, wwLambda);
00343 kmfit->AddTrack(1, wwLambdabar);
00344 kmfit->AddFourMomentum(0, ecms);
00345
00346 if(!kmfit->Fit()) return StatusCode::SUCCESS;
00347 m_chisq = kmfit->chisq();
00348 if(m_chisq > 50) return StatusCode::SUCCESS;
00349 HepLorentzVector kmf_pLambda = kmfit->pfit(0);
00350 HepLorentzVector kmf_pLambdabar = kmfit->pfit(1);
00351
00352 kmf_pLambda.boost(u_cms);
00353 kmf_pLambdabar.boost(u_cms);
00354 m_mLambda = kmf_pLambda.m();
00355 m_mLambdabar = kmf_pLambdabar.m();
00356 m_pLambda = kmf_pLambda.rho();
00357 m_pLambdabar = kmf_pLambdabar.rho();
00358
00359 if(fabs(m_mLambda-1.1157)>0.003||fabs(m_mLambdabar-1.1157)>0.003) return StatusCode::SUCCESS;
00360
00361
00362 m_tuple->write();
00363
00364
00366
00367
00368
00369
00370 (*(evtRecTrkCol->begin()+iplus[0]))->setPartId(pidp0);
00371 (*(evtRecTrkCol->begin()+iplus[1]))->setPartId(pidp1);
00372 (*(evtRecTrkCol->begin()+iminus[0]))->setPartId(pidm0);
00373 (*(evtRecTrkCol->begin()+iminus[1]))->setPartId(pidm1);
00374
00375
00376
00377
00378
00379 (*(evtRecTrkCol->begin()+iplus[0]))->setQuality(0);
00380 (*(evtRecTrkCol->begin()+iplus[1]))->setQuality(0);
00381 (*(evtRecTrkCol->begin()+iminus[0]))->setQuality(0);
00382 (*(evtRecTrkCol->begin()+iminus[1]))->setQuality(0);
00383
00384
00385
00386
00387 setFilterPassed(true);
00389
00390 return StatusCode::SUCCESS;
00391
00392 }
00393
00394
00395 StatusCode JsiLL::finalize() {
00396
00397 MsgStream log(msgSvc(), name());
00398 log << MSG::INFO << "in finalize()" << endmsg;
00399 return StatusCode::SUCCESS;
00400 }
00401
00402