| 1 | /*
|
|---|
| 2 | * ---- Delphes ----
|
|---|
| 3 | * A Fast Simulator for general purpose LHC detector
|
|---|
| 4 | * S. Ovyn ~~~~ severine.ovyn@uclouvain.be
|
|---|
| 5 | *
|
|---|
| 6 | * Center for Particle Physics and Phenomenology (CP3)
|
|---|
| 7 | * Universite Catholique de Louvain (UCL)
|
|---|
| 8 | * Louvain-la-Neuve, Belgium
|
|---|
| 9 | * */
|
|---|
| 10 |
|
|---|
| 11 | #include "interface/BFieldProp.h"
|
|---|
| 12 | #include<cmath>
|
|---|
| 13 | #include "TMath.h"
|
|---|
| 14 | using namespace std;
|
|---|
| 15 |
|
|---|
| 16 |
|
|---|
| 17 | //------------------------------------------------------------------------------
|
|---|
| 18 |
|
|---|
| 19 | TrackPropagation::TrackPropagation(const string DetDatacard):
|
|---|
| 20 | MAXITERATION(10000), q(-9999.), phi_0(-9999.), gammam(-9999.), omega(-9999.), r(-9999.),
|
|---|
| 21 | x_c(-9999.), y_c(-9999.), R_c(-9999.), Phi_c(-9999.),
|
|---|
| 22 | rr(-9999.), t(-9999.), t_z(-9999.), t_T(-9999.),
|
|---|
| 23 | x_t(-9999.), y_t(-9999.), z_t(-9999.),
|
|---|
| 24 | R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.),
|
|---|
| 25 | Px_t(-9999), Py_t(-9999), Pz_t(-9999), PT_t(-9999), p_t(-9999), E_t(-9999) {
|
|---|
| 26 |
|
|---|
| 27 | // if(DetDatacard="") { DET = new RESOLution(); }
|
|---|
| 28 | // else DET = new RESOLution(DetDatacard);
|
|---|
| 29 | DET = new RESOLution();
|
|---|
| 30 | DET->ReadDataCard(DetDatacard);
|
|---|
| 31 |
|
|---|
| 32 | // magnetic field parameters
|
|---|
| 33 | R_max = DET->TRACK_radius;
|
|---|
| 34 | z_max = DET->TRACK_length/2.;
|
|---|
| 35 | B_x = DET->TRACK_bfield_x;
|
|---|
| 36 | B_y = DET->TRACK_bfield_y;
|
|---|
| 37 | B_z = DET->TRACK_bfield_z;
|
|---|
| 38 |
|
|---|
| 39 | loop_overflow_counter=0;
|
|---|
| 40 | }
|
|---|
| 41 |
|
|---|
| 42 | void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
|
|---|
| 43 |
|
|---|
| 44 | q = Charge(Part->PID);
|
|---|
| 45 | if(q==0) return;
|
|---|
| 46 |
|
|---|
| 47 | if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
|
|---|
| 48 | if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
|
|---|
| 49 |
|
|---|
| 50 | if (B_x== 0 && B_y== 0) { // faster if only B_z
|
|---|
| 51 | if (B_z==0) return; // nothing to do
|
|---|
| 52 |
|
|---|
| 53 | // initial conditions:
|
|---|
| 54 | // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
|
|---|
| 55 | // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
|
|---|
| 56 |
|
|---|
| 57 | // 1. initial transverse momentum p_{T0} : Part->PT
|
|---|
| 58 | // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
|
|---|
| 59 | // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
|
|---|
| 60 | // giration frequency \omega = q/(gamma m) B_z
|
|---|
| 61 | // helix radius r = p_T0 / (omega gamma m)
|
|---|
| 62 | phi_0 = -atan2(Part->Px,Part->Py);
|
|---|
| 63 | gammam = Part->E; // here c==1
|
|---|
| 64 | //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
|
|---|
| 65 | omega = q * B_z /gammam;
|
|---|
| 66 | r = Part->PT / (omega * gammam);
|
|---|
| 67 |
|
|---|
| 68 | // 2. Helix parameters : center coordinates in transverse plane
|
|---|
| 69 | // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
|
|---|
| 70 | // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
|
|---|
| 71 | x_c = Part->X - r*cos(phi_0); /// TEST !!
|
|---|
| 72 | y_c = Part->Y - r*sin(phi_0);
|
|---|
| 73 | R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
|
|---|
| 74 | Phi_c = atan2(y_c,x_c);
|
|---|
| 75 |
|
|---|
| 76 | // 3. time evaluation t = min(t_T, t_z)
|
|---|
| 77 | // t_T : time to exit from the sides
|
|---|
| 78 | // t_T= [ Phi_c - phi_0 + atan( (R_max^2 - (R_c^2 + r^2))/(2rR_c) ) ]/omega
|
|---|
| 79 | // t_z : time to exit from the front or the back
|
|---|
| 80 | // t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
|
|---|
| 81 | rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
|
|---|
| 82 | t_T=0;
|
|---|
| 83 | t_z = gammam / Part->Pz * (-Part->Z + z_max* TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
|
|---|
| 84 | if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z;
|
|---|
| 85 | else {
|
|---|
| 86 | t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
|
|---|
| 87 | t = min(t_T,t_z);
|
|---|
| 88 | }
|
|---|
| 89 |
|
|---|
| 90 | // 4. position in terms of x(t), y(t), z(t)
|
|---|
| 91 | // x(t) = x_c + r cos (omega t + phi_0)
|
|---|
| 92 | // y(t) = y_c + r sin (omega t + phi_0)
|
|---|
| 93 | // z(t) = z_0 + (p_Z0/gammam) t
|
|---|
| 94 | x_t = x_c + r * cos(omega * t + phi_0);
|
|---|
| 95 | y_t = y_c + r * sin(omega * t + phi_0);
|
|---|
| 96 | z_t = Part->Z + Part->Pz / gammam * t;
|
|---|
| 97 |
|
|---|
| 98 | // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
|
|---|
| 99 | // R(t) = sqrt(x(t)² + y(t)²)
|
|---|
| 100 | // Phi(t) = atan(y(t)/x(t))
|
|---|
| 101 | // Theta(t) = atan(R(t)/z(t))
|
|---|
| 102 | // Eta(t) = -ln tan (Theta(t)/2)
|
|---|
| 103 | R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
|
|---|
| 104 | Phi_t = atan2( y_t, x_t);
|
|---|
| 105 | /* if(R_t>0) {
|
|---|
| 106 | Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
|
|---|
| 107 | Eta_t = - log(tan(Theta_t/2.));
|
|---|
| 108 | } else{
|
|---|
| 109 | Theta_t=0; Eta_t = 9999;
|
|---|
| 110 | }
|
|---|
| 111 | */
|
|---|
| 112 | Px_t = - Part->PT * sin(omega*t + phi_0);
|
|---|
| 113 | Py_t = Part->PT * cos(omega*t + phi_0);
|
|---|
| 114 | Pz_t = Part->Pz;
|
|---|
| 115 | PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
|
|---|
| 116 | p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
|
|---|
| 117 | E_t=sqrt(Part->M*Part->M +p_t);
|
|---|
| 118 | if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
|
|---|
| 119 | if(p_t>0) Theta_t = acos(Pz_t/p_t);
|
|---|
| 120 | momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
|
|---|
| 121 |
|
|---|
| 122 | // test zone ---
|
|---|
| 123 | /*
|
|---|
| 124 | cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
|
|---|
| 125 | double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
|
|---|
| 126 | double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
|
|---|
| 127 | if(0 && fabs(Eta_t -Eta_t2)>1e-310) {
|
|---|
| 128 | cout << "ERROR-BUG: Eta_t != Eta_t2\n";
|
|---|
| 129 | cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
|
|---|
| 130 | }
|
|---|
| 131 |
|
|---|
| 132 | double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
|
|---|
| 133 | if(fabs(R_t - R_t2) > 1e-7)
|
|---|
| 134 | cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
|
|---|
| 135 | if( fabs(E - gammam) > 1e-3 ) {
|
|---|
| 136 | cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
|
|---|
| 137 | cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
|
|---|
| 138 | if( fabs(PT_t - Part->PT) > 1e-10 ) {
|
|---|
| 139 | cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
|
|---|
| 140 | cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
|
|---|
| 141 | }
|
|---|
| 142 | if(momentum.Pz() != Pz_t)
|
|---|
| 143 | cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
|
|---|
| 144 |
|
|---|
| 145 | double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
|
|---|
| 146 | if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) {
|
|---|
| 147 | cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
|
|---|
| 148 | cout << temp_p << "\t" << temp_p0 << endl;
|
|---|
| 149 | }
|
|---|
| 150 |
|
|---|
| 151 | // if x_c == y_c ==0 (set it by hand!), easy cross-check
|
|---|
| 152 | //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
|
|---|
| 153 | */
|
|---|
| 154 |
|
|---|
| 155 | } else { // if B_x or B_y are non zero: longer computation
|
|---|
| 156 |
|
|---|
| 157 | float Xvertex1 = Part->X;
|
|---|
| 158 | float Yvertex1 = Part->Y;
|
|---|
| 159 | float Zvertex1 = Part->Z;
|
|---|
| 160 |
|
|---|
| 161 | //out of tracking coverage?
|
|---|
| 162 | if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
|
|---|
| 163 | if(fabs(Zvertex1) > z_max){return;}
|
|---|
| 164 |
|
|---|
| 165 | double px = Part->Px / 0.003;
|
|---|
| 166 | double py = Part->Py / 0.003;
|
|---|
| 167 | double pz = Part->Pz / 0.003;
|
|---|
| 168 | double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
|
|---|
| 169 | double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
|
|---|
| 170 |
|
|---|
| 171 | double M = Part->M;
|
|---|
| 172 | double vx = px/M;
|
|---|
| 173 | double vy = py/M;
|
|---|
| 174 | double vz = pz/M;
|
|---|
| 175 | double qm = q/M;
|
|---|
| 176 |
|
|---|
| 177 | double ax = qm*(B_z*vy - B_y*vz);
|
|---|
| 178 | double ay = qm*(B_x*vz - B_z*vx);
|
|---|
| 179 | double az = qm*(B_y*vx - B_x*vy);
|
|---|
| 180 | double dt = 1/p;
|
|---|
| 181 | if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
|
|---|
| 182 |
|
|---|
| 183 | double xold=Xvertex1; double x=xold;
|
|---|
| 184 | double yold=Yvertex1; double y=yold;
|
|---|
| 185 | double zold=Zvertex1; double z=zold;
|
|---|
| 186 |
|
|---|
| 187 | double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
|
|---|
| 188 |
|
|---|
| 189 | unsigned int k = 0;
|
|---|
| 190 | double VTratio=0;
|
|---|
| 191 | double R_max2 = R_max*R_max;
|
|---|
| 192 | double r2=0; // will be x*x+y*y
|
|---|
| 193 |
|
|---|
| 194 | while(k < MAXITERATION){
|
|---|
| 195 | k++;
|
|---|
| 196 |
|
|---|
| 197 | vx += ax*dt;
|
|---|
| 198 | vy += ay*dt;
|
|---|
| 199 | vz += az*dt;
|
|---|
| 200 |
|
|---|
| 201 | VTratio = VTold/sqrt(vx*vx+vy*vy);
|
|---|
| 202 | vx *= VTratio;
|
|---|
| 203 | vy *= VTratio;
|
|---|
| 204 |
|
|---|
| 205 | ax = qm*(B_z*vy - B_y*vz);
|
|---|
| 206 | ay = qm*(B_x*vz - B_z*vx);
|
|---|
| 207 | az = qm*(B_y*vx - B_x*vy);
|
|---|
| 208 |
|
|---|
| 209 | x += vx*dt;
|
|---|
| 210 | y += vy*dt;
|
|---|
| 211 | z += vz*dt;
|
|---|
| 212 | r2 = x*x + y*y;
|
|---|
| 213 |
|
|---|
| 214 | if( r2 > R_max2 ){
|
|---|
| 215 | x /= r2/R_max2;
|
|---|
| 216 | y /= r2/R_max2;
|
|---|
| 217 | break;
|
|---|
| 218 | }
|
|---|
| 219 | if( fabs(z)>z_max)break;
|
|---|
| 220 |
|
|---|
| 221 | xold = x;
|
|---|
| 222 | yold = y;
|
|---|
| 223 | zold = z;
|
|---|
| 224 | } // while loop
|
|---|
| 225 |
|
|---|
| 226 | if(k == MAXITERATION) loop_overflow_counter++;
|
|---|
| 227 | //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
|
|---|
| 228 |
|
|---|
| 229 | if(x!=0 && y!=0 && z!=0) {
|
|---|
| 230 | float Theta = atan2(sqrt(r2),z);
|
|---|
| 231 | double eta = -log(tan(Theta/2.));
|
|---|
| 232 | double phi = atan2(y,x);
|
|---|
| 233 | momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
|
|---|
| 234 | }
|
|---|
| 235 |
|
|---|
| 236 | } // if b_x or b_y non zero
|
|---|
| 237 | }
|
|---|
