void trd_elid_sim(Int_t nEvents = 20)
{
   Double_t momentum = 1.5;

   Double_t thetaMin, thetaMax;
   thetaMin = 2.5;
   thetaMax = 25.;

   TString inFile = "", parFile = "", outFile ="";
   TString caveGeom = "", targetGeom = "", pipeGeom   = "", magnetGeom = "", mvdGeom = "",
           mvdGeom = "",stsGeom = "", richGeom= "", trdGeom = "";
   TString fieldMap = "";

	TString outFile = "/d/cbm06/user/slebedev/trd/piel.0000.mc.root";
	TString parFile = "/d/cbm06/user/slebedev/trd/piel.0000.params.root";

	// -----  Geometries  -----------------------------------------------------
   caveGeom = "cave.geo";
   targetGeom = "target_au_250mu.geo";
   pipeGeom = "pipe_standard.geo";
   mvdGeom = "";//"mvd_v08a.geo";
   stsGeom = "sts/sts_v11a.geo";
   richGeom = "rich/rich_v08a.geo";
   trdGeom = "trd/trd_v10b.geo";
   fieldMap = "field_v10e";
   magnetGeom = "passive/magnet_v09e.geo";
	Double_t fieldZ = 50.; // field center z position
	Double_t fieldScale = 1.; // field scaling factor

   gDebug = 0;

	TStopwatch timer;
	timer.Start();

   // ----  Load libraries   -------------------------------------------------
   gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
   basiclibs();

   gROOT->LoadMacro("$VMCWORKDIR/macro/rich/cbmlibs.C");
   cbmlibs();

   FairRunSim* fRun = new FairRunSim();
   fRun->SetName("TGeant3"); // Transport engine
   fRun->SetOutputFile(outFile);
   FairRuntimeDb* rtdb = fRun->GetRuntimeDb();

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

   if ( caveGeom != "" ) {
      FairModule* cave = new CbmCave("CAVE");
      cave->SetGeometryFileName(caveGeom);
      fRun->AddModule(cave);
   }

   if ( pipeGeom != "" ) {
      FairModule* pipe = new CbmPipe("PIPE");
      pipe->SetGeometryFileName(pipeGeom);
      fRun->AddModule(pipe);
   }

   if ( targetGeom != "" ) {
      FairModule* target = new CbmTarget("Target");
      target->SetGeometryFileName(targetGeom);
      fRun->AddModule(target);
   }

   if ( magnetGeom != "" ) {
      FairModule* magnet = new CbmMagnet("MAGNET");
      magnet->SetGeometryFileName(magnetGeom);
      fRun->AddModule(magnet);
   }

   if ( mvdGeom != "" ) {
      FairDetector* mvd = new CbmMvd("MVD", kTRUE);
      mvd->SetGeometryFileName(mvdGeom);
      fRun->AddModule(mvd);
   }

   if ( stsGeom != "" ) {
      FairDetector* sts = new CbmSts("STS", kTRUE);
      sts->SetGeometryFileName(stsGeom);
      fRun->AddModule(sts);
   }

   if ( richGeom != "" ) {
      FairDetector* rich = new CbmRich("RICH", kTRUE);
      rich->SetGeometryFileName(richGeom);
      fRun->AddModule(rich);
   }

   if ( trdGeom != "" ) {
      FairDetector* trd = new CbmTrd("TRD",kTRUE );
      trd->SetGeometryFileName(trdGeom);
      fRun->AddModule(trd);
   }

   // -----   Create magnetic field   ----------------------------------------
   CbmFieldMap* magField = NULL;
   if (fieldMap == "field_electron_standard" || fieldMap == "field_v10e")
      magField = new CbmFieldMapSym2(fieldMap);
   else if (fieldMap == "field_muon_standard" )
      magField = new CbmFieldMapSym2(fieldMap);
   else if (fieldMap == "FieldMuonMagnet" )
      magField = new CbmFieldMapSym3(fieldMap);
   else {
      cout << "===> ERROR: Unknown field map " << fieldMap << endl;
      exit;
   }
   magField->SetPosition(0., 0., fieldZ);
   magField->SetScale(fieldScale);
   fRun->SetField(magField);


   FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

	Int_t kfCode1 = 11; // electrons
	Int_t kfCode2 = -11; // positrons
	Double_t minMomentum = momentum; //minimum momentum
	Double_t maxMomentum = momentum; //maximum momentum

	FairBoxGenerator* boxGen1 = new FairBoxGenerator(11, 60);
	boxGen1->SetPRange(minMomentum, maxMomentum);
//	boxGen1->SetXYZ(50., 50., 450.);
	boxGen1->SetPhiRange(0., 360.);
	boxGen1->SetThetaRange(thetaMin, thetaMax);
	boxGen1->SetCosTheta();
	boxGen1->Init();
	primGen->AddGenerator(boxGen1);

	FairBoxGenerator* boxGen2 = new FairBoxGenerator(-11, 60);
	boxGen2->SetPRange(minMomentum, maxMomentum);
//	boxGen2->SetXYZ(50., 50., 450.);
	boxGen2->SetPhiRange(0., 360.);
	boxGen2->SetThetaRange(thetaMin, thetaMax);
	boxGen1->SetCosTheta();
	boxGen2->Init();
	primGen->AddGenerator(boxGen2);

	FairBoxGenerator* boxGen3 = new FairBoxGenerator(211, 60);
	boxGen3->SetPRange(minMomentum, maxMomentum);
//	boxGen3->SetXYZ(50., 50., 450.);
	boxGen3->SetPhiRange(0., 360.);
	boxGen3->SetThetaRange(thetaMin, thetaMax);
	boxGen1->SetCosTheta();
	boxGen3->Init();
	primGen->AddGenerator(boxGen3);

	FairBoxGenerator* boxGen4 = new FairBoxGenerator(-211, 60);
	boxGen4->SetPRange(minMomentum, maxMomentum);
//	boxGen4->SetXYZ(50., 50., 450.);
	boxGen4->SetPhiRange(0., 360.);
	boxGen4->SetThetaRange(thetaMin, thetaMax);
	boxGen1->SetCosTheta();
	boxGen4->Init();
	primGen->AddGenerator(boxGen4);

	fRun->SetGenerator(primGen);


	fRun->Init();

	CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
	fieldPar->SetParameters(magField);
	fieldPar->setChanged();
	fieldPar->setInputVersion(fRun->GetRunId(), 1);
	Bool_t kParameterMerged = kTRUE;
	FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
	parOut->open(parFile.Data());
	rtdb->setOutput(parOut);
	rtdb->saveOutput();
	rtdb->print();

	fRun->Run(nEvents);

	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;


	exit(0);
}