00001 00002 // Class K0Event 00003 // 00004 // Author: Benno List, Jenny Boehme 00005 // Last update: $Date: 2005/01/12 10:11:45 $ 00006 // by: $Author: blist $ 00007 // 00008 // Description: class to generate and fit K0S decays 00009 // 00011 #include "jbltools/kinfit/K0Event.h" 00012 00013 #include "jbltools/kinfit/JetFitObject.h" 00014 #include "jbltools/kinfit/MassConstraint.h" 00015 #include "jbltools/kinfit/cernlib.h" 00016 #include "jbltools/kinfit/NeutralParticleTrack.h" 00017 #include "jbltools/kinfit/ChargedParticleTrack.h" 00018 #include "jbltools/kinfit/ThreeVector.h" 00019 #include "jbltools/kinfit/FourVector.h" 00020 #include "jbltools/kinfit/TTVertexConstraint.h" 00021 #include "jbltools/kinfit/VertexConstraint.h" 00022 #include "jbltools/kinfit/VertexFitObject.h" 00023 #include "jbltools/kinfit/TrackMomentumConstraint.h" 00024 00025 #include <iostream> // - cout 00026 using std::cout; 00027 using std::endl; 00028 #include <cmath> 00029 using std::pow; 00030 using std::atan; 00031 using std::exp; 00032 using std::sin; 00033 using std::cos; 00034 using std::log; 00035 using std::abs; 00036 00037 // constructor: 00038 K0Event::K0Event() 00039 { 00040 for (int i = 0; i < NFV; ++i) fv[i] = 0; 00041 for (int i = 0; i < NTFO; ++i) gentrack[i] = rectrack[i] = 0; 00042 for (int i = 0; i < NVER; ++i) genvert[i] = 0; 00043 }; 00044 00045 //destructor: 00046 K0Event::~K0Event() { 00047 for (int i = 0; i < NFV; ++i) delete fv[i]; 00048 for (int i = 0; i < NTFO; ++i) { 00049 delete gentrack[i]; 00050 delete rectrack[i]; 00051 } 00052 for (int i = 0; i < NVER; ++i) delete genvert[i]; 00053 }; 00054 00055 00056 // generate four vectors and tracks 00057 void K0Event::genEvent(){ 00059 FReal randoms[4]; 00060 rnorml (randoms, 3); 00061 ThreeVector *primvert = genvert[0] = new ThreeVector (-0.4+0.015*randoms[0], 00062 0.2+0.004*randoms[1], 00063 1.0+10.5*randoms[2]); 00064 // ThreeVector *primvert = genvert[0] = new ThreeVector (0.015*randoms[0], 00065 // 0.004*randoms[1], 00066 // 10.5*randoms[2]); 00067 // ThreeVector *primvert = genvert[0] = new ThreeVector (0, 0, 0); 00069 ranmar (randoms, 4); 00070 double phi = 6.2831853*randoms[0]; 00071 double ptmin = 1.0; 00072 double ptmax = 10.; 00073 int n = 4; // pt^-4 00074 double pt = pow (randoms[1]*(pow (ptmin, -n+1) - pow (ptmax, -n+1)) + pow (ptmax, -n+1), -1./(n-1)); 00075 double etamin = -1.5; 00076 double etamax = 1.5; 00077 double eta = randoms[1]*(etamax-etamin) + etamin; 00078 double theta = 2*atan (exp (-eta)); 00079 00080 // cout << "k0s phi = " << phi << ", pt = " << pt << ", eta = " << eta << ", theta = " << theta << endl; 00081 00082 ThreeVector k0smom (pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta)); 00083 FourVector *k0s = fv[0] = new FourVector (k0smom, 0.497648); 00084 00085 // cout << "k0s phi = " << k0s->getPhi() << ", pt = " << k0s->getPt() << ", eta = " 00086 // << k0s->getEta() << ", theta = " << k0s->getTheta() << endl; 00087 00088 00089 FourVector *pi1 = fv[1] = new FourVector (0.13957, 0, 0, 0); 00090 FourVector *pi2 = fv[2] = new FourVector (0.13957, 0, 0, 0); 00091 00092 k0s->decayto (*pi1, *pi2); 00093 00094 // cout << "k0s 4-vector: " << *k0s << endl; 00095 // cout << "pi1+p2 4-vector: " << *pi1+*pi2 << endl; 00096 // cout << "pi1 4-vector: " << *pi1 << endl; 00097 // cout << "pi2 4-vector: " << *pi2 << endl << endl; 00098 00099 00100 double slen = -sin(theta)*log (randoms[3])*2.6842*k0s->getBetaGamma().getMag(); 00101 // double slen = -sin(theta)*log (0.5)*2.6842*k0s->getBetaGamma().getMag(); 00102 00103 // cout << "slen = : " << slen << endl << endl; 00104 00105 NeutralParticleTrack *k0strack = new NeutralParticleTrack ("K0 gen", *primvert, k0smom, 0.497648); 00106 gentrack[0] = k0strack; 00107 00108 ThreeVector *decvert = genvert[1] = new ThreeVector; 00109 double sstop = k0strack->getArcLength(0) + slen; 00110 // k0strack->setArcLength(1, sstop); 00111 k0strack->getTrajectoryPointEx(sstop, *decvert); 00112 00113 ChargedParticleTrack *pi1track; 00114 ChargedParticleTrack *pi2track; 00115 gentrack[1] = pi1track = new ChargedParticleTrack ("pi1 gen", *decvert, pi1->getThreeVector(), 0.13957, +1); 00116 gentrack[2] = pi2track = new ChargedParticleTrack ("pi2 gen", *decvert, pi2->getThreeVector(), 0.13957, -1); 00117 00118 // cout << "k0s track: " << *gentrack[0] << endl; 00119 // cout << "pi1 track: " << *gentrack[1] << endl; 00120 // cout << "pi2 track: " << *gentrack[2] << endl << endl; 00121 // 00122 // cout << "primary vertex: " << *genvert[0] << endl; 00123 // cout << "k0s vertex: " << gentrack[0]->getVertex(0) << endl << endl; 00124 // 00125 // cout << "decay vertex: " << *genvert[1] << endl; 00126 // cout << "pi1 vertex: " << gentrack[1]->getVertex(0) << endl; 00127 // cout << "pi2 vertex: " << gentrack[2]->getVertex(0) << endl << endl; 00128 // 00129 // cout << "k0s 4-vector: " << *k0s << endl; 00130 // cout << "k0s track 4-vector: " << gentrack[0]->getMomentum(0) << endl; 00131 // cout << "pi1 4-vector: " << *pi1 << endl; 00132 // cout << "pi1 track 4-vector: " << gentrack[1]->getMomentum(0) << endl; 00133 // cout << "pi2 4-vector: " << *pi2 << endl; 00134 // cout << "pi2 track 4-vector: " << gentrack[2]->getMomentum(0) << endl << endl; 00135 00136 rectrack[1] = createSmearedChargedTrack ("pi1 rec", *pi1track); 00137 rectrack[2] = createSmearedChargedTrack ("pi2 rec", *pi2track); 00138 smtrack[1] = new ChargedParticleTrack (* (ChargedParticleTrack*)rectrack[1]); 00139 smtrack[2] = new ChargedParticleTrack (* (ChargedParticleTrack*)rectrack[2]); 00140 smtrack[1]->setName ("pi1 sm"); 00141 smtrack[2]->setName ("pi2 sm"); 00142 00143 ThreeVector nominalvertex; 00144 ThreeVector zeromomentum; 00145 rectrack[0] = new NeutralParticleTrack ("K0 rec", nominalvertex, zeromomentum, 0.497648); 00146 rectrack[0]->fixVertexParam(0); 00147 rectrack[0]->fixVertexParam(1, false); 00148 rectrack[0]->fixParam(5, true); 00149 00150 // cout << "After smearing:\n"; 00151 // cout << "k0s track 4-vector: " << rectrack[0]->getMomentum(0) << endl; 00152 // cout << "pi1 track 4-vector: " << rectrack[1]->getMomentum(0) << endl; 00153 // cout << "pi2 track 4-vector: " << rectrack[2]->getMomentum(0) << endl; 00154 00155 00156 }; 00157 00158 ChargedParticleTrack *K0Event::createSmearedChargedTrack (const char *name, const ChargedParticleTrack& in) { 00159 double kappa = in.getParam(0)*ChargedParticleTrack::parfact[0]; 00160 double phi0 = in.getParam(1)*ChargedParticleTrack::parfact[1]; 00161 double theta = in.getParam(2)*ChargedParticleTrack::parfact[2]; 00162 double dca = in.getParam(3)*ChargedParticleTrack::parfact[3]; 00163 double z0 = in.getParam(4)*ChargedParticleTrack::parfact[4]; 00164 00165 // pt=cBq/kappa => dpt = cBq/kappa^2 dkappa, dpt/pt = dkappa/kappa = pt * dkappa/cBq 00166 // => dkappa = dpt/pt^2*cBq = 0.01/GeV*cBq 00167 double dkappa = 0.01*abs(0.002998*TrackFitObject::bfield); 00168 double dphi0 = 0.0002; // 0.2mrad = 100um at 50cm 00169 double dtheta = 0.08; // = 2cm / 25cm 00170 double ddca = 0.05; // 0.5mm 00171 double dz0 = 2.; // 2cm 00172 00173 FReal randoms[5]; 00174 rnorml (randoms, 5); 00175 00176 kappa += dkappa*randoms[0]; 00177 phi0 += dphi0 *randoms[1]; 00178 theta += dtheta*randoms[2]; 00179 dca += ddca *randoms[3]; 00180 z0 += dz0 *randoms[4]; 00181 00182 ChargedParticleTrack *result = new ChargedParticleTrack (name, kappa, phi0, theta, dca, z0, 00183 in.getMass(), in.getCharge(), 0); 00184 00185 result->setError (0, dkappa/ChargedParticleTrack::parfact[0]); 00186 result->setError (1, dphi0/ChargedParticleTrack::parfact[1]); 00187 result->setError (2, dtheta/ChargedParticleTrack::parfact[2]); 00188 result->setError (3, ddca/ChargedParticleTrack::parfact[3]); 00189 result->setError (4, dz0/ChargedParticleTrack::parfact[4]); 00190 00191 return result; 00192 } 00193 00194 // fit it! 00195 int K0Event::fitEvent (BaseFitter& fitter){ 00196 00197 fitter.addFitObject (rectrack[0]); 00198 fitter.addFitObject (rectrack[1]); 00199 fitter.addFitObject (rectrack[2]); 00200 00201 VertexFitObject *k0decvertex = new VertexFitObject ("K0dec", 0, 0, 0); 00202 fitter.addFitObject (k0decvertex); 00203 00204 rectrack[0]->fixVertexParam (0, true); 00205 k0decvertex->addTrack (rectrack[0], true, false); 00206 k0decvertex->addTrack (rectrack[1], false, true); 00207 k0decvertex->addTrack (rectrack[2], false, true); 00208 00209 k0decvertex->addConstraints (fitter, VertexFitObject::VXYZ|VertexFitObject::PXYZ); 00210 00211 k0decvertex->initForFit(); 00212 00213 00214 // TTVertexConstraint vc_x1 (*rectrack[0], 0, *rectrack[1], 0, 0); 00215 // TTVertexConstraint vc_y1 (*rectrack[0], 0, *rectrack[1], 0, 1); 00216 // TTVertexConstraint vc_z1 (*rectrack[0], 0, *rectrack[1], 0, 2); 00217 // TTVertexConstraint vc_x2 (*rectrack[0], 0, *rectrack[2], 0, 0); 00218 // TTVertexConstraint vc_y2 (*rectrack[0], 0, *rectrack[2], 0, 1); 00219 // TTVertexConstraint vc_z2 (*rectrack[0], 0, *rectrack[2], 0, 2); 00220 // fitter.addConstraint (vc_x1); 00221 // fitter.addConstraint (vc_y1); 00222 // fitter.addConstraint (vc_z1); 00223 // fitter.addConstraint (vc_x2); 00224 // fitter.addConstraint (vc_y2); 00225 // fitter.addConstraint (vc_z2); 00226 00227 // VertexConstraint vc_x0 (*k0decvertex, *rectrack[0], 0, 0); 00228 // VertexConstraint vc_y0 (*k0decvertex, *rectrack[0], 0, 1); 00229 // VertexConstraint vc_z0 (*k0decvertex, *rectrack[0], 0, 2); 00230 // VertexConstraint vc_x1 (*k0decvertex, *rectrack[1], 0, 0); 00231 // VertexConstraint vc_y1 (*k0decvertex, *rectrack[1], 0, 1); 00232 // VertexConstraint vc_z1 (*k0decvertex, *rectrack[1], 0, 2); 00233 // VertexConstraint vc_x2 (*k0decvertex, *rectrack[2], 0, 0); 00234 // VertexConstraint vc_y2 (*k0decvertex, *rectrack[2], 0, 1); 00235 // VertexConstraint vc_z2 (*k0decvertex, *rectrack[2], 0, 2); 00236 00237 // cout << "vc_x0: " << vc_x0.getValue() << endl; 00238 // cout << "vc_y0: " << vc_y0.getValue() << endl; 00239 // cout << "vc_z0: " << vc_z0.getValue() << endl; 00240 // cout << "vc_x1: " << vc_x1.getValue() << endl; 00241 // cout << "vc_y1: " << vc_y1.getValue() << endl; 00242 // cout << "vc_z1: " << vc_z1.getValue() << endl; 00243 // cout << "vc_x2: " << vc_x2.getValue() << endl; 00244 // cout << "vc_y2: " << vc_y2.getValue() << endl; 00245 // cout << "vc_z2: " << vc_z2.getValue() << endl; 00246 00247 // fitter.addConstraint (vc_x0); 00248 // fitter.addConstraint (vc_y0); 00249 // fitter.addConstraint (vc_z0); 00250 // fitter.addConstraint (vc_x1); 00251 // fitter.addConstraint (vc_y1); 00252 // fitter.addConstraint (vc_z1); 00253 // fitter.addConstraint (vc_x2); 00254 // fitter.addConstraint (vc_y2); 00255 // fitter.addConstraint (vc_z2); 00256 // 00257 // TrackMomentumConstraint p_x (1, 0, 0); 00258 // p_x.addToFOList(*rectrack[0]); 00259 // p_x.addToFOList(*rectrack[1]); 00260 // p_x.addToFOList(*rectrack[2]); 00261 // TrackMomentumConstraint p_y (0, 1, 0); 00262 // p_y.addToFOList(*rectrack[0]); 00263 // p_y.addToFOList(*rectrack[1]); 00264 // p_y.addToFOList(*rectrack[2]); 00265 // TrackMomentumConstraint p_z (0, 0, 1); 00266 // p_z.addToFOList(*rectrack[0]); 00267 // p_z.addToFOList(*rectrack[1]); 00268 // p_z.addToFOList(*rectrack[2]); 00269 // fitter.addConstraint (p_x); 00270 // fitter.addConstraint (p_y); 00271 // fitter.addConstraint (p_z); 00272 // rectrack[0]->fixParam(0); 00273 // rectrack[0]->fixParam(1); 00274 // rectrack[0]->fixParam(2); 00275 00276 // cout << "Before fitting:\n"; 00277 // cout << "k0s track: " << *rectrack[0] << endl; 00278 // cout << "pi1 track: " << *rectrack[1] << endl; 00279 // cout << "pi2 track: " << *rectrack[2] << endl; 00280 // cout << "k0s track 4-vector: " << rectrack[0]->getMomentum(0) << endl; 00281 // cout << "pi1 track 4-vector: " << rectrack[1]->getMomentum(0) << endl; 00282 // cout << "pi2 track 4-vector: " << rectrack[2]->getMomentum(0) << endl; 00283 // cout << "k0s track vertex: " << rectrack[0]->getVertex(0) << endl; 00284 // cout << "pi1 track vertex: " << rectrack[1]->getVertex(0) << endl; 00285 // cout << "pi2 track vertex: " << rectrack[2]->getVertex(0) << endl; 00286 // cout << "Constraint vc_x1: " << vc_x1.getValue() << endl; 00287 // cout << "Constraint vc_y1: " << vc_y1.getValue() << endl; 00288 // cout << "Constraint vc_z1: " << vc_z1.getValue() << endl; 00289 // cout << "Constraint vc_x2: " << vc_x2.getValue() << endl; 00290 // cout << "Constraint vc_y2: " << vc_y2.getValue() << endl; 00291 // cout << "Constraint vc_z2: " << vc_z2.getValue() << endl; 00292 // cout << "Constraint p_x: " << p_x.getValue() << endl; 00293 // cout << "Constraint p_y: " << p_y.getValue() << endl; 00294 // cout << "Constraint p_z: " << p_z.getValue() << endl; 00295 00296 double prob = fitter.fit(); 00297 00298 // cout << "After fitting:\n"; 00299 // cout << "k0s track: " << *rectrack[0] << endl; 00300 // cout << "pi1 track: " << *rectrack[1] << endl; 00301 // cout << "pi2 track: " << *rectrack[2] << endl; 00302 // cout << "k0s track 4-vector: " << rectrack[0]->getMomentum(0) << endl; 00303 // cout << "pi1 track 4-vector: " << rectrack[1]->getMomentum(0) << endl; 00304 // cout << "pi2 track 4-vector: " << rectrack[2]->getMomentum(0) << endl; 00305 // cout << "k0s track vertex: " << rectrack[0]->getVertex(0) << endl; 00306 // cout << "pi1 track vertex: " << rectrack[1]->getVertex(0) << endl; 00307 // cout << "pi2 track vertex: " << rectrack[2]->getVertex(0) << endl; 00308 // cout << "Constraint vc_x1: " << vc_x1.getValue() << endl; 00309 // cout << "Constraint vc_y1: " << vc_y1.getValue() << endl; 00310 // cout << "Constraint vc_z1: " << vc_z1.getValue() << endl; 00311 // cout << "Constraint vc_x2: " << vc_x2.getValue() << endl; 00312 // cout << "Constraint vc_y2: " << vc_y2.getValue() << endl; 00313 // cout << "Constraint vc_z2: " << vc_z2.getValue() << endl; 00314 // cout << "Constraint p_x: " << p_x.getValue() << endl; 00315 // cout << "Constraint p_y: " << p_y.getValue() << endl; 00316 // cout << "Constraint p_z: " << p_z.getValue() << endl; 00317 // cout << "fit probability = " << prob << endl << endl; 00318 00319 return fitter.getError(); 00320 };
1.3.2