EXP1803

#include "Particle.h"

#include <sstream>

#include <string>

/////////////////////////////////////////////////////////

//

// Particle

//

// This class enables you to set

// a mass and an impulse of a particle,

// change values and get them back.

// fA and fZ variables are still ambiguous

// and need to be discussed.

// Each particle also stores information about its energy states

//

////////////////////////////////////////////////////////

ExcitationState::ExcitationState() {

        //Default constructor

        fMean = 0.;

        fWidth = 0.1;

        fShape = "gauss";

        fStrength = 1;

}

;

ExcitationState::ExcitationState(Double_t mean, Double_t width, TString shape,

                Int_t strength) :

                fMean(mean), fWidth(width), fShape(shape), fStrength(strength) {

        //Constructor initializing all necessary variables for exc. state

}

TString ExcitationState::CreateConfigString() {

        ConfigDictionary CD;

        CD.SetDouble("mean", fMean);

        CD.SetDouble("width", fWidth);

        CD.SetString("shape", fShape.Data());

        CD.SetInt("strength", fStrength);

        return TString(CD.ToString().c_str());

}

;

int ExcitationState::ReadConfigString(TString cs) {

        ConfigDictionary CD(cs.Data());

        try {

                fMean = CD.GetDouble("mean");

                fWidth = CD.GetDouble("width");

                fShape = CD.GetString("shape").c_str();

                fStrength = CD.GetInt("strength");

                return SUCCESS;

        } catch (std::string & e) {

                Error("ExcitationState::ReadConfigString",

                                "Couldn't find parameter: %s", e.c_str());

                return NOTFOUND;

        }

        if (fWidth <= 0.) {

                Error("ExcitationState::ReadConfigString", "Width parameter <= 0!");

                return FAILURE;

        }

        if (fMean <= 0.) {

                Error("ExcitationState::ReadConfigString", "Mean parameter <= 0!");

                return FAILURE;

        }

        if (fStrength <= 0) {

                Error("ExcitationState::ReadConfigString", "Strength parameter <= 0!");

                return FAILURE;

        }

        return SUCCESS;

}

;

/////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////////////

ClassImp(Particle);

Int_t Particle::numberOfParticles = 0;

Particle::Particle() {

// default constructor

// set Z = 0, A = 0

        Reset();

        fZ = 0;

        fA = 0;

        numberOfParticles++;

        std::stringstream ss;

        ss << "unnamed particle no. ";

        ss << numberOfParticles;

        SetName(ss.str().c_str());

//        cout << "kjadbhfkjasdf" << endl;

}

;

Particle::Particle(const char *name, Double_t mass, Int_t A, Int_t Z,

                Bool_t obs) {

        //name: name of the object

        //mass: rest-mass of the object, in MeV

        //A: mass number, number of nucleons

        //Z: charge, equals number in periodic table

        //set total energy equal to the mass of the particle

        fName = name;

        fMass = mass;

        fGroundStateMass = mass;

        Info("Particle::Particle", "Ground state mass was set to %f",

                        fGroundStateMass);

        fA = A;

        fZ = Z;

        fObservable = obs;

        fImpulse.SetE(mass);

}

Particle::Particle(TString cs) {

//config constructor

        ReadConfigString(cs);

        CreateStateMassFunctions();

}

Particle::~Particle() {

//destructor

        for (UInt_t i = 0; i < fState.size(); i++) {

                delete fState[i];

//                Info("DetManager::~DetManager", "Detector %s deleting", det->GetName());

//                delete det;

        }

        fState.clear();

}

//public functions

//_____________________________________________________________________________

void Particle::AddExcitationState(ExcitationState exstate) {

//add new excitation states to the particle

        fExcitationStates.push_back(exstate);

}

//_____________________________________________________________________________

void Particle::ClearExcitationStates() {

//clear excitation states

        fExcitationStates.clear();

}

//_____________________________________________________________________________

