00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/IDataProviderSvc.h"
00005 #include "GaudiKernel/PropertyMgr.h"
00006
00007 #include "EventModel/EventModel.h"
00008 #include "EventModel/Event.h"
00009 #include "EventModel/EventHeader.h"
00010 #include "EvtRecEvent/EvtRecEvent.h"
00011 #include "EvtRecEvent/EvtRecTrack.h"
00012 #include "EvtRecEvent/EvtRecDTag.h"
00013
00014 #include "EvtRecEvent/EvtRecPi0.h"
00015
00016 #include "TMath.h"
00017 #include "GaudiKernel/INTupleSvc.h"
00018 #include "GaudiKernel/NTuple.h"
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/IHistogramSvc.h"
00021 #include "CLHEP/Vector/ThreeVector.h"
00022 #include "CLHEP/Vector/LorentzVector.h"
00023 #include "CLHEP/Vector/TwoVector.h"
00024 #include "EmcRawEvent/EmcDigi.h"
00025 #include "RawEvent/RawDataUtil.h"
00026 #include "MdcRawEvent/MdcDigi.h"
00027
00028 #include "VertexFit/IVertexDbSvc.h"
00029 #include "VertexFit/KinematicFit.h"
00030 #include "VertexFit/VertexFit.h"
00031 #include "VertexFit/Helix.h"
00032
00033 #include "DstEvent/TofHitStatus.h"
00034
00035
00036 #include "GaudiKernel/Bootstrap.h"
00037 #include "GaudiKernel/ISvcLocator.h"
00038 using CLHEP::Hep3Vector;
00039 using CLHEP::Hep2Vector;
00040 using CLHEP::HepLorentzVector;
00041 #include "CLHEP/Geometry/Point3D.h"
00042 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00043 typedef HepGeom::Point3D<double> HepPoint3D;
00044 #endif
00045 #include "CalibEventSelect/CalibEventSelect.h"
00046
00047 #include <vector>
00048 #include "mysql.h"
00049
00050 typedef std::vector<int> Vint;
00051 typedef std::vector<HepLorentzVector> Vp4;
00052
00053 const double mpi = 0.13957;
00054 const double mkaon = 0.49367;
00055 const double mproton = 0.93827;
00056
00057
00058 double Px(RecMdcKalTrack *trk){
00059 double phi= trk->fi0();
00060 double kappa= trk->kappa();
00061 double tanl = trk->tanl();
00062
00063 double px=-sin(phi)/fabs(kappa);
00064 return px;
00065 }
00066
00067 double Py(RecMdcKalTrack *trk){
00068 double phi= trk->fi0();
00069 double kappa= trk->kappa();
00070 double tanl = trk->tanl();
00071
00072 double py=cos(phi)/fabs(kappa);
00073
00074 return py;
00075 }
00076
00077
00078 double Pz(RecMdcKalTrack *trk){
00079 double phi= trk->fi0();
00080 double kappa= trk->kappa();
00081 double tanl = trk->tanl();
00082
00083 double pz=tanl/fabs(kappa);
00084 return pz;
00085 }
00086
00087 double P(RecMdcKalTrack *trk){
00088 double phi= trk->fi0();
00089 double kappa= trk->kappa();
00090 double tanl = trk->tanl();
00091
00092 double px=-sin(phi)/fabs(kappa);
00093 double py=cos(phi)/fabs(kappa);
00094 double pz=tanl/fabs(kappa);
00095
00096 return sqrt(px*px+py*py+pz*pz);
00097 }
00098
00099 double Phi(RecMdcKalTrack *trk){
00100 double phi0=trk->fi0();
00101 double kappa = trk->kappa();
00102 double pxy=0;
00103 if(kappa!=0) pxy = 1.0/fabs(kappa);
00104
00105
00106 double px = pxy * (-sin(phi0));
00107 double py = pxy * cos(phi0);
00108
00109 return atan2(py,px);
00110
00111 }
00112
00113
00114 int CalibEventSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
00115 88,88,100,100,112,112,128,128,140,140,
00116 160,160,160,160,176,176,176,176,208,208,
00117 208,208,240,240,240,240,256,256,256,256,
00118 288,288,288};
00120
00121 CalibEventSelect::CalibEventSelect(const std::string& name, ISvcLocator* pSvcLocator) :
00122 Algorithm(name, pSvcLocator) {
00123
00124
00125 declareProperty("Output", m_output = false);
00126 declareProperty("Display", m_display = false);
00127 declareProperty("PrintMonitor", m_printmonitor=false);
00128 declareProperty("SelectRadBhabha", m_selectRadBhabha=false);
00129 declareProperty("SelectGGEE", m_selectGGEE=false);
00130 declareProperty("SelectGG4Pi", m_selectGG4Pi=false);
00131 declareProperty("SelectRadBhabhaT", m_selectRadBhabhaT=false);
00132 declareProperty("SelectRhoPi", m_selectRhoPi=false);
00133 declareProperty("SelectKstarK", m_selectKstarK=false);
00134 declareProperty("SelectPPPiPi", m_selectPPPiPi=false);
00135 declareProperty("SelectRecoJpsi", m_selectRecoJpsi=false);
00136 declareProperty("SelectBhabha", m_selectBhabha=false);
00137 declareProperty("SelectDimu", m_selectDimu=false);
00138 declareProperty("SelectCosmic", m_selectCosmic=false);
00139 declareProperty("SelectGenProton", m_selectGenProton=false);
00140 declareProperty("SelectPsipRhoPi", m_selectPsipRhoPi=false);
00141 declareProperty("SelectPsipKstarK", m_selectPsipKstarK=false);
00142 declareProperty("SelectPsipPPPiPi", m_selectPsipPPPiPi=false);
00143 declareProperty("SelectPsippCand", m_selectPsippCand=false);
00144
00145 declareProperty("BhabhaScale", m_radbhabha_scale=1000);
00146 declareProperty("RadBhabhaScale", m_bhabha_scale=1000);
00147 declareProperty("DimuScale", m_dimu_scale=10);
00148 declareProperty("CosmicScale", m_cosmic_scale=3);
00149 declareProperty("ProtonScale", m_proton_scale=100);
00150
00151 declareProperty("CosmicLowp", m_cosmic_lowp=1.0);
00152
00153 declareProperty("WriteDst", m_writeDst=false);
00154 declareProperty("WriteRec", m_writeRec=false);
00155 declareProperty("Ecm", m_ecm=3.1);
00156 declareProperty("Vr0cut", m_vr0cut=1.0);
00157 declareProperty("Vz0cut", m_vz0cut=10.0);
00158 declareProperty("Pt0HighCut", m_pt0HighCut=5.0);
00159 declareProperty("Pt0LowCut", m_pt0LowCut=0.05);
00160 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
00161 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00162 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00163
00164
00165 }
00166
00167
00168 StatusCode CalibEventSelect::initialize() {
00169 MsgStream log(msgSvc(), name());
00170
00171 log << MSG::INFO << "in initialize()" << endmsg;
00172
00173 m_run=0;
00174
00175
00176
00177 h_nGood= histoSvc()->book("1D/nGoodtracks", 1, "num of good chaged tracks", 20, 0, 20);
00178 h_nCharge= histoSvc()->book("1D/nCharge", 1, "net charge", 20, -10, 10);
00179 h_pmaxobeam= histoSvc()->book("1D/pmax", 1, "pmax over beam ratio", 100, 0, 3);
00180 h_eopmax= histoSvc()->book("1D/eopmax", 1, "e over pmax ratio", 100, 0, 3);
00181 h_eopsec= histoSvc()->book("1D/eopsec", 1, "e over psec ratio", 100, 0, 3);
00182 h_deltap= histoSvc()->book("1D/deltap", 1, "pmax minus psec", 100, -3, 3);
00183 h_esumoecm= histoSvc()->book("1D/esumoverecm", 1, "total energy over ecm ratio", 100, 0, 3);
00184 h_egmax= histoSvc()->book("1D/egmax", 1, "most energetic photon energy", 100, 0, 3);
00185 h_deltaphi= histoSvc()->book("1D/deltaphi", 1, "phi_pmax minus phi_sec", 150, -4, 4);
00186 h_Pt= histoSvc()->book("1D/Pt", 1, "total Transverse Momentum", 200,-1, 4);
00187
00188 StatusCode sc;
00189
00190 log << MSG::INFO << "creating sub-algorithms...." << endreq;
00191
00192
00193
00194
00195
00196 if(m_writeDst) {
00197 sc = createSubAlgorithm( "EventWriter", "WriteDst", m_subalg1);
00198 if( sc.isFailure() ) {
00199 log << MSG::ERROR << "Error creating Sub-Algorithm WriteDst" <<endreq;
00200 return sc;
00201 } else {
00202 log << MSG::INFO << "Success creating Sub-Algorithm WriteDst" <<endreq;
00203 }
00204 }
00205
00206 if(m_writeRec) {
00207 sc = createSubAlgorithm( "EventWriter", "WriteRec", m_subalg2);
00208 if( sc.isFailure() ) {
00209 log << MSG::ERROR << "Error creating Sub-Algorithm WriteRec" <<endreq;
00210 return sc;
00211 } else {
00212 log << MSG::INFO << "Success creating Sub-Algorithm WriteRec" <<endreq;
00213 }
00214 }
00215
00216
00217 if(m_selectRadBhabha) {
00218 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabha", m_subalg3);
00219 if( sc.isFailure() ) {
00220 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabha" <<endreq;
00221 return sc;
00222 } else {
00223 log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabha" <<endreq;
00224 }
00225 }
00226
00227 if(m_selectGGEE) {
00228 sc = createSubAlgorithm( "EventWriter", "SelectGGEE", m_subalg4);
00229 if( sc.isFailure() ) {
00230 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGGEE" <<endreq;
00231 return sc;
00232 } else {
00233 log << MSG::INFO << "Success creating Sub-Algorithm SelectGGEE" <<endreq;
00234 }
00235 }
00236
00237 if(m_selectGG4Pi) {
00238 sc = createSubAlgorithm( "EventWriter", "SelectGG4Pi", m_subalg5);
00239 if( sc.isFailure() ) {
00240 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGG4Pi" <<endreq;
00241 return sc;
00242 } else {
00243 log << MSG::INFO << "Success creating Sub-Algorithm SelectGG4Pi" <<endreq;
00244 }
00245 }
00246
00247
00248 if(m_selectRadBhabhaT) {
00249 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabhaT", m_subalg6);
00250 if( sc.isFailure() ) {
00251 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
00252 return sc;
00253 } else {
00254 log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
00255 }
00256 }
00257
00258
00259 if(m_selectRhoPi) {
00260 sc = createSubAlgorithm( "EventWriter", "SelectRhoPi", m_subalg7);
00261 if( sc.isFailure() ) {
00262 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRhoPi" <<endreq;
00263 return sc;
00264 } else {
00265 log << MSG::INFO << "Success creating Sub-Algorithm SelectRhoPi" <<endreq;
00266 }
00267 }
00268
00269
00270 if(m_selectKstarK) {
00271 sc = createSubAlgorithm( "EventWriter", "SelectKstarK", m_subalg8);
00272 if( sc.isFailure() ) {
00273 log << MSG::ERROR << "Error creating Sub-Algorithm SelectKstarK" <<endreq;
00274 return sc;
00275 } else {
00276 log << MSG::INFO << "Success creating Sub-Algorithm SelectKstarK" <<endreq;
00277 }
00278 }
00279
00280
00281
00282 if(m_selectPPPiPi) {
00283 sc = createSubAlgorithm( "EventWriter", "SelectPPPiPi", m_subalg9);
00284 if( sc.isFailure() ) {
00285 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPPPiPi" <<endreq;
00286 return sc;
00287 } else {
00288 log << MSG::INFO << "Success creating Sub-Algorithm SelectPPPiPi" <<endreq;
00289 }
00290 }
00291
00292
00293 if(m_selectRecoJpsi) {
00294 sc = createSubAlgorithm( "EventWriter", "SelectRecoJpsi", m_subalg10);
00295 if( sc.isFailure() ) {
00296 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRecoJpsi" <<endreq;
00297 return sc;
00298 } else {
00299 log << MSG::INFO << "Success creating Sub-Algorithm SelectRecoJpsi" <<endreq;
00300 }
00301 }
00302
00303
00304 if(m_selectBhabha) {
00305 sc = createSubAlgorithm( "EventWriter", "SelectBhabha", m_subalg11);
00306 if( sc.isFailure() ) {
00307 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabha" <<endreq;
00308 return sc;
00309 } else {
00310 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabha" <<endreq;
00311 }
00312 }
00313
00314 if(m_selectDimu) {
00315 sc = createSubAlgorithm( "EventWriter", "SelectDimu", m_subalg12);
00316 if( sc.isFailure() ) {
00317 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimu" <<endreq;
00318 return sc;
00319 } else {
00320 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimu" <<endreq;
00321 }
00322 }
00323
00324 if(m_selectCosmic) {
00325 sc = createSubAlgorithm( "EventWriter", "SelectCosmic", m_subalg13);
00326 if( sc.isFailure() ) {
00327 log << MSG::ERROR << "Error creating Sub-Algorithm SelectCosmic" <<endreq;
00328 return sc;
00329 } else {
00330 log << MSG::INFO << "Success creating Sub-Algorithm SelectCosmic" <<endreq;
00331 }
00332 }
00333
00334 if(m_selectGenProton) {
00335 sc = createSubAlgorithm( "EventWriter", "SelectGenProton", m_subalg14);
00336 if( sc.isFailure() ) {
00337 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGenProton" <<endreq;
00338 return sc;
00339 } else {
00340 log << MSG::INFO << "Success creating Sub-Algorithm SelectGenProton" <<endreq;
00341 }
00342 }
00343
00344
00345 if(m_selectPsipRhoPi) {
00346 sc = createSubAlgorithm( "EventWriter", "SelectPsipRhoPi", m_subalg15);
00347 if( sc.isFailure() ) {
00348 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
00349 return sc;
00350 } else {
00351 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
00352 }
00353 }
00354
00355
00356 if(m_selectPsipKstarK) {
00357 sc = createSubAlgorithm( "EventWriter", "SelectPsipKstarK", m_subalg16);
00358 if( sc.isFailure() ) {
00359 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipKstarK" <<endreq;
00360 return sc;
00361 } else {
00362 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipKstarK" <<endreq;
00363 }
00364 }
00365
00366
00367 if(m_selectPsipPPPiPi) {
00368 sc = createSubAlgorithm( "EventWriter", "SelectPsipPPPiPi", m_subalg17);
00369 if( sc.isFailure() ) {
00370 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
00371 return sc;
00372 } else {
00373 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
00374 }
00375 }
00376
00377
00378 if(m_selectPsippCand) {
00379 sc = createSubAlgorithm( "EventWriter", "SelectPsippCand", m_subalg18);
00380 if( sc.isFailure() ) {
00381 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsippCand" <<endreq;
00382 return sc;
00383 } else {
00384 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsippCand" <<endreq;
00385 }
00386 }
00387
00388
00389
00390
00391 if(m_output) {
00392 StatusCode status;
00393 NTuplePtr nt0(ntupleSvc(), "FILE1/hadron");
00394 if ( nt0 ) m_tuple0 = nt0;
00395 else {
00396 m_tuple0 = ntupleSvc()->book ("FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example");
00397 if ( m_tuple0 ) {
00398 status = m_tuple0->addItem ("esum", m_esum);
00399 status = m_tuple0->addItem ("eemc", m_eemc);
00400 status = m_tuple0->addItem ("etot", m_etot);
00401 status = m_tuple0->addItem ("nGood", m_nGood);
00402 status = m_tuple0->addItem ("nCharge", m_nCharge);
00403 status = m_tuple0->addItem ("nGam", m_nGam);
00404 status = m_tuple0->addItem ("ptot", m_ptot);
00405 status = m_tuple0->addItem ("pp", m_pp);
00406 status = m_tuple0->addItem ("pm", m_pm);
00407 status = m_tuple0->addItem ("run", m_runnb);
00408 status = m_tuple0->addItem ("event", m_evtnb);
00409 status = m_tuple0->addItem ("maxE", m_maxE);
00410 status = m_tuple0->addItem ("secE", m_secE);
00411 status = m_tuple0->addItem ("dthe", m_dthe);
00412 status = m_tuple0->addItem ("dphi", m_dphi);
00413 status = m_tuple0->addItem ("mdcHit1", m_mdcHit1);
00414 status = m_tuple0->addItem ("mdcHit2", m_mdcHit2);
00415 }
00416 else {
00417 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple0) << endmsg;
00418 return StatusCode::FAILURE;
00419 }
00420 }
00421
00422 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
00423 if ( nt1 ) m_tuple1 = nt1;
00424 else {
00425 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
00426 if ( m_tuple1 ) {
00427 status = m_tuple1->addItem ("vx0", m_vx0);
00428 status = m_tuple1->addItem ("vy0", m_vy0);
00429 status = m_tuple1->addItem ("vz0", m_vz0);
00430 status = m_tuple1->addItem ("vr0", m_vr0);
00431 status = m_tuple1->addItem ("theta0", m_theta0);
00432 status = m_tuple1->addItem ("p0", m_p0);
00433 status = m_tuple1->addItem ("pt0", m_pt0);
00434 }
00435 else {
00436 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00437 return StatusCode::FAILURE;
00438 }
00439 }
00440
00441 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
00442 if ( nt2 ) m_tuple2 = nt2;
00443 else {
00444 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
00445 if ( m_tuple2 ) {
00446 status = m_tuple2->addItem ("dthe", m_dthe);
00447 status = m_tuple2->addItem ("dphi", m_dphi);
00448 status = m_tuple2->addItem ("dang", m_dang);
00449 status = m_tuple2->addItem ("eraw", m_eraw);
00450 }
00451 else {
00452 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00453 return StatusCode::FAILURE;
00454 }
00455 }
00456 }
00457
00458
00459
00460
00461 m_events=0;
00462 m_radBhabhaNumber=0;
00463 m_GGEENumber=0;
00464 m_GG4PiNumber=0;
00465 m_radBhabhaTNumber=0;
00466 m_rhoPiNumber=0;
00467 m_kstarKNumber=0;
00468 m_ppPiPiNumber=0;
00469 m_recoJpsiNumber=0;
00470 m_bhabhaNumber=0;
00471 m_dimuNumber=0;
00472 m_cosmicNumber=0;
00473 m_genProtonNumber=0;
00474 m_psipRhoPiNumber=0;
00475 m_psipKstarKNumber=0;
00476 m_psipPPPiPiNumber=0;
00477
00478 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00479 return StatusCode::SUCCESS;
00480
00481 }
00482
00483
00484 StatusCode CalibEventSelect::execute() {
00485
00486
00487
00488 MsgStream log(msgSvc(), name());
00489 log << MSG::INFO << "in execute()" << endreq;
00490
00491 if( m_writeDst) m_subalg1->execute();
00492 if( m_writeRec) m_subalg2->execute();
00493
00494
00495 m_isRadBhabha = false;
00496 m_isGGEE = false;
00497 m_isGG4Pi = false;
00498 m_isRadBhabhaT = false;
00499 m_isRhoPi = false;
00500 m_isKstarK = false;
00501 m_isRecoJpsi = false;
00502 m_isPPPiPi = false;
00503 m_isBhabha = false;
00504 m_isDimu = false;
00505 m_isCosmic = false;
00506 m_isGenProton = false;
00507 m_isPsipRhoPi = false;
00508 m_isPsipKstarK = false;
00509 m_isPsipPPPiPi = false;
00510 m_isPsippCand = false;
00511
00512 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00513 if(!eventHeader)
00514 {
00515 cout<<" eventHeader "<<endl;
00516 return StatusCode::FAILURE;
00517 }
00518
00519 int run=eventHeader->runNumber();
00520 int event=eventHeader->eventNumber();
00521
00522
00523 if(m_run !=run){
00524 m_run=run;
00525 double beamE=0;
00526 readbeamEfromDb(run,beamE);
00527 cout<<"new run is:"<<m_run<<endl;
00528 cout<<"beamE is:"<<beamE<<endl;
00529 if(beamE>0 && beamE<3)
00530 m_ecm=2*beamE;
00531 }
00532
00533
00534
00535 if(m_display && m_events%1000==0){
00536 cout<< m_events << " -------- run,event: "<<run<<","<<event<<endl;
00537 cout<<"m_ecm="<<m_ecm<<endl;
00538 }
00539 m_events++;
00540
00541 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00542 if(!evtRecEvent ) {
00543 cout<<" evtRecEvent "<<endl;
00544 return StatusCode::FAILURE;
00545 }
00546
00547 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00548 << evtRecEvent->totalCharged() << " , "
00549 << evtRecEvent->totalNeutral() << " , "
00550 << evtRecEvent->totalTracks() <<endreq;
00551 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00552 if(!evtRecTrkCol){
00553 cout<<" evtRecTrkCol "<<endl;
00554 return StatusCode::FAILURE;
00555 }
00556
00557 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size()) return StatusCode::SUCCESS;
00558
00559
00560
00561 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
00562 if ( ! recPi0Col ) {
00563 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
00564 return StatusCode::FAILURE;
00565 }
00566
00567
00568 int nPi0 = recPi0Col->size();
00569 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
00570 if(nPi0==1){
00571 double mpi0 = (*itpi0)->unconMass();
00572 if ( fabs( mpi0 - 0.135 )> 0.015 )
00573 nPi0=0;
00574 }
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 Vint iGood;
00588 iGood.clear();
00589 vector<int> iCP, iCN;
00590 iCP.clear();
00591 iCN.clear();
00592
00593 Vint iCosmicGood;
00594 iCosmicGood.clear();
00595
00596 int nCharge = 0;
00597 int nCosmicCharge =0;
00598
00599 for(int i = 0;i < evtRecEvent->totalCharged(); i++)
00600 {
00601
00602 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00603
00604
00605
00606 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00607 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
00608
00609
00610
00611 double vx0 = mdcTrk->x();
00612 double vy0 = mdcTrk->y();
00613 double vz0 = mdcTrk->z();
00614 double vr0 = mdcTrk->r();
00615 double theta0 = mdcTrk->theta();
00616 double p0 = P(mdcTrk);
00617 double pt0 = sqrt( pow(Px(mdcTrk),2)+pow(Py(mdcTrk),2));
00618
00619
00620 if(m_output) {
00621 m_vx0 = vx0;
00622 m_vy0 = vy0;
00623 m_vz0 = vz0;
00624 m_vr0 = vr0;
00625 m_theta0 = theta0;
00626 m_p0 = p0;
00627 m_pt0 = pt0;
00628 m_tuple1->write();
00629 }
00630
00631
00632
00633 Hep3Vector xorigin(0,0,0);
00634 IVertexDbSvc* vtxsvc;
00635 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00636 if(vtxsvc->isVertexValid()){
00637 double* dbv = vtxsvc->PrimaryVertex();
00638 double* vv = vtxsvc->SigmaPrimaryVertex();
00639 xorigin.setX(dbv[0]);
00640 xorigin.setY(dbv[1]);
00641 xorigin.setZ(dbv[2]);
00642 }
00643 HepVector a = mdcTrk->getZHelix();
00644 HepSymMatrix Ea = mdcTrk->getZError();
00645 HepPoint3D point0(0.,0.,0.);
00646 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
00647 VFHelix helixip(point0,a,Ea);
00648 helixip.pivot(IP);
00649 HepVector vecipa = helixip.a();
00650 double db=fabs(vecipa[0]);
00651 double dz=vecipa[3];
00652
00653
00654
00655
00656 if(fabs(dz) >= m_vz0cut) continue;
00657 if(db >= m_vr0cut) continue;
00658
00659
00660 if(p0>m_cosmic_lowp && p0<20){
00661 iCosmicGood.push_back((*itTrk)->trackId());
00662 nCosmicCharge += mdcTrk->charge();
00663 }
00664
00665
00666
00667 if(pt0 >= m_pt0HighCut) continue;
00668 if(pt0 <= m_pt0LowCut) continue;
00669
00670 iGood.push_back((*itTrk)->trackId());
00671 nCharge += mdcTrk->charge();
00672 if(mdcTrk->charge()>0)
00673 iCP.push_back((*itTrk)->trackId());
00674 else if(mdcTrk->charge()<0)
00675 iCN.push_back((*itTrk)->trackId());
00676
00677
00678 }
00679 int nGood = iGood.size();
00680 int nCosmicGood = iCosmicGood.size();
00681
00682
00683 Vint iGam;
00684 iGam.clear();
00685 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00686
00687 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00688 if(!(*itTrk)->isEmcShowerValid()) continue;
00689 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00690 double eraw = emcTrk->energy();
00691 double time = emcTrk->time();
00692 double costh = cos(emcTrk->theta());
00693 if(fabs(costh)<0.83 && eraw<0.025) continue;
00694 if(fabs(costh)>=0.83 && eraw<0.05) continue;
00695 if(time<0 || time>14) continue;
00696
00697 iGam.push_back((*itTrk)->trackId());
00698 }
00699 int nGam = iGam.size();
00700
00701
00702 Vint ipip, ipim;
00703 ipip.clear();
00704 ipim.clear();
00705 Vp4 ppip, ppim;
00706 ppip.clear();
00707 ppim.clear();
00708
00709
00710 double echarge = 0.;
00711 double ptot = 0.;
00712 double etot = 0.;
00713 double eemc = 0.;
00714 double pp = 0.;
00715 double pm = 0.;
00716 double pmax=0.0;
00717 double psec=0.0;
00718 double eccmax=0.0;
00719 double eccsec=0.0;
00720 double phimax=0.0;
00721 double phisec=0.0;
00722 double eopmax=0.0;
00723 double eopsec=0.0;
00724 Hep3Vector Pt_charge(0,0,0);
00725
00726 for(int i = 0; i < nGood; i++) {
00727 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00728
00729 double p3=0;
00730
00731
00732
00733 if((*itTrk)->isMdcKalTrackValid()) {
00734 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
00735 mdcTrk->setPidType(RecMdcKalTrack::electron);
00736
00737 ptot += P(mdcTrk);
00738
00739
00740 double phi= Phi(mdcTrk);
00741
00742
00743 HepLorentzVector ptrk;
00744 ptrk.setPx(Px(mdcTrk));
00745 ptrk.setPy(Py(mdcTrk));
00746 ptrk.setPz(Pz(mdcTrk));
00747 p3 = fabs(ptrk.mag());
00748
00749
00750
00751
00752
00753
00754
00755 Hep3Vector ptemp(Px(mdcTrk),Py(mdcTrk),0);
00756 Pt_charge+=ptemp;
00757
00758 double ecc=0.0;
00759 if((*itTrk)->isEmcShowerValid()) {
00760 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00761 ecc = emcTrk->energy();
00762 }
00763
00764 if( p3 > pmax){
00765 psec=pmax;
00766 eccsec=eccmax;
00767 phisec=phimax;
00768 pmax=p3;
00769 eccmax=ecc;
00770 phimax=phi;
00771 }
00772 else if( p3 < pmax && p3> psec){
00773 psec=p3;
00774 eccsec=ecc;
00775 phisec=phi;
00776 }
00777
00778
00779
00780
00781
00782 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00783 ptrk = ptrk.boost(-0.011,0,0);
00784
00785
00786 echarge += ptrk.e();
00787 etot += ptrk.e();
00788
00789 if(mdcTrk->charge() >0) {
00790 ppip.push_back(ptrk);
00791 pp = ptrk.rho();
00792 } else {
00793 ppim.push_back(ptrk);
00794 pm = ptrk.rho();
00795 }
00796
00797 }
00798
00799 if((*itTrk)->isEmcShowerValid()) {
00800
00801 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00802 double eraw = emcTrk->energy();
00803 double phiemc = emcTrk->phi();
00804 double the = emcTrk->theta();
00805
00806 HepLorentzVector pemc;
00807 pemc.setPx(eraw*sin(the)*cos(phiemc));
00808 pemc.setPy(eraw*sin(the)*sin(phiemc));
00809 pemc.setPz(eraw*cos(the));
00810 pemc.setE(eraw);
00811 pemc = pemc.boost(-0.011,0,0);
00812
00813
00814 etot += pemc.e();
00815
00816 }
00817
00818
00819
00820
00821 }
00822
00823 if(pmax!=0) eopmax=eccmax/pmax;
00824 if(psec!=0) eopsec=eccsec/psec;
00825
00826 eemc=eccmax+eccsec;
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839 Vp4 pGam;
00840 pGam.clear();
00841 double eneu=0;
00842 double egmax=0;
00843 Hep3Vector Pt_neu(0,0,0);
00844
00845 for(int i = 0; i < nGam; i++) {
00846 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00847 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00848 double eraw = emcTrk->energy();
00849 double phi = emcTrk->phi();
00850 double the = emcTrk->theta();
00851 HepLorentzVector ptrk;
00852 ptrk.setPx(eraw*sin(the)*cos(phi));
00853 ptrk.setPy(eraw*sin(the)*sin(phi));
00854 ptrk.setPz(eraw*cos(the));
00855 ptrk.setE(eraw);
00856 ptrk = ptrk.boost(-0.011,0,0);
00857 pGam.push_back(ptrk);
00858 eneu += ptrk.e();
00859 etot += ptrk.e();
00860 eemc += eraw;
00861
00862 Hep3Vector ptemp(eraw*sin(the)*cos(phi), eraw*sin(the)*sin(phi),0);
00863 Pt_neu+=ptemp;
00864
00865 if(ptrk.e()>egmax)
00866 egmax=ptrk.e();
00867
00868 }
00869
00870 double esum = echarge + eneu;
00871 Hep3Vector Pt=Pt_charge+Pt_neu;
00872
00873
00874 double phid=phimax-phisec;
00875 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00876
00877
00878
00879
00880
00881
00882 if( nGood == 2 && nCharge==0 && (m_selectBhabha || m_selectDimu) ){
00883
00884
00885 if( m_events%m_bhabha_scale == 0 && m_selectBhabha && echarge/m_ecm>0.8 && eopmax>0.8 && eopsec>0.8){
00886 m_isBhabha=true;
00887 m_bhabhaNumber++;
00888 }
00889
00890
00891
00892 if( m_events%m_dimu_scale == 0 && m_selectDimu && eemc/m_ecm<0.3){
00893
00894 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCP[0];
00895 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCN[0];
00896
00897 double time1=-99,depth1=-99;
00898 double time2=-99,depth2=-99;
00899 if( (*itp)->isTofTrackValid() ){
00900 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
00901 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
00902
00903 for(;iter_tof!=tofTrkCol.end();iter_tof++){
00904 TofHitStatus* status =new TofHitStatus;
00905 status->setStatus( (*iter_tof)->status() );
00906 if(status->is_cluster()){
00907 time1=(*iter_tof)->tof();
00908 }
00909 delete status;
00910 }
00911 }
00912
00913 if( (*itp)->isMucTrackValid() ){
00914 RecMucTrack* mucTrk=(*itp)->mucTrack();
00915 depth1=mucTrk->depth();
00916 }
00917
00918 if( (*itn)->isTofTrackValid() ){
00919 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
00920 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
00921
00922 for(;iter_tof!=tofTrkCol.end();iter_tof++){
00923 TofHitStatus* status =new TofHitStatus;
00924 status->setStatus( (*iter_tof)->status() );
00925 if(status->is_cluster()){
00926 time2=(*iter_tof)->tof();
00927 }
00928 delete status;
00929 }
00930 }
00931
00932 if( (*itn)->isMucTrackValid() ){
00933 RecMucTrack* mucTrk=(*itn)->mucTrack();
00934 depth2=mucTrk->depth();
00935 }
00936
00937
00938
00939
00940
00941
00942 if( ((fabs(pmax)/0.5/m_ecm>0.15 && fabs(pmax)/0.5/m_ecm<.75) || (fabs(psec)/0.5/m_ecm>0.15 && fabs(psec)/0.5/m_ecm<.75)) && (eopmax<0.4 || eopsec<0.4)
00943 && (depth1>=3 || depth2>=3)){
00944 m_isDimu=true;
00945 m_dimuNumber++;
00946 }
00947 }
00948
00949
00950
00951 }
00952
00953
00954
00955 if(m_selectCosmic && nCosmicGood == 2 && eemc/m_ecm<0.3){
00956
00957 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCosmicGood[0];
00958 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCosmicGood[1];
00959
00960 double time1=-99,depth1=-99;
00961 double time2=-99,depth2=-99;
00962 if( (*itp)->isTofTrackValid() ){
00963 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
00964 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
00965
00966 for(;iter_tof!=tofTrkCol.end();iter_tof++){
00967 TofHitStatus* status =new TofHitStatus;
00968 status->setStatus( (*iter_tof)->status() );
00969 if(status->is_cluster()){
00970 time1=(*iter_tof)->tof();
00971 }
00972 delete status;
00973 }
00974 }
00975
00976 if( (*itp)->isMucTrackValid() ){
00977 RecMucTrack* mucTrk=(*itp)->mucTrack();
00978 depth1=mucTrk->depth();
00979 }
00980
00981 if( (*itn)->isTofTrackValid() ){
00982 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
00983 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
00984
00985 for(;iter_tof!=tofTrkCol.end();iter_tof++){
00986 TofHitStatus* status =new TofHitStatus;
00987 status->setStatus( (*iter_tof)->status() );
00988 if(status->is_cluster()){
00989 time2=(*iter_tof)->tof();
00990 }
00991 delete status;
00992 }
00993 }
00994
00995 if( (*itn)->isMucTrackValid() ){
00996 RecMucTrack* mucTrk=(*itn)->mucTrack();
00997 depth2=mucTrk->depth();
00998 }
00999
01000
01001
01002 if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 ){
01003 m_isCosmic=true;
01004 m_cosmicNumber++;
01005 }
01006
01007
01008 }
01009
01010
01011
01012
01013 if(m_events%m_proton_scale ==0 ){
01014
01015 bool protoncand=false;
01016 for(int i=0; i<nGood; i++){
01017
01018 EvtRecTrackIterator iter = evtRecTrkCol->begin() + iGood[i];
01019 RecMdcKalTrack *mdcTrk = (*iter)->mdcKalTrack();
01020 mdcTrk->setPidType(RecMdcKalTrack::proton);
01021
01022 double pp = P(mdcTrk);
01023 double chiP=-99;
01024 if((*iter)->isMdcDedxValid()){
01025 RecMdcDedx* dedxTrk=(*iter)->mdcDedx();
01026 chiP=dedxTrk->chiP();
01027 }
01028
01029 if( fabs(pp)<0.6 && fabs(chiP)<5){
01030 protoncand=true;
01031 break;
01032 }
01033 }
01034
01035 if( protoncand ){
01036 m_isGenProton=true;
01037 m_genProtonNumber++;
01038 }
01039
01040 }
01041
01042
01043
01044
01045
01046 if(m_printmonitor){
01047 h_nGood->fill(nGood);
01048 h_nCharge->fill(nCharge);
01049 h_pmaxobeam->fill(pmax/(m_ecm/2.0));
01050 h_eopmax->fill(eopmax);
01051 h_eopsec->fill(eopsec);
01052 h_deltap->fill(pmax-psec);
01053 h_esumoecm->fill(esum/m_ecm);
01054 h_egmax->fill(egmax);
01055 h_deltaphi->fill(phid);
01056 h_Pt->fill(Pt.mag());
01057 }
01058
01059
01060
01061
01062 if(nGood>1 && pmax/(m_ecm/2.0)>0.4 && eopmax>0.5 && esum/m_ecm>0.75 &&
01063 egmax>0.08 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi>2.86 )
01064 {
01065 m_isRadBhabha=true;
01066 m_radBhabhaNumber++;
01067 }
01068
01069 if( m_isRadBhabha )
01070 {
01071
01072 if(nGood==2 && nCharge==0 && eemc/m_ecm>=0.7 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 ){
01073
01074 if( fabs(fabs(pmax)-m_ecm/2.0)<0.1 && fabs(fabs(psec)-m_ecm/2.0)<0.1 ){
01075 if(m_events%1000==0){
01076 m_isRadBhabhaT=true;
01077 m_radBhabhaTNumber++;
01078 }
01079 }
01080 else{
01081 m_isRadBhabhaT=true;
01082 m_radBhabhaTNumber++;
01083 }
01084
01085 }
01086
01087
01088
01089 }
01090
01091
01092
01093 if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 && eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
01094 {
01095 m_isGGEE=true;
01096 m_GGEENumber++;
01097 }
01098
01099
01100 if(m_selectGG4Pi && nGood==4 && nCharge==0 && pmax/(m_ecm/2.0)>0.04 && pmax/(m_ecm/2.0)<0.9 && esum/m_ecm>0.04 && esum/m_ecm<0.75 && Pt.mag()<=0.2)
01101 {
01102 m_isGG4Pi=true;
01103 m_GG4PiNumber++;
01104 }
01105
01106
01107 if( (m_selectRhoPi || m_selectKstarK) && nGood == 2 && nCharge==0 && nPi0 == 1){
01108
01109 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iGood[0];
01110 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iGood[1];
01111
01112
01113 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
01114
01115 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
01116 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
01117
01118 HepLorentzVector p4trk1;
01119 p4trk1.setPx(Px(trk1));
01120 p4trk1.setPy(Py(trk1));
01121 p4trk1.setPz(Pz(trk1));
01122 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
01123
01124 HepLorentzVector p4trk2;
01125 p4trk2.setPx(Px(trk2));
01126 p4trk2.setPy(Py(trk2));
01127 p4trk2.setPz(Pz(trk2));
01128 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
01129
01130
01131 HepLorentzVector p4trk3;
01132 p4trk3.setPx(Px(trk1));
01133 p4trk3.setPy(Py(trk1));
01134 p4trk3.setPz(Pz(trk1));
01135 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01136
01137 HepLorentzVector p4trk4;
01138 p4trk4.setPx(Px(trk2));
01139 p4trk4.setPy(Py(trk2));
01140 p4trk4.setPz(Pz(trk2));
01141 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01142
01143
01144
01145 itpi0 = recPi0Col->begin();
01146 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
01147 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
01148 HepLorentzVector p4pi0 = p4gam1+p4gam2;
01149
01150 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
01151 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
01152
01153 double rhopimass=p4total2.m();
01154 double kstarkmass=p4total.m();
01155 if( (kstarkmass > 3.0 && kstarkmass < 3.15) || (rhopimass > 3.0 && rhopimass < 3.15) ){
01156
01157
01158
01159
01160 double eop1=0.0,eop2=0.0;
01161 if( (*itone)->isEmcShowerValid() ){
01162 RecEmcShower *emcTrk = (*itone)->emcShower();
01163 double etrk=emcTrk->energy();
01164
01165 if(P(trk1)!=0){
01166 eop1 = etrk/P(trk1);
01167
01168 }
01169 }
01170
01171 if( (*ittwo)->isEmcShowerValid() ){
01172 RecEmcShower *emcTrk = (*ittwo)->emcShower();
01173 double etrk=emcTrk->energy();
01174
01175 if(P(trk2)!=0){
01176 eop2 = etrk/P(trk2);
01177
01178 }
01179 }
01180
01181 if(eop1<0.8 && eop2<0.8){
01182
01183 if(rhopimass>3.0 && rhopimass<3.15){
01184
01185
01186 HepLorentzVector ecms(0.034,0,0,m_ecm);
01187
01188
01189 WTrackParameter wvpiTrk1, wvpiTrk2;
01190 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
01191 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01192
01193 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
01194 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
01195 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
01196 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
01197
01198 KinematicFit * kmfit = KinematicFit::instance();
01199 kmfit->init();
01200 kmfit->AddTrack(0, wvpiTrk1);
01201 kmfit->AddTrack(1, wvpiTrk2);
01202 kmfit->AddTrack(2, 0.0, gShr1);
01203 kmfit->AddTrack(3, 0.0, gShr2);
01204 kmfit->AddFourMomentum(0, ecms);
01205
01206 bool oksq1 = kmfit->Fit();
01207 double chi1 = kmfit->chisq();
01208
01209
01210 if(oksq1 && chi1<=60){
01211 m_isRhoPi = true;
01212 m_rhoPiNumber++;
01213 }
01214 }
01215
01216
01217 if(kstarkmass>3.0 && kstarkmass<3.15){
01218
01219
01220 double mkstarp=0, mkstarm=0;
01221 if( trk1->charge() >0 ){
01222 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
01223 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
01224 mkstarp = p4kstarp.m();
01225 mkstarm = p4kstarm.m();
01226 }
01227 else{
01228 HepLorentzVector p4kstarm = p4trk1 + p4pi0;
01229 HepLorentzVector p4kstarp = p4trk2 + p4pi0;
01230 mkstarp = p4kstarp.m();
01231 mkstarm = p4kstarm.m();
01232 }
01233
01234 if ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1) || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) ){
01235
01236
01237 HepLorentzVector ecms(0.034,0,0,m_ecm);
01238
01239 WTrackParameter wvpiTrk1, wvpiTrk2;
01240 wvpiTrk1 = WTrackParameter(mkaon, trk1->getZHelix(), trk1->getZError());
01241 wvpiTrk2 = WTrackParameter(mkaon, trk2->getZHelix(), trk2->getZError());
01242
01243 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
01244 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
01245 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
01246 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
01247
01248 KinematicFit * kmfit = KinematicFit::instance();
01249 kmfit->init();
01250 kmfit->AddTrack(0, wvpiTrk1);
01251 kmfit->AddTrack(1, wvpiTrk2);
01252 kmfit->AddTrack(2, 0.0, gShr1);
01253 kmfit->AddTrack(3, 0.0, gShr2);
01254 kmfit->AddFourMomentum(0, ecms);
01255
01256 bool oksq1 = kmfit->Fit();
01257 double chi1 = kmfit->chisq();
01258
01259
01260 if(oksq1 && chi1<=60){
01261 m_isKstarK = true;
01262 m_kstarKNumber++;
01263 }
01264
01265 }
01266
01267 }
01268
01269 }
01270
01271
01272
01273 }
01274 }
01275
01276
01277
01278 }
01279
01280
01281 if(m_selectPPPiPi && nGood ==4 && nCharge == 0 && nPi0==0){
01282
01283 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iCP[0];
01284 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iCP[1];
01285 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iCN[0];
01286 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iCN[1];
01287 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
01288 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
01289 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
01290 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
01291
01292 HepLorentzVector p4trkpp;
01293 HepLorentzVector p4trkpm;
01294 HepLorentzVector p4trkpip;
01295 HepLorentzVector p4trkpim;
01296
01297
01298
01299 trk1->setPidType(RecMdcKalTrack::proton);
01300 trk2->setPidType(RecMdcKalTrack::pion);
01301 trk3->setPidType(RecMdcKalTrack::proton);
01302 trk4->setPidType(RecMdcKalTrack::pion);
01303
01304
01305 p4trkpp.setPx(Px(trk1));
01306 p4trkpp.setPy(Py(trk1));
01307 p4trkpp.setPz(Pz(trk1));
01308 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
01309
01310 p4trkpm.setPx(Px(trk3));
01311 p4trkpm.setPy(Py(trk3));
01312 p4trkpm.setPz(Pz(trk3));
01313 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
01314
01315
01316 p4trkpip.setPx(Px(trk2));
01317 p4trkpip.setPy(Py(trk2));
01318 p4trkpip.setPz(Pz(trk2));
01319 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01320
01321 p4trkpim.setPx(Px(trk4));
01322 p4trkpim.setPy(Py(trk4));
01323 p4trkpim.setPz(Pz(trk4));
01324 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
01325
01326 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01327
01328
01329
01330
01331 trk1->setPidType(RecMdcKalTrack::proton);
01332 trk2->setPidType(RecMdcKalTrack::pion);
01333 trk3->setPidType(RecMdcKalTrack::pion);
01334 trk4->setPidType(RecMdcKalTrack::proton);
01335
01336
01337 p4trkpp.setPx(Px(trk1));
01338 p4trkpp.setPy(Py(trk1));
01339 p4trkpp.setPz(Pz(trk1));
01340 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
01341
01342 p4trkpm.setPx(Px(trk4));
01343 p4trkpm.setPy(Py(trk4));
01344 p4trkpm.setPz(Pz(trk4));
01345 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
01346
01347
01348 p4trkpip.setPx(Px(trk2));
01349 p4trkpip.setPy(Py(trk2));
01350 p4trkpip.setPz(Pz(trk2));
01351 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01352
01353 p4trkpim.setPx(Px(trk3));
01354 p4trkpim.setPy(Py(trk3));
01355 p4trkpim.setPz(Pz(trk3));
01356 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
01357
01358 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01359
01360
01361
01362
01363 trk1->setPidType(RecMdcKalTrack::pion);
01364 trk2->setPidType(RecMdcKalTrack::proton);
01365 trk3->setPidType(RecMdcKalTrack::proton);
01366 trk4->setPidType(RecMdcKalTrack::pion);
01367
01368
01369 p4trkpp.setPx(Px(trk2));
01370 p4trkpp.setPy(Py(trk2));
01371 p4trkpp.setPz(Pz(trk2));
01372 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
01373
01374 p4trkpm.setPx(Px(trk3));
01375 p4trkpm.setPy(Py(trk3));
01376 p4trkpm.setPz(Pz(trk3));
01377 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
01378
01379
01380 p4trkpip.setPx(Px(trk1));
01381 p4trkpip.setPy(Py(trk1));
01382 p4trkpip.setPz(Pz(trk1));
01383 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01384
01385 p4trkpim.setPx(Px(trk4));
01386 p4trkpim.setPy(Py(trk4));
01387 p4trkpim.setPz(Pz(trk4));
01388 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
01389
01390 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01391
01392
01393
01394
01395 trk1->setPidType(RecMdcKalTrack::pion);
01396 trk2->setPidType(RecMdcKalTrack::proton);
01397 trk3->setPidType(RecMdcKalTrack::pion);
01398 trk4->setPidType(RecMdcKalTrack::proton);
01399
01400
01401 p4trkpp.setPx(Px(trk2));
01402 p4trkpp.setPy(Py(trk2));
01403 p4trkpp.setPz(Pz(trk2));
01404 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
01405
01406 p4trkpm.setPx(Px(trk4));
01407 p4trkpm.setPy(Py(trk4));
01408 p4trkpm.setPz(Pz(trk4));
01409 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
01410
01411
01412 p4trkpip.setPx(Px(trk1));
01413 p4trkpip.setPy(Py(trk1));
01414 p4trkpip.setPz(Pz(trk1));
01415 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01416
01417 p4trkpim.setPx(Px(trk3));
01418 p4trkpim.setPy(Py(trk3));
01419 p4trkpim.setPz(Pz(trk3));
01420 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
01421
01422 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01423
01424
01425
01426
01427 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
01428 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
01429
01430
01431
01432
01433 double chi1, chi2, chi3, chi4;
01434 int type=0;
01435 double chimin=-999;
01436 HepLorentzVector ecms(0.034,0,0,m_ecm);
01437
01438 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 ;
01439
01440 {
01441 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
01442 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01443 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
01444 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
01445
01446
01447 KinematicFit * kmfit = KinematicFit::instance();
01448 kmfit->init();
01449 kmfit->AddTrack(0, wvpiTrk1);
01450 kmfit->AddTrack(1, wvpiTrk2);
01451 kmfit->AddTrack(2, wvpiTrk3);
01452 kmfit->AddTrack(3, wvpiTrk4);
01453 kmfit->AddFourMomentum(0, ecms);
01454
01455 bool oksq1 = kmfit->Fit();
01456 chi1 = kmfit->chisq();
01457 if(oksq1) {
01458 chimin=chi1;
01459 type=1;
01460 }
01461
01462 }
01463
01464
01465 {
01466 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
01467 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01468 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
01469 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
01470
01471
01472 KinematicFit * kmfit = KinematicFit::instance();
01473 kmfit->init();
01474 kmfit->AddTrack(0, wvpiTrk1);
01475 kmfit->AddTrack(1, wvpiTrk2);
01476 kmfit->AddTrack(2, wvpiTrk3);
01477 kmfit->AddTrack(3, wvpiTrk4);
01478 kmfit->AddFourMomentum(0, ecms);
01479
01480 bool oksq1 = kmfit->Fit();
01481 chi2 = kmfit->chisq();
01482 if(oksq1){
01483 if(type==0){
01484 chimin=chi2;
01485 type=2;
01486 }
01487 else if(chi2<chimin){
01488 chimin=chi2;
01489 type=2;
01490 }
01491 }
01492
01493 }
01494
01495 {
01496 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
01497 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
01498 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
01499 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
01500
01501 KinematicFit * kmfit = KinematicFit::instance();
01502 kmfit->init();
01503 kmfit->AddTrack(0, wvpiTrk1);
01504 kmfit->AddTrack(1, wvpiTrk2);
01505 kmfit->AddTrack(2, wvpiTrk3);
01506 kmfit->AddTrack(3, wvpiTrk4);
01507
01508 kmfit->AddFourMomentum(0, ecms);
01509
01510 bool oksq1 = kmfit->Fit();
01511 chi3 = kmfit->chisq();
01512 if(oksq1){
01513 if(type==0){
01514 chimin=chi3;
01515 type=3;
01516 }
01517 else if(chi3<chimin){
01518 chimin=chi3;
01519 type=3;
01520 }
01521 }
01522
01523
01524 }
01525
01526
01527 {
01528 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
01529 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
01530 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
01531 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
01532
01533 KinematicFit * kmfit = KinematicFit::instance();
01534 kmfit->init();
01535 kmfit->AddTrack(0, wvpiTrk1);
01536 kmfit->AddTrack(1, wvpiTrk2);
01537 kmfit->AddTrack(2, wvpiTrk3);
01538 kmfit->AddTrack(3, wvpiTrk4);
01539
01540 kmfit->AddFourMomentum(0, ecms);
01541
01542 bool oksq1 = kmfit->Fit();
01543 chi4 = kmfit->chisq();
01544
01545 if(oksq1){
01546 if(type==0){
01547 chimin=chi4;
01548 type=4;
01549 }
01550 else if(chi4<chimin){
01551 chimin=chi4;
01552 type=4;
01553 }
01554 }
01555
01556
01557 }
01558
01559 if(type!=0 && chimin<100){
01560 m_isPPPiPi = true;
01561 m_ppPiPiNumber++;
01562 }
01563
01564
01565
01566
01567 }
01568
01569 }
01570
01571
01572
01573 if( (m_selectPsipRhoPi || m_selectPsipKstarK || m_selectPsipPPPiPi) && (nGood==4 || nGood==6) && nCharge==0 ){
01574 EvtRecTrackIterator pione, pitwo;
01575 RecMdcKalTrack *pitrk1;
01576 RecMdcKalTrack *pitrk2;
01577
01578 double bestmass=1.0;
01579 int pindex,nindex;
01580 vector<int> iJPsiP,iJPsiN;
01581 for(int ip=0; ip<iCP.size(); ip++){
01582 for(int in=0; in<iCN.size();in++){
01583 pione = evtRecTrkCol->begin() + iCP[ip];
01584 pitwo = evtRecTrkCol->begin() + iCN[in];
01585 pitrk1 = (*pione)->mdcKalTrack();
01586 pitrk2 = (*pitwo)->mdcKalTrack();
01587 Hep3Vector p1(Px(pitrk1), Py(pitrk1), Pz(pitrk1));
01588 Hep3Vector p2(Px(pitrk2), Py(pitrk2), Pz(pitrk2));
01589 double E1=sqrt( pow(P(pitrk1),2)+mpi*mpi);
01590 double E2=sqrt( pow(P(pitrk2),2)+mpi*mpi);
01591 double recomass = sqrt(pow(3.686-E1-E2,2)- (-(p1+p2)).mag2() );
01592
01593 if(fabs(recomass-3.096)< fabs(bestmass-3.096)){
01594 bestmass=recomass;
01595 pindex=ip;
01596 nindex=in;
01597 }
01598 }
01599 }
01600
01601
01602
01603 pione = evtRecTrkCol->begin() + iCP[pindex];
01604 pitwo = evtRecTrkCol->begin() + iCN[nindex];
01605 pitrk1 = (*pione)->mdcKalTrack();
01606 pitrk2 = (*pitwo)->mdcKalTrack();
01607
01608
01609
01610 double jpsicharge=0;
01611 for(int ip=0; ip<iCP.size(); ip++){
01612 if(ip!=pindex){
01613 iJPsiP.push_back(iCP[ip]);
01614 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCP[ip];
01615 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
01616 jpsicharge+=trktmp->charge();
01617 }
01618 }
01619
01620 for(int in=0; in<iCN.size(); in++){
01621 if(in!=nindex){
01622 iJPsiN.push_back(iCN[in]);
01623 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCN[in];
01624 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
01625 jpsicharge+=trktmp->charge();
01626 }
01627 }
01628
01629
01630 HepLorentzVector ecms(0.034,0,0,m_ecm);
01631
01632
01633 if( (m_selectPsipRhoPi || m_selectPsipKstarK) && (bestmass>3.075 && bestmass<3.120) && nGood==4 && jpsicharge==0 && nPi0==1){
01634
01635 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
01636 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiN[0];
01637
01638
01639 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
01640
01641 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
01642 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
01643
01644 HepLorentzVector p4trk1;
01645 p4trk1.setPx(Px(trk1));
01646 p4trk1.setPy(Py(trk1));
01647 p4trk1.setPz(Pz(trk1));
01648 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
01649
01650 HepLorentzVector p4trk2;
01651 p4trk2.setPx(Px(trk2));
01652 p4trk2.setPy(Py(trk2));
01653 p4trk2.setPz(Pz(trk2));
01654 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
01655
01656
01657 HepLorentzVector p4trk3;
01658 p4trk3.setPx(Px(trk1));
01659 p4trk3.setPy(Py(trk1));
01660 p4trk3.setPz(Pz(trk1));
01661 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01662
01663 HepLorentzVector p4trk4;
01664 p4trk4.setPx(Px(trk2));
01665 p4trk4.setPy(Py(trk2));
01666 p4trk4.setPz(Pz(trk2));
01667 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01668
01669
01670 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
01671 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
01672 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
01673 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
01674 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
01675 RecEmcShower *gShr3 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
01676 RecEmcShower *gShr4 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
01677
01678
01679 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
01680 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
01681 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
01682 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
01683 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
01684
01685
01686 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
01687 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
01688 double mkstarp = p4kstarp.m();
01689 double mkstarm = p4kstarm.m();
01690
01691 if( (p4total.m() > 3.0 && p4total.m() < 3.15) || (p4total2.m() > 3.0 && p4total2.m() < 3.15) ){
01692
01693
01694
01695
01696
01697
01698
01699 WTrackParameter wvpiTrk1, wvpiTrk2;
01700 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
01701 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01702
01703 WTrackParameter wvkTrk1, wvkTrk2;
01704 wvkTrk1 = WTrackParameter(mkaon, trk1->getZHelixK(), trk1->getZErrorK());
01705 wvkTrk2 = WTrackParameter(mkaon, trk2->getZHelixK(), trk2->getZErrorK());
01706
01707
01708 WTrackParameter wvpiTrk3, wvpiTrk4;
01709 wvpiTrk3 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
01710 wvpiTrk4 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
01711
01712 if(m_selectPsipRhoPi){
01713 KinematicFit * kmfit = KinematicFit::instance();
01714 kmfit->init();
01715 kmfit->AddTrack(0, wvpiTrk1);
01716 kmfit->AddTrack(1, wvpiTrk2);
01717 kmfit->AddTrack(2, 0.0, gShr1);
01718 kmfit->AddTrack(3, 0.0, gShr2);
01719 kmfit->AddTrack(4, wvpiTrk3);
01720 kmfit->AddTrack(5, wvpiTrk4);
01721 kmfit->AddFourMomentum(0, ecms);
01722
01723 bool oksq1 = kmfit->Fit();
01724 double chi1 = kmfit->chisq();
01725
01726
01727 if(oksq1 && chi1>0 && chi1<100){
01728 m_isPsipRhoPi = true;
01729 m_psipRhoPiNumber++;
01730 }
01731 }
01732 if(m_selectPsipKstarK){
01733 KinematicFit * kmfit2 = KinematicFit::instance();
01734 kmfit2->init();
01735 kmfit2->AddTrack(0, wvkTrk1);
01736 kmfit2->AddTrack(1, wvkTrk2);
01737 kmfit2->AddTrack(2, 0.0, gShr3);
01738 kmfit2->AddTrack(3, 0.0, gShr4);
01739 kmfit2->AddTrack(4, wvpiTrk3);
01740 kmfit2->AddTrack(5, wvpiTrk4);
01741 kmfit2->AddFourMomentum(0, ecms);
01742
01743
01744 bool oksq2 = kmfit2->Fit();
01745 double chi2 = kmfit2->chisq();
01746
01747 if(oksq2 && chi2>0 && chi2<200 &&
01748 ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1)
01749 || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) )){
01750 m_isPsipKstarK = true;
01751 m_psipKstarKNumber++;
01752 }
01753
01754 }
01755
01756
01757 }
01758
01759 }
01760
01761 }
01762
01763
01764
01765
01766 if(m_selectPsipPPPiPi && (bestmass>3.075 && bestmass<3.120) && nGood==6 && jpsicharge==0 && nPi0==0){
01767
01768
01769 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
01770 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiP[1];
01771 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iJPsiN[0];
01772 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iJPsiN[1];
01773 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
01774 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
01775 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
01776 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
01777
01778 HepLorentzVector p4trkpp;
01779 HepLorentzVector p4trkpm;
01780 HepLorentzVector p4trkpip;
01781 HepLorentzVector p4trkpim;
01782
01783
01784
01785 trk1->setPidType(RecMdcKalTrack::proton);
01786 trk2->setPidType(RecMdcKalTrack::pion);
01787 trk3->setPidType(RecMdcKalTrack::proton);
01788 trk4->setPidType(RecMdcKalTrack::pion);
01789
01790
01791 p4trkpp.setPx(Px(trk1));
01792 p4trkpp.setPy(Py(trk1));
01793 p4trkpp.setPz(Pz(trk1));
01794 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
01795
01796 p4trkpm.setPx(Px(trk3));
01797 p4trkpm.setPy(Py(trk3));
01798 p4trkpm.setPz(Pz(trk3));
01799 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
01800
01801
01802 p4trkpip.setPx(Px(trk2));
01803 p4trkpip.setPy(Py(trk2));
01804 p4trkpip.setPz(Pz(trk2));
01805 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01806
01807 p4trkpim.setPx(Px(trk4));
01808 p4trkpim.setPy(Py(trk4));
01809 p4trkpim.setPz(Pz(trk4));
01810 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
01811
01812 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01813
01814
01815
01816
01817 trk1->setPidType(RecMdcKalTrack::proton);
01818 trk2->setPidType(RecMdcKalTrack::pion);
01819 trk3->setPidType(RecMdcKalTrack::pion);
01820 trk4->setPidType(RecMdcKalTrack::proton);
01821
01822
01823 p4trkpp.setPx(Px(trk1));
01824 p4trkpp.setPy(Py(trk1));
01825 p4trkpp.setPz(Pz(trk1));
01826 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
01827
01828 p4trkpm.setPx(Px(trk4));
01829 p4trkpm.setPy(Py(trk4));
01830 p4trkpm.setPz(Pz(trk4));
01831 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
01832
01833
01834 p4trkpip.setPx(Px(trk2));
01835 p4trkpip.setPy(Py(trk2));
01836 p4trkpip.setPz(Pz(trk2));
01837 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
01838
01839 p4trkpim.setPx(Px(trk3));
01840 p4trkpim.setPy(Py(trk3));
01841 p4trkpim.setPz(Pz(trk3));
01842 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
01843
01844 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01845
01846
01847
01848
01849 trk1->setPidType(RecMdcKalTrack::pion);
01850 trk2->setPidType(RecMdcKalTrack::proton);
01851 trk3->setPidType(RecMdcKalTrack::proton);
01852 trk4->setPidType(RecMdcKalTrack::pion);
01853
01854
01855 p4trkpp.setPx(Px(trk2));
01856 p4trkpp.setPy(Py(trk2));
01857 p4trkpp.setPz(Pz(trk2));
01858 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
01859
01860 p4trkpm.setPx(Px(trk3));
01861 p4trkpm.setPy(Py(trk3));
01862 p4trkpm.setPz(Pz(trk3));
01863 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
01864
01865
01866 p4trkpip.setPx(Px(trk1));
01867 p4trkpip.setPy(Py(trk1));
01868 p4trkpip.setPz(Pz(trk1));
01869 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01870
01871 p4trkpim.setPx(Px(trk4));
01872 p4trkpim.setPy(Py(trk4));
01873 p4trkpim.setPz(Pz(trk4));
01874 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
01875
01876 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01877
01878
01879
01880
01881 trk1->setPidType(RecMdcKalTrack::pion);
01882 trk2->setPidType(RecMdcKalTrack::proton);
01883 trk3->setPidType(RecMdcKalTrack::pion);
01884 trk4->setPidType(RecMdcKalTrack::proton);
01885
01886
01887 p4trkpp.setPx(Px(trk2));
01888 p4trkpp.setPy(Py(trk2));
01889 p4trkpp.setPz(Pz(trk2));
01890 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
01891
01892 p4trkpm.setPx(Px(trk4));
01893 p4trkpm.setPy(Py(trk4));
01894 p4trkpm.setPz(Pz(trk4));
01895 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
01896
01897
01898 p4trkpip.setPx(Px(trk1));
01899 p4trkpip.setPy(Py(trk1));
01900 p4trkpip.setPz(Pz(trk1));
01901 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
01902
01903 p4trkpim.setPx(Px(trk3));
01904 p4trkpim.setPy(Py(trk3));
01905 p4trkpim.setPz(Pz(trk3));
01906 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
01907
01908 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
01909
01910 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
01911 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
01912
01913
01914
01915 double chi1, chi2, chi3, chi4;
01916 int type=0;
01917 double chimin=-999;
01918
01919
01920 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 , wvpiTrk5, wvpiTrk6 ;
01921
01922 {
01923 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
01924 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01925 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
01926 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
01927 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
01928 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
01929
01930
01931
01932 KinematicFit * kmfit = KinematicFit::instance();
01933 kmfit->init();
01934 kmfit->AddTrack(0, wvpiTrk1);
01935 kmfit->AddTrack(1, wvpiTrk2);
01936 kmfit->AddTrack(2, wvpiTrk3);
01937 kmfit->AddTrack(3, wvpiTrk4);
01938 kmfit->AddTrack(4, wvpiTrk5);
01939 kmfit->AddTrack(5, wvpiTrk6);
01940 kmfit->AddFourMomentum(0, ecms);
01941
01942 bool oksq1 = kmfit->Fit();
01943 chi1 = kmfit->chisq();
01944 if(oksq1){
01945 chimin=chi1;
01946 type=1;
01947 }
01948
01949 }
01950
01951
01952 {
01953 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
01954 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
01955 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
01956 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
01957 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
01958 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
01959
01960
01961
01962 KinematicFit * kmfit = KinematicFit::instance();
01963 kmfit->init();
01964 kmfit->AddTrack(0, wvpiTrk1);
01965 kmfit->AddTrack(1, wvpiTrk2);
01966 kmfit->AddTrack(2, wvpiTrk3);
01967 kmfit->AddTrack(3, wvpiTrk4);
01968 kmfit->AddTrack(4, wvpiTrk5);
01969 kmfit->AddTrack(5, wvpiTrk6);
01970 kmfit->AddFourMomentum(0, ecms);
01971
01972 bool oksq1 = kmfit->Fit();
01973 chi2 = kmfit->chisq();
01974 if(oksq1){
01975 if(type==0){
01976 chimin=chi2;
01977 type=2;
01978 }
01979 else if(chi2<chimin){
01980 chimin=chi2;
01981 type=2;
01982 }
01983
01984 }
01985
01986 }
01987
01988 {
01989 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
01990 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
01991 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
01992 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
01993 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
01994 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
01995
01996
01997
01998 KinematicFit * kmfit = KinematicFit::instance();
01999 kmfit->init();
02000 kmfit->AddTrack(0, wvpiTrk1);
02001 kmfit->AddTrack(1, wvpiTrk2);
02002 kmfit->AddTrack(2, wvpiTrk3);
02003 kmfit->AddTrack(3, wvpiTrk4);
02004 kmfit->AddTrack(4, wvpiTrk5);
02005 kmfit->AddTrack(5, wvpiTrk6);
02006 kmfit->AddFourMomentum(0, ecms);
02007
02008 bool oksq1 = kmfit->Fit();
02009 chi3 = kmfit->chisq();
02010 if(oksq1){
02011 if(type==0){
02012 chimin=chi3;
02013 type=3;
02014 }
02015 else if(chi3<chimin){
02016 chimin=chi3;
02017 type=3;
02018 }
02019 }
02020
02021 }
02022
02023 {
02024 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
02025 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
02026 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
02027 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
02028 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
02029 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
02030
02031
02032 KinematicFit * kmfit = KinematicFit::instance();
02033 kmfit->init();
02034 kmfit->AddTrack(0, wvpiTrk1);
02035 kmfit->AddTrack(1, wvpiTrk2);
02036 kmfit->AddTrack(2, wvpiTrk3);
02037 kmfit->AddTrack(3, wvpiTrk4);
02038 kmfit->AddTrack(4, wvpiTrk5);
02039 kmfit->AddTrack(5, wvpiTrk6);
02040 kmfit->AddFourMomentum(0, ecms);
02041
02042 bool oksq1 = kmfit->Fit();
02043 chi4 = kmfit->chisq();
02044 if(oksq1){
02045 if(type==0){
02046 chimin=chi4;
02047 type=4;
02048 }
02049 else if(chi4<chimin){
02050 chimin=chi4;
02051 type=4;
02052 }
02053 }
02054
02055 }
02056
02057
02058 if(chimin>0 && chimin<200){
02059 m_isPsipPPPiPi = true;
02060 m_psipPPPiPiNumber++;
02061 }
02062
02063 }
02064
02065 }
02066
02067
02068
02069
02070 }
02071
02072
02073
02074
02075
02076 if(m_selectPsippCand){
02077
02078 SmartDataPtr<EvtRecDTagCol> evtRecDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
02079 if ( ! evtRecDTagCol ) {
02080 log << MSG::FATAL << "Could not find EvtRecDTagCol" << endreq;
02081 return StatusCode::FAILURE;
02082 }
02083 if(evtRecDTagCol->size()>0){
02084
02085 EvtRecDTagCol::iterator m_iterbegin=evtRecDTagCol->begin();
02086 EvtRecDTagCol::iterator m_iterend=evtRecDTagCol->end();
02087 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
02088
02089 if( ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
02090 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0Pi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
02091 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
02092 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
02093 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
02094 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
02095 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPi && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
02096 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
02097 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
02098 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPiPi&& fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ) {
02099 m_isPsippCand = true;
02100 m_psippCandNumber++;
02101 }
02102
02103
02104 }
02105
02106
02107 }
02108
02109 }
02110
02111
02112
02113 if( m_selectRadBhabha && m_isRadBhabha ) m_subalg3->execute();
02114 if( m_selectGGEE && m_isGGEE ) m_subalg4->execute();
02115 if( m_selectGG4Pi && m_isGG4Pi ) m_subalg5->execute();
02116 if( m_selectRadBhabhaT && m_isRadBhabhaT ) m_subalg6->execute();
02117 if( m_selectRhoPi && m_isRhoPi ) m_subalg7->execute();
02118 if( m_selectKstarK && m_isKstarK ) m_subalg8->execute();
02119 if( m_selectPPPiPi && m_isPPPiPi ) m_subalg9->execute();
02120 if( m_selectRecoJpsi && m_isRecoJpsi ) m_subalg10->execute();
02121 if( m_selectBhabha && m_isBhabha ) m_subalg11->execute();
02122 if( m_selectDimu && m_isDimu ) m_subalg12->execute();
02123 if( m_selectCosmic && m_isCosmic ) m_subalg13->execute();
02124 if( m_selectGenProton && m_isGenProton ) m_subalg14->execute();
02125 if( m_selectPsipRhoPi && m_isPsipRhoPi ) m_subalg15->execute();
02126 if( m_selectPsipKstarK && m_isPsipKstarK ) m_subalg16->execute();
02127 if( m_selectPsipPPPiPi && m_isPsipPPPiPi ) m_subalg17->execute();
02128 if( m_selectPsippCand && m_isPsippCand ) m_subalg18->execute();
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152 return StatusCode::SUCCESS;
02153
02154 }
02155
02156
02157 StatusCode CalibEventSelect::finalize() {
02158
02159 MsgStream log(msgSvc(), name());
02160 log << MSG::INFO << "in finalize()" << endmsg;
02161
02162 cout<<"Total events: "<<m_events<<endl;
02163
02164
02165 if(m_selectRadBhabha) {
02166 cout << "Selected rad bhabha: " << m_radBhabhaNumber << endl;
02167 }
02168
02169
02170 if(m_selectGGEE) {
02171 cout << "Selected ggee events: " << m_GGEENumber << endl;
02172 }
02173
02174 if(m_selectGG4Pi) {
02175 cout << "Selected gg4pi events: " << m_GG4PiNumber << endl;
02176 }
02177
02178 if(m_selectRadBhabhaT) {
02179 cout << "Selected rad bhabha tight: " << m_radBhabhaTNumber << endl;
02180 }
02181
02182 if(m_selectRhoPi) {
02183 cout << "Selected rhopi events: " << m_rhoPiNumber << endl;
02184 }
02185
02186 if(m_selectKstarK) {
02187 cout << "Selected kstark events: " << m_kstarKNumber << endl;
02188 }
02189
02190 if(m_selectPPPiPi) {
02191 cout << "Selected pppipi events: " << m_ppPiPiNumber << endl;
02192 }
02193
02194 if(m_selectRecoJpsi) {
02195 cout << "Selected recoil jsi events: " << m_recoJpsiNumber << endl;
02196 }
02197
02198
02199 if(m_selectBhabha) {
02200 cout << "Selected bhabha events: " << m_bhabhaNumber << endl;
02201 }
02202
02203
02204 if(m_selectDimu) {
02205 cout << "Selected dimu events: " << m_dimuNumber << endl;
02206 }
02207
02208
02209 if(m_selectCosmic) {
02210 cout << "Selected cosmic events: " << m_cosmicNumber << endl;
02211 }
02212
02213 if(m_selectGenProton) {
02214 cout << "Selected generic proton events: " << m_genProtonNumber << endl;
02215 }
02216
02217 if(m_selectPsipRhoPi) {
02218 cout << "Selected recoil rhopi events: " << m_psipRhoPiNumber << endl;
02219 }
02220
02221 if(m_selectPsipKstarK) {
02222 cout << "Selected recoil kstark events: " << m_psipKstarKNumber << endl;
02223 }
02224
02225 if(m_selectPsipPPPiPi) {
02226 cout << "Selected recoil pppipi events: " << m_psipPPPiPiNumber << endl;
02227 }
02228
02229 if(m_selectPsippCand) {
02230 cout << "Selected psi'' candi events: " << m_psippCandNumber << endl;
02231 }
02232
02233 return StatusCode::SUCCESS;
02234 }
02235
02236 bool CalibEventSelect::WhetherSector(double ph,double ph1,double ph2) {
02237 double phi1=min(ph1,ph2);
02238 double phi2=max(ph1,ph2);
02239 double delta=0.0610865;
02240 if((phi2-phi1)<CLHEP::pi){
02241 phi1 -=delta;
02242 phi2 +=delta;
02243 if(phi1<0.) phi1 += CLHEP::twopi;
02244 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
02245 double tmp1=min(phi1,phi2);
02246 double tmp2=max(phi1,phi2);
02247 phi1=tmp1;
02248 phi2=tmp2;
02249 }
02250 else{
02251 phi1 +=delta;
02252 phi2 -=delta;
02253 }
02254 if((phi2-phi1)<CLHEP::pi){
02255 if(ph<=phi2&&ph>=phi1) return true;
02256 else return false;
02257 }
02258 else{
02259 if(ph>=phi2||ph<=phi1) return true;
02260 else return false;
02261 }
02262 }
02263
02264
02265
02266 void CalibEventSelect::readbeamEfromDb(int runNo, double & beamE){
02267
02268 const char host[] = "202.122.37.69";
02269 const char user[] = "guest";
02270 const char passwd[] = "guestpass";
02271 const char db[] = "run";
02272 unsigned int port_num = 3306;
02273
02274
02275 MYSQL* mysql = mysql_init(NULL);
02276 mysql = mysql_real_connect(mysql, host, user, passwd, db, port_num,
02277 NULL,
02278 0);
02279
02280 if (mysql == NULL) {
02281 fprintf(stderr, "can not open database: RunInfo for run %d lum\n", runNo);
02282 }
02283
02284 char stmt[1024];
02285
02286 snprintf(stmt, 1024,
02287 "select BER_PRB, BPR_PRB "
02288 "from RunParams where run_number = %d", runNo);
02289 if (mysql_real_query(mysql, stmt, strlen(stmt))) {
02290 fprintf(stderr, "query error\n");
02291 }
02292
02293 MYSQL_RES* result_set = mysql_store_result(mysql);
02294 MYSQL_ROW row = mysql_fetch_row(result_set);
02295 if (!row) {
02296 fprintf(stderr, "cannot find data for RunNo %d\n", runNo);
02297 return ;
02298 }
02299
02300 double E_E=0, E_P=0;
02301 sscanf(row[0], "%lf", &E_E);
02302 sscanf(row[1], "%lf", &E_P);
02303
02304 beamE=(E_E+E_P)/2.0;
02305 }