//Author: Alex Sage
//June 2013
//Outline:
//  need to figure out how to get cherenkov radiation (angle, intensity)
//  need to figure out how to get geometry of lqbar
//Progress:
//  Can do basic arithmetic on LVectors and update them 
//  still don't understand what include statements are necessary
//Problems:
//  no idea how to get formula for any of this stuff
//     can't find formula for wavength based on speed of proton
//     n is based on lambda
//     there is a formula for dE but it is based on n (intensity)
//     angle given off is based on Beta (don't know how to use yet)
//  I couldn't find a funtion in root to simply call 
//  there has got to be an easier way

#include <iostream>
#include "TMath.h"
#include "Math/Vector4D.h"
#include "Math/Vector3D.h"
#include "TLorentzVector.h"

using namespace ROOT::Math;


void mytry()
{
//after talking to rasmus lets try the TLorentzVectors again
//so v1 can be for position. and v2 can be for velocity
TLorentzVector v1(0.0, 2., 0.0, 2.);
TLorentzVector v2(0.0, 4., 0.0, 1.);
Double_t s;  //mag of vel
Double_t n; //index of refraction
Double_t lambda; // wavelength
Double_t bet;

lambda = 50.0 * pow(10., -9.);     //picked 50nm wavelgnth out of guess since cherenkov 
                             //radiation is usually in hard UV
  //dispertion fomula is used for n based on fused silica in a range of 
  //lambda = .21 - 3.71 um; formula also under assumption of room temperature 
n = sqrt( 1 + 0.6961663*pow(lambda,2.)/(pow(lambda,2.)-pow(0.0684043,2.)) + 0.4079426*pow(lambda,2.)/(pow(lambda,2.)-pow(0.1162414,2.)) + 0.8974794*pow(lambda,2.)/(pow(lambda,2.)-pow(9.896161,2.)) );

printf("\nInitialize the vector");
printf("\nLorentzVector = (%.2f, %.2f, %.2f, %.2f) \n", v1.X(), v1.Y(), v1.Z(), v1.T());
cout << "\nafter arithmetic" << endl;
TVector3 p = v2.Vect(); s = p.Dot(p);
printf("Tryting to get the mag of velocity %.2f \n", s);
bet = v1.Beta();
printf("Beta? = %.2f \n", bet);

/*
while (v1.Y() < 50)
{
v1 += v2;
printf("LorentzVector = (%.2f, %.2f, %.2f, %.2f) \n", v1.X(), v1.Y(), v1.Z(), v1.T());
}
*/
}