EXP_10He

void sim_digi (Int_t nEvents = 100000) {

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

  Double_t BeamDetLToF = 1232.0;     // [cm] 12348

  Double_t BeamDetPosZToF = -95.3;  // [cm] 

  Double_t BeamDetPosZ1MWPC = -81.5;     // [cm]

  Double_t BeamDetPosZ2MWPC = -27.;  // [cm]  

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

  Double_t beamStartPosition = -1600.;  // [cm]

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

  Double_t targetD2Thickness = 0.42;  // [cm] this parameter should coincide with target H2 thickness in /macro/geo/create_target_D2_geo.C

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

  TString outFile= "sim_digi.root";

  TString parFile= "par.root";

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

 // TString paramFileQTelescope = workDirPath

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

  TString paramFileQTelescope = workDirPath 

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

  TString paramFileBeamDet = workDirPath

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

  TString targetGeoFileName = workDirPath + "/geometry/target.T.gas.root";

  // TString targetGeoFileName = workDirPath + "/geometry/target.D2.gas.root";

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

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

  TStopwatch timer; 

  timer.Start();

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

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

  ERRunSim* run = new ERRunSim();

  /** Select transport engine

  * TGeant3

  * TGeant4exp1904_sim_digi.C

  **/

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

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

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

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

  FairRuntimeDb* rtdb = run->GetRuntimeDb();

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

  ER10Heto8HeEventHeader* decayMCheader = new ER10Heto8HeEventHeader();

  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", BeamDetPosZ1MWPC);   //  of Z-coordinate of part.

  setupBeamDet->SetMWPCnumberingInvOrderX(); 

  setupBeamDet->AddMWPC("MWPC1", BeamDetPosZ2MWPC);            //должно быть MWPC2???

  setupBeamDet->SetMWPCnumberingInvOrderX(); 

  //setupBeamDet->SetSensitiveTarget();

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

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

  target->SetGeometryFileName(targetGeoFileName);

  run->AddModule(target);

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

  ERQTelescopeSetup* setupQTelescope = ERQTelescopeSetup::Instance();

  setupQTelescope->SetXMLParametersFile(paramFileQTelescope);

  setupQTelescope->SetGeoName("QTelescopeTmp");

  // ----- instead of annular parameters ----------------------------------------------------

  Double_t xPos, yPos, zPos;

  xPos = 0.;

  yPos = 0.;

  zPos = -10.05;

  TVector3 rotationC(0., 0., 0.);

  ERGeoSubAssembly* assembly_proton= new ERGeoSubAssembly("Telescope_proton", TVector3(xPos, yPos, zPos), rotationC);

  ERQTelescopeGeoComponentDoubleSi* det_proton = new ERQTelescopeGeoComponentDoubleSi("DoubleSi", "DoubleSi_DSD", 

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

  assembly_proton->AddComponent(det_proton);

  setupQTelescope->AddSubAssembly(assembly_proton);

  // ----- square telescope parameters ----------------------------------------------------

  Double_t x,y,z_d1,z_d2,z;

  x = 0.;

  y = 0.;

  z = 25.15;

  z_d1 = 0.75;

  z_d2 = 1.05;

  TVector3 fZeroRotation(0., 0., 0.);

  ERGeoSubAssembly* assembly_he8 = new ERGeoSubAssembly("Telescope_he8", TVector3(x, y, z), fZeroRotation);

  ERQTelescopeGeoComponentSingleSi* det_he8_X = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20", 

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

  ERQTelescopeGeoComponentSingleSi* det_he8_Y1 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20_1", 

                                                                                  TVector3(0., 0., z_d1), TVector3(), "Y");

  ERQTelescopeGeoComponentSingleSi* det_he8_Y2 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20_2", 

                                                                                  TVector3(0., 0., z_d1+z_d2), TVector3(), "Y");

  ERQTelescopeGeoComponentSingleSi* det_he8_Y3 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20_3", 

                                                                                  TVector3(0., 0., z_d1+z_d2*2), TVector3(), "Y");

  ERQTelescopeGeoComponentSingleSi* det_he8_Y4 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20_4", 

                                                                                  TVector3(0., 0., z_d1+z_d2*3), TVector3(), "Y");

  ERQTelescopeGeoComponentSingleSi* det_he8_Y5 = new ERQTelescopeGeoComponentSingleSi("SingleSi", "SingleSi_SSD20_5", 

                                                                                  TVector3(0., 0., z_d1+z_d2*4), TVector3(), "Y");                                                                                                                                                                  

  assembly_he8->AddComponent(det_he8_X);

  assembly_he8->AddComponent(det_he8_Y1);

  assembly_he8->AddComponent(det_he8_Y2);

  assembly_he8->AddComponent(det_he8_Y3);

  assembly_he8->AddComponent(det_he8_Y4);

  assembly_he8->AddComponent(det_he8_Y5);

  setupQTelescope->AddSubAssembly(assembly_he8);

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

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

  run->AddModule(qtelescope);

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

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

  run->AddModule(beamdet);

  // ------ND ---------------------------------------------------------------

  ERND* nd= new ERND("ERND", kTRUE,verbose);

  nd->SetGeometryFileName(ndGeoFileName);

  run->AddModule(nd);

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

  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

  Double_t  kinE_MevPerNucleon = 21.5;

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

  TString ionName = "8He";

  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->SetThetaSigma(0, sigmaTheta);

  generator->SetPhiRange(0, 360);

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

  generator->SpreadingOnTarget(); 

  primGen->AddGenerator(generator);

  run->SetGenerator(primGen);

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

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

  Double_t massHe10 = 9.36268;//7.5061760;  // [GeV] 10He

  ERDecayer* decayer = new ERDecayer();

  ERDecay10Heto8He* targetDecay = new ERDecay10Heto8He();

  targetDecay->SetHe10Mass(massHe10);

  // // ERDecayEXP1811* targetDecay = new ERDecayEXP1811();

  // // targetDecay->SetH7Mass(massHe10);

  targetDecay->SetInteractionVolumeName("tubeD2");

  targetDecay->SetNuclearInteractionLength(20.);

  targetDecay->SetAngularDistribution("cos_tetta_cross.txt");

  // targetDecay->SetDecayFile("pmom-pv-1_short.dat", 0.0005 /*excitation in file [GeV]*/);

  targetDecay->SetDecayFile("10he_int-3030_short.dat", 0.001 /*excitation in file [GeV]*/);

  //targetDecay->SetH7Exitation(0.0004, 0.00002355, 1);

  //targetDecay->SetH7Exitation(0.0012, 0.0002355, 1);

  targetDecay->SetMinStep(1e-1);

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

  decayer->AddDecay(targetDecay);

  run->SetDecayer(decayer);

  // ------- 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);

  // -----  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);

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

  ERNDDigitizer* ndDigitizer = new ERNDDigitizer(1);

  ndDigitizer->SetEdepError(0.0,0.04,0.02);

  ndDigitizer->SetLYError(0.0,0.04,0.02);

  ndDigitizer->SetTimeError(0.001);

  ndDigitizer->SetQuenchThreshold(0.005);

  ndDigitizer->SetLYThreshold(0.004);

  ndDigitizer->SetProbabilityB(0.1);

  ndDigitizer->SetProbabilityC(0.3);

  run->AddTask(ndDigitizer);

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

  //run->SetStoreTraj(kTRUE);

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

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

  // -----   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();

  run->CreateGeometryFile("setup_exp1904.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 " << outFile << endl;

  cout << "Parameter file is " << parFile << endl;

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

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

  TDatabasePDG *tdb = TDatabasePDG::Instance();

//  tdb->Print();

}