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 #include "SD0TagAlg/Kpipi0pi0.h"
00026 #include "SD0TagAlg/SingleBase.h"
00027
00028
00029 Kpipi0pi0::Kpipi0pi0()
00030 {}
00031
00032 Kpipi0pi0::~Kpipi0pi0()
00033 {}
00034
00035
00036 void Kpipi0pi0::MTotal(double event,SmartDataPtr<EvtRecTrackCol> evtRecTrkCol, Vint iGood,Vint
00037 iGam, double Ebeam, int PID_flag, int Charge_candidate_D)
00038 {
00039
00040 int nGood=iGood.size();
00041 int nGam=iGam.size();
00042
00043 iGoodtag.clear();
00044 iGamtag.clear();
00045
00046 double mass_bcgg,delE_tag_temp;
00047 int m_chargetag, m_chargek,m_chargepi;
00048 int ika_temp,ipi_temp, iGam1_temp, iGam2_temp, iGam3_temp, iGam4_temp;
00049 HepLorentzVector pddd;
00050 HepLorentzVector pddd_temp;
00051
00052 int cqtm_temp;
00053 IDataProviderSvc* eventSvc = NULL;
00054 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00055 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc,EventModel::EvtRec::EvtRecEvent);
00056 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
00057
00058 int runNo=eventHeader->runNumber();
00059 int rec=eventHeader->eventNumber();
00060
00061 double xecm=2*Ebeam;
00062
00063 kpipi0pi0md=false;
00064 double tagmode=0;
00065
00066 if((evtRecEvent->totalCharged() < 2||nGam<4)){ return; }
00067
00068 double ecms = xecm;
00069
00070
00071 ISimplePIDSvc* simple_pid;
00072 Gaudi::svcLocator()->service("SimplePIDSvc", simple_pid);
00073
00074 double deltaE_tem = 0.20;
00075 int ncount1 = 0;
00076
00077 Hep3Vector xorigin(0,0,0);
00078 IVertexDbSvc* vtxsvc;
00079 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00080 if(vtxsvc->isVertexValid())
00081 {
00082 double* dbv = vtxsvc->PrimaryVertex();
00083 double* vv = vtxsvc->SigmaPrimaryVertex();
00084 xorigin.setX(dbv[0]);
00085 xorigin.setY(dbv[1]);
00086 xorigin.setZ(dbv[2]);
00087 }
00088
00089 double xv=xorigin.x();
00090 double yv=xorigin.y();
00091 double zv=xorigin.z();
00092
00093 HepPoint3D point0(0.,0.,0.);
00094 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
00096 HepLorentzVector p2gfit1, p2gfit2;
00097 HepLorentzVector p2gg;
00098
00099 HepLorentzVector ptrk1_temp, ptrk2_temp, ptrk3_temp, ptrk4_temp, ptrk5_temp, ptrk6_temp, ptrk7_temp, ptrk8_temp;
00100 for(int i = 0; i < evtRecEvent->totalCharged(); i++) {
00101 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
00102
00103 int ika= (*itTrk)->trackId();
00104
00105 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00106 RecMdcKalTrack* mdcKalTrk1 = (*itTrk)->mdcKalTrack();
00107 RecMdcKalTrack::setPidType(RecMdcKalTrack::kaon);
00109 m_chargek=mdcKalTrk1->charge();
00110 if(Charge_candidate_D != 0) {
00111 if(m_chargek != -Charge_candidate_D) continue;
00112 }
00113 if(Charge_candidate_D == 0) {
00114 if(abs(m_chargek) != 1) continue;
00115 }
00117 HepVector a1 = mdcKalTrk1->getZHelixK();
00118 HepSymMatrix Ea1 = mdcKalTrk1->getZErrorK();
00119 VFHelix helixip3_1(point0,a1,Ea1);
00120 helixip3_1.pivot(IP);
00121 HepVector vecipa1 = helixip3_1.a();
00122
00123 double dr1 = fabs(vecipa1[0]);
00124 double dz1 = fabs(vecipa1[3]);
00125 double costheta1 = cos(mdcKalTrk1->theta());
00126 if ( dr1 >= 1.0) continue;
00127 if ( dz1 >= 10.0) continue;
00128 if ( fabs(costheta1) >= 0.93) continue;
00130 if(PID_flag == 5) {
00131 simple_pid->preparePID(*itTrk);
00132 if(simple_pid->probKaon() < 0.0 || simple_pid->probKaon() < simple_pid->probPion()) continue;
00133 }
00135
00136
00137 WTrackParameter kam(xmass[3],mdcKalTrk1->getZHelixK(),mdcKalTrk1->getZErrorK() );
00138
00139
00140
00141
00142 for(int j = 0; j< evtRecEvent->totalCharged(); j++) {
00143 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + j;
00144
00145 int ipi= (*itTrk)->trackId();
00146 if(ipi==ika) continue;
00147
00148 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00149 RecMdcKalTrack* mdcKalTrk2 = (*itTrk)->mdcKalTrack();
00150 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00152 m_chargepi=mdcKalTrk2->charge();
00153 if((m_chargek + m_chargepi) != 0) continue;
00155 HepVector a2 = mdcKalTrk2->getZHelix();
00156 HepSymMatrix Ea2 = mdcKalTrk2->getZError();
00157 VFHelix helixip3_2(point0,a2,Ea2);
00158 helixip3_2.pivot(IP);
00159 HepVector vecipa2 = helixip3_2.a();
00160
00161 double dr2 = fabs(vecipa2[0]);
00162 double dz2 = fabs(vecipa2[3]);
00163 double costheta2 = cos(mdcKalTrk2->theta());
00164 if ( dr2 >= 1.0) continue;
00165 if ( dz2 >= 10.0) continue;
00166 if ( fabs(costheta2) >= 0.93) continue;
00168 if(PID_flag == 5) {
00169 simple_pid->preparePID(*itTrk);
00170 if(simple_pid->probPion() < 0.0 || simple_pid->probPion() < simple_pid->probKaon()) continue;
00171 }
00173
00174 WTrackParameter pip(xmass[2],mdcKalTrk2->getZHelix(),mdcKalTrk2->getZError() );
00175
00176 for(int m = 0; m < nGam-1; m++) {
00177 if(iGam[m]==-1) continue;
00178 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[m]))->emcShower();
00179 double eraw1 = g1Trk->energy();
00180 double phi1 = g1Trk->phi();
00181 double the1 = g1Trk->theta();
00182 HepLorentzVector ptrkg1,ptrkg10,ptrkg12;
00183 ptrkg1.setPx(eraw1*sin(the1)*cos(phi1));
00184 ptrkg1.setPy(eraw1*sin(the1)*sin(phi1));
00185 ptrkg1.setPz(eraw1*cos(the1));
00186 ptrkg1.setE(eraw1);
00187 ptrkg10 = ptrkg1;
00188 ptrkg12 = ptrkg1.boost(-0.011,0,0);
00189
00190 for(int n = m+1; n < nGam; n++) {
00191 if(iGam[n]==-1) continue;
00192 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[n]))->emcShower();
00193 double eraw2 = g2Trk->energy();
00194 double phi2 = g2Trk->phi();
00195 double the2 = g2Trk->theta();
00196 HepLorentzVector ptrkg2,ptrkg20,ptrkg22;
00197 ptrkg2.setPx(eraw2*sin(the2)*cos(phi2));
00198 ptrkg2.setPy(eraw2*sin(the2)*sin(phi2));
00199 ptrkg2.setPz(eraw2*cos(the2));
00200 ptrkg2.setE(eraw2);
00201 ptrkg20 = ptrkg2;
00202 ptrkg22 = ptrkg2.boost(-0.011,0,0);
00203
00205 HepLorentzVector ptrkpi0;
00206 ptrkpi0 = ptrkg12+ptrkg22;
00207 double m_xmpi0_tem = ptrkpi0.mag();
00208 if(m_xmpi0_tem>0.150||m_xmpi0_tem<0.115) continue;
00210 bool IsEndcap1 = false; bool IsEndcap2 = false;
00211 if(fabs(cos(the1)) > 0.86 && fabs(cos(the1)) < 0.92) IsEndcap1 = true;
00212 if(fabs(cos(the2)) > 0.86 && fabs(cos(the2)) < 0.92) IsEndcap2 = true;
00213 if(IsEndcap1 && IsEndcap2) continue;
00215
00216 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
00217 kmfit->init();
00218 kmfit->setChisqCut(2500);
00219 kmfit->AddTrack(0, 0.0 , g1Trk);
00220 kmfit->AddTrack(1, 0.0 , g2Trk);
00221 kmfit->AddResonance(0, mpi0, 0, 1);
00222 kmfit->Fit(0);
00223 kmfit->BuildVirtualParticle(0);
00224
00225 double pi0_chisq = kmfit->chisq(0);
00226 if ( pi0_chisq >= 2500) continue;
00227 HepLorentzVector p2gfit1 = kmfit->pfit(0) + kmfit->pfit(1);
00228 p2gfit1.boost(-0.011,0,0);
00229
00230 for(int s = 0; s < nGam-1; s++) {
00231 if(iGam[s]==-1) continue;
00232 RecEmcShower *g3Trk = (*(evtRecTrkCol->begin()+iGam[s]))->emcShower();
00233 if(iGam[s] == iGam[m] || iGam[s] == iGam[n] ) continue;
00234 double eraw3 = g3Trk->energy();
00235 double phi3 = g3Trk->phi();
00236 double the3 = g3Trk->theta();
00237 HepLorentzVector ptrkg3,ptrkg30,ptrkg32;
00238 ptrkg3.setPx(eraw3*sin(the3)*cos(phi3));
00239 ptrkg3.setPy(eraw3*sin(the3)*sin(phi3));
00240 ptrkg3.setPz(eraw3*cos(the3));
00241 ptrkg3.setE(eraw3);
00242 ptrkg30 = ptrkg3;
00243 ptrkg32 = ptrkg3.boost(-0.011,0,0);
00244
00245 for(int r = s+1; r < nGam; r++) {
00246 if(iGam[r]==-1) continue;
00247 RecEmcShower *g4Trk = (*(evtRecTrkCol->begin()+iGam[r]))->emcShower();
00248 if(iGam[r] == iGam[m] || iGam[r] == iGam[n] ) continue;
00249 double eraw4 = g4Trk->energy();
00250 double phi4 = g4Trk->phi();
00251 double the4 = g4Trk->theta();
00252 HepLorentzVector ptrkg4,ptrkg40,ptrkg42;
00253 ptrkg4.setPx(eraw4*sin(the4)*cos(phi4));
00254 ptrkg4.setPy(eraw4*sin(the4)*sin(phi4));
00255 ptrkg4.setPz(eraw4*cos(the4));
00256 ptrkg4.setE(eraw4);
00257 ptrkg40 = ptrkg4;
00258 ptrkg42 = ptrkg4.boost(-0.011,0,0);
00259
00261 HepLorentzVector ptrkpi0_2;
00262 ptrkpi0_2 = ptrkg32+ptrkg42;
00263 double m_xmpi0_2_tem = ptrkpi0_2.mag();
00264 if(m_xmpi0_2_tem>0.150||m_xmpi0_2_tem<0.115) continue;
00266 bool IsEndcap3 = false; bool IsEndcap4 = false;
00267 if(fabs(cos(the3)) > 0.86 && fabs(cos(the3)) < 0.92) IsEndcap3 = true;
00268 if(fabs(cos(the4)) > 0.86 && fabs(cos(the4)) < 0.92) IsEndcap4 = true;
00269 if(IsEndcap3 && IsEndcap4) continue;
00271
00272 KalmanKinematicFit * kmfit2 = KalmanKinematicFit::instance();
00273 kmfit2->init();
00274 kmfit2->setChisqCut(2500);
00275 kmfit2->AddTrack(0, 0.0 , g3Trk);
00276 kmfit2->AddTrack(1, 0.0 , g4Trk);
00277 kmfit2->AddResonance(0, mpi0, 0, 1);
00278
00279 kmfit2->Fit(0);
00280 kmfit2->BuildVirtualParticle(0);
00281
00282 double pi0_2_chisq = kmfit2->chisq(0);
00283 if ( pi0_2_chisq >= 2500) continue;
00284 HepLorentzVector p2gfit2 = kmfit2->pfit(0) + kmfit2->pfit(1);
00285 p2gfit2.boost(-0.011,0,0);
00286
00287
00289 HepPoint3D vx(xorigin.x(), xorigin.y(), xorigin.z());
00290 HepSymMatrix Evx(3, 0);
00291 double bx = 1E+6; Evx[0][0] = bx*bx;
00292 double by = 1E+6; Evx[1][1] = by*by;
00293 double bz = 1E+6; Evx[2][2] = bz*bz;
00294 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
00296
00297 VertexFit* vtxfit = VertexFit::instance();
00298 vtxfit->init();
00299 vtxfit->AddTrack(0, kam);
00300 vtxfit->AddTrack(1, pip);
00301 vtxfit->AddVertex(0, vxpar, 0, 1);
00302 if(!vtxfit->Fit(0)) continue;
00303 vtxfit->Swim(0);
00304
00305 WTrackParameter wkam = vtxfit->wtrk(0);
00306 WTrackParameter wpip = vtxfit->wtrk(1);
00307
00308 HepVector kam_val = HepVector(7,0);
00309 kam_val = wkam.w();
00310 HepVector pip_val = HepVector(7,0);
00311 pip_val = wpip.w();
00312
00313 HepLorentzVector P_KAM(kam_val[0],kam_val[1],kam_val[2],kam_val[3]);
00314 HepLorentzVector P_PIP(pip_val[0],pip_val[1],pip_val[2],pip_val[3]);
00315
00316 P_KAM.boost(-0.011,0,0);
00317 P_PIP.boost(-0.011,0,0);
00318 pddd = P_KAM + P_PIP + p2gfit1 + p2gfit2;
00319
00320 double pkpipi0pi0=pddd.rho();
00321
00322 double temp1 = (ecms/2)*(ecms/2)-pkpipi0pi0*pkpipi0pi0 ;
00323 if(temp1<0) temp1 =0;
00324 double mass_bc_tem = sqrt(temp1);
00325 if(mass_bc_tem < 1.82 || mass_bc_tem > 1.89) continue;
00326
00327 double delE_tag_tag = ecms/2-pddd.e();
00328
00329 if(fabs(delE_tag_tag)<deltaE_tem) {
00330 deltaE_tem = fabs(delE_tag_tag);
00331 delE_tag_temp = delE_tag_tag;
00332 mass_bcgg = mass_bc_tem;
00333
00334 pddd_temp = pddd;
00335 cqtm_temp = m_chargek;
00336
00337 ika_temp=ika;
00338 ipi_temp=ipi;
00339
00340 iGam1_temp = iGam[m];
00341 iGam2_temp = iGam[n];
00342 iGam3_temp = iGam[s];
00343 iGam4_temp = iGam[r];
00344 ncount1 = 1;
00345
00346 }
00347 }
00348 }
00349 }
00350 }
00351 }
00352 }
00353
00354 if(ncount1 == 1){
00355 tagmode=25;
00356 if(cqtm_temp<0) tagmode=-25;
00357 tagmd=tagmode;
00358 mass_bc = mass_bcgg;
00359 delE_tag = delE_tag_temp;
00360
00361 cqtm = -1.0*cqtm_temp;
00362
00363 iGoodtag.push_back(ipi_temp);
00364 iGoodtag.push_back(ika_temp);
00365
00366 iGamtag.push_back(iGam1_temp);
00367 iGamtag.push_back(iGam2_temp);
00368 iGamtag.push_back(iGam3_temp);
00369 iGamtag.push_back(iGam4_temp);
00370
00371 ptag = pddd_temp;
00372
00373 kpipi0pi0md = true;
00374
00375 }
00376
00377
00378 }
00379
00380
00381