Fork me on GitHub

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

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

Remove datacard bug + CaloTowers OK

File size: 2.8 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(string DetDatacard) {
28
29 DET = new RESOLution();
30 DET->ReadDataCard(DetDatacard);
31 MAXITERATION = 10000;
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(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > DET->TRACK_radius){return;}
48 if(fabs(Zvertex1) > DET->TRACK_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(px*px+py*py);
58 double p = sqrt(px*px+py*py+pz*pz);
59
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;
86
87 while(k < MAXITERATION){
88 k++;
89
90 vx += ax*dt;
91 vy += ay*dt;
92 vz += az*dt;
93
94 double VTratio = VTold/sqrt(vx*vx+vy*vy);
95 vx *= VTratio;
96 vy *= VTratio;
97
98 ax = (q/M)*(Bz*vy - By*vz);
99 ay = (q/M)*(Bx*vz - Bz*vx);
100 az = (q/M)*(By*vx - Bx*vy);
101
102 x += vx*dt;
103 y += vy*dt;
104 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;
108
109 xold = x;
110 yold = y;
111 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));
119 double phi = atan2(y,x);
120 genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
121 }
122 }
123}
Note: See TracBrowser for help on using the repository browser.