|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
{
//do before in test directory: ln -s ../../../GeneratorInterface GeneratorInterface
cout << "macro e_loss: determine energy loss of muons from surface to detector" << endl;
gROOT->Reset();
gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CosmicMuonGenerator.cc");
gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CMSCGEN.cc");
gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CMSCGENnorm.cc");
gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/SingleParticleEvent.cc");
SingleParticleEvent OneMuoEvt;
OneMuoEvt.PlugVx = 0.;
OneMuoEvt.PlugVz = -14000.;
//OneMuoEvt.PlugVz = 5000.;
double Px_in = 0.;
double Py_in = -30.0;
double Pz_in = 5.;
double Vx_in = 0.;
double Vy_in = 0.;
double Vz_in = 0.;
//determine vertex of muon at Surface (+PlugWidth)
double P_in = sqrt(Px_in*Px_in + Py_in*Py_in + Pz_in*Pz_in);
double dy = Vy_in - (SurfaceOfEarth+PlugWidth);
double Vy_sf = Vy_in - dy;
double Vx_sf = Vx_in - dy*Px_in/Py_in;
double Vz_sf = Vz_in - dy*Pz_in/Py_in;
cout << "Vx_sf=" << Vx_sf << " Vy_sf=" << Vy_sf << " Vz_sf=" << Vz_sf << endl;
double E_in = sqrt(P_in*P_in + MuonMass*MuonMass);
int id=13; //dummy
double T0=0.; //dummy
//waterEquivalents are determined without knowledge of absolute values of Px,y,z's
//just needed for unaltered muon direction
OneMuoEvt.create(id, Px_in, Py_in, Pz_in, E_in, MuonMass, Vx_sf, Vy_sf, Vz_sf, T0);
double ElossScaleFactor = 1.0;
double RadiusOfTarget = 8000.;
double ZDistOfTarget = 150000.;
double ZCentrOfTarget = 0.;
bool TrackerOnly = false;
bool MTCCHalf = false;
OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
double waterEquivalents = OneMuoEvt.WaterEquivalents();
cout << " waterEquivalents=" << waterEquivalents << endl;
OneMuoEvt.setEug(E_in);
//double dEguess = 25.;
//double E_sf = E_in + dEguess;
//try now different initial energies to match the energy of the muon at CMS
//cout << "Check: dE = " << OneMuoEvt.Eloss(waterEquivalents,E_sf) << endl;
//double minEdiff = OneMuoEvt.deltaEmin(E_sf); //E_ug - E_ug_guess
//simple Newton root finder
double x0 = E_in; //minimal possible energy at surface (E_sf) to start with
double eps = 0.001; //accuracy in GeV
double err = 1 + eps;
double dx = 0.01;
double x = x0;
while (err > eps) {
double deriv = (OneMuoEvt.deltaEmin(x+dx)-OneMuoEvt.deltaEmin(x))/dx;
double y = x - OneMuoEvt.deltaEmin(x)/deriv;
err = fabs(y - x);
x = y;
}
double E_sf = x;
cout << "Newton: x=E_sf=" << x << endl;
double minEdiff = OneMuoEvt.deltaEmin(x);
cout << "Newton: minEdiff=" << minEdiff << endl;
double E_ug_guess = E_in - minEdiff;
double dE = E_sf - E_ug_guess;
double P_sf = sqrt(E_sf*E_sf - MuonMass*MuonMass);
double Px_sf = Px_in *P_sf/P_in;
double Py_sf = Py_in *P_sf/P_in;
double Pz_sf = Pz_in *P_sf/P_in;
double Vx_ug = OneMuoEvt.vx();
double Vy_ug = OneMuoEvt.vy();
double Vz_ug = OneMuoEvt.vz();
cout << "E_ug_guess=" << E_ug_guess << " E_sf=" << E_sf << " dE=" << dE << endl;
cout << "E_ug_orig=" << E_in << endl;
cout << "Vx_ug=" << Vx_ug << " Vy_ug=" << Vy_ug
<< " R_ug=" << sqrt(Vx_ug*Vx_ug + Vy_ug*Vy_ug) << " Vz_ug=" << Vz_ug << endl;
}
|