/**
 * \file global_sim.C
 *
 * \brief Macro runs simulation for "electron" or "muon" setups of CBM.
 *
 * \author Andrey Lebedev <andrey.lebedev@gsi.de>
 * \date 2010
 **/

void global_sim(Int_t nEvents = 10)
{
	TString script = TString(gSystem->Getenv("SCRIPT"));

	// Specify "electron" or "muon" setup of CBM
//	TString setup = "muon";
	TString setup = "electron";

	// Event parameters
	Int_t NMUONSPLUS = 5; // number of embedded muons
	Int_t NMUONSMINUS = 5; // number of embedded muons
	Int_t NELECTRONS = 5; // number of embedded electrons
	Int_t NPOSITRONS = 5; // number of embedded positrons
	Int_t NPIONSPLUS = 5; // number of embedded pions
	Int_t NPIONSMINUS = 5; // number of embedded pions
	Int_t NPLUTO = 10; // number of embedded particles from pluto
	TString urqmd = "yes"; // If "yes" than UrQMD will be used as background
	TString muons = "no"; // If "yes" than primary muons will be generated
	TString electrons = "yes"; // If "yes" than primary electrons will be generated
	TString pions = "no"; // If "yes" than primary pions will be generated
	TString pluto = "no"; // If "yes" PLUTO particles will be embedded
   TString useUnigen = "no"; // If "yes" than CbmUnigenGenerator will be used instead of FairUrqmdGenerator

	// Files
	TString inFile  = "/d/cbm03/urqmd/auau/25gev/centr/urqmd.auau.25gev.centr.0000.ftn14"; // input UrQMD file
	TString dir = "/data.local1/andrey/events/trd_v11a/"; //directory for output simulation files
	TString mcFile = dir + "mc.0000.root"; //MC file name
	TString parFile = dir + "param.0000.root"; //Parameter file name

	// Pluto files
	TString plutoDir = "/home/d/andrey/pluto/"; // if necessary specify input pluto file to embed signal particles
	std::vector<TString> plutoFile;
	if (script != "yes") {
		plutoFile.resize(NPLUTO);
		for (Int_t i = 0; i < NPLUTO; i++) {
			char name[256];
			if (i < 10) sprintf(name, "Jpsi.000%d", i);
			else if (i < 100) sprintf(name, "Jpsi.00%d", i);
			else sprintf(name, "Jpsi.0%d", i);
			plutoFile[0] = plutoDir + TString(name) + ".root";
		}
	}

	// Geometries
	TString caveGeom = "", targetGeom = "", pipeGeom = "", shieldGeom = "",
			mvdGeom = "", stsGeom = "", muchGeom = "", richGeom = "", trdGeom = "",
			tofGeom = "", ecalGeom = "", fieldMap = "", magnetGeom = "";
	if (setup == "muon") {
		caveGeom   = "cave.geo";
		targetGeom = "target_au_250mu.geo";
		pipeGeom   = "pipe_much.geo";
		shieldGeom = "shield_standard.geo";
		stsGeom    = "sts/sts_v11a.geo";
		muchGeom   = "much/much_v11a.geo";
		trdGeom    = "";//"trd_muon_setup_new.geo";
		tofGeom    = "tof/tof_v07a.geo";
		fieldMap   = "field_v10e";
		magnetGeom = "passive/magnet_v09m.geo";
	} else if (setup == "electron") {
		caveGeom   = "cave.geo";
		targetGeom = "target_au_250mu.geo";
		pipeGeom   = "pipe_standard.geo";
		mvdGeom    = "";//"mvd_v07a.geo";
		stsGeom    = "sts/sts_v11a.geo";
		richGeom   = "rich/rich_v08a.geo";
		trdGeom    = "trd/trd_v11a.geo";//"trd/trd_v10b.geo";
		tofGeom    = "tof/tof_v07a.geo";
		ecalGeom   = "";//"ecal_FastMC.geo";
		fieldMap   = "field_v10e";
		magnetGeom = "passive/magnet_v09e.geo";
	}

	//if SCRIPT environment variable is set to "yes", i.e. macro is run via script
	if (script == "yes") {
		inFile  = TString(gSystem->Getenv("INFILE"));
		mcFile = TString(gSystem->Getenv("MCFILE"));
		parFile = TString(gSystem->Getenv("PARFILE"));

		NMUONSPLUS = TString(gSystem->Getenv("NMUONSPLUS")).Atoi();
		NMUONSMINUS = TString(gSystem->Getenv("NMUONSMINUS")).Atoi();
		NELECTRONS = TString(gSystem->Getenv("NELECTRONS")).Atoi();
		NPOSITRONS = TString(gSystem->Getenv("NPOSITRONS")).Atoi();
		NPIONSPLUS = TString(gSystem->Getenv("NPIONSPLUS")).Atoi();
		NPIONSMINUS = TString(gSystem->Getenv("NPIONSMINUS")).Atoi();
		NPLUTO = TString(gSystem->Getenv("NPLUTO")).Atoi();
		urqmd = TString(gSystem->Getenv("URQMD"));
		muons = TString(gSystem->Getenv("MUONS"));
		electrons = TString(gSystem->Getenv("ELECTRONS"));
		pions = TString(gSystem->Getenv("PIONS"));
		pluto = TString(gSystem->Getenv("PLUTO"));
      useUnigen = TString(gSystem->Getenv("USEUNIGEN"));

		plutoFile.resize(NPLUTO);
		for (Int_t i = 0; i < NPLUTO; i++) {
			char name[256];
			if (i < 10) sprintf(name, "PLUTOFILE%d", i);
			plutoFile[i] = TString(gSystem->Getenv(name));
		}

		caveGeom = TString(gSystem->Getenv("CAVEGEOM"));
		targetGeom = TString(gSystem->Getenv("TARGETGEOM"));
		pipeGeom = TString(gSystem->Getenv("PIPEGEOM"));
		shieldGeom = TString(gSystem->Getenv("SHIELDGEOM"));
		mvdGeom = TString(gSystem->Getenv("MVDGEOM"));
		stsGeom = TString(gSystem->Getenv("STSGEOM"));
		muchGeom = TString(gSystem->Getenv("MUCHGEOM"));
		richGeom = TString(gSystem->Getenv("RICHGEOM"));
		trdGeom = TString(gSystem->Getenv("TRDGEOM"));
		tofGeom = TString(gSystem->Getenv("TOFGEOM"));
		ecalGeom = TString(gSystem->Getenv("ECALGEOM"));
		fieldMap = TString(gSystem->Getenv("FIELDMAP"));
		magnetGeom = TString(gSystem->Getenv("MAGNETGEOM"));
	}
	// -----   Magnetic field   -----------------------------------------------
	Double_t fieldZ = 50.; // field center z position
	Double_t fieldScale = 1.; // field scaling factor

	TStopwatch timer;
	timer.Start();

	gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
	basiclibs();
	gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/cbmrootlibs.C");
	cbmrootlibs();

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

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

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

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

	if ( shieldGeom != "" ) {
		FairModule* shield = new CbmShield("SHIELD");
		shield->SetGeometryFileName(shieldGeom);
		fRun->AddModule(shield);
		cout << "    --- " << shieldGeom << endl;
	}

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

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

	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);
		cout << "    --- " << stsGeom << endl;
	}

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

	if ( muchGeom != "" ) {
		FairDetector* much = new CbmMuch("MUCH", kTRUE);
		much->SetGeometryFileName(muchGeom);
		fRun->AddModule(much);
		cout << "    --- " << muchGeom << endl;
	}

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

	if ( tofGeom != "" ) {
		FairDetector* tof = new CbmTof("TOF", kTRUE);
		tof->SetGeometryFileName(tofGeom);
		fRun->AddModule(tof);
		cout << "    --- " << tofGeom << endl;
	}

	if ( ecalGeom != "" ) {
		FairDetector* ecal = new CbmEcal("ECAL", kTRUE, ecalGeom.Data());
		fRun->AddModule(ecal);
	}
	// ------------------------------------------------------------------------

	// -----   Create magnetic field   ----------------------------------------
	std::cout << "FIELD MAP " << fieldMap << std::endl;
	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();

	if (urqmd == "yes" && useUnigen == "yes") {
		CbmUnigenGenerator*  urqmdGen = new CbmUnigenGenerator(inFile);
		primGen->AddGenerator(urqmdGen);
	}

	if (urqmd == "yes" && useUnigen != "yes") {
	    FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(inFile);
            primGen->AddGenerator(urqmdGen);
	}

	if (pluto == "yes") {
		for (Int_t i = 0; i < NPLUTO; i++) {
			std::cout << "B PLUTO:" << plutoFile[i] << std::endl;
			FairPlutoGenerator *plutoGen= new FairPlutoGenerator(plutoFile[i]);
			primGen->AddGenerator(plutoGen);
			std::cout << "A PLUTO:" << plutoFile[i] << std::endl;
		}
	}

	if (muons == "yes" || electrons == "yes") {
      CbmLitPolarizedGenerator *polGen;
      if (muons == "yes") polGen = new CbmLitPolarizedGenerator(443, NMUONSPLUS);
      else if (electrons == "yes") polGen = new CbmLitPolarizedGenerator(443, NELECTRONS);
      polGen->SetDistributionPt(0.176);        // 25 GeV
      polGen->SetDistributionY(1.9875,0.228);  // 25 GeV
      polGen->SetRangePt(0.,3.);
      polGen->SetRangeY(1.,3.);
      polGen->SetBox(0);
      polGen->SetRefFrame(CbmLitPolarizedGenerator::kHelicity);
      if (muons == "yes") polGen->SetDecayMode(CbmLitPolarizedGenerator::kDiMuon);
      else if (electrons == "yes") polGen->SetDecayMode(CbmLitPolarizedGenerator::kDiElectron);
      polGen->SetAlpha(0);
      polGen->Init();
      primGen->AddGenerator(polGen);
	}