Int_t Particle::GetNumberOfExStates() {

//return number of excitation states

        return fExcitationStates.size();

}

//_____________________________________________________________________________

ExcitationState Particle::GetExcitationState(Int_t index) {

//return "index"th excitation state

        return fExcitationStates[index];

}

void Particle::Reset() {

        // Set mass and impulse of particle as zero value

        fMass = 0.;

        fGroundStateMass = 0.;

        fImpulse.SetPxPyPzE(0., 0., 0., 0.);

        return;

}

int Particle::CopyValues(Particle * other) {

        //Copies kinetic energy and direction angles from other particle

        //Dont check if other == null to increase speed, be careful!

        fImpulse.SetPxPyPzE(other->GetPx(), other->GetPy(), other->GetPz(),

                        other->GetE());

        return SUCCESS;

}

void Particle::SetMPxPyPz(Double_t m, Double_t px, Double_t py, Double_t pz) {

        // m: set mass in the MeV

        // px: x cordinate of impulse

        // py: y cordinate of impulse

        // pz: z cordinate of impulse

        fImpulse.SetPx(px);

        fImpulse.SetPy(py);

        fImpulse.SetPz(pz);

        fImpulse.SetE(CalcT(m, px, py, pz) + m);

        fMass = m;

}

void Particle::SetMTNxNyNz(Double_t m, Double_t T, Double_t nx, Double_t ny,

                Double_t nz) {

// m: set mass in MeV

// T: set kinetic energy in MeV

// set direction vector of impulse with (nx,ny,nz) components

//todo: function restored, check if it works properly

        fMass = m;

        TVector3 vect(nx, ny, nz);

        vect.SetMag(CalcP(m, T));

        fImpulse.SetVectMag(vect, m);

        fImpulse.SetE(Sqrt(vect.Mag2() + m * m));

        return;

}

;

void Particle::SetMTDir(Double_t m, Double_t T, TVector3 dir) {

        // dir: set direction vector of impulse

        // m: set mass in MeV

        // T: set kinetic energy in MeV

        // imp: need to use TVector3 dir(dir(0),dir(1),dir(2))

        // where dir(0), dir(1), dir(2) are coordinates of direction vector of impulse

        Double_t p = CalcP(m, T); // impulse magnitude

        fMass = m;

        fImpulse.SetE(T + m);

        fImpulse.SetPx(dir.Unit().x() * p);

        fImpulse.SetPy(dir.Unit().y() * p);

        fImpulse.SetPz(dir.Unit().z() * p);

}

void Particle::SetP(TVector3 P) {

        // set impulse of the particle

        fImpulse.SetXYZM(P(0), P(1), P(2), fMass);

}

void Particle::SetPx(Double_t px) {

// set cartesian x coordinate of impulse in MeV

        fImpulse.SetPx(px);

        //fImpulse.SetE(0.);

        fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));

        return;

}

;

void Particle::SetPy(Double_t py) {

// set cartesian y coordinate of impulse in MeV

        fImpulse.SetPy(py);

        //fImpulse.SetE(0.);

        fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));

}

;

void Particle::SetPz(Double_t pz) {

// set cartesian z coordinate of impulse in MeV

        fImpulse.SetPz(pz);

        fImpulse.SetE(0.);

        fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));

}

;

void Particle::SetE(Double_t E) {

        // set total energy in MeV

        fImpulse.SetE(E);

}

;

void Particle::SetObservable(Bool_t obs) {

//whether object is observable or not

        fObservable = obs;

}

;

void Particle::SetMass(Double_t mass) {

        // Set particle mass.

        // Impulse and Kinetic energy are remained

        Double_t T = GetT();

        fMass = mass;

        fImpulse.SetE(T + mass);

}

