Fork me on GitHub

source: svn/trunk/src/BFieldProp.cc@ 82

Last change on this file since 82 was 62, checked in by severine ovyn, 16 years ago

Read BField info from datacard

File size: 3.1 KB
Line 
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 "TRandom.h"
13
14#include <iostream>
15#include <sstream>
16#include <fstream>
17#include <iomanip>
18
19#include<cmath>
20
21
22using namespace std;
23
24
25//------------------------------------------------------------------------------
26
27TrackPropagation::TrackPropagation() {
28
29 MAXITERATION = 20000;
30 MINSEGLENGTH = 70;
31
32}
33
34void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum)
35{
36
37 double q = Charge(Part->PID);
38 if(q==0)return;
39
40 float Xvertex1 = Part->X;
41 float Yvertex1 = Part->Y;
42 float Zvertex1 = Part->Z;
43
44 //out of trackibg coverage?
45 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACKING_RADIUS){return;}
46 if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
47
48 float Px = Part->Px;
49 float Py = Part->Py;
50 float Pz = Part->Pz;
51
52 double px = Px / 0.003;
53 double py = Py / 0.003;
54 double pz = Pz / 0.003;
55 double pt = sqrt(px*px+py*py);
56 double p = sqrt(px*px+py*py+pz*pz);
57
58 if(q!=0){
59 double e = Part->E / 0.003;
60// double M = sqrt(e*e - (px*px + py*py + pz*pz) );
61 double M = Part->M;
62/*if(fabs(Part->PID)==11)
63{
64 cout<<"genMomentum.M() "<<genMomentum.M()<<" "<<Part->M<<endl;
65 cout<<"e*e - (px*px + py*py + pz*pz) "<<e*e - (px*px + py*py + pz*pz)<<endl;
66}*/
67 double vx = px/M;
68 double vy = py/M;
69 double vz = pz/M;
70
71 double Bx = BFIELD_X;
72 double By = BFIELD_Y;
73 double Bz = BFIELD_Z;
74
75 double ax = (q/M)*(Bz*vy - By*vz);
76 double ay = (q/M)*(Bx*vz - Bz*vx);
77 double az = (q/M)*(By*vx - Bx*vy);
78/*
79cout<<"M "<<M<<" et le pid "<<Part->PID<<endl;
80cout<<"ax "<<ax<<" ay "<<ay<<" az "<<az<<endl;
81*/
82 double dt = 1/p;
83 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz);
84
85 double xold=Xvertex1; double x=xold;
86 double yold=Yvertex1; double y=yold;
87 double zold=Zvertex1; double z=zold;
88
89 double VTold = sqrt(vx*vx+vy*vy);
90
91 int k = 0;
92
93 while(k < MAXITERATION){
94 k++;
95
96 vx += ax*dt;
97 vy += ay*dt;
98 vz += az*dt;
99
100 double VTratio = VTold/sqrt(vx*vx+vy*vy);
101 vx *= VTratio;
102 vy *= VTratio;
103
104 ax = (q/M)*(Bz*vy - By*vz);
105 ay = (q/M)*(Bx*vz - Bz*vx);
106 az = (q/M)*(By*vx - Bx*vy);
107
108 x += vx*dt;
109 y += vy*dt;
110 z += vz*dt;
111
112 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;}
113 if( fabs(z)>TRACKING_LENGTH)break;
114
115 xold = x;
116 yold = y;
117 zold = z;
118 }
119
120 if(x!=0 && y!=0 && z!=0)
121 {
122 double theta = fabs(atan(y/z));
123 double eta;
124 if(z > 0)eta = -log(tan(theta/2));
125 else eta = -(-log(tan(theta/2)));
126 double phi = atan2(y,x);
127 genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
128 }
129 // cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
130 }
131}
Note: See TracBrowser for help on using the repository browser.