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
00008 #include "EventModel/EventModel.h"
00009 #include "EventModel/Event.h"
00010
00011 #include "EvtRecEvent/EvtRecEvent.h"
00012 #include "EvtRecEvent/EvtRecTrack.h"
00013 #include "DstEvent/TofHitStatus.h"
00014 #include "EventModel/EventHeader.h"
00015 #include "McTruth/McParticle.h"
00016 #include "HltEvent/HltInf.h"
00017 #include "HltEvent/DstHltInf.h"
00018
00019 #include "TMath.h"
00020 #include "GaudiKernel/INTupleSvc.h"
00021 #include "GaudiKernel/NTuple.h"
00022 #include "GaudiKernel/Bootstrap.h"
00023 #include "GaudiKernel/IHistogramSvc.h"
00024 #include "CLHEP/Vector/ThreeVector.h"
00025 #include "CLHEP/Vector/LorentzVector.h"
00026 #include "CLHEP/Vector/TwoVector.h"
00027 #include "VertexFit/IVertexDbSvc.h"
00028 #include "GaudiKernel/Bootstrap.h"
00029 #include "GaudiKernel/ISvcLocator.h"
00030
00031
00032 using CLHEP::Hep3Vector;
00033 using CLHEP::Hep2Vector;
00034 using CLHEP::HepLorentzVector;
00035 #include "CLHEP/Geometry/Point3D.h"
00036 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00037 typedef HepGeom::Point3D<double> HepPoint3D;
00038 #endif
00039 #include "Gam4pikpAlg/Gam4pikp.h"
00040
00041 #include "VertexFit/KalmanKinematicFit.h"
00042 #include "VertexFit/VertexFit.h"
00043 #include "VertexFit/SecondVertexFit.h"
00044 #include "VertexFit/Helix.h"
00045 #include "ParticleID/ParticleID.h"
00046 #include <vector>
00047 const double mpi = 0.13957;
00048 const double mk = 0.493677;
00049 const double mpro = 0.938272;
00050 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00051 const double velc = 299.792458;
00052 typedef std::vector<int> Vint;
00053 typedef std::vector<HepLorentzVector> Vp4;
00054 typedef std::vector<double> Vdouble;
00055 typedef std::vector<WTrackParameter> Vw;
00056 typedef std::vector<VertexParameter> Vv;
00057 static int Ncut[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
00058 static int Mcut[10]={0,0,0,0,0,0,0,0,0,0};
00060
00061 Gam4pikp::Gam4pikp(const std::string& name, ISvcLocator* pSvcLocator) :
00062 Algorithm(name, pSvcLocator) {
00063
00064
00065 declareProperty("skim4pi", m_skim4pi=false);
00066 declareProperty("skim4k", m_skim4k=false);
00067 declareProperty("skim2pi2pr", m_skim2pi2pr=false);
00068 declareProperty("rootput", m_rootput=false);
00069 declareProperty("Ecms", m_ecms=3.6862);
00070 declareProperty("Vr0cut", m_vr0cut=1.0);
00071 declareProperty("Vz0cut", m_vz0cut=5.0);
00072 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
00073 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00074 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00075 declareProperty("GammaDangCut", m_gammadangcut=20.0);
00076
00077
00078
00079
00080 }
00081
00082
00083 StatusCode Gam4pikp::initialize(){
00084 MsgStream log(msgSvc(), name());
00085
00086 log << MSG::INFO << "in initialize()" << endmsg;
00087
00088 StatusCode status;
00089
00090 if(m_rootput)
00091 {
00092 NTuplePtr nt1(ntupleSvc(), "FILE1/total4c");
00093 if ( nt1 ) m_tuple1 = nt1;
00094 else {
00095 m_tuple1 = ntupleSvc()->book ("FILE1/total4c", CLID_ColumnWiseTuple, "ks N-Tuple example");
00096 if ( m_tuple1 ) {
00097
00098 status = m_tuple1->addItem ("run", m_run );
00099 status = m_tuple1->addItem ("rec", m_rec );
00100 status = m_tuple1->addItem ("mpprecall", m_mpprecall );
00101 status = m_tuple1->addItem ("meeall", m_meeall );
00102 status = m_tuple1->addItem ("ncgjs", m_ncgjs );
00103 status = m_tuple1->addItem ("cla2kpi", m_cla2kpi );
00104 status = m_tuple1->addItem("indexmc", m_idxmc, 0, 100);
00105 status = m_tuple1->addIndexedItem("pdgid", m_idxmc, m_pdgid);
00106
00107 status = m_tuple1->addIndexedItem("motheridx", m_idxmc, m_motheridx);
00108 status = m_tuple1->addItem("indexmdc", m_idxmdc, 0, 5000);
00109 status = m_tuple1->addIndexedItem ("x0js", m_idxmdc, m_x0js);
00110 status = m_tuple1->addIndexedItem ("y0js", m_idxmdc, m_y0js);
00111 status = m_tuple1->addIndexedItem ("z0js",m_idxmdc, m_z0js);
00112 status = m_tuple1->addIndexedItem ("r0js",m_idxmdc, m_r0js);
00113 status = m_tuple1->addIndexedItem ("Rxyjs",m_idxmdc, m_Rxyjs);
00114 status = m_tuple1->addIndexedItem ("Rzjs",m_idxmdc, m_Rzjs);
00115 status = m_tuple1->addIndexedItem ("Rnxyjs",m_idxmdc, m_Rnxyjs);
00116 status = m_tuple1->addIndexedItem ("phinjs",m_idxmdc, m_phinjs);
00117 status = m_tuple1->addIndexedItem ("Rnzjs",m_idxmdc, m_Rnzjs);
00118 status = m_tuple1->addItem ("ncy20", m_ncy20);
00119 status = m_tuple1->addItem ("ncy30", m_ncy30);
00120 status = m_tuple1->addIndexedItem("angjs5", m_idxmdc, m_angjs5);
00121 status = m_tuple1->addIndexedItem("nearjs5", m_idxmdc, m_nearjs5);
00122 status = m_tuple1->addIndexedItem("angjs6", m_idxmdc, m_angjs6);
00123 status = m_tuple1->addIndexedItem("nearjs6", m_idxmdc, m_nearjs6);
00124 status = m_tuple1->addIndexedItem("ang4pi5", m_idxmdc, m_ang4pi5);
00125 status = m_tuple1->addIndexedItem("near4pi5", m_idxmdc, m_near4pi5);
00126 status = m_tuple1->addIndexedItem("ang4pi6", m_idxmdc, m_ang4pi6);
00127 status = m_tuple1->addIndexedItem("near4pi6", m_idxmdc, m_near4pi6);
00128 status = m_tuple1->addIndexedItem("ppmdcjs", m_idxmdc, m_ppmdcjs);
00129 status = m_tuple1->addIndexedItem("pxmdcjs", m_idxmdc, m_pxmdcjs);
00130 status = m_tuple1->addIndexedItem("pymdcjs", m_idxmdc, m_pymdcjs);
00131 status = m_tuple1->addIndexedItem("pzmdcjs", m_idxmdc, m_pzmdcjs);
00132 status = m_tuple1->addIndexedItem("ppkaljs", m_idxmdc, m_ppkaljs);
00133 status = m_tuple1->addIndexedItem("ptmdcjs", m_idxmdc, m_ptmdcjs);
00134 status = m_tuple1->addIndexedItem("ptkaljs", m_idxmdc, m_ptkaljs);
00135 status = m_tuple1->addIndexedItem("ppmdc2kpi", m_idxmdc, m_ppmdc2kpi);
00136 status = m_tuple1->addIndexedItem("pxmdc2kpi", m_idxmdc, m_pxmdc2kpi);
00137 status = m_tuple1->addIndexedItem("pymdc2kpi", m_idxmdc, m_pymdc2kpi);
00138 status = m_tuple1->addIndexedItem("pzmdc2kpi", m_idxmdc, m_pzmdc2kpi);
00139 status = m_tuple1->addIndexedItem("ppkal2kpi", m_idxmdc, m_ppkal2kpi);
00140 status = m_tuple1->addIndexedItem("ptmdc2kpi", m_idxmdc, m_ptmdc2kpi);
00141 status = m_tuple1->addIndexedItem("charge2kpi", m_idxmdc, m_charge2kpi);
00142 status = m_tuple1->addIndexedItem("ptkal2kpi", m_idxmdc, m_ptkal2kpi);
00143 status = m_tuple1->addItem ("cy2pi", m_cy2kpi, 0, 100 );
00144 status = m_tuple1->addIndexedItem("comcs2kpi", m_cy2kpi, m_comcs2kpi);
00145 status = m_tuple1->addItem ("chiejs", m_idxmdc, m_chiejs);
00146 status = m_tuple1->addItem ("chimujs", m_idxmdc, m_chimujs);
00147 status = m_tuple1->addItem ("chipijs", m_idxmdc, m_chipijs);
00148 status = m_tuple1->addItem ("chikjs", m_idxmdc, m_chikjs);
00149 status = m_tuple1->addItem ("chipjs", m_idxmdc, m_chipjs);
00150 status = m_tuple1->addItem ("ghitjs", m_idxmdc, m_ghitjs);
00151 status = m_tuple1->addItem ("thitjs", m_idxmdc, m_thitjs);
00152 status = m_tuple1->addIndexedItem("probphjs", m_idxmdc, m_probphjs);
00153 status = m_tuple1->addIndexedItem("normphjs", m_idxmdc, m_normphjs);
00154 status = m_tuple1->addItem ("pdg", m_idxmdc, m_pdg);
00155 status = m_tuple1->addItem ("cbmc", m_idxmdc, m_cbmc);
00156 status = m_tuple1->addIndexedItem("sigmaetof2kpi", m_idxmdc, m_sigmaetof2kpi);
00157 status = m_tuple1->addIndexedItem("sigmamutof2kpi", m_idxmdc, m_sigmamutof2kpi);
00158 status = m_tuple1->addIndexedItem("sigmapitof2kpi", m_idxmdc, m_sigmapitof2kpi);
00159 status = m_tuple1->addIndexedItem("sigmaktof2kpi", m_idxmdc, m_sigmaktof2kpi);
00160 status = m_tuple1->addIndexedItem("sigmaprtof2kpi", m_idxmdc, m_sigmaprtof2kpi);
00161 status = m_tuple1->addIndexedItem("t0tof2kpi", m_idxmdc, m_t0tof2kpi);
00162 status = m_tuple1->addIndexedItem("errt0tof2kpi", m_idxmdc, m_errt0tof2kpi);
00163
00164 status = m_tuple1->addItem ("chie2kpi", m_idxmdc, m_chie2kpi);
00165 status = m_tuple1->addItem ("chimu2kpi", m_idxmdc, m_chimu2kpi);
00166 status = m_tuple1->addItem ("chipi2kpi", m_idxmdc, m_chipi2kpi);
00167 status = m_tuple1->addItem ("chik2kpi", m_idxmdc, m_chik2kpi);
00168 status = m_tuple1->addItem ("chip2kpi", m_idxmdc, m_chip2kpi);
00169 status = m_tuple1->addItem ("ghit2kpi", m_idxmdc, m_ghit2kpi);
00170 status = m_tuple1->addItem ("thit2kpi", m_idxmdc, m_thit2kpi);
00171 status = m_tuple1->addIndexedItem("probph2kpi", m_idxmdc, m_probph2kpi);
00172 status = m_tuple1->addIndexedItem("normph2kpi", m_idxmdc, m_normph2kpi);
00173 status = m_tuple1->addIndexedItem("pidnum2kpi", m_idxmdc, m_pidnum2kpi);
00174 status = m_tuple1->addIndexedItem("bjmucjs", m_idxmdc, m_bjmucjs);
00175 status = m_tuple1->addIndexedItem("bjmuc2kpi", m_idxmdc, m_bjmuc2kpi);
00176 status = m_tuple1->addIndexedItem("bjemcjs", m_idxmdc, m_bjemcjs);
00177 status = m_tuple1->addIndexedItem("bjemc2kpi", m_idxmdc, m_bjemc2kpi);
00178 status = m_tuple1->addIndexedItem("bjtofjs", m_idxmdc, m_bjtofjs);
00179 status = m_tuple1->addIndexedItem("bjtof2kpi", m_idxmdc, m_bjtof2kpi);
00180 status = m_tuple1->addIndexedItem("bjtofvaljs", m_idxmdc, m_bjtofvaljs);
00181 status = m_tuple1->addIndexedItem("bjtofval2kpi", m_idxmdc, m_bjtofval2kpi);
00182
00183 status = m_tuple1->addIndexedItem("emcjs", m_idxmdc, m_emcjs);
00184 status = m_tuple1->addIndexedItem("evpjs", m_idxmdc, m_evpjs);
00185 status = m_tuple1->addIndexedItem("timecgjs", m_idxmdc, m_timecgjs);
00186 status = m_tuple1->addIndexedItem("depthjs", m_idxmdc, m_depthmucjs);
00187 status = m_tuple1->addIndexedItem("layermucjs", m_idxmdc, m_layermucjs);
00188
00189 status = m_tuple1->addIndexedItem("emc2kpi", m_idxmdc, m_emc2kpi);
00190 status = m_tuple1->addIndexedItem("evp2kpi", m_idxmdc, m_evp2kpi);
00191 status = m_tuple1->addIndexedItem("timecg2kpi", m_idxmdc, m_timecg2kpi);
00192 status = m_tuple1->addIndexedItem("depth2kpi", m_idxmdc, m_depthmuc2kpi);
00193 status = m_tuple1->addIndexedItem("layermuc2kpi", m_idxmdc, m_layermuc2kpi);
00194
00195 status = m_tuple1->addIndexedItem("cotof1js", m_idxmdc, m_cotof1js);
00196 status = m_tuple1->addIndexedItem("cotof2js", m_idxmdc, m_cotof2js);
00197 status = m_tuple1->addIndexedItem("counterjs", m_idxmdc, m_counterjs);
00198 status = m_tuple1->addIndexedItem("barreljs", m_idxmdc, m_barreljs);
00199 status = m_tuple1->addIndexedItem("layertofjs", m_idxmdc, m_layertofjs);
00200 status = m_tuple1->addIndexedItem("readoutjs", m_idxmdc, m_readoutjs);
00201 status = m_tuple1->addIndexedItem("clusterjs", m_idxmdc, m_clusterjs);
00202 status = m_tuple1->addIndexedItem("betajs", m_idxmdc, m_betajs);
00203 status = m_tuple1->addIndexedItem("tofjs", m_idxmdc, m_tofjs);
00204 status = m_tuple1->addIndexedItem("tofpathjs", m_idxmdc, m_tofpathjs);
00205 status = m_tuple1->addIndexedItem("zhitjs", m_idxmdc, m_zhitjs);
00206 status = m_tuple1->addIndexedItem("tofIDjs", m_idxmdc, m_tofIDjs);
00207 status = m_tuple1->addIndexedItem("clusterIDjs", m_idxmdc, m_clusterIDjs);
00208 status = m_tuple1->addIndexedItem("texejs", m_idxmdc, m_texejs);
00209 status = m_tuple1->addIndexedItem("texmujs", m_idxmdc, m_texmujs);
00210 status = m_tuple1->addIndexedItem("texpijs", m_idxmdc, m_texpijs);
00211 status = m_tuple1->addIndexedItem("texkjs", m_idxmdc, m_texkjs);
00212 status = m_tuple1->addIndexedItem("texprjs", m_idxmdc, m_texprjs);
00213 status = m_tuple1->addIndexedItem("dtejs", m_idxmdc, m_dtejs);
00214 status = m_tuple1->addIndexedItem("dtmujs", m_idxmdc, m_dtmujs);
00215 status = m_tuple1->addIndexedItem("dtpijs", m_idxmdc, m_dtpijs);
00216 status = m_tuple1->addIndexedItem("dtkjs", m_idxmdc, m_dtkjs);
00217 status = m_tuple1->addIndexedItem("dtprjs", m_idxmdc, m_dtprjs);
00218 status = m_tuple1->addIndexedItem("sigmaetofjs", m_idxmdc, m_sigmaetofjs);
00219 status = m_tuple1->addIndexedItem("sigmamutofjs", m_idxmdc, m_sigmamutofjs);
00220 status = m_tuple1->addIndexedItem("sigmapitofjs", m_idxmdc, m_sigmapitofjs);
00221 status = m_tuple1->addIndexedItem("sigmaktofjs", m_idxmdc, m_sigmaktofjs);
00222 status = m_tuple1->addIndexedItem("sigmaprtofjs", m_idxmdc, m_sigmaprtofjs);
00223 status = m_tuple1->addIndexedItem("t0tofjs", m_idxmdc,m_t0tofjs);
00224 status = m_tuple1->addIndexedItem("errt0tofjs", m_idxmdc,m_errt0tofjs);
00225 status = m_tuple1->addIndexedItem("cotof12kpi", m_idxmdc, m_cotof12kpi);
00226 status = m_tuple1->addIndexedItem("cotof22kpi", m_idxmdc, m_cotof22kpi);
00227 status = m_tuple1->addIndexedItem("counter2kpi", m_idxmdc, m_counter2kpi);
00228 status = m_tuple1->addIndexedItem("barrel2kpi", m_idxmdc, m_barrel2kpi);
00229 status = m_tuple1->addIndexedItem("layertof2kpi", m_idxmdc, m_layertof2kpi);
00230 status = m_tuple1->addIndexedItem("readout2kpi", m_idxmdc, m_readout2kpi);
00231 status = m_tuple1->addIndexedItem("cluster2kpi", m_idxmdc, m_cluster2kpi);
00232 status = m_tuple1->addIndexedItem("beta2kpi", m_idxmdc, m_beta2kpi);
00233 status = m_tuple1->addIndexedItem("tof2kpi", m_idxmdc, m_tof2kpi);
00234 status = m_tuple1->addIndexedItem("tofpath2kpi", m_idxmdc, m_tofpath2kpi);
00235 status = m_tuple1->addIndexedItem("zhit2kpi", m_idxmdc, m_zhit2kpi);
00236 status = m_tuple1->addIndexedItem("tofID2kpi", m_idxmdc, m_tofID2kpi);
00237 status = m_tuple1->addIndexedItem("clusterID2kpi", m_idxmdc, m_clusterID2kpi);
00238 status = m_tuple1->addIndexedItem("texe2kpi", m_idxmdc, m_texe2kpi);
00239 status = m_tuple1->addIndexedItem("texmu2kpi", m_idxmdc, m_texmu2kpi);
00240 status = m_tuple1->addIndexedItem("texpi2kpi", m_idxmdc, m_texpi2kpi);
00241 status = m_tuple1->addIndexedItem("texk2kpi", m_idxmdc, m_texk2kpi);
00242 status = m_tuple1->addIndexedItem("texpr2kpi", m_idxmdc, m_texpr2kpi);
00243 status = m_tuple1->addIndexedItem("dte2kpi", m_idxmdc, m_dte2kpi);
00244 status = m_tuple1->addIndexedItem("dtmu2kpi", m_idxmdc, m_dtmu2kpi);
00245 status = m_tuple1->addIndexedItem("dtpi2kpi", m_idxmdc, m_dtpi2kpi);
00246 status = m_tuple1->addIndexedItem("dtk2kpi", m_idxmdc, m_dtk2kpi);
00247 status = m_tuple1->addIndexedItem("dtpr2kpi", m_idxmdc, m_dtpr2kpi);
00248 status = m_tuple1->addIndexedItem("costpid2kpi", m_idxmdc, m_costpid2kpi);
00249 status = m_tuple1->addIndexedItem("dedxpid2kpi", m_idxmdc, m_dedxpid2kpi);
00250 status = m_tuple1->addIndexedItem("tof1pid2kpi", m_idxmdc, m_tof1pid2kpi);
00251 status = m_tuple1->addIndexedItem("tof2pid2kpi", m_idxmdc, m_tof2pid2kpi);
00252 status = m_tuple1->addIndexedItem("probe2kpi", m_idxmdc, m_probe2kpi);
00253 status = m_tuple1->addIndexedItem("probmu2kpi", m_idxmdc, m_probmu2kpi);
00254 status = m_tuple1->addIndexedItem("probpi2kpi", m_idxmdc, m_probpi2kpi);
00255 status = m_tuple1->addIndexedItem("probk2kpi", m_idxmdc, m_probk2kpi);
00256 status = m_tuple1->addIndexedItem("probpr2kpi", m_idxmdc, m_probpr2kpi);
00257
00258 status = m_tuple1->addIndexedItem("chipidxpid2kpi", m_idxmdc, m_chipidxpid2kpi);
00259 status = m_tuple1->addIndexedItem("chipitof1pid2kpi", m_idxmdc, m_chipitof1pid2kpi);
00260 status = m_tuple1->addIndexedItem("chipitof2pid2kpi", m_idxmdc, m_chipitof2pid2kpi);
00261 status = m_tuple1->addIndexedItem("chipitofpid2kpi", m_idxmdc, m_chipitofpid2kpi);
00262 status = m_tuple1->addIndexedItem("chipitofepid2kpi", m_idxmdc, m_chipitofepid2kpi);
00263 status = m_tuple1->addIndexedItem("chipitofqpid2kpi", m_idxmdc, m_chipitofqpid2kpi);
00264 status = m_tuple1->addIndexedItem("probpidxpid2kpi", m_idxmdc, m_probpidxpid2kpi);
00265 status = m_tuple1->addIndexedItem("probpitofpid2kpi", m_idxmdc, m_probpitofpid2kpi);
00266 status = m_tuple1->addIndexedItem("chikdxpid2kpi", m_idxmdc, m_chikdxpid2kpi);
00267 status = m_tuple1->addIndexedItem("chiktof1pid2kpi", m_idxmdc, m_chiktof1pid2kpi);
00268 status = m_tuple1->addIndexedItem("chiktof2pid2kpi", m_idxmdc, m_chiktof2pid2kpi);
00269 status = m_tuple1->addIndexedItem("chiktofpid2kpi", m_idxmdc, m_chiktofpid2kpi);
00270 status = m_tuple1->addIndexedItem("chiktofepid2kpi", m_idxmdc, m_chiktofepid2kpi);
00271 status = m_tuple1->addIndexedItem("chiktofqpid2kpi", m_idxmdc, m_chiktofqpid2kpi);
00272 status = m_tuple1->addIndexedItem("probkdxpid2kpi", m_idxmdc, m_probkdxpid2kpi);
00273 status = m_tuple1->addIndexedItem("probktofpid2kpi", m_idxmdc, m_probktofpid2kpi);
00274
00275 status = m_tuple1->addIndexedItem("chiprdxpid2kpi", m_idxmdc, m_chiprdxpid2kpi);
00276 status = m_tuple1->addIndexedItem("chiprtof1pid2kpi", m_idxmdc, m_chiprtof1pid2kpi);
00277 status = m_tuple1->addIndexedItem("chiprtof2pid2kpi", m_idxmdc, m_chiprtof2pid2kpi);
00278 status = m_tuple1->addIndexedItem("chiprtofpid2kpi", m_idxmdc, m_chiprtofpid2kpi);
00279 status = m_tuple1->addIndexedItem("chiprtofepid2kpi", m_idxmdc, m_chiprtofepid2kpi);
00280 status = m_tuple1->addIndexedItem("chiprtofqpid2kpi", m_idxmdc, m_chiprtofqpid2kpi);
00281 status = m_tuple1->addIndexedItem("probprdxpid2kpi", m_idxmdc, m_probprdxpid2kpi);
00282 status = m_tuple1->addIndexedItem("probprtofpid2kpi", m_idxmdc, m_probprtofpid2kpi);
00283
00284 status = m_tuple1->addIndexedItem("cosmdcjs", m_idxmdc, m_cosmdcjs);
00285 status = m_tuple1->addIndexedItem("phimdcjs", m_idxmdc, m_phimdcjs);
00286 status = m_tuple1->addIndexedItem("cosmdc2kpi", m_idxmdc, m_cosmdc2kpi);
00287 status = m_tuple1->addIndexedItem("phimdc2kpi", m_idxmdc, m_phimdc2kpi);
00288
00289 status = m_tuple1->addIndexedItem("dedxpidjs", m_idxmdc, m_dedxpidjs);
00290 status = m_tuple1->addIndexedItem("tof1pidjs", m_idxmdc, m_tof1pidjs);
00291 status = m_tuple1->addIndexedItem("tof2pidjs", m_idxmdc, m_tof2pidjs);
00292 status = m_tuple1->addIndexedItem("probejs", m_idxmdc, m_probejs);
00293 status = m_tuple1->addIndexedItem("probmujs", m_idxmdc, m_probmujs);
00294 status = m_tuple1->addIndexedItem("probpijs", m_idxmdc, m_probpijs);
00295 status = m_tuple1->addIndexedItem("probkjs", m_idxmdc, m_probkjs);
00296 status = m_tuple1->addIndexedItem("probprjs", m_idxmdc, m_probprjs);
00297 status = m_tuple1->addItem ("mchic2kpi", m_mchic2kpi);
00298 status = m_tuple1->addItem ("mpsip2kpi", m_mpsip2kpi);
00299 status = m_tuple1->addItem ("chis2kpi", m_chis2kpi);
00300 status = m_tuple1->addItem ("mchic4c2kpi", m_mchic4c2kpi);
00301 status = m_tuple1->addItem ("mpsip4c2kpi", m_mpsip4c2kpi);
00302 status = m_tuple1->addItem ("chis4c2kpi", m_chis4c2kpi);
00303
00304 status = m_tuple1->addItem("indexemc", m_idxemc, 0, 5000);
00305 status = m_tuple1->addIndexedItem("numHits", m_idxemc, m_numHits);
00306 status = m_tuple1->addIndexedItem("secmom", m_idxemc, m_secondmoment);
00307 status = m_tuple1->addIndexedItem("latmom", m_idxemc, m_latmoment);
00308 status = m_tuple1->addIndexedItem("timegm", m_idxemc, m_timegm);
00309 status = m_tuple1->addIndexedItem("cellId", m_idxemc, m_cellId);
00310 status = m_tuple1->addIndexedItem("module", m_idxemc, m_module);
00311 status = m_tuple1->addIndexedItem("a20Moment", m_idxemc, m_a20Moment);
00312 status = m_tuple1->addIndexedItem("a42Moment", m_idxemc, m_a42Moment);
00313 status = m_tuple1->addIndexedItem("getEAll", m_idxemc, m_getEAll);
00314 status = m_tuple1->addIndexedItem("getShowerId", m_idxemc, m_getShowerId);
00315 status = m_tuple1->addIndexedItem("getClusterId", m_idxemc, m_getClusterId);
00316 status = m_tuple1->addIndexedItem("x", m_idxemc, m_x);
00317 status = m_tuple1->addIndexedItem("y", m_idxemc, m_y);
00318 status = m_tuple1->addIndexedItem("z", m_idxemc, m_z);
00319 status = m_tuple1->addIndexedItem("cosemc", m_idxemc, m_cosemc);
00320 status = m_tuple1->addIndexedItem("phiemc", m_idxemc, m_phiemc);
00321 status = m_tuple1->addIndexedItem("energy", m_idxemc, m_energy);
00322 status = m_tuple1->addIndexedItem("e1", m_idxemc, m_eSeed);
00323 status = m_tuple1->addIndexedItem("e9", m_idxemc, m_e3x3);
00324 status = m_tuple1->addIndexedItem("e25", m_idxemc, m_e5x5);
00325 status = m_tuple1->addIndexedItem("dang4c", m_idxemc, m_dang4c);
00326 status = m_tuple1->addIndexedItem("dthe4c", m_idxemc, m_dthe4c);
00327 status = m_tuple1->addIndexedItem("dphi4c", m_idxemc, m_dphi4c);
00328 status = m_tuple1->addIndexedItem("dang4crt", m_idxemc, m_dang4crt);
00329 status = m_tuple1->addIndexedItem("dthe4crt", m_idxemc, m_dthe4crt);
00330 status = m_tuple1->addIndexedItem("dphi4crt", m_idxemc, m_dphi4crt);
00331 status = m_tuple1->addIndexedItem("phtof", m_idxemc, 3, m_phgmtof,0.0,10000.0);
00332 status = m_tuple1->addIndexedItem("phgmtof0", m_idxemc, m_phgmtof0);
00333 status = m_tuple1->addIndexedItem("phgmtof1", m_idxemc, m_phgmtof1);
00334 status = m_tuple1->addIndexedItem("phgmtof2", m_idxemc, m_phgmtof2);
00335
00336 }
00337 else {
00338 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00339 return StatusCode::FAILURE;
00340 }
00341 }
00342 }
00343
00344 StatusCode sc;
00345
00346 if(m_skim4pi)
00347 {
00348 sc = createSubAlgorithm( "EventWriter", "Selectgam4pi", m_subalg1);
00349 if( sc.isFailure() ) {
00350 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4pi" <<endreq;
00351 return sc;
00352 } else {
00353 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4pi" <<endreq;
00354 }
00355 }
00356
00357
00358 if(m_skim4k)
00359 {
00360 sc = createSubAlgorithm( "EventWriter", "Selectgam4k", m_subalg2);
00361 if( sc.isFailure() ) {
00362 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4k" <<endreq;
00363 return sc;
00364 } else {
00365 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4k" <<endreq;
00366 }
00367 }
00368
00369 if(m_skim2pi2pr)
00370 {
00371 sc = createSubAlgorithm( "EventWriter", "Selectgam2pi2pr", m_subalg3);
00372 if( sc.isFailure() ) {
00373 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
00374 return sc;
00375 } else {
00376 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
00377 }
00378 }
00379
00380
00381
00382
00383
00384
00385
00386
00387 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00388 return StatusCode::SUCCESS;
00389
00390 }
00391
00392
00393 StatusCode Gam4pikp::execute() {
00394
00395
00396 setFilterPassed(false);
00397
00398
00399 MsgStream log(msgSvc(), name());
00400 log << MSG::INFO << "in execute()" << endreq;
00401
00402 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00403 int run=eventHeader->runNumber();
00404 int event=eventHeader->eventNumber();
00405 log << MSG::DEBUG <<"run, evtnum = "
00406 << run << " , "
00407 << event <<endreq;
00408
00409
00410 if(m_rootput)
00411 {
00412 m_run=run;
00413 m_rec=event;
00414 }
00415 Ncut[0]++;
00416 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00417 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00418 << evtRecEvent->totalCharged() << " , "
00419 << evtRecEvent->totalNeutral() << " , "
00420 << evtRecEvent->totalTracks() <<endreq;
00421
00422 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00423
00424 if(m_rootput)
00425 {
00426 Gam4pikp::InitVar();
00427 }
00428
00429
00430 if(m_rootput)
00431 {
00432 if (eventHeader->runNumber()<0)
00433 {
00434 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
00435 int m_numParticle = 0;
00436 if (!mcParticleCol)
00437 {
00438
00439 return StatusCode::FAILURE;
00440 }
00441 else
00442 {
00443 bool psipDecay = false;
00444 int rootIndex = -1;
00445 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
00446 for (; iter_mc != mcParticleCol->end(); iter_mc++)
00447 {
00448 if ((*iter_mc)->primaryParticle()) continue;
00449 if (!(*iter_mc)->decayFromGenerator()) continue;
00450 if ((*iter_mc)->particleProperty()==100443)
00451 {
00452 psipDecay = true;
00453 rootIndex = (*iter_mc)->trackIndex();
00454 }
00455 if (!psipDecay) continue;
00456 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
00457 int pdgid = (*iter_mc)->particleProperty();
00458 m_pdgid[m_numParticle] = pdgid;
00459 m_motheridx[m_numParticle] = mcidx;
00460 m_numParticle += 1;
00461 }
00462 }
00463 m_idxmc = m_numParticle;
00464 }
00465 }
00466
00467 Vint iGood, ipip, ipim;
00468 iGood.clear();
00469 ipip.clear();
00470 ipim.clear();
00471 Vp4 ppip, ppim;
00472 ppip.clear();
00473 ppim.clear();
00474 int nCharge = 0;
00475 int ngdcgx=0;
00476 int ncgx=0;
00477 Vdouble pTrkCh;
00478 pTrkCh.clear();
00479 Hep3Vector v(0,0,0);
00480 Hep3Vector vv(0,0,0);
00481
00482 IVertexDbSvc* vtxsvc;
00483 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00484 if(vtxsvc->isVertexValid())
00485 {
00486 double* vertex = vtxsvc->PrimaryVertex();
00487 double* vertexsigma = vtxsvc->SigmaPrimaryVertex();
00488 v.setX(vertex[0]);
00489 v.setY(vertex[1]);
00490 v.setZ(vertex[2]);
00491 vv.setX(vertexsigma[0]);
00492 vv.setY(vertexsigma[1]);
00493 vv.setZ(vertexsigma[2]);
00494 }
00495
00496
00497 if(run<0)
00498 {
00499 v[0]=0.;
00500 v[1]=0.;
00501 v[2]=0.;
00502 vv[0]=0.;
00503 vv[1]=0.;
00504 vv[2]=0.;
00505
00506 }
00507
00508
00509 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00510 if(i>=evtRecTrkCol->size()) break;
00511 ncgx++;
00512 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00513 if(!(*itTrk)->isMdcTrackValid()) continue;
00514 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00515 ngdcgx++;
00516 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00517
00518 HepVector a = mdcTrk->helix();
00519 HepSymMatrix Ea = mdcTrk->err();
00520 HepPoint3D pivot(0.,0.,0.);
00521 HepPoint3D IP(v[0],v[1],v[2]);
00522 VFHelix helixp(pivot,a,Ea);
00523 helixp.pivot(IP);
00524 HepVector vec = helixp.a();
00525 pTrkCh.push_back(mdcTrk->p()*mdcTrk->charge());
00526 iGood.push_back((*itTrk)->trackId());
00527 nCharge += mdcTrk->charge();
00528 }
00529
00530
00531 int nGood = iGood.size();
00532 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00533 if((nGood < 4))
00534 {
00535 return StatusCode::SUCCESS;
00536 }
00537
00538 Ncut[1]++;
00539 Vint iGam;
00540 int ngamx=0;
00541 int ngdgamx=0;
00542 iGam.clear();
00543 Vint iGamnofit;
00544 iGamnofit.clear();
00545 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00546 ngamx++;
00547 if(i>=evtRecTrkCol->size()) break;
00548 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00549 if(!(*itTrk)->isEmcShowerValid()) continue;
00550 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00551 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00552
00553
00554
00555
00556 double dthe = 200.;
00557 double dphi = 200.;
00558 double dang = 200.;
00559
00560 ngdgamx++;
00561 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
00562 if(j>=evtRecTrkCol->size()) break;
00563 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
00564 if(!(*jtTrk)->isExtTrackValid()) continue;
00565 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00566 if(extTrk->emcVolumeNumber() == -1) continue;
00567 Hep3Vector extpos = extTrk->emcPosition();
00568
00569 double angd = extpos.angle(emcpos);
00570 double thed = extpos.theta() - emcpos.theta();
00571 double phid = extpos.deltaPhi(emcpos);
00572 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00573 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00574
00575
00576
00577
00578
00579 if(angd < dang){ dang=angd;
00580 dthe=thed;
00581 dphi=phid;
00582 }
00583
00584 }
00585
00586 double eraw = emcTrk->energy();
00587 dthe = dthe * 180 / (CLHEP::pi);
00588 dphi = dphi * 180 / (CLHEP::pi);
00589 dang = dang * 180 / (CLHEP::pi);
00590
00591
00592
00593
00594 iGam.push_back((*itTrk)->trackId());
00595 if(eraw < m_energyThreshold) continue;
00596 if(dang < 20.0) continue;
00597 iGamnofit.push_back((*itTrk)->trackId());
00598
00599 }
00600
00601
00602 int nGam = iGam.size();
00603 int nGamnofit = iGamnofit.size();
00604
00605 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
00606 if(nGam<1){
00607 return StatusCode::SUCCESS;
00608 }
00609 Ncut[2]++;
00610
00611 if(nGood>20||nGam>30) return StatusCode::SUCCESS;
00612
00613 if(m_rootput)
00614 {
00615 m_idxmdc = nGood;
00616 m_idxemc=nGam;
00617 }
00618
00619 Vp4 pGam;
00620 pGam.clear();
00621 for(int i = 0; i < nGam; i++) {
00622 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00623 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00624 double eraw = emcTrk->energy();
00625 double phi = emcTrk->phi();
00626 double the = emcTrk->theta();
00627 HepLorentzVector ptrk;
00628 ptrk.setPx(eraw*sin(the)*cos(phi));
00629 ptrk.setPy(eraw*sin(the)*sin(phi));
00630 ptrk.setPz(eraw*cos(the));
00631 ptrk.setE(eraw);
00632
00633
00634
00635 pGam.push_back(ptrk);
00636 }
00637
00638
00639
00640 ParticleID *pid = ParticleID::instance();
00641
00642
00643 for(int i = 0; i < nGood; i++) {
00644 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00645 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00646
00647 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00648 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
00649
00650 if(mdcKalTrk->charge() >0) {
00651 ipip.push_back(iGood[i]);
00652 HepLorentzVector ptrk;
00653 ptrk.setPx(mdcKalTrk->px());
00654 ptrk.setPy(mdcKalTrk->py());
00655 ptrk.setPz(mdcKalTrk->pz());
00656 double p3 = ptrk.mag();
00657 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00658
00659
00660 ppip.push_back(ptrk);
00661 } else {
00662 ipim.push_back(iGood[i]);
00663 HepLorentzVector ptrk;
00664 ptrk.setPx(mdcKalTrk->px());
00665 ptrk.setPy(mdcKalTrk->py());
00666 ptrk.setPz(mdcKalTrk->pz());
00667 double p3 = ptrk.mag();
00668 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00669
00670
00671
00672 ppim.push_back(ptrk);
00673 }
00674 }
00675 int npip = ipip.size();
00676 int npim = ipim.size();
00677
00678 if((npip < 2)||(npim < 2)) return SUCCESS;
00679
00680 Ncut[3]++;
00681
00682
00683
00684
00685
00686 int icgp1=-1;
00687 int icgp2=-1;
00688 int icgm1=-1;
00689 int icgm2=-1;
00690 int igam=-1;
00691 double chisqtrk = 9999.;
00692 WTrackParameter wcgp1;
00693 WTrackParameter wcgp2;
00694 WTrackParameter wcgm1;
00695 WTrackParameter wcgm2;
00696 int ipip1js=-1;
00697 int ipip2js=-1;
00698 int ipim1js=-1;
00699 int ipim2js=-1;
00700 double mcompall=9999;
00701 double mppreclst=9999;
00702 double meelst=9999;;
00703 double mchic2kpilst=9999;
00704 double chis4c2kpilst=9999;
00705 int type2kpilst=0;
00706 double dtpr2kpilst[4]={9999,9999,9999,9999};
00707
00708 double mc1=mpi,mc2=mpi,mc3=mpi,mc4=mpi;
00709 double mchic01=0.0;
00710 HepLorentzVector pgam=(0,0,0,0);
00711 HepPoint3D xorigin(0.,0.,0.);
00712 xorigin[0]=9999.;
00713 xorigin[1]=9999.;
00714 xorigin[2]=9999.;
00715 HepSymMatrix xem(3,0);
00716
00717 int cl4pi=0;
00718 int clajs=0;
00719 HepLorentzVector p4psipjs(0.011*m_ecms, 0.0, 0.0, m_ecms);
00720 double psipBetajs = (p4psipjs.vect()).mag()/(p4psipjs.e());
00721 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.0, m_ecms);
00722 double psipBeta = (p4psip.vect()).mag()/(p4psip.e());
00723
00724 for(int ii = 0; ii < npip-1; ii++) {
00725 RecMdcKalTrack *pip1Trk = (*(evtRecTrkCol->begin()+ipip[ii]))->mdcKalTrack();
00726 for(int ij = ii+1; ij < npip; ij++) {
00727 RecMdcKalTrack *pip2Trk = (*(evtRecTrkCol->begin()+ipip[ij]))->mdcKalTrack();
00728 for(int ik = 0; ik < npim-1; ik++) {
00729 RecMdcKalTrack *pim1Trk = (*(evtRecTrkCol->begin()+ipim[ik]))->mdcKalTrack();
00730 for(int il = ik+1; il < npim; il++) {
00731 RecMdcKalTrack *pim2Trk = (*(evtRecTrkCol->begin()+ipim[il]))->mdcKalTrack();
00732 double squar[3]={9999.,9999.,9999.};
00733 double squarkpi[6]={9999.,9999.,9999.,9999.,9999.,9999.};
00734 WTrackParameter wvpip1Trk, wvpim1Trk, wvpip2Trk, wvpim2Trk;
00735
00736 Vint iGoodjs;
00737 iGoodjs.clear();
00738 Vdouble pTrkjs;
00739 pTrkjs.clear();
00740
00741
00742
00743 pTrkjs.push_back(pip1Trk->p()*pip1Trk->charge());
00744 pTrkjs.push_back(pip2Trk->p()*pip2Trk->charge());
00745 pTrkjs.push_back(pim1Trk->p()*pim1Trk->charge());
00746 pTrkjs.push_back(pim2Trk->p()*pim2Trk->charge());
00747 iGoodjs.push_back(ipip[ii]);
00748 iGoodjs.push_back(ipip[ij]);
00749 iGoodjs.push_back(ipim[ik]);
00750 iGoodjs.push_back(ipim[il]);
00751
00752 Gam4pikp::BubbleSort(pTrkjs, iGoodjs);
00753 Vint jGoodjs;
00754 jGoodjs.clear();
00755 Vp4 p4chTrkjs;
00756 p4chTrkjs.clear();
00757 Vint i4cpip1js, i4cpip2js, i4cpim1js, i4cpim2js;
00758 i4cpip1js.clear();
00759 i4cpip2js.clear();
00760 i4cpim1js.clear();
00761 i4cpim2js.clear();
00762 i4cpip1js.push_back(iGoodjs[2]);
00763 i4cpip2js.push_back(iGoodjs[3]);
00764 i4cpim1js.push_back(iGoodjs[1]);
00765 i4cpim2js.push_back(iGoodjs[0]);
00766 jGoodjs.push_back(i4cpip1js[0]);
00767 jGoodjs.push_back(i4cpim1js[0]);
00768 jGoodjs.push_back(i4cpip2js[0]);
00769 jGoodjs.push_back(i4cpim2js[0]);
00770
00771 for (int i = 0; i < 4; i++)
00772 {
00773 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGoodjs[i];
00774 if (!(*itTrk)->isMdcTrackValid()) continue;
00775 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00776 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00777 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00778 double ptrk;
00779 HepLorentzVector p4trk;
00780 if (i<2)
00781 {
00782 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00783 ptrk = mdcKalTrk->p();
00784 p4trk.setPx(mdcKalTrk->px());
00785 p4trk.setPy(mdcKalTrk->py());
00786 p4trk.setPz(mdcKalTrk->pz());
00787 p4trk.setE(sqrt(ptrk*ptrk+mpi*mpi));
00788 }
00789 else{
00790 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
00791 ptrk = mdcKalTrk->p();
00792 p4trk.setPx(mdcKalTrk->px());
00793 p4trk.setPy(mdcKalTrk->py());
00794 p4trk.setPz(mdcKalTrk->pz());
00795 p4trk.setE(sqrt(ptrk*ptrk+xmass[0]*xmass[0]));
00796 }
00797 p4trk.boost(-1.0*psipBetajs, 0.0, 0.0);
00798 p4chTrkjs.push_back(p4trk);
00799 }
00800 p4psipjs.boost(-1.0*psipBetajs, 0.0, 0.0);
00801
00802 HepLorentzVector p4pipijs = p4chTrkjs[0] + p4chTrkjs[1];
00803 HepLorentzVector p4eejs = p4chTrkjs[2] + p4chTrkjs[3];
00804 HepLorentzVector p4pipirecjs = p4psipjs - p4pipijs;
00805 double mpprecjs=p4pipirecjs.m();
00806 double mpipijs=p4pipijs.m();
00807 double meejs=p4eejs.m();
00808 double mcomp=sqrt((mpprecjs-3.097)*(mpprecjs-3.097)+(meejs-3.097)*(meejs-3.097));
00809 if(mcomp<mcompall)
00810 {
00811 mcompall=mcomp;
00812 ipip1js=i4cpip1js[0];
00813 ipim1js=i4cpim1js[0];
00814 ipip2js=i4cpip2js[0];
00815 ipim2js=i4cpim2js[0];
00816 mppreclst=mpprecjs;
00817 meelst=meejs;
00818 }
00819
00820 if(m_rootput)
00821 {
00822 m_mpprecall=mppreclst;
00823 m_meeall=meelst;
00824 }
00825 }
00826
00827 }
00828 }
00829 }
00830
00831
00832
00833
00834 Vint iGood4c;
00835 iGood4c.clear();
00836 Vdouble pTrk4c;
00837 pTrk4c.clear();
00838 Vint jGood;
00839 jGood.clear();
00840 Vint jsGood;
00841 jsGood.clear();
00842
00843 if(mcompall>9997)
00844 return StatusCode::SUCCESS;
00845
00846 jsGood.push_back(ipip1js);
00847 jsGood.push_back(ipim1js);
00848 jsGood.push_back(ipip2js);
00849 jsGood.push_back(ipim2js);
00850
00851 for(int i = 0; i < evtRecEvent->totalCharged(); i++)
00852 {
00853 if(i>=evtRecTrkCol->size()) break;
00854 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00855 if(!(*itTrk)->isMdcTrackValid()) continue;
00856 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00857 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00858 if((i!=ipip1js)&&(i!=ipim1js)&&(i!=ipip2js)&&(i!=ipim2js))
00859 {
00860 jsGood.push_back(i);
00861 }
00862 }
00863
00864 int njsGood=jsGood.size();
00865
00866 if(m_rootput)
00867 {
00868 for (int i = 0; i < njsGood; i++)
00869 {
00870 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jsGood[i];
00871 if (!(*itTrk)->isMdcTrackValid()) continue;
00872 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00873 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00874 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00875 double ptrk;
00876 ptrk = mdcKalTrk->p();
00877 m_x0js[i] = mdcTrk->x();
00878 m_y0js[i] = mdcTrk->y();
00879 m_z0js[i] = mdcTrk->z();
00880 m_r0js[i] = mdcTrk->r();
00881 m_ppmdcjs[i] = mdcTrk->p();
00882 m_pxmdcjs[i] = mdcTrk->px();
00883 m_pymdcjs[i] = mdcTrk->py();
00884 m_pzmdcjs[i] = mdcTrk->pz();
00885 m_ppkaljs[i] = mdcKalTrk->p();
00886 Hep3Vector p3jsi(mdcTrk->px(),mdcTrk->py(),mdcTrk->pz());
00887 if(njsGood>4){
00888 EvtRecTrackIterator itTrk5 = evtRecTrkCol->begin() + jsGood[4];
00889 RecMdcTrack *mdcTrk5 = (*itTrk5)->mdcTrack();
00890 Hep3Vector p3js5(mdcTrk5->px(),mdcTrk5->py(),mdcTrk5->pz());
00891 m_angjs5[i]=p3jsi.angle(p3js5);
00892 m_nearjs5[i]=p3jsi.howNear(p3js5);
00893 }
00894 if(njsGood>5){
00895 EvtRecTrackIterator itTrk6 = evtRecTrkCol->begin() + jsGood[5];
00896 RecMdcTrack *mdcTrk6 = (*itTrk6)->mdcTrack();
00897 Hep3Vector p3js6(mdcTrk6->px(),mdcTrk6->py(),mdcTrk6->pz());
00898 m_angjs6[i]=p3jsi.angle(p3js6);
00899 m_nearjs6[i]=p3jsi.howNear(p3js6);
00900 }
00901
00902
00903 m_ptmdcjs[i] = mdcTrk->pxy();
00904 m_ptkaljs[i] = mdcKalTrk->pxy();
00905 double x0=mdcTrk->x();
00906 double y0=mdcTrk->y();
00907 double z0=mdcTrk->z();
00908 double phi0=mdcTrk->helix(1);
00909 double xv=v[0];
00910 double yv=v[1];
00911 double zv=v[2];
00912 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00913 double Rz=z0-zv;
00914 m_Rxyjs[i]=Rxy;
00915 m_Rzjs[i]=Rz;
00916 HepVector a = mdcTrk->helix();
00917 HepSymMatrix Ea = mdcTrk->err();
00918 HepPoint3D pivot(0.,0.,0.);
00919 HepPoint3D IP(v[0],v[1],v[2]);
00920 VFHelix helixp(pivot,a,Ea);
00921 helixp.pivot(IP);
00922 HepVector vec = helixp.a();
00923 m_Rnxyjs[i]=vec[0];
00924 m_Rnzjs[i]=vec[3];
00925 m_phinjs[i]=vec[1];
00926
00927 if ((*itTrk)->isTofTrackValid())
00928 { m_bjtofvaljs[i]=1;
00929 m_cotof1js[i]=1;
00930 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00931 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
00932 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
00933 TofHitStatus *status = new TofHitStatus;
00934 status->setStatus((*iter_tof)->status());
00935 if(status->is_cluster()){
00936 m_bjtofjs[i]=1;
00937 m_counterjs[i] = status->is_counter();
00938 m_barreljs[i] = status->is_barrel();
00939 m_layertofjs[i] = status->layer();
00940 m_readoutjs[i] = status->is_readout();
00941 m_clusterjs[i] = status->is_cluster();
00942 m_cotof2js[i]=2;
00943 m_betajs[i]= (*iter_tof)->beta();
00944 m_tofjs[i] = (*iter_tof)->tof();
00945 m_tofpathjs[i] = (*iter_tof)->path();
00946 m_zhitjs[i]= (*iter_tof)->zrhit();
00947 m_texejs[i] = (*iter_tof)->texpElectron();
00948 m_texmujs[i] = (*iter_tof)->texpMuon();
00949 m_texpijs[i] = (*iter_tof)->texpPion();
00950 m_texkjs[i] = (*iter_tof)->texpKaon();
00951 m_texprjs[i] = (*iter_tof)->texpProton();
00952 m_dtejs[i] = m_tofjs[i] - m_texejs[i];
00953 m_dtmujs[i] = m_tofjs[i] - m_texmujs[i];
00954 m_dtpijs[i] = m_tofjs[i] - m_texpijs[i];
00955 m_dtkjs[i] = m_tofjs[i] - m_texkjs[i];
00956 m_dtprjs[i] = m_tofjs[i] - m_texprjs[i];
00957
00958 m_sigmaetofjs[i]=(*iter_tof)->sigma(0);
00959 m_sigmamutofjs[i]=(*iter_tof)->sigma(1);
00960 m_sigmapitofjs[i]=(*iter_tof)->sigma(2);
00961 m_sigmaktofjs[i]=(*iter_tof)->sigma(3);
00962 m_sigmaprtofjs[i]=(*iter_tof)->sigma(4);
00963 m_t0tofjs[i]=(*iter_tof)->t0();
00964 m_errt0tofjs[i]=(*iter_tof)->errt0();
00965
00966 m_tofIDjs[i]=(*iter_tof)->tofID();
00967 m_clusterIDjs[i]=(*iter_tof)->tofTrackID();
00968 }
00969 delete status;
00970 }
00971 }
00972
00973 if ((*itTrk)->isMdcDedxValid())
00974 {
00975
00976 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00977 m_chiejs[i] = dedxTrk->chiE();
00978 m_chimujs[i] = dedxTrk->chiMu();
00979 m_chipijs[i] = dedxTrk->chiPi();
00980 m_chikjs[i] = dedxTrk->chiK();
00981 m_chipjs[i] = dedxTrk->chiP();
00982 m_ghitjs[i] = dedxTrk->numGoodHits();
00983 m_thitjs[i] = dedxTrk->numTotalHits();
00984 m_probphjs[i] = dedxTrk->probPH();
00985 m_normphjs[i] = dedxTrk->normPH();
00986 }
00987
00988
00989
00990
00991 if ( (*itTrk)->isEmcShowerValid() )
00992 {
00993 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00994 m_bjemcjs[i]=1;
00995 m_emcjs[i] = emcTrk->energy();
00996 m_evpjs[i] = emcTrk->energy()/ptrk;
00997 m_timecgjs[i] = emcTrk->time();
00998
00999 }
01000
01001
01002 if ( (*itTrk)->isMucTrackValid() )
01003 {
01004 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
01005 double dpthp = mucTrk->depth()/25.0;
01006
01007 m_bjmucjs[i]=1;
01008 m_depthmucjs[i]=mucTrk->depth();
01009 m_layermucjs[i] = mucTrk->numLayers();
01010 }
01011
01012 m_cosmdcjs[i] = cos(mdcTrk->theta());
01013 m_phimdcjs[i] = mdcTrk->phi();
01014
01015 pid->init();
01016 pid->setMethod(pid->methodProbability());
01017 pid->setChiMinCut(4);
01018 pid->setRecTrack(*itTrk);
01019 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
01020 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
01021 pid->identify(pid->onlyElectron() | pid->onlyMuon());
01022 pid->calculate();
01023 if(!(pid->IsPidInfoValid())) continue;
01024
01025 m_dedxpidjs[i] = pid->chiDedx(2);
01026 m_tof1pidjs[i] = pid->chiTof1(2);
01027 m_tof2pidjs[i] = pid->chiTof2(2);
01028 m_probejs[i] = pid->probElectron();
01029 m_probmujs[i] = pid->probMuon();
01030 m_probpijs[i] = pid->probPion();
01031 m_probkjs[i] = pid->probKaon();
01032 m_probprjs[i] = pid->probProton();
01033
01034
01035
01036 }
01037
01038 }
01039
01040 Vint jGam2kpi;
01041 jGam2kpi.clear();
01042
01043 Vint iGood2kpi, ipip2kpi, ipim2kpi;
01044 iGood2kpi.clear();
01045 ipip2kpi.clear();
01046 ipim2kpi.clear();
01047 Vp4 ppip2kpi, ppim2kpi;
01048 ppip2kpi.clear();
01049 ppim2kpi.clear();
01050
01051 Vint ipipnofit, ipimnofit, ikpnofit, ikmnofit, ipropnofit, ipromnofit;
01052 ipipnofit.clear();
01053 ipimnofit.clear();
01054 ikpnofit.clear();
01055 ikmnofit.clear();
01056 ipropnofit.clear();
01057 ipromnofit.clear();
01058 Vp4 ppipnofit, ppimnofit, pkpnofit, pkmnofit, ppropnofit, ppromnofit;
01059 ppipnofit.clear();
01060 ppimnofit.clear();
01061 pkpnofit.clear();
01062 pkmnofit.clear();
01063 ppropnofit.clear();
01064 ppromnofit.clear();
01065
01066 Vdouble p3pip2kpi;
01067 p3pip2kpi.clear();
01068 Vdouble p3pim2kpi;
01069 p3pim2kpi.clear();
01070
01071 Vint itrak2kpi;
01072 itrak2kpi.clear();
01073 int cls2kpi;
01074 double chis2kpi=9999.;
01075 double m4piall=9999.;
01076 double m4kall=9999.;
01077 double mgam4piall=9999.;
01078 double mgam4kall=9999.;
01079 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
01080 if(i>=evtRecTrkCol->size()) break;
01081 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
01082 if(!(*itTrk)->isMdcTrackValid()) continue;
01083 if (!(*itTrk)->isMdcKalTrackValid()) continue;
01084 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
01085 double z02kpi = mdcTrk->z();
01086 double r02kpi = mdcTrk->r();
01087 HepVector a = mdcTrk->helix();
01088 HepSymMatrix Ea = mdcTrk->err();
01089 HepPoint3D pivot(0.,0.,0.);
01090 HepPoint3D IP(v[0],v[1],v[2]);
01091 VFHelix helixp(pivot,a,Ea);
01092 helixp.pivot(IP);
01093 HepVector vec = helixp.a();
01094 double Rnxy02kpi=vec[0];
01095 double Rnz02kpi=vec[3];
01096
01097
01098 iGood2kpi.push_back((*itTrk)->trackId());
01099 }
01100 int nGood2kpi = iGood2kpi.size();
01101 if(nGood2kpi==4)
01102 {
01103 for(int i = 0; i < nGood2kpi; i++) {
01104 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood2kpi[i];
01105 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
01106 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
01107 if(mdcKalTrk->charge() >0) {
01108 ipip2kpi.push_back(iGood2kpi[i]);
01109 p3pip2kpi.push_back(mdcKalTrk->p());
01110 HepLorentzVector ptrk;
01111 ptrk.setPx(mdcKalTrk->px());
01112 ptrk.setPy(mdcKalTrk->py());
01113 ptrk.setPz(mdcKalTrk->pz());
01114 double p3 = ptrk.mag();
01115 ptrk.setE(sqrt(p3*p3+mpi*mpi));
01116 ppip2kpi.push_back(ptrk);
01117 } else {
01118 ipim2kpi.push_back(iGood2kpi[i]);
01119 p3pim2kpi.push_back(mdcKalTrk->p());
01120 HepLorentzVector ptrk;
01121 ptrk.setPx(mdcKalTrk->px());
01122 ptrk.setPy(mdcKalTrk->py());
01123 ptrk.setPz(mdcKalTrk->pz());
01124 double p3 = ptrk.mag();
01125 ptrk.setE(sqrt(p3*p3+mpi*mpi));
01126 ppim2kpi.push_back(ptrk);
01127 }
01128 }
01129 int npip2kpi = ipip2kpi.size();
01130 int npim2kpi = ipim2kpi.size();
01131
01132
01133
01134 ParticleID *pid = ParticleID::instance();
01135
01136 if(m_rootput)
01137 {
01138 m_cy2kpi=6;
01139 }
01140
01141 if((npip2kpi == 2)&&(npim2kpi == 2))
01142 {
01143
01144
01145
01146 Ncut[4]++;
01147
01148 HepPoint3D xorigin2kpi(0.,0.,0.);
01149 xorigin2kpi[0]=9999.;
01150 xorigin2kpi[1]=9999.;
01151 xorigin2kpi[2]=9999.;
01152 HepSymMatrix xem2kpi(3,0);
01153
01154 Gam4pikp::BubbleSort(p3pip2kpi, ipip2kpi);
01155 Gam4pikp::BubbleSort(p3pim2kpi, ipim2kpi);
01156
01157 Mcut[1]++;
01158 RecMdcKalTrack *pip12kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[0]))->mdcKalTrack();
01159 RecMdcKalTrack *pip22kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[1]))->mdcKalTrack();
01160 RecMdcKalTrack *pim12kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[0]))->mdcKalTrack();
01161 RecMdcKalTrack *pim22kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[1]))->mdcKalTrack();
01162 double squar2kpi[10]={9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
01163 double mc12kpi,mc22kpi,mc32kpi,mc42kpi;
01164
01165 WTrackParameter wcgp12kpi;
01166 WTrackParameter wcgp22kpi;
01167 WTrackParameter wcgm12kpi;
01168 WTrackParameter wcgm22kpi;
01169 int icgp12kpi=-1;
01170 int icgp22kpi=-1;
01171 int icgm12kpi=-1;
01172 int icgm22kpi=-1;
01173 int igam2kpi=-1;
01174 double m2kpi=9999;
01175
01176 int n20=0;
01177 int n30=0;
01178 WTrackParameter wvpip12kpiTrk, wvpim12kpiTrk, wvpip22kpiTrk, wvpim22kpiTrk;
01179 for(int k=0;k<6;k++)
01180 {
01181 if(k==0){mc12kpi=mpi;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpi;}
01182 if(k==1){mc12kpi=mk;mc22kpi=mk;mc32kpi=mk;mc42kpi=mk;}
01183 if(k==2){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpi;mc42kpi=mpro;}
01184 if(k==3){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpro;mc42kpi=mpi;}
01185 if(k==4){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpro;mc42kpi=mpi; }
01186 if(k==5){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpro;}
01187
01188
01189
01190 wvpip12kpiTrk = WTrackParameter(mc12kpi, pip12kpiTrk->getZHelix(), pip12kpiTrk->getZError());
01191 wvpip22kpiTrk = WTrackParameter(mc22kpi, pip22kpiTrk->getZHelix(), pip22kpiTrk->getZError());
01192 wvpim12kpiTrk = WTrackParameter(mc32kpi, pim12kpiTrk->getZHelix(), pim12kpiTrk->getZError());
01193 wvpim22kpiTrk = WTrackParameter(mc42kpi, pim22kpiTrk->getZHelix(), pim22kpiTrk->getZError());
01194 HepPoint3D vx(0., 0., 0.);
01195 HepSymMatrix Evx(3, 0);
01196 double bx = 1E+6;
01197 double by = 1E+6;
01198 double bz = 1E+6;
01199 Evx[0][0] = bx*bx;
01200 Evx[1][1] = by*by;
01201 Evx[2][2] = bz*bz;
01202 VertexParameter vxpar;
01203 vxpar.setVx(vx);
01204 vxpar.setEvx(Evx);
01205 VertexFit* vtxfit = VertexFit::instance();
01206 vtxfit->init();
01207 vtxfit->AddTrack(0, wvpip12kpiTrk);
01208 vtxfit->AddTrack(1, wvpim12kpiTrk);
01209 vtxfit->AddTrack(2, wvpip22kpiTrk);
01210 vtxfit->AddTrack(3, wvpim22kpiTrk);
01211 vtxfit->AddVertex(0, vxpar,0, 1, 2, 3);
01212 if(!vtxfit->Fit(0)) continue;
01213 vtxfit->Swim(0);
01214 WTrackParameter wpip12kpi = vtxfit->wtrk(0);
01215 WTrackParameter wpim12kpi = vtxfit->wtrk(1);
01216 WTrackParameter wpip22kpi = vtxfit->wtrk(2);
01217 WTrackParameter wpim22kpi = vtxfit->wtrk(3);
01218 WTrackParameter wpip12kpiys = vtxfit->wtrk(0);
01219 WTrackParameter wpim12kpiys = vtxfit->wtrk(1);
01220 WTrackParameter wpip22kpiys = vtxfit->wtrk(2);
01221 WTrackParameter wpim22kpiys = vtxfit->wtrk(3);
01222 xorigin2kpi = vtxfit->vx(0);
01223 xem2kpi = vtxfit->Evx(0);
01224 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
01225
01226 HepLorentzVector ecms(0.040547,0,0,3.68632);
01227
01228 double chisq = 9999.;
01229 int ig1 = -1;
01230 for(int i = 0; i < nGam; i++) {
01231 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
01232 kmfit->init();
01233 kmfit->setBeamPosition(xorigin2kpi);
01234 kmfit->setVBeamPosition(xem2kpi);
01235 kmfit->AddTrack(0, wpip12kpi);
01236 kmfit->AddTrack(1, wpim12kpi);
01237 kmfit->AddTrack(2, wpip22kpi);
01238 kmfit->AddTrack(3, wpim22kpi);
01239 kmfit->AddTrack(4, 0.0, g1Trk);
01240 kmfit->AddFourMomentum(0, p4psip);
01241
01242 bool oksq = kmfit->Fit();
01243 if(oksq) {
01244 double chi2 = kmfit->chisq();
01245 if(chi2 < chisq) {
01246 chisq = chi2;
01247 squar2kpi[k]=chi2;
01248 ig1 = iGam[i];
01249
01250 }
01251 }
01252 }
01253 if(squar2kpi[k]<200&&squar2kpi[k]<chis2kpi)
01254 {
01255
01256 chis2kpi=squar2kpi[k];
01257 if(squar2kpi[k]<20) n20=n20+1;
01258 if(squar2kpi[k]<30) n30=n30+1;
01259
01260 icgp12kpi=ipip2kpi[0];
01261 icgp22kpi=ipip2kpi[1];
01262 icgm12kpi=ipim2kpi[0];
01263 icgm22kpi=ipim2kpi[1];
01264 igam2kpi=ig1;
01265 wcgp12kpi=wpip12kpiys;
01266 wcgp22kpi=wpip22kpiys;
01267 wcgm12kpi=wpim12kpiys;
01268 wcgm22kpi=wpim22kpiys;
01269 cls2kpi=k;
01270
01271 if(k==0){
01272 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
01273 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
01274 }
01275
01276 if(k==1){
01277 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
01278 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
01279 }
01280
01281
01282 if(k==2){
01283 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
01284 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
01285 }
01286
01287 if(k==3){
01288 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
01289 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
01290 }
01291
01292 if(k==4){
01293 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
01294 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
01295 }
01296
01297 if(k==5){
01298 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
01299 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
01300 }
01301
01302
01303
01304
01305
01306 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+igam2kpi))->emcShower();
01307 kmfit->init();
01308 kmfit->setBeamPosition(xorigin2kpi);
01309 kmfit->setVBeamPosition(xem2kpi);
01310 kmfit->AddTrack(0, wpip12kpi);
01311 kmfit->AddTrack(1, wpim12kpi);
01312 kmfit->AddTrack(2, wpip22kpi);
01313 kmfit->AddTrack(3, wpim22kpi);
01314 kmfit->AddTrack(4, 0.0, g1Trk);
01315 kmfit->AddFourMomentum(0, p4psip);
01316 bool oksq = kmfit->Fit();
01317 if(oksq){
01318 HepLorentzVector pchic2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3);
01319 HepLorentzVector ppsip2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3) + kmfit->pfit(4);
01320 mchic2kpilst=pchic2kpi.m();
01321 chis4c2kpilst=kmfit->chisq();
01322 if(m_rootput)
01323 {
01324 m_mchic2kpi = pchic2kpi.m();
01325 m_chis2kpi = kmfit->chisq();
01326 m_mpsip2kpi=ppsip2kpi.m();
01327 }
01328
01329
01330
01331 }
01332
01333
01334 }
01335 }
01336
01337
01338 if(chis2kpi<200)
01339 {
01340 Ncut[5]++;
01341 if(m_rootput)
01342 {
01343 m_ncy20=n20;
01344 m_ncy30=n30;
01345 m_cla2kpi=cls2kpi;
01346 }
01347 type2kpilst=cls2kpi;
01348
01349 Vp4 p2kpi;
01350 p2kpi.clear();
01351
01352 for (int i = 0; i < 4; i++)
01353 {
01354 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
01355 if (!(*itTrk)->isMdcTrackValid()) continue;
01356 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
01357 if (!(*itTrk)->isMdcKalTrackValid()) continue;
01358 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
01359 double ptrk2kpi;
01360 HepLorentzVector p4trk2kpi;
01361 if(cls2kpi==1)
01362 {
01363 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
01364 ptrk2kpi = mdcKalTrk->p();
01365 p4trk2kpi.setPx(mdcKalTrk->px());
01366 p4trk2kpi.setPy(mdcKalTrk->py());
01367 p4trk2kpi.setPz(mdcKalTrk->pz());
01368 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mk*mk));
01369 p2kpi.push_back(p4trk2kpi);
01370 }
01371
01372 if(cls2kpi==2)
01373 {
01374 if (i<2)
01375 {
01376 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
01377 ptrk2kpi = mdcKalTrk->p();
01378 p4trk2kpi.setPx(mdcKalTrk->px());
01379 p4trk2kpi.setPy(mdcKalTrk->py());
01380 p4trk2kpi.setPz(mdcKalTrk->pz());
01381 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
01382 p2kpi.push_back(p4trk2kpi);
01383
01384 }
01385 else{
01386 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
01387 ptrk2kpi = mdcKalTrk->p();
01388 p4trk2kpi.setPx(mdcKalTrk->px());
01389 p4trk2kpi.setPy(mdcKalTrk->py());
01390 p4trk2kpi.setPz(mdcKalTrk->pz());
01391 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpro*mpro));
01392 p2kpi.push_back(p4trk2kpi);
01393
01394 }
01395
01396 }
01397 if(cls2kpi!=1&&cls2kpi!=2)
01398 {
01399 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
01400 ptrk2kpi = mdcKalTrk->p();
01401 p4trk2kpi.setPx(mdcKalTrk->px());
01402 p4trk2kpi.setPy(mdcKalTrk->py());
01403 p4trk2kpi.setPz(mdcKalTrk->pz());
01404 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
01405 p2kpi.push_back(p4trk2kpi);
01406 }
01407
01408
01409
01410
01411 }
01412
01413 for (int i = 0; i < 4; i++)
01414 {
01415 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
01416 if (!(*itTrk)->isMdcTrackValid()) continue;
01417 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
01418 if (!(*itTrk)->isMdcKalTrackValid()) continue;
01419 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
01420 if ((*itTrk)->isTofTrackValid())
01421 {
01422 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
01423 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
01424 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
01425 TofHitStatus *status = new TofHitStatus;
01426 status->setStatus((*iter_tof)->status());
01427 if(status->is_cluster()){ dtpr2kpilst[i]=((*iter_tof)->tof()-(*iter_tof)->texpProton());
01428 }
01429 delete status;
01430 }
01431 }
01432 }
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466 Mcut[2]++;
01467 if(m_rootput)
01468 {
01469 for (int i = 0; i < 4; i++)
01470 {
01471 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
01472 if (!(*itTrk)->isMdcTrackValid()) continue;
01473 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
01474 if (!(*itTrk)->isMdcKalTrackValid()) continue;
01475 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
01476 m_ppmdc2kpi[i]=mdcTrk->p();
01477 m_pxmdc2kpi[i]=mdcTrk->px();
01478 m_pymdc2kpi[i]=mdcTrk->py();
01479 m_pzmdc2kpi[i]=mdcTrk->pz();
01480 m_ppkal2kpi[i]=mdcKalTrk->p();
01481 m_charge2kpi[i]=mdcKalTrk->charge();
01482 double ptrk;
01483 ptrk=mdcKalTrk->p();
01484
01485 if (eventHeader->runNumber()<0)
01486 {double mcall=9999;
01487 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
01488 int m_numParticle = 0;
01489 if (!mcParticleCol)
01490 {
01491
01492 return StatusCode::FAILURE;
01493 }
01494 else
01495 {
01496 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
01497 for (; iter_mc != mcParticleCol->end(); iter_mc++)
01498 { double commc;
01499 int pdgcode = (*iter_mc)->particleProperty();
01500 float px,py,pz;
01501 px=(*iter_mc)->initialFourMomentum().x();
01502 py=(*iter_mc)->initialFourMomentum().y();
01503 pz=(*iter_mc)->initialFourMomentum().z();
01504
01505 commc=ptrk*ptrk-px*px-py*py-pz*pz;
01506 if(fabs(commc)<fabs(mcall))
01507 {
01508 mcall=commc;
01509 m_pdg[i]=pdgcode;
01510 m_cbmc[i]=commc;
01511 }
01512
01513 }
01514
01515 }
01516
01517 }
01518
01519 if ((*itTrk)->isTofTrackValid())
01520 { m_bjtofval2kpi[i]=1;
01521 m_cotof12kpi[i]=1;
01522 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
01523 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
01524 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
01525 TofHitStatus *status = new TofHitStatus;
01526 status->setStatus((*iter_tof)->status());
01527
01528 if(status->is_cluster()){
01529
01530
01531 m_bjtof2kpi[i]=1;
01532 m_counter2kpi[i] = status->is_counter();
01533 m_barrel2kpi[i] = status->is_barrel();
01534 m_layertof2kpi[i] = status->layer();
01535 m_readout2kpi[i] = status->is_readout();
01536 m_cluster2kpi[i] = status->is_cluster();
01537 m_cotof22kpi[i]=2;
01538 m_beta2kpi[i]= (*iter_tof)->beta();
01539 m_tof2kpi[i] = (*iter_tof)->tof();
01540 m_tofpath2kpi[i] = (*iter_tof)->path();
01541 m_zhit2kpi[i]= (*iter_tof)->zrhit();
01542 m_texe2kpi[i] = (*iter_tof)->texpElectron();
01543 m_texmu2kpi[i] = (*iter_tof)->texpMuon();
01544 m_texpi2kpi[i] = (*iter_tof)->texpPion();
01545 m_texk2kpi[i] = (*iter_tof)->texpKaon();
01546 m_texpr2kpi[i] = (*iter_tof)->texpProton();
01547 m_dte2kpi[i] = m_tof2kpi[i] - m_texe2kpi[i];
01548 m_dtmu2kpi[i] = m_tof2kpi[i] - m_texmu2kpi[i];
01549 m_dtpi2kpi[i] = m_tof2kpi[i] - m_texpi2kpi[i];
01550 m_dtk2kpi[i] = m_tof2kpi[i] - m_texk2kpi[i];
01551 m_dtpr2kpi[i] = m_tof2kpi[i] - m_texpr2kpi[i];
01552 m_tofID2kpi[i]=(*iter_tof)->tofID();
01553 m_clusterID2kpi[i]=(*iter_tof)->tofTrackID();
01554 m_sigmaetof2kpi[i]=(*iter_tof)->sigma(0);
01555 m_sigmamutof2kpi[i]=(*iter_tof)->sigma(1);
01556 m_sigmapitof2kpi[i]=(*iter_tof)->sigma(2);
01557 m_sigmaktof2kpi[i]=(*iter_tof)->sigma(3);
01558 m_sigmaprtof2kpi[i]=(*iter_tof)->sigma(4);
01559 m_t0tof2kpi[i]=(*iter_tof)->t0();
01560 m_errt0tof2kpi[i]=(*iter_tof)->errt0();
01561
01562 }
01563 delete status;
01564 }
01565 }
01566
01567
01568
01569
01570
01571
01572 if ((*itTrk)->isMdcDedxValid())
01573 {
01574 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
01575 m_chie2kpi[i] = dedxTrk->chiE();
01576 m_chimu2kpi[i] = dedxTrk->chiMu();
01577 m_chipi2kpi[i] = dedxTrk->chiPi();
01578 m_chik2kpi[i] = dedxTrk->chiK();
01579 m_chip2kpi[i] = dedxTrk->chiP();
01580 m_ghit2kpi[i] = dedxTrk->numGoodHits();
01581 m_thit2kpi[i] = dedxTrk->numTotalHits();
01582 m_probph2kpi[i] = dedxTrk->probPH();
01583 m_normph2kpi[i] = dedxTrk->normPH();
01584 }
01585
01586
01587 if ( (*itTrk)->isEmcShowerValid() )
01588 {
01589 RecEmcShower* emcTrk = (*itTrk)->emcShower();
01590 m_bjemc2kpi[i]=1;
01591 m_emc2kpi[i] = emcTrk->energy();
01592 m_evp2kpi[i] = emcTrk->energy()/ptrk;
01593 m_timecg2kpi[i] = emcTrk->time();
01594 }
01595
01596
01597
01598 if ( (*itTrk)->isMucTrackValid() )
01599 {
01600 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
01601 double dpthp = mucTrk->depth()/25.0;
01602 m_bjmuc2kpi[i]=1;
01603 m_depthmuc2kpi[i]=mucTrk->depth();
01604 m_layermuc2kpi[i] = mucTrk->numLayers();
01605 }
01606
01607 m_cosmdc2kpi[i] = cos(mdcTrk->theta());
01608 m_phimdc2kpi[i] = mdcTrk->phi();
01609
01610 m_pidnum2kpi[i]=(*itTrk)->partId();
01611
01612 if(m_skim4pi)
01613 {
01614 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
01615 {
01616 if(i<4) (*itTrk)->setPartId(3);
01617 }
01618
01619 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1&&dtpr2kpilst[0]<-0.4&&dtpr2kpilst[1]<-0.4&&dtpr2kpilst[2]<-0.4&&dtpr2kpilst[3]<-0.4)
01620 {
01621 if(i<4) (*itTrk)->setPartId(4);
01622 }
01623
01624
01625 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
01626 {
01627 if(i==0||i==1) (*itTrk)->setPartId(3);
01628 if(i==2||i==3) (*itTrk)->setPartId(5);
01629 }
01630
01631 }
01632
01633 ParticleID *pid = ParticleID::instance();
01634 pid->init();
01635 pid->setMethod(pid->methodProbability());
01636 pid->setChiMinCut(4);
01637 pid->setRecTrack(*itTrk);
01638 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
01639 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
01640 pid->identify(pid->onlyElectron() | pid->onlyMuon());
01641 pid->calculate();
01642 if(!(pid->IsPidInfoValid())) continue;
01643 m_costpid2kpi[i] = cos(mdcTrk->theta());
01644
01645
01646 m_probe2kpi[i] = pid->probElectron();
01647 m_probmu2kpi[i] = pid->probMuon();
01648 m_probpi2kpi[i] = pid->probPion();
01649 m_probk2kpi[i] = pid->probKaon();
01650 m_probpr2kpi[i] = pid->probProton();
01651
01652
01653
01654
01655
01656
01657 }
01658 }
01659
01660 }
01661
01662
01663 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
01664 {
01665 jGam2kpi.push_back(igam2kpi);
01666 }
01667
01668 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
01669 if(i>=evtRecTrkCol->size()) break;
01670 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
01671 if(!(*itTrk)->isEmcShowerValid()) continue;
01672 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
01673 {
01674 if(i!=igam2kpi) jGam2kpi.push_back((*itTrk)->trackId());
01675 }
01676 else{
01677 jGam2kpi.push_back((*itTrk)->trackId());
01678 }
01679
01680 }
01681
01682 int ngam2kpi=jGam2kpi.size();
01683
01684 if(m_rootput)
01685 {
01686 for(int i = 0; i< ngam2kpi; i++) {
01687 if(i>=evtRecTrkCol->size()) break;
01688 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGam2kpi[i];
01689 if(!(*itTrk)->isEmcShowerValid()) continue;
01690 RecEmcShower *emcTrk = (*(evtRecTrkCol->begin()+jGam2kpi[i]))->emcShower();
01691 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
01692 double dthe = 200.;
01693 double dphi = 200.;
01694 double dang = 200.;
01695
01696
01697
01698 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
01699 if(j>=evtRecTrkCol->size()) break;
01700 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
01701 if(!(*jtTrk)->isExtTrackValid()) continue;
01702 RecExtTrack *extTrk = (*jtTrk)->extTrack();
01703 if(extTrk->emcVolumeNumber() == -1) continue;
01704 Hep3Vector extpos = extTrk->emcPosition();
01705 double angd = extpos.angle(emcpos);
01706 double thed = extpos.theta() - emcpos.theta();
01707 double phid = extpos.deltaPhi(emcpos);
01708 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
01709 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
01710
01711
01712
01713 if(angd < dang) { dang = angd;
01714 dthe = thed;
01715 dphi = phid;
01716 }
01717
01718 }
01719 if(dang>=200) continue;
01720 double eraw = emcTrk->energy();
01721 dthe = dthe * 180 / (CLHEP::pi);
01722 dphi = dphi * 180 / (CLHEP::pi);
01723 dang = dang * 180 / (CLHEP::pi);
01724
01725
01726
01727 m_numHits[i]= emcTrk->numHits();
01728 m_secondmoment[i] = emcTrk->secondMoment();
01729 m_latmoment[i] = emcTrk->latMoment();
01730 m_timegm[i] = emcTrk->time();
01731 m_cellId[i]=emcTrk->cellId();
01732 m_module[i]=emcTrk->module();
01733 m_a20Moment[i]=emcTrk->a20Moment();
01734 m_a42Moment[i]=emcTrk->a42Moment();
01735 m_getShowerId[i]=emcTrk->getShowerId();
01736 m_getClusterId[i]=emcTrk->getClusterId();
01737 m_getEAll[i]=emcTrk->getEAll();
01738 m_x[i]= emcTrk->x();
01739 m_y[i]= emcTrk->y();
01740 m_z[i]= emcTrk->z();
01741 m_cosemc[i] = cos(emcTrk->theta());
01742 m_phiemc[i] = emcTrk->phi();
01743 m_energy[i] = emcTrk->energy();
01744 m_eSeed[i] = emcTrk->eSeed();
01745 m_e3x3[i] = emcTrk->e3x3();
01746 m_e5x5[i] = emcTrk->e5x5();
01747 m_dang4c[i] = dang;
01748 m_dthe4c[i] = dthe;
01749 m_dphi4c[i] = dphi;
01750
01751
01752
01753
01754 }
01755 }
01756 }
01757
01758 }
01759
01760
01761
01762
01763
01764
01765 if(m_skim4pi)
01766 {
01767
01768 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
01769 {m_subalg1->execute();
01770 Ncut[6]++;
01771 }
01772 }
01773
01774
01775 if(m_skim4k)
01776 {
01777
01778 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1&&dtpr2kpilst[0]<-0.4&&dtpr2kpilst[1]<-0.4&&dtpr2kpilst[2]<-0.4&&dtpr2kpilst[3]<-0.4)
01779 {m_subalg2->execute();
01780 Ncut[7]++;
01781 }
01782 }
01783
01784 if(m_skim2pi2pr)
01785 {
01786 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
01787 {m_subalg3->execute();
01788
01789 Ncut[8]++;
01790 }
01791 }
01792
01793
01794
01795
01796 if(m_rootput)
01797 {
01798
01799 if((mppreclst<3.06&&chis4c2kpilst<40)||((meelst>3.06&&meelst<3.12)&&(fabs(mppreclst-3.097)<0.01)))
01800 { Ncut[9]++;
01801 m_tuple1->write();
01802
01803 }
01804 }
01805
01806 return StatusCode::SUCCESS;
01807 }
01808
01809
01810
01811 StatusCode Gam4pikp::finalize() {
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822 MsgStream log(msgSvc(), name());
01823 log << MSG::INFO << "in finalize()" << endmsg;
01824 return StatusCode::SUCCESS;
01825 }
01826
01827 void Gam4pikp::InitVar()
01828 {
01829
01830
01831 m_chis4c2kpi=9999.;
01832 m_chis2kpi=9999.;
01833 m_cla2kpi=9999.;
01834 m_ncy20=9999;
01835 m_ncy30=9999;
01836 m_meeall=9999;
01837 m_mpprecall=9999;
01838 for (int i=0; i<100; i++)
01839 {
01840
01841 m_angjs5[i]=9999;
01842 m_nearjs5[i]=9999;
01843 m_angjs6[i]=9999;
01844 m_nearjs6[i]=9999;
01845 m_ang4pi5[i]=9999;
01846 m_near4pi5[i]=9999;
01847 m_ang4pi6[i]=9999;
01848 m_near4pi6[i]=9999;
01849 m_probe2kpi[i]=-1;
01850 m_probmu2kpi[i]=-1;
01851 m_probpi2kpi[i]=-1;
01852 m_probk2kpi[i]=-1;
01853 m_probpr2kpi[i]=-1;
01854
01855 m_probejs[i]=-1;
01856 m_probmujs[i]=-1;
01857 m_probpijs[i]=-1;
01858 m_probkjs[i]=-1;
01859 m_probprjs[i]=-1;
01860
01861 m_cbmc[i]=9999;
01862 m_bjemcjs[i]=0;
01863 m_bjemc2kpi[i]=0;
01864 m_bjtofjs[i]=0;
01865 m_bjtof2kpi[i]=0;
01866 m_bjmucjs[i]=0;
01867 m_bjmuc2kpi[i]=0;
01868 m_charge2kpi[i]=9999;
01869 m_ghitjs[i] = 9999;
01870 m_thitjs[i] = 9999;
01871 m_probphjs[i] = 99999;
01872 m_normphjs[i] = 99999;
01873
01874 m_ghit2kpi[i] = 9999;
01875 m_thit2kpi[i] = 9999;
01876 m_probph2kpi[i] = 99999;
01877 m_normph2kpi[i] = 99999;
01878
01879
01880
01881 m_counterjs[i] = 9999;
01882 m_barreljs[i] = 9999;
01883 m_layertofjs[i] = 9999;
01884 m_readoutjs[i] = 9999;
01885 m_clusterjs[i] = 9999;
01886
01887 m_counter2kpi[i] = 9999;
01888 m_barrel2kpi[i] = 9999;
01889 m_layertof2kpi[i] = 9999;
01890 m_readout2kpi[i] = 9999;
01891 m_cluster2kpi[i] = 9999;
01892
01893 m_comcs2kpi[i]=9999;
01894 m_dte2kpi[i]=9999;
01895 m_dtmu2kpi[i]=9999;
01896 m_dtpi2kpi[i]=9999;
01897 m_dtk2kpi[i]=9999;
01898 m_dtpr2kpi[i]=9999;
01899 m_sigmaetof2kpi[i]=9999;
01900 m_sigmamutof2kpi[i]=9999;
01901 m_sigmapitof2kpi[i]=9999;
01902 m_sigmaktof2kpi[i]=9999;
01903 m_sigmaprtof2kpi[i]=9999;
01904 m_t0tof2kpi[i]=9999;
01905 m_errt0tof2kpi[i]=9999;
01906
01907 m_sigmaetofjs[i]=9999;
01908 m_sigmamutofjs[i]=9999;
01909 m_sigmapitofjs[i]=9999;
01910 m_sigmaktofjs[i]=9999;
01911 m_sigmaprtofjs[i]=9999;
01912 m_t0tofjs[i]=9999;
01913 m_errt0tofjs[i]=9999;
01914
01915 m_dtejs[i]=9999;
01916 m_dtmujs[i]=9999;
01917 m_dtpijs[i]=9999;
01918 m_dtkjs[i]=9999;
01919 m_dtprjs[i]=9999;
01920
01921 m_chie2kpi[i]=9999;
01922 m_chimu2kpi[i]=9999;
01923 m_chipi2kpi[i]=9999;
01924 m_chik2kpi[i]=9999;
01925 m_chip2kpi[i]=9999;
01926 m_pidnum2kpi[i]=9999;
01927 m_chiejs[i]=9999;
01928 m_chimujs[i]=9999;
01929 m_chipijs[i]=9999;
01930 m_chikjs[i]=9999;
01931 m_chipjs[i]=9999;
01932
01933
01934
01935
01936 }
01937
01938
01939 }
01940
01941 void Gam4pikp::BubbleSort(std::vector<double> &p, std::vector<int> &trkid)
01942 {
01943 if ( (int) p.size() != (int) trkid.size() )
01944 {
01945 std::cout << "the two vector have different length!" << std::endl;
01946 std::cout << "the size of p is " << (int) p.size() << std::endl;
01947 std::cout << "the size of trkid is " << (int) trkid.size() << std::endl;
01948 std::cout << "Please check your code!" << std::endl;
01949 exit(1);
01950 }
01951 unsigned int vsize = p.size();
01952 double ptemp;
01953 int idtemp;
01954 for ( unsigned int i=0; i<vsize-1; i++ )
01955 {
01956 for ( unsigned int j=i+1; j<vsize; j++ )
01957 {
01958 if ( p[i] > p[j] )
01959 {
01960 ptemp = p[i]; idtemp = trkid[i];
01961 p[i] = p[j]; trkid[i] = trkid[j];
01962 p[j] = ptemp; trkid[j] = idtemp;
01963 }
01964 }
01965 }
01966
01967 }
01968