//	if (muons == "yes") {
//		FairBoxGenerator* boxGen1 = new FairBoxGenerator(13, NMUONSMINUS);
//		boxGen1->SetPRange(1.5, 16.);
//		boxGen1->SetPhiRange(0.,360.);
//		boxGen1->SetThetaRange(2.5, 25.);
//		boxGen1->SetCosTheta();
//		boxGen1->Init();
//		primGen->AddGenerator(boxGen1);
//
//		FairBoxGenerator* boxGen2 = new FairBoxGenerator(-13, NMUONSPLUS);
//		boxGen2->SetPRange(1.5, 16.);
//		boxGen2->SetPhiRange(0.,360.);
//		boxGen2->SetThetaRange(2.5, 25.);
//		boxGen2->SetCosTheta();
//		boxGen2->Init();
//		primGen->AddGenerator(boxGen2);
//	}
//
//	if (electrons == "yes") {
//		FairBoxGenerator* boxGen1 = new FairBoxGenerator(11, NELECTRONS);
//		boxGen1->SetPtRange(0.,3.);
////		boxGen1->SetPRange(30, 30);
//		boxGen1->SetPhiRange(0.,360.);
//		boxGen1->SetThetaRange(2.5,25.);
//		boxGen1->SetCosTheta();
//		boxGen1->Init();
//		primGen->AddGenerator(boxGen1);
//
//		FairBoxGenerator* boxGen2 = new FairBoxGenerator(-11, NPOSITRONS);
//		boxGen2->SetPtRange(0.,3.);
////		boxGen2->SetPRange(30, 30);
//		boxGen2->SetPhiRange(0.,360.);
//		boxGen2->SetThetaRange(2.5,25.);
//		boxGen2->SetCosTheta();
//		boxGen2->Init();
//		primGen->AddGenerator(boxGen2);
//	}

	if (pions == "yes") {
		FairBoxGenerator* boxGen1 = new FairBoxGenerator(211, NPIONSPLUS);
//		boxGen1->SetPtRange(0.,3.);
		boxGen1->SetPRange(0.3, 1.);
		boxGen1->SetPhiRange(0.,360.);
		boxGen1->SetThetaRange(2.5,25.);
		boxGen1->SetCosTheta();
		boxGen1->Init();
		primGen->AddGenerator(boxGen1);

		FairBoxGenerator* boxGen2 = new FairBoxGenerator(-211, NPIONSMINUS);
//		boxGen2->SetPtRange(0.,3.);
		boxGen2->SetPRange(0.3, 1.);
		boxGen2->SetPhiRange(0.,360.);
		boxGen2->SetThetaRange(2.5,25.);
		boxGen2->SetCosTheta();
		boxGen2->Init();
		primGen->AddGenerator(boxGen2);
	}

	fRun->SetGenerator(primGen);
	fRun->Init();

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

	// -----   Start run   ----------------------------------------------------
	fRun->Run(nEvents);
	// ------------------------------------------------------------------------

	// -----   Finish   -------------------------------------------------------
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	cout << endl << endl;
	cout << "Macro finished successfully." << endl;
	cout << "Output file is " << mcFile << endl;
	cout << "Real time used: " << rtime << "s " << endl;
	cout << "CPU time used : " << ctime << "s " << endl << endl << endl;
	// ------------------------------------------------------------------------
}