- Timestamp:
- Nov 27, 2008, 8:49:44 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/BFieldProp.cc
r53 r59 45 45 46 46 //out of trackibg coverage? 47 if(sqrt( pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}47 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACKING_RADIUS){return;} 48 48 if(fabs(Zvertex1) > TRACKING_LENGTH){return;} 49 49 … … 55 55 double py = Py / 0.003; 56 56 double pz = Pz / 0.003; 57 double pt = sqrt(pow(px,2)+pow(py,2)); 57 double pt = sqrt(px*px+py*py); 58 double p = sqrt(px*px+py*py+pz*pz); 59 60 if(q!=0){ 61 double e = Part->E / 0.003; 62 // double M = sqrt(e*e - (px*px + py*py + pz*pz) ); 63 double M = Part->M; 64 /*if(fabs(Part->PID)==11) 65 { 66 cout<<"genMomentum.M() "<<genMomentum.M()<<" "<<Part->M<<endl; 67 cout<<"e*e - (px*px + py*py + pz*pz) "<<e*e - (px*px + py*py + pz*pz)<<endl; 68 }*/ 69 double vx = px/M; 70 double vy = py/M; 71 double vz = pz/M; 58 72 59 double p = sqrt(pow(px,2)+pow(py,2)+pow(pz,2)); 60 61 double m = sqrt(pow(Part->E,2) - pow(Part->Px,2)+ pow(Part->Py,2)+ pow(Part->Pz,2)); 62 63 if(q!=0) 64 { 65 double e = Part->E / 0.003; 66 double M = sqrt(e*e - (px*px + py*py + pz*pz) ); 67 double vx = px/M; 68 double vy = py/M; 69 double vz = pz/M; 70 71 /*double Bx = BField_[0]; 72 double By = BField_[1]; 73 double Bz = BField_[2]; 74 */ 75 76 double Bx = 0; 77 double By = 0; 78 double Bz = 3.8; 79 80 double ax = (q/M)*(Bz*vy - By*vz); 81 double ay = (q/M)*(Bx*vz - Bz*vx); 82 double az = (q/M)*(By*vx - Bx*vy); 83 84 double xold=Xvertex1; 85 double yold=Yvertex1; 86 double zold=Zvertex1; 87 88 double dt = 1/p; 89 90 double x = xold; 91 double y = yold; 92 double z = zold; 93 94 double VTold = sqrt(pow(vx,2)+pow(vy,2)); 95 96 int k = 0; 97 98 while(k < MAXITERATION) 99 { 100 k++; 101 102 // 103 vx += ax*dt; 104 vy += ay*dt; 105 vz += az*dt; 106 107 double VTratio = VTold/sqrt(pow(vx,2)+pow(vy,2)); 108 vx *= VTratio; 109 vy *= VTratio; 110 111 ax = (q/M)*(Bz*vy - By*vz); 112 ay = (q/M)*(Bx*vz - Bz*vx); 113 az = (q/M)*(By*vx - Bx*vy); 114 115 x += vx*dt; 116 y += vy*dt; 117 z += vz*dt; 118 119 if( (x*x+y*y) > (TRACKING_RADIUS*TRACKING_RADIUS) ) 120 { 121 x /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS); 122 y /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS); 123 break; 124 } 125 if( fabs(z) > TRACKING_LENGTH)break; 126 if( fabs((pow(x,2)+pow(y,2)+pow(z,2)) - (pow(xold,2)+pow(yold,2)+pow(zold,2))) < MINSEGLENGTH )continue; 127 128 xold = x; 129 yold = y; 130 zold = z; 131 } 73 double Bx = 0; 74 double By = 0; 75 double Bz = 3.8; 76 77 double ax = (q/M)*(Bz*vy - By*vz); 78 double ay = (q/M)*(Bx*vz - Bz*vx); 79 double az = (q/M)*(By*vx - Bx*vy); 80 /* 81 cout<<"M "<<M<<" et le pid "<<Part->PID<<endl; 82 cout<<"ax "<<ax<<" ay "<<ay<<" az "<<az<<endl; 83 */ 84 double dt = 1/p; 85 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); 86 87 double xold=Xvertex1; double x=xold; 88 double yold=Yvertex1; double y=yold; 89 double zold=Zvertex1; double z=zold; 90 91 double VTold = sqrt(vx*vx+vy*vy); 92 93 int k = 0; 94 95 while(k < MAXITERATION){ 96 k++; 97 98 vx += ax*dt; 99 vy += ay*dt; 100 vz += az*dt; 101 102 double VTratio = VTold/sqrt(vx*vx+vy*vy); 103 vx *= VTratio; 104 vy *= VTratio; 105 106 ax = (q/M)*(Bz*vy - By*vz); 107 ay = (q/M)*(Bx*vz - Bz*vx); 108 az = (q/M)*(By*vx - Bx*vy); 109 110 x += vx*dt; 111 y += vy*dt; 112 z += vz*dt; 113 114 if( (x*x+y*y) > TRACKING_RADIUS*TRACKING_RADIUS ){ x /= (x*x+y*y)/(TRACKING_RADIUS*TRACKING_RADIUS); y /= (x*x+y*y)/(TRACKING_RADIUS*TRACKING_RADIUS); break;} 115 if( fabs(z)>TRACKING_LENGTH)break; 116 117 xold = x; 118 yold = y; 119 zold = z; 120 } 121 132 122 if(x!=0 && y!=0 && z!=0) 133 123 { … … 139 129 genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E); 140 130 } 141 cout<<"genMomentum final "<<genMomentum.Pt()<<endl;131 // cout<<"genMomentum final "<<genMomentum.Pt()<<endl; 142 132 } 143 133 }
Note:
See TracChangeset
for help on using the changeset viewer.