EXPERT ROOT

//#include "FairBoxGenerator.h"

//#include "FairIonGenerator.h"

using std::cout;

using std::endl;

void RTelescope_sim_test(Int_t nEvents = 10){

        //---------------------Files-----------------------------------------------

        TString outFile= "sim.root";

        TString parFile= "par.root";

        // ------------------------------------------------------------------------

        // -----   Timer   --------------------------------------------------------

        TStopwatch timer;

        timer.Start();

        // ------------------------------------------------------------------------

        // -----   Create simulation run   ----------------------------------------

        FairRunSim* run = new FairRunSim();

        /** Select transport engine

         * TGeant3

         * TGeant4

         **/

        run->SetName("TGeant4");              // Transport engine

        run->SetOutputFile(outFile.Data());          // Output file

        // ------------------------------------------------------------------------

        // -----   Runtime database   ---------------------------------------------

        FairRuntimeDb* rtdb = run->GetRuntimeDb();

        // ------------------------------------------------------------------------

        // -----   Create media   -------------------------------------------------

        run->SetMaterials("media.geo");       // Materials

        // ------------------------------------------------------------------------

        // -----   Create detectors  ----------------------------------------------

        FairModule* cave= new ERCave("CAVE");

        cave->SetGeometryFileName("cave.geo");

        run->AddModule(cave);

        // ER NeuRad definition

        /* Select verbosity level

         * 1 - only standard logs

         * 2 - Print points after each event

         * 3 - - GEANT Step information

         */

        Int_t verbose = 1;

        ERRTelescope* RTelescope= new ERRTelescope("ERRTelescope", kTRUE,verbose);

        RTelescope->SetGeometryFileName("RTelescope.v3.geo.root");

        run->AddModule(RTelescope);

        // ------------------------------------------------------------------------

        // -----   Create PrimaryGenerator   --------------------------------------

        FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

        //proton

        Int_t pdgId = 2212; // proton  beam

        Double32_t theta1 = 0;  // polar angle distribution

        Double32_t theta2 = 0;

        Double32_t kin_energy = 0.01; //GeV

        Double_t mass = TDatabasePDG::Instance()->GetParticle(pdgId)->Mass();

        Double32_t momentum = TMath::Sqrt(kin_energy*kin_energy + 2.*kin_energy*mass); //GeV

        FairBoxGenerator* boxGen = new FairBoxGenerator(pdgId, 1);

        boxGen->SetThetaRange(theta1, theta2);

        boxGen->SetPRange(momentum, momentum);

        boxGen->SetPhiRange(0, 360);

        boxGen->SetBoxXYZ(2.,0.,2.,0.,-2.);

        //proton 2

//        Int_t pdgId = 2212; // proton  beam

//        Double32_t theta1 = 0;  // polar angle distribution

//        Double32_t theta2 = 0;

        Double32_t kin_energy_p2 = 0.01; //GeV

        Double_t mass_p2 = TDatabasePDG::Instance()->GetParticle(pdgId)->Mass();

        Double32_t momentum_p2 = TMath::Sqrt(kin_energy_p2*kin_energy_p2 + 2.*kin_energy_p2*mass_p2); //GeV

        FairBoxGenerator* boxGen_p2 = new FairBoxGenerator(pdgId, 1);

        boxGen_p2->SetThetaRange(theta1, theta2);

        boxGen_p2->SetPRange(momentum, momentum);

        boxGen_p2->SetPhiRange(0, 360);

        boxGen_p2->SetBoxXYZ(1.,-3.8,1.,-3.8,-2.);

        //alpha

        Int_t A = 4;

        Int_t Z = 2;

        Int_t Q = 2;

        Double32_t theta1A = 0;  // polar angle distribution

        Double32_t theta2A = 0;

        Double32_t kin_energy_alpha = 0.07; //GeV

        Double_t massA = 3.727;

        Double32_t momentumAlpha = TMath::Sqrt(kin_energy_alpha*kin_energy_alpha + 2.*kin_energy_alpha*massA); //GeV

        Double_t xA = 1.;

        Double_t yA = 3.;

        Double_t zA = -2.;

        FairIonGenerator *ionGen = new FairIonGenerator(Z, A, Q, 1, 0., 0., momentumAlpha, xA, yA, zA);

        primGen->AddGenerator(boxGen);

//        primGen->AddGenerator(boxGen_p2);

        primGen->AddGenerator(ionGen);

        run->SetGenerator(primGen);

        // ------------------------------------------------------------------------

        //-------Set visualisation flag to true------------------------------------

        run->SetStoreTraj(kTRUE);

        //-------Set LOG verbosity  -----------------------------------------------

        FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");

        // -----   Initialize simulation run   ------------------------------------

        run->Init();

        Int_t nSteps = -15000;

        //gMC->SetMaxNStep(nSteps);

        // -----   Runtime database   ---------------------------------------------

        Bool_t kParameterMerged = kTRUE;

        FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);

        parOut->open(parFile.Data());

        rtdb->setOutput(parOut);

        rtdb->saveOutput();

        rtdb->print();

        // ---------------------------------------------------------

        // -----   Run simulation  ------------------------------------------------

        run->Run(nEvents);

        // -----   Finish   -------------------------------------------------------

        //neuRad->WriteHistos();

        timer.Stop();

        Double_t rtime = timer.RealTime();

        Double_t ctime = timer.CpuTime();

        cout << endl << endl;

        cout << "Macro finished succesfully." << endl;

        cout << "Output file is sim.root" << endl;

        cout << "Parameter file is par.root" << endl;

        cout << "Real time " << rtime << " s, CPU time " << ctime

                        << "s" << endl << endl;

}