EXPERT ROOT

void exp1803_sim_digi(Int_t nEvents = 1000) {

  // --------------- Telescope T1 -------------------------------------------

  Double_t T1Dl = 0.5;         // [cm]      

  Double_t T1PosZ = 10.;       // [cm] 

  Double_t T1D1Thick = 0.002;  // [cm]  

  Double_t T1D2Thick = 0.1;    // [cm] 

  Double_t T1Side = 6.2;       // [cm]  this parameter should coincide with SD1 side length from /db/QTelescope/QTelescopeParts.xml

  Double_t T1Aperture = 3.1;   // [cm]  

  // --------------- Telescope D1 -------------------------------------------

  Double_t D1PosZ = 20.;       // [cm]

  Double_t D1Thick = 0.03;     // [cm]

  // --------------- BeamDet ------------------------------------------------

  Double_t BeamDetLToF = 1500.;     // [cm] 

  Double_t BeamDetPosZToF = -50;  // [cm] 

  Double_t BeamDetLMWPC = 32.;     // [cm]

  Double_t BeamDetPosZMWPC = -8;  // [cm]  

  // --------------- Beam start position ------------------------------------

  Double_t beamStartPosition = -1600;  // [cm]

  // --------------- Target -------------------------------------------------

  Double_t targetH2Thickness = 0.4;  // [cm] this parameter should coincide with target H2 thickness in /macro/geo/create_target_h2_geo.C

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

  TString outFile= "sim_digi.root";

  TString parFile= "par.root";

  TString workDirPath = gSystem->Getenv("VMCWORKDIR");

  TString paramFileQTelescope = workDirPath

                         + "/db/QTelescope/QTelescopeParts.xml";

  TString paramFileBeamDet = workDirPath

                         + "/db/BeamDet/BeamDetParts.xml";

  TString targetGeoFileName = workDirPath + "/geometry/target_CD2_geo.root";

  TString gadastGeoFileName = workDirPath + "/geometry/partOfGadast.v1.geo.root";

  TString ndGeoFileName = workDirPath + "/geometry/ND.geo.root";

  TString magnetGeoFileName = workDirPath + "/geometry/magnet.geo.root";

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

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

  TStopwatch timer; 

  timer.Start();

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

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

  ERRunSim* run = new ERRunSim();

  /** Select transport engine

  * TGeant3

  * TGeant4

  **/

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

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

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

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

  FairRuntimeDb* rtdb = run->GetRuntimeDb();

  //-------- Set MC event header --------------------------------------------

  ERDecayMCEventHeader* decayMCheader = new ERDecayMCEventHeader();

  run->SetMCEventHeader(decayMCheader);

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

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

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

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

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

  run->AddModule(cave);

  Int_t verbose = 0;

  // -----  BeamDet Setup ---------------------------------------------------

  ERBeamDetSetup* setupBeamDet = ERBeamDetSetup::Instance();

  setupBeamDet->SetXmlParametersFile(paramFileBeamDet);

  // -----  BeamDet parameters ----------------------------------------------

  setupBeamDet->AddToF("ToF1", BeamDetPosZToF - BeamDetLToF);       // 

  setupBeamDet->AddToF("ToF1", BeamDetPosZToF);                     //  BeamDet parts should be added in ascending order   

  setupBeamDet->AddMWPC("MWPC1", BeamDetPosZMWPC - BeamDetLMWPC);   //  of Z-coordinate of part.

  setupBeamDet->AddMWPC("MWPC1", BeamDetPosZMWPC);                  // 

  // setupBeamDet->SetSensitiveTarget();

  // -----   Create target  -------------------------------------------------

  FairModule* target = new ERTarget("targetH2", kTRUE, 1);

  target->SetGeometryFileName(targetGeoFileName);

  run->AddModule(target);

  // -----   Create Part of Gadast ------------------------------------------

  ERGadast* gadast = new ERGadast("PartofGadast", kTRUE, 1);

  gadast->SetGeometryFileName(gadastGeoFileName);

  //run->AddModule(gadast);

  // -----   Create Stilbene wall -------------------------------------------

  ERND* neutronDetector = new ERND("StilbeneWall", kTRUE, 1);

  neutronDetector->SetGeometryFileName(ndGeoFileName);

  //run->AddModule(neutronDetector);

  // -----   Create Magnet geometry -----------------------------------------

  FairModule* magnet = new ERTarget("magnet", 1, kTRUE);

  magnet->SetGeometryFileName(magnetGeoFileName);

  //run->AddModule(magnet);

  // -----  QTelescope Setup ------------------------------------------------

  ERQTelescopeSetup* setupQTelescope = ERQTelescopeSetup::Instance();

  setupQTelescope->SetXMLParametersFile(paramFileQTelescope);

  setupQTelescope->SetGeoName("QTelescopeTmp");

  // ----- T1 parameters ----------------------------------------------------

  Double_t xPos, yPos, zPos;

  Double_t radius = 23.;

  TVector3 rotationT1(0., 12., 0.);

  xPos = radius * TMath::Sin(rotationT1.Y() * TMath::DegToRad());

  yPos = 0.;

  zPos = radius * TMath::Cos(rotationT1.Y() * TMath::DegToRad());

  ERGeoSubAssembly* assemblyT1 = new ERGeoSubAssembly("T1", TVector3(xPos, yPos, zPos), rotationT1);

  ERQTelescopeGeoComponentSingleSi* thinT1 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_1", 

                                                                                  TVector3(0., 0., -5.36), TVector3(), "X");

  ERQTelescopeGeoComponentDoubleSi* thickT1 = new ERQTelescopeGeoComponentDoubleSi("DoubleSi", "DoubleSi_SD2", 

                                                                                  TVector3(0, 0, 0.), TVector3(), "X");

  ERQTelescopeGeoComponentCsI* csi1 = new ERQTelescopeGeoComponentCsI("CsI", "CsI_1", TVector3(0, 0, 5.), TVector3());

  assemblyT1->AddComponent(thinT1);

  assemblyT1->AddComponent(thickT1);

  assemblyT1->AddComponent(csi1);

  setupQTelescope->AddSubAssembly(assemblyT1);

  // ----- T2 parameters ----------------------------------------------------

  radius = 29.;

  TVector3 rotationT2(0., -8.27, 0.);

  xPos = radius * TMath::Sin(rotationT2.Y() * TMath::DegToRad());

  yPos = 0.;

  zPos = radius * TMath::Cos(rotationT2.Y() * TMath::DegToRad());

  ERGeoSubAssembly* assemblyT2 = new ERGeoSubAssembly("T2", TVector3(xPos, yPos, zPos), rotationT2);

  /*ERQTelescopeGeoComponentSingleSi* thinT2 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_1", 

                                                                                  TVector3(0, 0, 0.), TVector3(), "X");*/

  ERQTelescopeGeoComponentDoubleSi* thickT2 = new ERQTelescopeGeoComponentDoubleSi("DoubleSi", "DoubleSi_SD2", 

                                                                                  TVector3(0, 0, 0.), TVector3(), "X");

  ERQTelescopeGeoComponentCsI* csi2 = new ERQTelescopeGeoComponentCsI("CsI", "CsI_1", TVector3(0, 0, 5.), TVector3());

  //assemblyT2->AddComponent(thinT2);

  assemblyT2->AddComponent(thickT2);

  assemblyT2->AddComponent(csi2);

  setupQTelescope->AddSubAssembly(assemblyT2);

  // ------QTelescope -------------------------------------------------------

  ERQTelescope* qtelescope= new ERQTelescope("ERQTelescope", kTRUE,verbose);

  run->AddModule(qtelescope);

  // ------BeamDet ----------------------------------------------------------

  ERBeamDet* beamdet= new ERBeamDet("ERBeamDet", kTRUE,verbose);

  run->AddModule(beamdet);

  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

  Double_t  kinE_MevPerNucleon = 40.;

  // Int_t     Z = 1, A = 3, Q = 1;

  // TString   ionName = "3H";

    Int_t Z = 2, A = 6, Q = 2;

  TString ionName = "6He";

  ERIonMixGenerator* generator = new ERIonMixGenerator(ionName, Z, A, Q, 1);

  Double32_t kin_energy = kinE_MevPerNucleon * 1e-3 * A; //GeV

  generator->SetKinE(kin_energy);

  generator->SetPSigmaOverP(0);

  Double32_t sigmaTheta = 0.004*TMath::RadToDeg();

  // generator->SetKinERange(0,kin_energy);

  generator->SetThetaSigma(0, sigmaTheta);

  generator->SetPhiRange(0, 360);

  generator->SetBoxXYZ(0, 0, 0., 0., beamStartPosition);

  generator->SpreadingOnTarget(); 

  /* List of internal FairRoot names of particles that should coincide with 

     user defined ones (character case is important):

      * Alpha and AntiAlpha (2, 4, 2),

      * Deuteron and AntiDeuteron (1, 2, 1),

      * Triton and AntiTriton (1, 3, 1),

      * HE3 and AntiHE3 (2, 3, 2).

  */

  generator->AddBackgroundIon("Deuteron", 1, 2, 1, 1. / 13.);

  generator->AddBackgroundIon("Triton", 1, 3, 1, 1. / 13.);

  generator->AddBackgroundIon("HE3", 2, 3, 2, 1. / 13.);

  generator->AddBackgroundIon("Alpha", 2, 4, 2, 1. / 13.);

  generator->AddBackgroundIon("6Li", 3, 6, 3, 1. / 13.);

  generator->AddBackgroundIon("7Li", 3, 7, 3, 1. / 13.);

  generator->AddBackgroundIon("8Li", 3, 8, 3, 1. / 13.);

  generator->AddBackgroundIon("9Li", 3, 9, 3, 1. / 13.);

  generator->AddBackgroundIon("9Be", 4, 9, 4, 1. / 13.);

  generator->AddBackgroundIon("10Be", 4, 10, 4, 1. / 13.);

  generator->AddBackgroundIon("11Be", 4, 11, 4, 1. / 13.);

  generator->AddBackgroundIon("12Be", 4, 12, 4, 1. / 13.);

  primGen->AddGenerator(generator);

  run->SetGenerator(primGen);

  // ------- Decayer --------------------------------------------------------

  Double_t massH5 = 4.69036244;  // [GeV]

  ERDecayer* decayer = new ERDecayer();

  ERDecayEXP1803* targetDecay = new ERDecayEXP1803();

  targetDecay->SetInteractionVolumeName("boxCD");

  targetDecay->SetNuclearInteractionLength(1e-3);

  targetDecay->SetAngularDistribution("Cs_6He_d_3He_5H_35-25AMeV.txt");

  targetDecay->SetTargetThickness(targetH2Thickness);

  targetDecay->SetH5Mass(massH5);

  targetDecay->SetH5Exitation(0.0004, 0.00002355, 1);

  targetDecay->SetH5Exitation(0.0012, 0.0002355, 1);

  targetDecay->SetMinStep(1e-5);

  targetDecay->SetMaxPathLength(2e-4 * 10 * 1.1);

  decayer->AddDecay(targetDecay);

  run->SetDecayer(decayer);

  // ------- Magnetic field -------------------------------------------------

  // ERFieldMap* magField = new ERFieldMap("exp1803Field","A"); //exp1803Field, testField

  // magField->SetVolume("magnet");

  // run->SetField(magField);

    // ------- QTelescope Digitizer -------------------------------------------

  ERQTelescopeDigitizer* qtelescopeDigitizer = new ERQTelescopeDigitizer(verbose);

  qtelescopeDigitizer->SetSiElossThreshold(0);

  qtelescopeDigitizer->SetSiElossSigma(0);

  qtelescopeDigitizer->SetSiTimeSigma(0);

  qtelescopeDigitizer->SetCsIElossThreshold(0);

  qtelescopeDigitizer->SetCsIElossSigma(0);

  qtelescopeDigitizer->SetCsITimeSigma(0);

  run->AddTask(qtelescopeDigitizer);

  // ------  Gadast Digitizer -----------------------------------------------

  ERGadastDigitizer* gadastDigitizer = new ERGadastDigitizer(verbose);

  //run->AddTask(gadastDigitizer);

  // -----  BeamDet Digitizer ----------------------------------------------

  ERBeamDetDigitizer* beamDetDigitizer = new ERBeamDetDigitizer(verbose);

  // beamDetDigitizer->SetMWPCElossThreshold(0.006);

  // beamDetDigitizer->SetToFElossThreshold(0.006);

  // beamDetDigitizer->SetToFElossSigmaOverEloss(0);

  // beamDetDigitizer->SetToFTimeSigma(1e-10);

  run->AddTask(beamDetDigitizer);

  ERBeamDetTrackFinder* trackFinder = new ERBeamDetTrackFinder(verbose);

  trackFinder->SetTargetVolume("boxCD");

  run->AddTask(trackFinder);

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

  run->SetStoreTraj(kTRUE);

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

  FairLogger::GetLogger()->SetLogScreenLevel("DEBUG");

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

  run->Init();

  Int_t nSteps = -15000;

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

  Bool_t kParameterMerged = kTRUE;

  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);

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

  rtdb->setOutput(parOut);

  rtdb->saveOutput();

  rtdb->print();

  //gMC->SetMaxNStep(nSteps);

  run->CreateGeometryFile("setup.root");

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

  run->Run(nEvents);

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

  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;

}