// Example for usage of event filter in combination with dpm direct
// Martin Galuska, Katja Kleeberg


sim_filter_ex2(Int_t nEvents = 100, TString  SimEngine ="TGeant3", Float_t mom = 6.231552)
{
	//-----User Settings:-----------------------------------------------
	TString  OutputFile     ="sim_filter_ex1.root";
	TString  ParOutputfile  ="simparams.root";
	TString  MediaFile      ="media_pnd.geo";
	gDebug                  = 0;
	TString digiFile        = "all.par"; //The emc run the hit producer directly
	// choose your event generator
	Bool_t UseEvtGen	      =kFALSE;
	Bool_t UseEvtGenDirect      =kFALSE;
	Bool_t UseDpm 	      =kTRUE;
	Bool_t UseBoxGenerator      =kFALSE;

	Double_t BeamMomentum = 0.; // beam momentum ONLY for the scaling of the dipole field.
	if (UseBoxGenerator)
	{
		BeamMomentum   =15.0; // ** change HERE if you run Box generator
	}
	else
	{
		BeamMomentum = mom;  // for DPM/EvtGen BeamMomentum is always = mom
	}

	//------------------------------------------------------------------
	TStopwatch timer;
	timer.Start();
	gRandom->SetSeed();

	// Create the Simulation run manager--------------------------------
	FairRunSim *fRun = new FairRunSim();
	fRun->SetName(SimEngine.Data() );
	fRun->SetOutputFile(OutputFile.Data());
	fRun->SetGenerateRunInfo(kFALSE);
	fRun->SetBeamMom(BeamMomentum);
	fRun->SetMaterials(MediaFile.Data());
	FairRuntimeDb *rtdb=fRun->GetRuntimeDb();

	// Set the parameters
	//-------------------------------
	TString allDigiFile = gSystem->Getenv("VMCWORKDIR");
	allDigiFile += "/macro/params/";
	allDigiFile += digiFile;


	//-------Set the parameter output --------------------
	FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
	parIo1->open(allDigiFile.Data(),"in");
	rtdb->setFirstInput(parIo1);

	//---------------------Set Parameter output      ----------
	Bool_t kParameterMerged=kTRUE;
	FairParRootFileIo* output=new FairParRootFileIo(kParameterMerged);
	output->open(ParOutputfile.Data());
	rtdb->setOutput(output);

	// Create and add detectors

	//-------------------------  CAVE      -----------------

	FairModule *Cave= new PndCave("CAVE");
	Cave->SetGeometryFileName("pndcave.geo");
	fRun->AddModule(Cave);
	//-------------------------  Magnet   -----------------
	FairModule *Magnet= new PndMagnet("MAGNET");
	//Magnet->SetGeometryFileName("FullSolenoid_V842.root");
	Magnet->SetGeometryFileName("FullSuperconductingSolenoid_v831.root");
	fRun->AddModule(Magnet);
	FairModule *Dipole= new PndMagnet("MAGNET");
	Dipole->SetGeometryFileName("dipole.geo");
	fRun->AddModule(Dipole);
	//-------------------------  Pipe     -----------------
	FairModule *Pipe= new PndPipe("PIPE");
	Pipe->SetGeometryFileName("beampipe_201309.root");
	fRun->AddModule(Pipe);
//	//-------------------------  STT       -----------------
//	FairDetector *Stt= new PndStt("STT", kTRUE);
//	Stt->SetGeometryFileName("straws_skewed_blocks_35cm_pipe.geo");
//	fRun->AddModule(Stt);
//	//-------------------------  MVD       -----------------
//	FairDetector *Mvd = new PndMvdDetector("MVD", kTRUE);
//	Mvd->SetGeometryFileName("Mvd-2.1_FullVersion.root");
//	fRun->AddModule(Mvd);
//	//-------------------------  GEM       -----------------
//	FairDetector *Gem = new PndGemDetector("GEM", kTRUE);
//	Gem->SetGeometryFileName("gem_3Stations.root");
//	fRun->AddModule(Gem);
//	//-------------------------  EMC       -----------------
//	PndEmc *Emc = new PndEmc("EMC",kTRUE);
//	Emc->SetGeometryVersion(1);
//	Emc->SetStorageOfData(kFALSE);
//	fRun->AddModule(Emc);
//	//-------------------------  SCITIL    -----------------
//	FairDetector *SciT = new PndSciT("SCIT",kTRUE);
//	SciT->SetGeometryFileName("barrel-SciTil_07022013.root");
//	fRun->AddModule(SciT);
//	//-------------------------  DRC       -----------------
//	PndDrc *Drc = new PndDrc("DIRC", kTRUE);
//	Drc->SetGeometryFileName("dirc_l0_p0_updated.root");
//	Drc->SetRunCherenkov(kFALSE);
//	fRun->AddModule(Drc);
//	//-------------------------  DISC      -----------------
//	PndDsk* Dsk = new PndDsk("DSK", kTRUE);
//	Dsk->SetStoreCerenkovs(kFALSE);
//	Dsk->SetStoreTrackPoints(kFALSE);
//	fRun->AddModule(Dsk);
//	//-------------------------  MDT       -----------------
//	PndMdt *Muo = new PndMdt("MDT",kTRUE);
//	Muo->SetBarrel("fast");
//	Muo->SetEndcap("fast");
//	Muo->SetMuonFilter("fast");
//	Muo->SetForward("fast");
//	Muo->SetMdtMagnet(kTRUE);
//	Muo->SetMdtMFIron(kTRUE);
//	fRun->AddModule(Muo);
//	//-------------------------  FTS       -----------------
//	FairDetector *Fts= new PndFts("FTS", kTRUE);
//	Fts->SetGeometryFileName("fts.geo");
//	fRun->AddModule(Fts);
//	//-------------------------  FTOF      -----------------
//	FairDetector *FTof = new PndFtof("FTOF",kTRUE);
//	FTof->SetGeometryFileName("ftofwall.root");
//	fRun->AddModule(FTof);
//	//-------------------------  RICH       ----------------
//	FairDetector *Rich= new PndRich("RICH",kFALSE);
//	Rich->SetGeometryFileName("rich_v2_shift.geo");
//	fRun->AddModule(Rich);

	// Create and Set Event Generator
	//-------------------------------
	FairFilteredPrimaryGenerator* primGen = new FairFilteredPrimaryGenerator();
	fRun->SetGenerator(primGen);

	if(UseBoxGenerator){	// Box Generator
		FairBoxGenerator* boxGen = new FairBoxGenerator(22, 5); // 13 = muon; 1 = multipl.
		boxGen->SetPRange(mom,mom); // GeV/c
		boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
		boxGen->SetThetaRange(0., 90.); // Polar angle in lab system range [degree]
		boxGen->SetXYZ(0., 0., 0.); // cm
		primGen->AddGenerator(boxGen);
	}
	if(UseDpm){
		PndDpmDirect *Dpm= new PndDpmDirect(mom,1);
		primGen->AddGenerator(Dpm);
	}
	if(UseEvtGen){
		TString  EvtInput =gSystem->Getenv("VMCWORKDIR");
		EvtInput+="/input/psi2s_jpsi2pi_1k.evt";
		FairEvtGenGenerator* evtGen = new FairEvtGenGenerator(EvtInput.Data());
		primGen->AddGenerator(evtGen);
	}
	if(UseEvtGenDirect){
		TString  EvtInput =gSystem->Getenv("VMCWORKDIR");
		EvtInput+="/macro/run/psi2s_Jpsi2pi_Jpsi_mumu.dec";
		PndEvtGenDirect *EvtGen = new PndEvtGenDirect("pbarpSystem", EvtInput.Data(), mom);
		EvtGen->SetStoreTree(kTRUE);
		primGen->AddGenerator(EvtGen);
	}


	// now add the event filter(s)

	// set the maximal number of attempts to find an event with matching properties
	primGen->SetFilterMaxTries(1000000); // for testing small number, for real production set usually to 9999999 or something very big
	primGen->SetVerbose(); // highest commenting level of the FairPrimaryGenerator

	// filter on events with pi+ pi- only
	FairEvtFilterOnSingleParticleCounts* pionFilter = new FairEvtFilterOnSingleParticleCounts("pionFilter");
//	pionFilter->SetVerbose();// highest commenting level of the FairEvtFilterOnCounts
	pionFilter->AndMinMaxPdgCodes(1,1,211);// min. 1 and max. 1 pi+
	pionFilter->AndMinMaxPdgCodes(1,1,-211);// min. 1 and max. 1 pi-
	pionFilter->AndMinMaxAllParticles(2,2);// min. 2 and max. 2 particles in total
	primGen->AndFilter(pionFilter);// add filter to fFilterList with logical operation AND (this is the first filter, so AND and OR would give the same result)






	//---------------------Create and Set the Field(s)----------
	PndMultiField *fField= new PndMultiField("AUTO");
	fRun->SetField(fField);

	// EMC Hit producer
	//-------------------------------
	PndEmcHitProducer* emcHitProd = new PndEmcHitProducer();
	fRun->AddTask(emcHitProd);

	//-------------------------  Initialize the RUN  -----------------
	fRun->Init();
	//-------------------------  Run the Simulation  -----------------
	fRun->Run(nEvents);
	primGen->WriteEvtFilterStatsToRootFile();
	//-------------------------  Save the parameters -----------------
	rtdb->saveOutput();
	//------------------------Print some info and exit----------------
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

	cout << " Test passed" << endl;
	cout << " All ok " << endl;

	//exit(0);

}