EXPERT ROOT

// Global declarations

TGraph* gr, igr;

// C++ defined function

Double_t f(Double_t *x, Double_t *par)

{

        return igr->Eval(x[0]);

}

ofstream out("CMom_B11_N15.txt");

void test2_p()

{

        // Read the text file into the tree

        // TFile file("myfile.root","RECREATE");

        TTree* MyTree = new TTree("MyTree", "MyTree");

        Double_t tetta, crossSection;

        Double_t CMmom, shorty,phi, px,py,pz,pz_B,px_B,py_B;

        MyTree->ReadFile("cos_tetta_cross.txt", "tetta/D:crossSection/D");

 TBranch *bpt = MyTree->Branch("px",&px,"px/F");

        MyTree->SetBranchAddress("tetta", &tetta);

        MyTree->SetBranchAddress("crossSection", &crossSection);

        /*

        MyTree->SetBranchAddress("pz", &pz);

        MyTree->SetBranchAddress("px", &px);

        MyTree->SetBranchAddress("py", &py);

        MyTree->SetBranchAddress("pz_B", &pz_B);

        MyTree->SetBranchAddress("px_B", &px_B);

        MyTree->SetBranchAddress("py_B", &py_B);

        */

        // Create two vectors to be filled from the tree

        const Int_t nEntries = MyTree->GetEntries(); // number of entries in the input file

        TVectorD tet(nEntries);

        TVectorD sigma(nEntries);

        // Fill two vectors from the tree

        for (int i=0; i<nEntries; i++) {

                MyTree->GetEntry(i);

                tet(i) = tetta;

                sigma(i) = crossSection;

        }

        // Create the graph

        gr = new TGraph(tet, sigma);

        igr = new TGraph(sigma,tet);

        // Draw the graph

        TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);

        c1->cd();

        igr->Draw();

        igr->SetLineColor(kRed);

        igr->SetLineWidth(3);

        gPad->SetGrid(1, 1);

        // ROOT-TF1 object with the reference to the C++ defined function

        TF1* func = new TF1("func", f, 0, 1, 0);

        // Draw the func

        TCanvas *c2 = new TCanvas("c2", "c2", 800, 600);

        c2->cd();

        func->Draw();

        func->SetLineColor(kGreen);

        func->SetLineWidth(3);

        gPad->SetGrid(1, 1);

        // Uniform random number generator

        TRandom3* rndGen = new TRandom3();

        TH1D* histo1 = new TH1D("histo1", "histo1", 1800, 0., 1.);

        TH1D* histo2 = new TH1D("histo2", "histo2", 1800, 0., 180.);

        TH1D* histo3 = new TH1D("histo3", "histo3", 1800, 0., 180.);

        TH1D* histo4 = new TH1D("histo4", "histo4", 1800, 0., 180.);

        TH1D* histo5 = new TH1D("histo5", "histo5", 1800, 0., 180.);

        UInt_t nIterations = 500000;

        // Loop over all iterations

        for (UInt_t iIter=0; iIter<nIterations; iIter++)

        {

                if (iIter%100000 == 0) cout << "Iteration " << iIter << endl;

                Double_t x = rndGen->Rndm();        // from 0.0 to 1.0

                Double_t y = func->Eval(x);

   //___________fill px py pz -----------------

Double_t         invariant;

                Double_t AMU=931.5 ;

                Double_t m11B=AMU*11.;

                Double_t m15N=AMU*15;

                Double_t T_N15=43;

Double_t CMmom, shorty,phi, px,py,pz;

                invariant = (m11B+m15N)*(m11B+m15N)+2*m11B*T_N15;

Double_t shorty=(invariant-m11B*m11B-m15N*m15N)*(invariant-m11B*m11B-m15N*m15N);

CMmom = sqrt((shorty-4*m15N*m15N*m11B*m11B)/(4*invariant));

                //phi= CLHEP::RandFlat::shoot(0.0,2*double(CLHEP::pi));

                Double_t PI = 3.14159265358;

                phi = rand()% 2*PI;

                pz=CMmom*cos(y);

                px=CMmom*sin(y)*cos(phi*180/PI) ;//cos(phi);

                py=CMmom*sin(y)*sin(phi*180/PI); //sin(phi);

                pz_B=-pz;

                px_B=-px;

                py_B=-py;

                out<<pz<<"\t"<<px<<"\t"<< py<<"\t"<<-pz<<"\t"<<-px<<"\t"<< -py<<endl;

                // Histograms filling

                histo1->Fill(x);

                histo2->Fill(y);

                histo3->Fill(pz);

                histo4->Fill(px);

                histo5->Fill(py);

                //file.Write();

        }

        // Draw

        TCanvas* canv1 = new TCanvas("canv1", "canv1");

        canv1->cd();

        histo1->Draw();

        gStyle->SetOptStat(111111);

        gPad->SetGrid(1, 1);

        // Draw

        TCanvas* canv2 = new TCanvas("canv2", "canv2");

        canv2->cd();

        histo2->Draw();

        gStyle->SetOptStat(111111);

        gPad->SetGrid(1, 1);

                TCanvas* canv3 = new TCanvas("canv3", "canv3");

        canv3->cd();

        histo3->Draw();

        gStyle->SetOptStat(111111);

        gPad->SetGrid(1, 1);

                TCanvas* canv4 = new TCanvas("canv4", "canv4");

        canv3->cd();

        histo4->Draw();

        gStyle->SetOptStat(111111);

        gPad->SetGrid(1, 1);

                TCanvas* canv5 = new TCanvas("canv5", "canv5");

        canv3->cd();

        histo5->Draw();

        gStyle->SetOptStat(111111);

        gPad->SetGrid(1, 1);

}