GeneratorInterface/CosmicMuonGenerator/test/e_loss - master - githubmirror/cmssw

GeneratorInterface/CosmicMuonGenerator/test/e_loss
102 строки3 KB
Giulio Eulisse
Snapshot of CMSSW_6_2_0_pre8.
26 июн 2013, 18:12
26 июн 2013, 18:12214ab73
{
  //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;	

}