00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "SD0TagAlg/SD0Tag.h"
00047 #include "SD0TagAlg/Sing.h"
00048 #include "SD0TagAlg/SingleBase.h"
00049 #include "DTagTool/DTagTool.h"
00050
00051 #include <iostream>
00052 #include <fstream>
00053 int nD0 = 0;
00054 int n1 = 0;
00055 int n2 = 0;
00056 int ND0 = 0;
00057 int NDp = 0;
00058
00059
00060
00062
00063 SD0Tag::SD0Tag(const std::string& name, ISvcLocator* pSvcLocator) :
00064 Algorithm(name, pSvcLocator) {
00065
00066 declareProperty("PID_FLAG", PID_flag = 1);
00067 declareProperty("MC_sample", m_MC_sample = 1);
00068
00069 declareProperty("m_Seperate_Charge", Seperate_Charge = 1);
00070 declareProperty("m_Charge_default", Charge_default = 1);
00071
00072
00073 }
00074
00075
00076
00077 StatusCode SD0Tag::initialize(){
00078
00079
00080 MsgStream log(msgSvc(), name());
00081
00082 log << MSG::INFO << "in initialize()" << endmsg;
00083
00084 StatusCode status;
00085
00086 ifstream readrunEcm;
00087
00088
00089 readrunEcm.open("../share/psipp_ecms_calrunbyrun_runno_3648runs.dat");
00090 cout<<"readrunEcm = "<<readrunEcm.is_open()<<endl;
00091 for(int i = 0; i < 3467; i++) {
00092 readrunEcm>>p_run[i]; readrunEcm>>p_Ecm[i];
00093 }
00094
00095
00096 NTuplePtr nt1(ntupleSvc(), "FILE1/tag");
00097 if ( nt1 ) m_tuple1 = nt1;
00098 else {
00099 m_tuple1 = ntupleSvc()->book ("FILE1/tag", CLID_ColumnWiseTuple, "tag N-Tuple example");
00100 if ( m_tuple1 ) {
00101 status = m_tuple1->addItem ("tagmode", m_tagmode);
00102 status = m_tuple1->addItem ("mass_bc", m_mass_bc);
00103 status = m_tuple1->addItem ("delE_tag", m_delE_tag);
00104
00105 status = m_tuple1->addItem ("cosmic_ok", m_cosmic_ok);
00106 status = m_tuple1->addItem ("EGam_max_0", m_EGam_max_0);
00107 status = m_tuple1->addItem ("nGood", m_nGood);
00108 status = m_tuple1->addItem ("nGam", m_nGam);
00109 status = m_tuple1->addItem ("runNo", m_runNo);
00110 status = m_tuple1->addItem ("event", m_event);
00111 }
00112 else {
00113 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) <<endmsg;
00114 return StatusCode::FAILURE;
00115 }
00116 }
00117
00118 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00119 return StatusCode::SUCCESS;
00120 }
00121
00122
00123
00124
00125
00126
00127 StatusCode SD0Tag::execute() {
00128
00129 MsgStream log(msgSvc(), name());
00130
00131 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00132 int runNo=eventHeader->runNumber();
00133 int event=eventHeader->eventNumber();
00134
00135 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00136 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00137
00138 nD0++;
00139
00140
00141
00142
00143 double Ebeam = 3.773/2.0;
00144
00145 if(runNo >=11414 && runNo <= 23500) {
00146 for(int i = 0; i < 3467; i++) {
00147 if(runNo == p_run[i]) {Ebeam = p_Ecm[i]/2.0;}
00148 }
00149 }
00150
00151 if(runNo < 0)
00152 {
00153 int irun = abs(runNo);
00154
00155
00156
00157 }
00158
00159 if(nD0%1000==0) cout<<"SD0TagAlg-13-11-15pp = "<<nD0<<" "<<runNo<<" "<<event<<" "<<Ebeam*2<<endl;
00160
00161 Hep3Vector xorigin(0,0,0);
00162 IVertexDbSvc* vtxsvc;
00163 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00164 if(vtxsvc->isVertexValid()){
00165 double* dbv = vtxsvc->PrimaryVertex();
00166 double* vv = vtxsvc->SigmaPrimaryVertex();
00167 xorigin.setX(dbv[0]);
00168 xorigin.setY(dbv[1]);
00169 xorigin.setZ(dbv[2]);
00170 }
00171
00173 Vint iCharge_good; iCharge_good.clear();
00174 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00175 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00176 if(!(*itTrk)->isMdcTrackValid()) continue;
00177 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00178
00179 HepVector a = mdcTrk->helix();
00180 HepSymMatrix Ea = mdcTrk->err();
00181 HepPoint3D point0(0.,0.,0.);
00182 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
00183 VFHelix helixip(point0,a,Ea);
00184 helixip.pivot(IP);
00185 HepVector vecipa = helixip.a();
00186 double Rvxy0=fabs(vecipa[0]);
00187 double Rvz0=vecipa[3];
00188 double costheta = cos(mdcTrk->theta());
00189
00190 if(fabs(Rvxy0) >= 1.0) continue;
00191 if(fabs(Rvz0) >= 10.0) continue;
00192 if(fabs(costheta) >= 0.930) continue;
00193 iCharge_good.push_back(i);
00194 }
00195
00197 Vint iGam; iGam.clear();
00198 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00199 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00200 if(!(*itTrk)->isEmcShowerValid()) continue;
00201 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00202
00203
00204
00205 if(abs(runNo)>=11615&&abs(runNo)<=11617&&emcTrk->cellId()==805438481) continue;
00206 if(abs(runNo)>=11615&&abs(runNo)<=11655&&emcTrk->cellId()==805376008) continue;
00207 if(abs(runNo)>=13629&&abs(runNo)<=13669&&emcTrk->cellId()==805374754) continue;
00208 if(abs(runNo)>=21190&&abs(runNo)<=21231&&emcTrk->cellId()==805375262) continue;
00209
00210 int emctdc = emcTrk->time();
00211 if(emctdc>14||emctdc<0) continue;
00212
00213 double eraw = emcTrk->energy();
00214 double theta = emcTrk->theta();
00215 int Module = 0;
00216 if(fabs(cos(theta)) < 0.80 && eraw > 0.025){ Module = 1; }
00217 if(fabs(cos(theta)) > 0.86 && fabs(cos(theta)) < 0.92 && eraw > 0.05) { Module = 2; }
00218 if(Module == 0) continue;
00219 iGam.push_back((*itTrk)->trackId());
00220 }
00222 DTagTool trk;
00223 bool cosmic_ok = trk.cosmicandleptonVeto();
00225
00226 int ntk;
00227 double tagmass,ebeam,tagmode,ksmass,dlte;
00228
00229 Vint tagtrk; tagtrk.clear();
00230 Vint tagGam; tagGam.clear();
00231
00232 HepLorentzVector tagp;
00233
00234 double mass_bc_temp, mass_kf_temp, kf_chi2_temp, mks_temp, mpi0_temp, ptag_temp;
00235 int Charge_candidate_D = 0;
00236 double EGam_max_0 = 0;
00237
00238 for(int i1 = 0; i1 < 2; i1++) {
00239 if(Seperate_Charge == 2 ) { Charge_candidate_D = Charge_default; i1 = 2;}
00240 if(Seperate_Charge == 1 ) { Charge_candidate_D = pow(-1,i1); }
00241 if(Seperate_Charge == 0 ) { Charge_candidate_D = 0; i1 = 2; }
00242
00243 for(int i = 0; i < 15; i++) {
00244 EGam_max_0 = 0;
00245 int mdslct=pow(2.,i);
00246 Sing sing;
00247 sing.Mdset(event,evtRecTrkCol,iCharge_good,iGam,mdslct,Ebeam, PID_flag,Charge_candidate_D);
00248 bool oktag=sing.Getoktg();
00249 if(oktag) {
00250
00251
00252
00253
00254 tagtrk=sing.Gettagtrk1();
00255 tagmode = sing.Gettagmd();
00256
00257 if((abs(tagmode) == 11 || abs(tagmode) == 15 || abs(tagmode) == 19)) {
00258
00259 for(int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++) {
00260 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00261 int itrk = (*itTrk)->trackId();
00262 if(!(*itTrk)->isEmcShowerValid()) continue;
00263 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00264
00265 int emctdc = emcTrk->time();
00266 if(emctdc>14||emctdc<0) continue;
00267
00268 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00269 double dang = 200.;
00270 for(int j = 0; j < tagtrk.size(); j++) {
00271 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + tagtrk[j];
00272 if(!(*jtTrk)->isExtTrackValid()) continue;
00273 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00274 if(extTrk->emcVolumeNumber() == -1) continue;
00275 Hep3Vector extpos = extTrk->emcPosition();
00276 double angd = extpos.angle(emcpos);
00277 if(angd < dang) dang = angd;
00278 }
00279 dang = dang * 180 / (CLHEP::pi);
00280 if(fabs(dang)<20.) continue;
00281
00282 double eraw = emcTrk->energy();
00283 if(eraw > EGam_max_0) EGam_max_0 = eraw;
00284 }
00285 }
00286
00287 m_cosmic_ok = cosmic_ok;
00288 m_EGam_max_0 = EGam_max_0;
00289 m_nGood = iCharge_good.size();
00290 m_nGam = iGam.size();
00291 m_runNo = runNo;
00292 m_event = event;
00293
00294 m_tagmode = sing.Gettagmd();
00295 m_mass_bc = sing.Getmass_bc();
00296 m_delE_tag = sing.GetdelE_tag();
00297
00298 m_tuple1->write();
00299
00300
00301
00302
00303
00304 }
00305
00306 }
00307 }
00308
00309 return StatusCode::SUCCESS;
00310
00311 }
00312
00313
00314
00315
00316
00317
00318 StatusCode SD0Tag::finalize() {
00319 MsgStream log(msgSvc(), name());
00320 cout<<"nD0 = "<<nD0<<endl;
00321 cout<<"n1 = "<<n1<<endl;
00322 cout<<"n2 = "<<n2<<endl;
00323 cout<<"ND0 ==== "<<ND0<<endl;
00324 cout<<"NDp ==== "<<NDp<<endl;
00325 log << MSG::INFO << "in finalize()" << endmsg;
00326 return StatusCode::SUCCESS;
00327 }
00328