Fork me on GitHub

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

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

1 file -> several
File size: 3.2 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 TRACKING_RADIUS = 129;
30 TRACKING_LENGTH = 300;
31 MAXITERATION = 20000;
32 MINSEGLENGTH = 70;
33
34}
35
36void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum)
37{
38
39 double q = Charge(Part->PID);
40 if(q==0)return;
41
42 float Xvertex1 = Part->X;
43 float Yvertex1 = Part->Y;
44 float Zvertex1 = Part->Z;
45
46 //out of trackibg coverage?
47 if(sqrt(pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}
48 if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
49
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(pow(px,2)+pow(py,2));
58
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 }
132 if(x!=0 && y!=0 && z!=0)
133 {
134 double theta = fabs(atan(y/z));
135 double eta;
136 if(z > 0)eta = -log(tan(theta/2));
137 else eta = -(-log(tan(theta/2)));
138 double phi = atan2(y,x);
139 genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
140 }
141 cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
142 }
143}
Note: See TracBrowser for help on using the repository browser.