void Particle::SetT(Double_t T) {

//set kinetic energy in MeV

//NB! impulse should be nonzero else you'll receive message "zero vector can't be stretched"

        if (T <= 0)

                fImpulse = TLorentzVector(0, 0, 0, fMass);

        else {

                fImpulse.SetRho(CalcP(fMass, T));

                fImpulse.SetE(fMass + T);

        }

}

void Particle::SetImpulse(TLorentzVector *P) {

        //set TLorenzVector

                fImpulse.SetPx(P->Px()*1000.);

                fImpulse.SetPy(P->Py()*1000.);

                fImpulse.SetPz(P->Pz()*1000.);

                fImpulse.SetE(P->E()*1000.);

}

void Particle::SetTPhiTheta(Double_t T, Double_t phi, Double_t theta) {

//change kinetic energy and direction

//T: kinetic energy in MeV

//phi: azimuthal angle of impulse in rad

//theta: polar angle of impulse in rad

//NB! impulse should be nonzero else you'll receive message "zero vector can't be stretched"

        fImpulse.SetRho(CalcP(fMass, T));

        fImpulse.SetE(fMass + T);

        fImpulse.SetPhi(phi);

        fImpulse.SetTheta(theta);

}

;

void Particle::Print(Option_t * option){

        // Print values:

        //                particle name, whether particle is observable,

        //                mass number and charge (in units of electron charge) of the particle,

        //                rest-mass and total energy in MeV, kinetic energy in MeV,

        //                coordinates of impulse and value of impulse in MeV,

        //                phi and theta angles in rad

        if (option != NULL) {

                //just checking

        }

        cout << "Name: " << fName;

        if (fObservable == 0)

                cout << " is not observable" << endl;

        else

                cout << " is observable" << endl;

        cout << "A, Z = " << fA << ",  " << fZ << endl << "Mass = " << fMass << endl

                        << "Energy = " << fImpulse.E() << endl << "Kinetic Energy = "

                        << GetT() << endl << "Impulse (Px;Py;Pz) = " << "(" << fImpulse.Px()

                        << " " << fImpulse.Py() << " " << fImpulse.Pz() << ")" << endl

                        << "P = " << fImpulse.Rho() << endl << "Phi = " << fImpulse.Phi()

                        << endl << "Theta = " << fImpulse.Theta() << endl;

}

;

TString Particle::CreateConfigString() {

        //Creates string containing config info about particle

        //so it can be saved to file or used anywhere else (gui?)

        //It saves direction of impulse as phi and theta angles

        //Its format looks like: (for example proton) #mass is stored in MeV I presume?

        //"name"="proton" "mass"="938" "observable"="true" "a_number"="1" "z_number"="1" "phi"="0.0" "theta"="0.0"

        ConfigDictionary CD;

        CD.SetString("name", GetName());

        CD.SetDouble("mass", GetM());

        CD.SetBool("observable", fObservable);

        CD.SetInt("a_number", GetA());

        CD.SetInt("z_number", GetZ());

        CD.SetDouble("phi", GetPhi());

        CD.SetDouble("theta", GetTheta());

        CD.SetDouble("k_energy", GetT()); //kinetic energy = GetT()

        CD.SetInt("exNumber", GetNumberOfExStates());

        TString ret(CD.ToString());

        return ret;

}

void Particle::ReadConfigString(TString conf) {

        //conf - string formatted as list of "key=value" pairs

        //Reads and sets parameters from config string

        //Parameters for impulse direction are phi and theta angles

        ConfigDictionary CD(conf.Data());

        try {

                SetName(CD.GetString("name").c_str());

                SetMass(CD.GetDouble("mass"));

                fGroundStateMass = GetM();

                Info("Particle::ReadConfigString", "Ground state mass was set to %f",

                                fGroundStateMass);

                SetObservable(CD.GetBool("observable"));

                SetA(CD.GetInt("a_number"));

                SetZ(CD.GetInt("z_number"));

                SetPz(1.);

                SetTPhiTheta(CD.GetDouble("k_energy"), CD.GetDouble("phi"),

                                CD.GetDouble("theta"));

        } catch (...) {

                Error("Particle::ReadConfigString", "Couldn't find all parameters!");

                return;

        }

}

