00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/PropertyMgr.h"
00007 #include "GaudiKernel/Bootstrap.h"
00008 #include "GaudiKernel/IJobOptionsSvc.h"
00009 #include "EventModel/EventModel.h"
00010 #include "EventModel/Event.h"
00011 #include "EventModel/EventHeader.h"
00012 #include "EvtRecEvent/EvtRecEvent.h"
00013 #include "EvtRecEvent/EvtRecTrack.h"
00014
00015 #include "MdcRecEvent/RecMdcTrack.h"
00016 #include "EmcRecEventModel/RecEmcEventModel.h"
00017 #include "TofRecEvent/RecTofTrack.h"
00018 #include "DstEvent/TofHitStatus.h"
00019 #include "MucRawEvent/MucDigi.h"
00020 #include "MucRecEvent/RecMucTrack.h"
00021
00022 #include "EventPreSelect/DimuPreSelect.h"
00023
00024 DimuPreSelect::DimuPreSelect()
00025 {
00026 m_output=false;
00027 m_totevent = m_currun = m_curevent = 0;
00028 for(int i = 0; i < 20; i++) m_cutpass[i] = 0;
00029 m_subpass[0] = m_subpass[1] = m_subpass[2] = m_subpass[3] = 0;
00030 m_totalpass = 0;
00031
00032 m_propMgr.declareProperty("CmsEnergy", m_ecms = 3.686);
00033 m_propMgr.declareProperty("Vr0Cut", m_vr0cut = 1.0);
00034 m_propMgr.declareProperty("Vz0Cut", m_vz0cut = 5.0);
00035 m_propMgr.declareProperty("PUpCut", m_pcut_up = 2.0);
00036 m_propMgr.declareProperty("PDownCut", m_pcut_down = 0.5);
00037 m_propMgr.declareProperty("PSymCut", m_psymcut = 0.5);
00038 m_propMgr.declareProperty("TCut", m_tcut = 4);
00039 m_propMgr.declareProperty("EUpCut", m_ecut_up = 1.0);
00040 m_propMgr.declareProperty("EDownCut", m_ecut_down = 0.1);
00041 m_propMgr.declareProperty("DThetaCut", m_dthetacut = 0.05);
00042 m_propMgr.declareProperty("DPhiCut", m_dphicut = 0.4);
00043 m_propMgr.declareProperty("PartSelect", m_partselect= 0);
00044 m_propMgr.declareProperty("MuDigiCut", m_mudigicut = 6);
00045 m_propMgr.declareProperty("MuTrkCut", m_mutrkcut = 1);
00046
00047 IJobOptionsSvc* jobSvc;
00048 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
00049 jobSvc->setMyProperties("DimuSelect", &m_propMgr);
00050 }
00051
00052 void DimuPreSelect::BookNtuple(NTuple::Tuple*& tuple)
00053 {
00054 m_output=true;
00055 m_passtuple=tuple;
00056 if(!m_passtuple) {
00057 std::cerr << "invalid ntuple in DimuPreSelect!" << std::endl;
00058 } else {
00059 m_passtuple->addItem ("run", m_run);
00060 m_passtuple->addItem ("event", m_event);
00061 m_passtuple->addItem ("part", m_part);
00062 m_passtuple->addItem ("c1", m_c1);
00063 m_passtuple->addItem ("c2", m_c2);
00064 m_passtuple->addItem ("r1", m_r1);
00065 m_passtuple->addItem ("r2", m_r2);
00066 m_passtuple->addItem ("z1", m_z1);
00067 m_passtuple->addItem ("z2", m_z2);
00068 m_passtuple->addItem ("p1", m_p1);
00069 m_passtuple->addItem ("p2", m_p2);
00070 m_passtuple->addItem ("t1", m_t1);
00071 m_passtuple->addItem ("t2", m_t2);
00072 m_passtuple->addItem ("e1", m_e1);
00073 m_passtuple->addItem ("e2", m_e2);
00074 m_passtuple->addItem ("theta1",m_theta1);
00075 m_passtuple->addItem ("theta2",m_theta2);
00076 m_passtuple->addItem ("phi1", m_phi1);
00077 m_passtuple->addItem ("phi2", m_phi2);
00078 m_passtuple->addItem ("mudigi",m_mudigi);
00079 m_passtuple->addItem ("mutrk", m_mutrk);
00080
00081 m_passtuple->addItem ("zeroC", m_zeroCFlag);
00082 m_passtuple->addItem ("vtRZ", m_vtRZFlag);
00083 m_passtuple->addItem ("pLim", m_pLimFlag);
00084 m_passtuple->addItem ("pSym", m_pSymFlag);
00085 m_passtuple->addItem ("tLim", m_tLimFlag);
00086 m_passtuple->addItem ("eLim", m_eLimFlag);
00087 m_passtuple->addItem ("eBB", m_eBBFlag);
00088 m_passtuple->addItem ("partSlt",m_partFlag);
00089 m_passtuple->addItem ("muDigiN",m_mudigiFlag);
00090 m_passtuple->addItem ("muTrkN", m_mutrkFlag);
00091
00092 m_passtuple->addItem ("mdc", m_mdcFlag);
00093 m_passtuple->addItem ("tof", m_tofFlag);
00094 m_passtuple->addItem ("emc", m_emcFlag);
00095 m_passtuple->addItem ("muc", m_mucFlag);
00096 }
00097 }
00098
00099 int DimuPreSelect::IsDimu()
00100 {
00101 IMessageSvc* msgSvc;
00102 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00103 MsgStream log(msgSvc, "DimuPreSelect");
00104
00105 IDataProviderSvc* eventSvc;
00106 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00107
00108 m_mdcPass = false;
00109 m_tofPass = false;
00110 m_emcPass = false;
00111 m_mucPass = false;
00112
00113
00114 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
00115 if(!mdcTrackCol) {
00116 log << MSG::ERROR << "Could not find RecMdcTrackCol!" << endreq;
00117 return false;
00118 }
00119 log << MSG::INFO << "MDC tracks:\t" << mdcTrackCol->size() << endreq;
00120
00121 if( mdcTrackCol->size() != 2 ) return 0;
00122
00123 m_cutpass[0] += 1;
00124
00125 double c1, c2, r1, r2, z1, z2, p1, p2;
00126 c1 = c2 = r1 = r2 = z1 = z2 = p1 = p2 = 0.;
00127
00128 c1 = (*mdcTrackCol)[0]->charge(); c2 = (*mdcTrackCol)[1]->charge();
00129 r1 = (*mdcTrackCol)[0]->r(); r2 = (*mdcTrackCol)[1]->r();
00130 z1 = (*mdcTrackCol)[0]->z(); z2 = (*mdcTrackCol)[1]->z();
00131 p1 = (*mdcTrackCol)[0]->p(); p2 = (*mdcTrackCol)[1]->p();
00132
00133
00134 bool mdcflag1 = c1 + c2 == 0;
00135 bool mdcflag2 = fabs(r1)<=m_vr0cut && fabs(z1)<m_vz0cut;
00136 bool mdcflag3 = fabs(r2)<=m_vr0cut && fabs(z2)<m_vz0cut;
00137 bool mdcflag4 = p1<m_pcut_up && p2<m_pcut_up;
00138
00139
00140 bool mdcflag5 = fabs( p1-p2 )/( p1+p2 ) < m_psymcut;
00141
00142 log << MSG::INFO << "r1:\t" << r1 << "\tz1:" << z1 << endreq;
00143 log << MSG::INFO << "r2:\t" << r2 << "\tz2:" << z2 << endreq;
00144 log << MSG::INFO << "p1:\t" << p1 << "\tp2:" << p2 << endreq;
00145
00146 double zeroCFlag,vtRZFlag,pLimFlag,pSymFlag;
00147 if( mdcflag1 ) { m_cutpass[1] += 1; zeroCFlag = 1; }
00148 if( mdcflag2 && mdcflag3 ) { m_cutpass[2] += 1; vtRZFlag = 1; }
00149 if( mdcflag4 ) { m_cutpass[3] += 1; pLimFlag = 1; }
00150 if( mdcflag5 ) { m_cutpass[4] += 1; pSymFlag = 1; }
00151
00152 if( mdcflag1 && mdcflag2 && mdcflag3 && mdcflag4 && mdcflag5 ) { m_mdcPass = true; m_subpass[0] += 1; }
00153 log << MSG::INFO << "MDC selection done!" << endreq;
00154
00155
00156 SmartDataPtr<RecTofTrackCol> tofTrackCol(eventSvc,"/Event/Recon/RecTofTrackCol");
00157 if(!tofTrackCol) {
00158 log << MSG::ERROR << "Could not find RecTofTrackCol!" << endreq;
00159 return 0;
00160 }
00161 log << MSG::INFO << "TOF tracks:\t" << tofTrackCol->size() << endreq;
00162
00163 double t1, t2;
00164 t1 = 0., t2 = 1000;
00165
00166 if(tofTrackCol->size() > 7 && tofTrackCol->size() < 21)
00167 {
00168 int goodtof = 0;
00169 for(int itof = 0; itof < tofTrackCol->size(); itof++)
00170 {
00171 TofHitStatus *status = new TofHitStatus;
00172 status->setStatus((*tofTrackCol)[itof]->status());
00173
00174 if( !(status->is_cluster()) ) { delete status; continue; }
00175 if(goodtof==0) t1 = (*tofTrackCol)[itof]->tof();
00176 if(goodtof==1) t2 = (*tofTrackCol)[itof]->tof();
00177
00178 goodtof++;
00179 delete status;
00180 }
00181 }
00182
00183 double tLimFlag;
00184 bool tofflag1 = fabs( t1-t2 ) < m_tcut;
00185 log << MSG::INFO << "t1:\t" << t1 << "\tt2:" << t2 << "dt:\t" << fabs(t1-t2) << endreq;
00186 if( tofflag1 ) { m_cutpass[5] += 1; tLimFlag = 1; m_tofPass = true; m_subpass[1] += 1; }
00187 log << MSG::INFO << "TOF selection done!" << endreq;
00188
00189
00190 SmartDataPtr<RecEmcShowerCol> emcShowerCol(eventSvc,"/Event/Recon/RecEmcShowerCol");
00191 if (!emcShowerCol) {
00192 log << MSG::ERROR << "Could not find RecEmcShowerCol!" << endreq;
00193 return 0;
00194 }
00195 log << MSG::INFO << "EMC showers:\t" << emcShowerCol->size() << endreq;
00196
00197 if( emcShowerCol->size() < 2 ) return 0;
00198
00199 double e1, e2, theta1, theta2, phi1, phi2;
00200 int part;
00201
00202 e1 = (*emcShowerCol)[0]->energy(); e2 = (*emcShowerCol)[1]->energy();
00203 theta1 = (*emcShowerCol)[0]->theta(); theta2 = (*emcShowerCol)[1]->theta();
00204 phi1 = (*emcShowerCol)[0]->phi(); phi2 = (*emcShowerCol)[1]->phi();
00205 part = (*emcShowerCol)[0]->module();
00206
00207 bool emcFlag1 = e1 < m_ecut_up && e1 > m_ecut_down;
00208 bool emcFlag2 = e2 < m_ecut_up && e2 > m_ecut_down;
00209 bool emcFlag3 = fabs(theta1 + theta2 - CLHEP::pi) < m_dthetacut;
00210 bool emcFlag4 = fabs(phi1 - phi2) - CLHEP::pi < m_dphicut;
00211 bool emcFlag5 = !m_partselect || (m_partselect==1&&part==1) || (m_partselect==2&&part!=1);
00212
00213 log << MSG::INFO << "e1:\t" << e1 << "\te2:\t" << e2 << endreq;
00214 log << MSG::INFO << "theta1:\t" << theta1 << "\ttheta2:\t" << theta2 << endreq;
00215 log << MSG::INFO << "phi1:\t" << phi1 << "\tphi2:\t" << phi2 << endreq;
00216 log << MSG::INFO << "part:\t" << part << "\tpartFlag:\t" << emcFlag5 << endreq;
00217
00218 double eLimFlag,eBBFlag,partFlag;
00219 if( emcFlag1 && emcFlag2 ) { m_cutpass[6] += 1; eLimFlag = 1; }
00220 if( emcFlag3 && emcFlag4 ) { m_cutpass[7] += 1; eBBFlag = 1; }
00221 if( emcFlag5 ) { m_cutpass[8] += 1; partFlag = 1; }
00222 if( emcFlag1 && emcFlag2 && emcFlag3 && emcFlag4 && emcFlag5 ) { m_emcPass = true; m_subpass[2] += 1; }
00223 log << MSG::INFO << "EMC selection done!" << endreq;
00224
00225
00226 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc,"/Event/Digi/MucDigiCol");
00227 if(!mucDigiCol) {
00228 log << MSG::ERROR << "Could not find MucDigiCol!" << endreq;
00229 return 0;
00230 }
00231 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc,"/Event/Recon/RecMucTrackCol");
00232 if (!mucTrackCol) {
00233 log << MSG::ERROR << "Could not find RecMucTrackCol" << endreq;
00234 return 0;
00235 }
00236
00237 int mudigiNum, mutrkNum;
00238 mudigiNum = mutrkNum = 0;
00239 mudigiNum = mucDigiCol->size(); mutrkNum = mucTrackCol->size();
00240
00241 bool mucflag1 = mudigiNum >= m_mudigicut;
00242 bool mucflag2 = mutrkNum >= m_mutrkcut;
00243
00244 log << MSG::INFO << "MUC digi:\t" << mudigiNum << "\tMUC track:\t" << mutrkNum << endreq;
00245
00246 double mudigiFlag,mutrkFlag;
00247 if( mucflag1 ) { m_cutpass[9] += 1; mudigiFlag = 1; }
00248 if( mucflag2 ) { m_cutpass[10] += 1; mutrkFlag = 1; }
00249 if( mucflag1 && mucflag2 ) { m_mucPass = true; m_subpass[3] += 1; }
00250 log << MSG::INFO << "MUC selection done!" << endreq;
00251
00252 if(m_output) {
00253 m_c1 = c1; m_c2 = c2; m_r1 = r1; m_r2 = r2; m_z1 = z1; m_z2 = z2; m_p1 = p1; m_p2 = p2;
00254 m_t1 = t1; m_t2 = t2;
00255 m_e1 = e1; m_e2 = e2; m_theta1 = theta1; m_theta2 = theta2; m_phi1 = phi1; m_phi2 = phi2;
00256 m_part = part;
00257 m_mudigi = mudigiNum; m_mutrk = mutrkNum;
00258 m_zeroCFlag = zeroCFlag; m_vtRZFlag = vtRZFlag; m_pLimFlag = pLimFlag; m_pSymFlag = pSymFlag;
00259 m_tLimFlag = tLimFlag; m_eLimFlag = eLimFlag; m_eBBFlag = eBBFlag; m_partFlag = partFlag;
00260 m_mudigiFlag = mudigiFlag; m_mutrkFlag = mutrkFlag;
00261
00262 m_mdcFlag = m_mdcPass; m_tofFlag = m_tofPass; m_emcFlag = m_emcPass; m_mucFlag = m_mucPass;
00263
00264 m_passtuple->write();
00265 }
00266
00267
00268 if( m_mdcPass && m_tofPass && m_emcPass && m_mucPass )
00269 {
00270 m_totalpass += 1;
00271 if( part==1 ) return 1;
00272 else return 2;
00273
00274 }
00275
00276 return 0;
00277
00278 }
00279
00280 void DimuPreSelect::Print()
00281 {
00282 const string str[3] = {"all","barrel","endcap"};
00283
00284 cout << "pass 0 - 2 MDC tracks : " << m_cutpass[0] << endl;
00285 cout << "pass 1 - 0 charge : " << m_cutpass[1] << endl;
00286 cout << "pass 2 - |r|<"<<m_vr0cut<<" && |z|<"<<m_vz0cut<<" : " << m_cutpass[2] << endl;
00287 cout << "pass 3 - p < "<<m_pcut_up<<" GeV/c : " << m_cutpass[3] << endl;
00288 cout << "pass 4 - p_sym < "<<m_psymcut<<" : " << m_cutpass[4] << endl;
00289 cout << "pass 5 - |t1-t1| < "<<m_tcut<<" ns : " << m_cutpass[5] << endl;
00290 cout << "pass 6 - "<<m_ecut_down<<" < e < "<<m_ecut_up<<" : " << m_cutpass[6] << endl;
00291 cout << "pass 7 - |dth|<"<<m_dthetacut<<" && |dphi|<"<<m_dphicut<<": " << m_cutpass[7] << endl;
00292 cout << "pass 8 - "<< str[(int)m_partselect] <<" part is selected : " << m_cutpass[8] << endl;
00293 cout << "pass 9 - mudigi >= "<<m_mudigicut<<" : " << m_cutpass[9] << endl;
00294 cout << "pass 10- mutrk >= "<<m_mutrkcut<<" : " << m_cutpass[10]<< endl;
00295
00296 cout << "pass MDC : " << m_subpass[0] << endl;
00297 cout << "pass TOF : " << m_subpass[1] << endl;
00298 cout << "pass EMC : " << m_subpass[2] << endl;
00299 cout << "pass MUC : " << m_subpass[3] << endl;
00300
00301 cout<<"Total event: " << m_totevent <<endl;
00302 cout<<"Dimu event: " << m_totalpass <<endl;
00303 }