void exp1803_sim_digi(Int_t nEvents = 100) {
	// --------------- 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]
//	Double_t beamStartPosition = -150;
	// --------------- 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.h2.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);

	// ----- T1 parameters ----------------------------------------------------
	// ----- T1.1--------------------------------------------------------------
	setupQTelescope->AddSi("DoubleSi_SD1", TVector3( T1Side/2 + T1Aperture/2,
			T1Side/2 - T1Aperture/2,
			T1PosZ + T1D1Thick/2), "X");
	setupQTelescope->AddSi("DoubleSi_SD2", TVector3( T1Side/2 + T1Aperture/2,
			T1Side/2 - T1Aperture/2,
			T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
	// ----- T1.2--------------------------------------------------------------
	setupQTelescope->AddSi("DoubleSi_SD1", TVector3( T1Side/2 - T1Aperture/2,
			-T1Side/2 - T1Aperture/2,
			T1PosZ + T1D1Thick/2), "X");
	setupQTelescope->AddSi("DoubleSi_SD2", TVector3( T1Side/2 - T1Aperture/2,
			-T1Side/2 - T1Aperture/2,
			T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
	// ----- T1.3 -------------------------------------------------------------
	setupQTelescope->AddSi("DoubleSi_SD1", TVector3(-T1Side/2 - T1Aperture/2,
			-T1Side/2 + T1Aperture/2,
			T1PosZ + T1D1Thick/2), "X");
	setupQTelescope->AddSi("DoubleSi_SD2", TVector3(-T1Side/2 - T1Aperture/2,
			-T1Side/2 + T1Aperture/2,
			T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");
	// ----- T1.4--------------------------------------------------------------
	setupQTelescope->AddSi("DoubleSi_SD1", TVector3(-T1Side/2 + T1Aperture/2,
			T1Side/2 + T1Aperture/2,
			T1PosZ + T1D1Thick/2), "X");
	setupQTelescope->AddSi("DoubleSi_SD2", TVector3(-T1Side/2 + T1Aperture/2,
			T1Side/2 + T1Aperture/2,
			T1PosZ + T1D1Thick +T1Dl + T1D2Thick/2), "X");

	// ----- D1 parameters ----------------------------------------------------
	setupQTelescope->AddSi("DoubleSi_D1", TVector3( 0,
			0,
			D1PosZ + D1Thick/2), "X");

	// ------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 = 30.;
//	Double_t  kinE_MevPerNucleonMax = 32.;
	Double_t  kinE_MevPerNucleon = 1.;
	Double_t  kinE_MevPerNucleonMax = 32.;
	// 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);
//	ERIonMixGenerator* generator = new ERIonMixGenerator("Alpha", 2, 4, 2, 1);

//	generator->AddBackgroundIon("Triton", 1, 3, 1, 0.2 / 0.20);
//	generator->AddBackgroundIon("Alpha", 1, 4, 1, 0.4 / 0.60);
//	generator->AddBackgroundIon("HE3", 2, 3, 2, 0.4 / 0.20);
//	generator->AddBackgroundIon("Deuteron", 1, 2, 1, 0.2 / 0.20);
//	generator->AddBackgroundIon("proton", 1, 1, 1, 0.2 / 0.20);

//	generator->AddBackgroundIon("proton", 1, 1, 1, 0.2 / 0.20);
//	generator->AddBackgroundIon("8He", 2, 8, 2, 0.2 / 0.40);
//	generator->AddBackgroundIon("11Li", 3, 11, 3, 0.2 / 0.40);
//	generator->AddBackgroundIon("9Li", 3, 9, 3, 0.1 / 0.40);

//	generator->AddBackgroundIon("4He", 2, 4, 2, 20);
//	generator->AddBackgroundIon("8He", 2, 8, 2, 0.2/0.8);
//	generator->AddBackgroundIon("9Li", 3, 9, 3, 0.05/0.8);


	Double32_t kin_energy = kinE_MevPerNucleon * 1e-3 * A; //GeV
	Double32_t kin_energy_max = kinE_MevPerNucleonMax * 1e-3 * A; //GeV
//	generator->SetKinE(kin_energy);
//	generator->SetPSigmaOverP(0);
	generator->SetKinERange(kin_energy, kin_energy_max);
	Double32_t sigmaTheta = 0.004*TMath::RadToDeg();
	generator->SetThetaSigma(0, 0);
//	generator->SetThetaRange(1.3, 1.3);
//	generator->SetPhiRange(0, 360);
//	generator->SetPhiRange(0, 270);
//	generator->SetPhiRange(45, 45);
//	generator->SetBoxXYZ(0, 0, 0, 0, beamStartPosition);
//	generator->SetBoxXYZ(-1, -1, -1, -1, -45.);
//	generator->SetBoxXYZ(0, 0, 0, 0, -45.);
//	generator->SetBoxXYZ(-2, -2, -2, -2, -45.);
//	generator->SetBoxXYZ(-1, -1, -1, -1, beamStartPosition);
//	generator->SetBoxXYZ(5, 5, 5, 5, beamStartPosition);
//	generator->SetBoxXYZ(5, 5, 5, 5, 0.);
	generator->SetBoxXYZ(5, 0., 5.5, 1., 2.);

//	generator->SpreadingOnTarget();

	primGen->AddGenerator(generator);
	run->SetGenerator(primGen);
	// ------- Decayer --------------------------------------------------------
//	Double_t massH5 = 4.8;  // [GeV]
	Double_t massH5 = (2*939.565 + 2809.431)*0.001;  // [GeV]

	ERDecayer* decayer = new ERDecayer();
	ERDecayEXP1803* targetDecay = new ERDecayEXP1803();
	targetDecay->SetTargetVolumeName("tubeH2");
	targetDecay->SetTargetThickness(targetH2Thickness);
	// targetDecay->SetH5Mass(massH5);
	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->SetSiElossThreshold(0.0005);
	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.0000012);
//	beamDetDigitizer->SetMWPCElossThreshold(0.00000001);
	// beamDetDigitizer->SetToFElossThreshold(0.006);
	// beamDetDigitizer->SetToFElossSigmaOverEloss(0);
//	beamDetDigitizer->SetToFTimeSigma(1e-10);
	beamDetDigitizer->SetToFTimeSigma(0.4/2.35);
	run->AddTask(beamDetDigitizer);
	// ------- BeamDet TrackFinder -------------------------------------------
	ERBeamDetTrackFinder* trackFinder = new ERBeamDetTrackFinder(verbose);
	run->AddTask(trackFinder);
	// -----------------------BeamDetTrackPID----------------------------------
	ERBeamDetPID* pid = new ERBeamDetPID(verbose);
	pid->SetIonMassNumber(A);
	pid->SetBoxPID(0., 1000., 0., 1000.);
	pid->SetOffsetToF(0.);
	pid->SetProbabilityThreshold(0);
	run->AddTask(pid);
	// ------- QTelescope TrackFinder -------------------------------------------
	ERQTelescopeTrackFinder* qtelescopeTrackFinder = new ERQTelescopeTrackFinder(verbose);
	qtelescopeTrackFinder->SetHitStation("DoubleSi_SD2_XY_1");
	qtelescopeTrackFinder->SetHitStation("DoubleSi_SD2_XY_3");
	qtelescopeTrackFinder->SetStripEdepRange(0., 100.);          // [GeV]
	qtelescopeTrackFinder->SetTargetPoint(0., 0., 0.);
	// qtelescopeTrackFinder->SetStripEdepRange(0.0097, 100.);   // [GeV]
	// qtelescopeTrackFinder->SetEdepDiffXY(5.);                 // [GeV]
	qtelescopeTrackFinder->SetEdepMaxDiffXY(0.5);
	run->AddTask(qtelescopeTrackFinder);
	//-------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();
	//gMC->SetMaxNStep(nSteps);
	// -----   Run simulation  ------------------------------------------------
	run->Run(nEvents);

//	TDatabasePDG *tdb = TDatabasePDG::Instance();
//	tdb->Print();

	// -----   Finish   -------------------------------------------------------
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	std::cout << std::endl << std::endl;
	std::cout << "Macro finished succesfully." << std::endl;
	std::cout << "Output file is sim.root" << std::endl;
	std::cout << "Parameter file is par.root" << std::endl;
	std::cout << "Real time " << rtime << " s, CPU time " << ctime
			<< "s" << std::endl << std::endl;
}