Changeset 193 in svn
- Timestamp:
- Jan 26, 2009, 3:43:42 PM (16 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/interface/BFieldProp.h ¶
r100 r193 12 12 * */ 13 13 14 #include <vector>15 14 #include "TLorentzVector.h" 16 17 15 #include "interface/SmearUtil.h" 18 19 16 #include "Utilities/ExRootAnalysis/interface/BlockClasses.h" 20 17 #include "Utilities/ExRootAnalysis/interface/TSimpleArray.h" … … 26 23 public: 27 24 // Constructor 28 TrackPropagation( string DetDatacard);25 TrackPropagation(const string DetDatacard); 29 26 30 27 void Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum); 31 28 32 int MAXITERATION;33 int MINSEGLENGTH;34 29 35 30 private: 31 unsigned int MAXITERATION; 32 RESOLution *DET; 36 33 37 RESOLution *DET; 34 35 /// radial/longitudinal extensions of magnetic field volume 36 double R_max, z_max; 37 /// magnetic field components 38 double B_x, B_y, B_z; 39 /// particle charge 40 double q; 41 /// initial coordinate 42 double phi_0; 43 /// relativistic gamma \times m 44 double gammam; 45 /// giration frequency and radius 46 double omega, r; 47 /// center of the helix in the transverse plane 48 double x_c, y_c, R_c, Phi_c; 49 // variable for an intermediate computing 50 double rr; 51 /// times for exiting the magnetic field volume 52 double t, t_z, t_T; 53 /// coordinates of the exit point 54 double x_t, y_t, z_t, R_t, Phi_t, Theta_t, Eta_t; 55 56 unsigned int loop_overflow_counter; 38 57 }; 39 58 -
TabularUnified trunk/src/BFieldProp.cc ¶
r189 r193 10 10 11 11 #include "interface/BFieldProp.h" 12 #include "TRandom.h"13 14 #include <iostream>15 #include <sstream>16 #include <fstream>17 #include <iomanip>18 19 12 #include<cmath> 20 21 13 #include "TMath.h" 22 14 using namespace std; 23 15 … … 25 17 //------------------------------------------------------------------------------ 26 18 27 TrackPropagation::TrackPropagation(string DetDatacard) { 28 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 26 // if(DetDatacard="") { DET = new RESOLution(); } 27 // else DET = new RESOLution(DetDatacard); 29 28 DET = new RESOLution(); 30 29 DET->ReadDataCard(DetDatacard); 31 MAXITERATION = 10000; 32 MINSEGLENGTH = 70; 33 30 31 // magnetic field parameters 32 R_max = DET->TRACK_radius; 33 z_max = DET->TRACK_length/2.; 34 B_x = DET->TRACK_bfield_x; 35 B_y = DET->TRACK_bfield_y; 36 B_z = DET->TRACK_bfield_z; 37 38 loop_overflow_counter=0; 34 39 } 35 40 36 void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum) 37 { 38 39 double q = Charge(Part->PID); 40 if(q==0)return; 41 41 void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) { 42 43 q = Charge(Part->PID); 44 if(q==0) return; 45 46 if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;} 47 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;} 48 49 if (0 && B_x== 0 && B_y== 0) { // faster if only B_z 50 if (B_z==0) return; // nothing to do 51 52 // initial conditions: 53 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT; 54 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z; 55 56 // 1. initial transverse momentum p_{T0} : Part->PT 57 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) 58 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m 59 // giration frequency \omega = q/(gamma m) B_z 60 // helix radius r = p_T0 / (omega gamma m) 61 phi_0 = -atan2(Part->Px,Part->Py); 62 gammam = Part->E; // here c==1 63 //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl; 64 omega = q * B_z /gammam; 65 r = Part->PT / (omega * gammam); 66 67 // 2. Helix parameters : center coordinates in transverse plane 68 // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0) 69 // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c} 70 x_c = Part->X - r*cos(phi_0); 71 y_c = Part->Y - r*sin(phi_0); 72 R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) ); 73 Phi_c = atan2(y_c,x_c); 74 75 // 3. time evaluation t = min(t_T, t_z) 76 // t_T : time to exit from the sides 77 // t_z : time to exit from the front or the back 78 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable 79 t_T=0; 80 t_z = gammam / Part->Pz * (-Part->Z + DET-> TRACK_length/2. * TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ; 81 if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z; 82 else { 83 t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega; 84 t = min(t_T,t_z); 85 } 86 87 // 4. position in terms of x(t), y(t), z(t) 88 // x(t) = x_c + r cos (omega t + phi_0) 89 // y(t) = y_c + r sin (omega t + phi_0) 90 // z(t) = z_0 + (p_Z0/gammam) 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 = Part->Z + Part->Pz / gammam * t; 94 95 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t) 96 // R(t) = sqrt(x(t)² + y(t)²) 97 // Phi(t) = atan(y(t)/x(t)) 98 // Theta(t) = atan(R(t)/z(t)) 99 // Eta(t) = -ln tan (Theta(t)/2) 100 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) ); 101 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check 102 if(fabs(R_t - R_t2) > 1e-11) 103 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl; 104 Phi_t = atan2( y_t, x_t); 105 Theta_t = atan2( R_t , z_t); 106 Eta_t = - log(tan(Theta_t/2.)); 107 108 double Px_t = - Part->PT * sin(omega*t + phi_0); 109 double Py_t = Part->PT * cos(omega*t + phi_0); 110 double Pz_t = Part->Pz; 111 double PT_t = sqrt(Px_t*Px_t + Py_t*Py_t); 112 double E=sqrt(Part->M*Part->M + PT_t*PT_t + Pz_t*Pz_t); // cross-check 113 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E); 114 TLorentzVector mom_temp; 115 mom_temp.SetPtEtaPhiM(PT_t,Eta_t,Phi_t+PI,Part->M); 116 117 if (momentum != mom_temp) { 118 // cout << "momentum et mom_temp different\n"; 119 /* if (momentum.Px() != mom_temp.Px() ) cout << " px: " << momentum.Px() << " & " << mom_temp.Px() << endl; 120 if (momentum.Py() != mom_temp.Py() ) cout << " py: " << momentum.Py() << " & " << mom_temp.Py() << endl; 121 if (momentum.Pz() != mom_temp.Pz() ) cout << " pz: " << momentum.Pz() << " & " << mom_temp.Pz() << endl; 122 if (momentum.Pt() != mom_temp.Pt() ) cout << " pt: " << momentum.Pt() << " & " << mom_temp.Pt() << endl;*/ 123 // cout << "eta: " << momentum.Eta() << " " << mom_temp.Eta() << endl; 124 // cout << "the: " << momentum.Theta() << " " << mom_temp.Theta() << endl; 125 // cout << "phi: " << momentum.Phi() << " " << mom_temp.Phi() << endl; 126 } 127 //cout << Part->PT << " or " << momentum.Pt() << " or " << mom_temp.Pt() << endl; 128 129 if( fabs(E - gammam) > 1e-4 ) { 130 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n"; 131 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; } 132 if( fabs(PT_t - Part->PT) > 1e-10 ) { 133 cout << "ERROR-BUG: PT is not conserved in src/BFieldProp.cc. "; 134 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n"; 135 } 136 if(momentum.Pz() != Pz_t) 137 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n"; 138 139 //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E); 140 141 } else { // if B_x or B_y are non zero: longer computation 142 42 143 float Xvertex1 = Part->X; 43 144 float Yvertex1 = Part->Y; 44 145 float Zvertex1 = Part->Z; 45 146 46 //out of tracki bg coverage?47 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > DET->TRACK_radius){return;}48 if(fabs(Zvertex1) > DET->TRACK_length/2.){return;}147 //out of tracking coverage? 148 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;} 149 if(fabs(Zvertex1) > z_max){return;} 49 150 50 float Px = Part->Px; 51 float Py = Part->Py; 52 float Pz = Part->Pz; 53 54 double px = Px / 0.003; 55 double py = Py / 0.003; 56 double pz = Pz / 0.003; 57 double pt = sqrt(px*px+py*py); 151 double px = Part->Px / 0.003; 152 double py = Part->Py / 0.003; 153 double pz = Part->Pz / 0.003; 154 double pt = Part->PT / 0.003; // sqrt(px*px+py*py); 58 155 double p = sqrt(px*px+py*py+pz*pz); 59 156 60 if(q!=0){ 61 double M = Part->M; 62 double vx = px/M; 63 double vy = py/M; 64 double vz = pz/M; 65 66 double Bx = DET->TRACK_bfield_x; 67 double By = DET->TRACK_bfield_y; 68 double Bz = DET->TRACK_bfield_z; 69 70 double ax = (q/M)*(Bz*vy - By*vz); 71 double ay = (q/M)*(Bx*vz - Bz*vx); 72 double az = (q/M)*(By*vx - Bx*vy); 73 double dt = 1/p; 74 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); 75 76 double xold=Xvertex1; double x=xold; 77 double yold=Yvertex1; double y=yold; 78 double zold=Zvertex1; double z=zold; 79 80 double VTold = sqrt(vx*vx+vy*vy); 81 82 int k = 0; 83 84 double Radius=DET->TRACK_radius; 85 double Length=DET->TRACK_length/2.; 86 87 while(k < MAXITERATION){ 157 double M = Part->M; 158 double vx = px/M; 159 double vy = py/M; 160 double vz = pz/M; 161 double qm = q/M; 162 163 double ax = qm*(B_z*vy - B_y*vz); 164 double ay = qm*(B_x*vz - B_z*vx); 165 double az = qm*(B_y*vx - B_x*vy); 166 double dt = 1/p; 167 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ????? 168 169 double xold=Xvertex1; double x=xold; 170 double yold=Yvertex1; double y=yold; 171 double zold=Zvertex1; double z=zold; 172 173 double VTold = pt/M; //=sqrt(vx*vx+vy*vy); 174 175 unsigned int k = 0; 176 double VTratio=0; 177 double R_max2 = R_max*R_max; 178 double r2=0; // will be x*x+y*y 179 180 while(k < MAXITERATION){ 88 181 k++; 89 182 … … 92 185 vz += az*dt; 93 186 94 doubleVTratio = VTold/sqrt(vx*vx+vy*vy);187 VTratio = VTold/sqrt(vx*vx+vy*vy); 95 188 vx *= VTratio; 96 189 vy *= VTratio; 97 190 98 ax = (q/M)*(Bz*vy - By*vz);99 ay = (q/M)*(Bx*vz - Bz*vx);100 az = (q/M)*(By*vx - Bx*vy);191 ax = qm*(B_z*vy - B_y*vz); 192 ay = qm*(B_x*vz - B_z*vx); 193 az = qm*(B_y*vx - B_x*vy); 101 194 102 195 x += vx*dt; 103 196 y += vy*dt; 104 197 z += vz*dt; 105 106 if( (x*x+y*y) > Radius*Radius ){ x /= (x*x+y*y)/(Radius*Radius); y /= (x*x+y*y)/(Radius*Radius); break;} 107 if( fabs(z)>Length)break; 198 r2 = x*x + y*y; 199 200 if( r2 > R_max2 ){ 201 x /= r2/R_max2; 202 y /= r2/R_max2; 203 break; 204 } 205 if( fabs(z)>z_max)break; 108 206 109 207 xold = x; 110 208 yold = y; 111 209 zold = z; 112 } 113 114 if(x!=0 && y!=0 && z!=0) 115 { 116 double eta; 117 float Theta = atan2(sqrt(x*x+y*y),z); 118 eta = -log(tan(Theta/2)); 210 } // while loop 211 212 if(k == MAXITERATION) loop_overflow_counter++; 213 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl; 214 215 if(x!=0 && y!=0 && z!=0) { 216 float Theta = atan2(sqrt(r2),z); 217 double eta = -log(tan(Theta/2.)); 119 218 double phi = atan2(y,x); 120 genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E); 121 } 122 } 219 momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E); 220 } 221 222 } // if b_x or b_y non zero 123 223 }
Note:
See TracChangeset
for help on using the changeset viewer.