TVector3 Particle::GetBoost() {

        // return vector beta (impulse components divided by the time component)

        return fImpulse.BoostVector();

}

;

void Particle::BoostTransform(TVector3 beta) {

//perform a boost transformation

//to correct working,we input (-beta), see TLorentzVector documentation

//this bug will be corrected soon

        fImpulse.Boost(-beta);

        return;

}

;

void Particle::CreateStatesWeigthsHist() {

        fStatesWeigths.SetBins(GetNumberOfExStates(), 0, GetNumberOfExStates());

//        fStatesWeigths.SetBins(8, 1, 8+1);

        for (Int_t i = 0; i < GetNumberOfExStates(); i++) {

//                cout << "\t\tjhagvdjhavsdasvd" << endl;

//                cout << i << endl;

//                cout << GetExcitationState(i).GetStrength() << endl;

//                cout << fStatesWeigths.GetNbinsX() << endl;

                fStatesWeigths.AddBinContent(i + 1,

                                GetExcitationState(i).GetStrength());

        }

}

void Particle::DrawWeigths() {

        fStatesWeigths.Draw();

        return;

}

void Particle::DrawMassDistribution(Int_t i, Option_t *option) {

        //"i" is number of excitation state

        //"option" is draw option in TF1

        if (fState.empty()) {

                Error("Particle::DrawMassDistribution",

                                "There are no defined states for particle %s", GetName());

                return;

        }

        if (i > (Int_t) fState.size()) {

                Error("Particle::DrawMassDistribution", "Maximum possible index is %d",

                                (Int_t) (fState.size()));

                return;

        }

        fState[i]->Draw(option);

        return;

}

void Particle::CreateStateMassFunctions() {

        //method creates TF1 function for each state as defined e.g.

        //in configuration file

        //

        //Two distributions are realized for the present in the code:

        //Gauss and Lorentz distributions.

        //In future releases other distributions will be added.

        //

        //In order to optimize the calculation process there exist

        //limits of distribution:

        //(-3*sigma ; 4*sigma) for Gauss;

        //(-3*FWHM ; 4*FWHM) for Lorentz.

        //This limits provide 98% of integral area.

        CreateStatesWeigthsHist();

        TF1 *massSpectrum;

        massSpectrum = NULL;

        Double_t xmin = 0.; //lower and

        Double_t xmax = 0.; //upper range for each state

        for (Int_t i = 0; i < GetNumberOfExStates(); i++) {

                // cout << GetNumberOfExStates()<< " WTF!!!!! " << GetName() << endl; 

                //gauss distribution of the excited state

                if (GetExcitationState(i).GetShape() == "gauss") {

                        const Double_t mean = GetExcitationState(i).GetMean();

                        const Double_t sigma = GetExcitationState(i).GetWidth() / 2.35;

                        Info(

                                        "Particle::CreateStateMassFunctions",

                                        "Creating gauss distribution with mean %f and FWHM of %f MeV",

                                        mean, sigma);

//                        const Double_t mean = GetExcitationState(i).GetMean();

//                        const Double_t sigma = GetExcitationState(i).GetWidth()/2.35;

                        xmin = mean - 4 * sigma; // important cuz for high mean and little sigma it doesnt work cuz of fNpx (see TRandom)

                        if (mean - 3 * sigma < 0. || sigma < 0.) {

                                Warning("Particle::CreateStateMassFunctions",

                                                "There is something weird in this state");

                                xmin = mean - 3 * sigma;

                        } //if

                        xmax = mean + 4 * sigma;

                        Info("Particle::CreateStateMassFunctions",

                                        "Range for this state was set from %f to %f", xmin, xmax);

                        //        (-3*sigma, +4*sigma)

                        massSpectrum = new TF1("massGauss", "TMath::Gaus(x, [0], [1])",

                                        xmin, xmax);

                        massSpectrum->SetParameters(mean, sigma); //1st parameter is mean value

                                                                                                          //2nd parameter is sigma

                } //if gauss

                if (GetExcitationState(i).GetShape() == "lorentzian") {

                        const Double_t mean = GetExcitationState(i).GetMean();

                        const Double_t fwhm = GetExcitationState(i).GetWidth();

                        Info(

                                        "Particle::CreateStateMassFunctions",

                                        "Creating lorentz distribution with mean %f and FWHM of %f MeV",

                                        mean, fwhm);

                        //                        const Double_t mean = GetExcitationState(i).GetMean();

                        //                        const Double_t sigma = GetExcitationState(i).GetWidth()/2.35;

                        if (-3 * fwhm < 0.) {

                                Warning("Particle::CreateStateMassFunctions",

                                                "There is something weird in this state");

                                xmin = mean - 3 * fwhm;

                        } //if

                        xmax = mean + 4 * fwhm;

                        Info("Particle::CreateStateMassFunctions",

                                        "Range for this state was set from %f to %f", xmin, xmax);

                        //        (-3*sigma, +4*sigma)

                        massSpectrum = new TF1("massLorentz",

                                        "TMath::BreitWigner(x, [0], [1])", xmin, xmax);

                        massSpectrum->SetParameters(mean, fwhm); //1st parameter is mean value

                } //if lorentz

                if (massSpectrum) {

                        fState.push_back(massSpectrum);

                        Info("Particle::CreateStateMassFunctions",

                                        "Function %s was created", fState[i]->GetName());

                } else {

                        Warning("Particle::CreateStateMassFunctions",

                                        "Function for e.s. number %d of type %s was not created", i,

                                        GetExcitationState(i).GetShape().Data());

                }

        }

        return;

}

