Warning, /GeneratorInterface/CosmicMuonGenerator/test/e_loss is written in an unsupported language. File is not indexed.
0001 {
0002 //do before in test directory: ln -s ../../../GeneratorInterface GeneratorInterface
0003 cout << "macro e_loss: determine energy loss of muons from surface to detector" << endl;
0004 gROOT->Reset();
0005 gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CosmicMuonGenerator.cc");
0006 gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CMSCGEN.cc");
0007 gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/CMSCGENnorm.cc");
0008 gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/SingleParticleEvent.cc");
0009
0010
0011
0012 SingleParticleEvent OneMuoEvt;
0013 OneMuoEvt.PlugVx = 0.;
0014 OneMuoEvt.PlugVz = -14000.;
0015 //OneMuoEvt.PlugVz = 5000.;
0016
0017 double Px_in = 0.;
0018 double Py_in = -30.0;
0019 double Pz_in = 5.;
0020 double Vx_in = 0.;
0021 double Vy_in = 0.;
0022 double Vz_in = 0.;
0023
0024
0025 //determine vertex of muon at Surface (+PlugWidth)
0026 double P_in = sqrt(Px_in*Px_in + Py_in*Py_in + Pz_in*Pz_in);
0027 double dy = Vy_in - (SurfaceOfEarth+PlugWidth);
0028 double Vy_sf = Vy_in - dy;
0029 double Vx_sf = Vx_in - dy*Px_in/Py_in;
0030 double Vz_sf = Vz_in - dy*Pz_in/Py_in;
0031 cout << "Vx_sf=" << Vx_sf << " Vy_sf=" << Vy_sf << " Vz_sf=" << Vz_sf << endl;
0032
0033 double E_in = sqrt(P_in*P_in + MuonMass*MuonMass);
0034
0035
0036 int id=13; //dummy
0037 double T0=0.; //dummy
0038 //waterEquivalents are determined without knowledge of absolute values of Px,y,z's
0039 //just needed for unaltered muon direction
0040 OneMuoEvt.create(id, Px_in, Py_in, Pz_in, E_in, MuonMass, Vx_sf, Vy_sf, Vz_sf, T0);
0041
0042 double ElossScaleFactor = 1.0;
0043 double RadiusOfTarget = 8000.;
0044 double ZDistOfTarget = 150000.;
0045 double ZCentrOfTarget = 0.;
0046 bool TrackerOnly = false;
0047 bool MTCCHalf = false;
0048
0049 OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
0050 double waterEquivalents = OneMuoEvt.WaterEquivalents();
0051 cout << " waterEquivalents=" << waterEquivalents << endl;
0052
0053
0054 OneMuoEvt.setEug(E_in);
0055
0056 //double dEguess = 25.;
0057 //double E_sf = E_in + dEguess;
0058
0059 //try now different initial energies to match the energy of the muon at CMS
0060 //cout << "Check: dE = " << OneMuoEvt.Eloss(waterEquivalents,E_sf) << endl;
0061 //double minEdiff = OneMuoEvt.deltaEmin(E_sf); //E_ug - E_ug_guess
0062
0063
0064 //simple Newton root finder
0065 double x0 = E_in; //minimal possible energy at surface (E_sf) to start with
0066 double eps = 0.001; //accuracy in GeV
0067 double err = 1 + eps;
0068 double dx = 0.01;
0069 double x = x0;
0070 while (err > eps) {
0071 double deriv = (OneMuoEvt.deltaEmin(x+dx)-OneMuoEvt.deltaEmin(x))/dx;
0072 double y = x - OneMuoEvt.deltaEmin(x)/deriv;
0073 err = fabs(y - x);
0074 x = y;
0075 }
0076
0077 double E_sf = x;
0078 cout << "Newton: x=E_sf=" << x << endl;
0079 double minEdiff = OneMuoEvt.deltaEmin(x);
0080 cout << "Newton: minEdiff=" << minEdiff << endl;
0081
0082
0083 double E_ug_guess = E_in - minEdiff;
0084 double dE = E_sf - E_ug_guess;
0085
0086 double P_sf = sqrt(E_sf*E_sf - MuonMass*MuonMass);
0087 double Px_sf = Px_in *P_sf/P_in;
0088 double Py_sf = Py_in *P_sf/P_in;
0089 double Pz_sf = Pz_in *P_sf/P_in;
0090
0091 double Vx_ug = OneMuoEvt.vx();
0092 double Vy_ug = OneMuoEvt.vy();
0093 double Vz_ug = OneMuoEvt.vz();
0094
0095
0096
0097 cout << "E_ug_guess=" << E_ug_guess << " E_sf=" << E_sf << " dE=" << dE << endl;
0098 cout << "E_ug_orig=" << E_in << endl;
0099 cout << "Vx_ug=" << Vx_ug << " Vy_ug=" << Vy_ug
0100 << " R_ug=" << sqrt(Vx_ug*Vx_ug + Vy_ug*Vy_ug) << " Vz_ug=" << Vz_ug << endl;
0101
0102 }