void Particle::GenerateMass() {

        //calculate and set mass according to mass distribution functions

        if (!GetNumberOfExStates())

                return;

        //choose state:

        Int_t state = (Int_t) fStatesWeigths.GetRandom();

        //fixme Vratislav: random mass

        //for a few reaction in chain the same random numbers are obtained in each event

        //same situation as for angles

        Double_t mass = fGroundStateMass + fState[state]->GetRandom();

        SetMass(mass);

        return;

}

Bool_t Particle::IsObservable() {

//whether particle is observable (1) or not (0)

        return fObservable;

        /*

         if(fObservable==0)

         {

         cout << fName << " is not observable" << endl;

         return fObservable;

         }

         else

         {

         cout << fName << " is observable" << endl;

         return fObservable;

         }*/

}

//private functions

Double_t Particle::CalcT(Double_t m, Double_t px, Double_t py, Double_t pz) {

// Calculate kinetic energy

        TVector3 p(px, py, pz);

        return (Sqrt(p.Mag2() + m * m) - m);

}

;

Double_t Particle::CalcP(Double_t m, Double_t T) {

// Calculate value of impulse

        return (Sqrt(T * T + 2 * T * m));

}

;

Double_t Particle::CalcE() {

        return GetE();

}

;

Bool_t Particle::CheckEnergyConservation() {

// Check energy conservation

        Double_t ls = Power(fImpulse.Energy(), 2);

        Double_t ps = Power(fImpulse.Rho(), 2) + fMass * fMass;

        printf("%2.18f\t\n", ls);

        printf("%2.18f\t\n", ps);

        printf("%2.18f\t\n", TMath::Abs(ls / ps - 1.0));

        printf("%2.18f\t\n", 1.0e-15);

        printf("%f\t%f\n", Power(fImpulse.Rho(), 2), fMass * fMass);

        if (TMath::Abs(ls / ps - 1.0) < 1.0e-14) {

                cout << "Energy is conserved" << endl;

                return true;

        } else {

                cout << "Energy's not conserved" << endl;

                return false;

        }

}