Primary generator with added event filtering capabilities. More...

#include <PndFilteredPrimaryGenerator.h>

Inheritance diagram for PndFilteredPrimaryGenerator:

Public Member Functions
	PndFilteredPrimaryGenerator ()
	Default constructor. More...

	PndFilteredPrimaryGenerator (TString filterset)
	Constructor with filter-set string. More...

virtual	~PndFilteredPrimaryGenerator ()
	Destructor. More...

virtual Bool_t	Init ()
	Initialize the event generator(s) and the event (veto) filter(s). More...

void	AddFilter (TString filterStr)
	Registers a filter as a string to be parsed Each filter set consist of one filter definition or a number of filters connected with logical && (AND) (= a certain event type). All filter sets are combined with a logical \|\| (OR) (all event types to be accepted). Syntax: "Filt_1.1 && Filt_1.2 \|\| Filt_2.1 \|\| Filt_3.1 && Filt_3.2 && Filt_3.2". More...

virtual Bool_t	GenerateEvent (FairGenericStack *pStack)
	Calls event generators and the event filters. More...

void	SetFilterMaxTries (Int_t maxTries=99999)
	Define the maximum number of times that this object should try to find an event which suits all event filters. More...

Int_t	GetNumberOfFilterMaxTries ()
	returns the maximum number of times that this object should try to find an event which suits all event filters. More...

Int_t	GetNumberOfGeneratedEvents ()
	returns the total (accepted + rejected) number of events generated by the event generators. If no event filters are used this number is equal to the number of simulated events. More...

void	SetEventPrintFrequency (int freq)
	Sets the frequency (accepted events) for printout (verbose>0) of accepted and generated events. More...

Int_t	GetNumberOfFilterFailedEvents ()
	Returns the number of cases in which no matching event was found within the set max. tries. More...

void	WriteEvtFilterStatsToRootFile (TFile *outputFile=NULL)
	Writes all relevant event filter information to the output root file. More...

void	SetVerbose (Int_t verbose=12)
	Set the level of commenting output. More...

Protected Member Functions
PndSmpCandList	CombineList (int pdg, PndSmpCandList l0, PndSmpCandList l1, PndSmpCandList l2=0, PndSmpCandList l3=0, PndSmpCandList *l4=0)
	Combines upt to five particle lists of PndSmpCand with overlap and double counting prevention. More...

void	PrintSmpCandList (PndSmpCandList l, TString name="")
	Prints a candidate lits. More...

StrVec	SplitString (TString s, TString delim=" ")
	Splits a TString to substrings. More...

bool	CheckKinematic (const PndSmpFilt &f, const TLorentzVector &p4)
	Checks whether P4 kinematics match the criteria of a PndSmpFilt. More...

int	AntiPdgCode (int pdg)
	Gets anti-pdg code, if exists. If not returns the code itself (particle is its anti-particle) More...

void	GetRangeDouble (TString s, double &a, double &b, TString delim=",", bool forceset=false)
	Turns a string of the form <float><delim><float> (e.g. "124.2,178.3") to a range a,b. More...

void	GetRangeInt (TString s, int &a, int &b, TString delim="..")
	Turns a string of the form <int><delim><int> (e.g. "3..8") to a range a,b. More...

Protected Attributes
std::vector< PndSmpFilterSet >	fFilterSets
	Contains the filter-sets. Each filter set consist of one filter definition or a number of filters connected with logical AND (= a certain event type). All filter sets are combined with a logical OR (all event types to be accepted). Structure is like: (Filt_1.1 AND Filt_1.2) OR (Filt_2.1) OR (Filt_3.1 AND Filt_3.2 AND Filt_3.2) More...

FairEvtFilterParams	fEvtFilterStat
	Contains the statistics of the event filtering process. More...

Int_t	fVerbose
	Level of commenting output, 0 means no output, higher gives more output. More...

Int_t	fEventNrFiltered
	Event number (Set by the filtered primary generator. More...

Int_t	fEventPrintFreq
	Print frequency for filtered events. More...

TDatabasePDG *	fdbPdg
	Shortcut to TDatabasePDG. More...

std::vector< TString >	fPartNames
	particle names for particles to count (with and w/o charged specification, also simply tracks and neutrals) More...

std::vector< int >	fNamePdg
	particle codes for particles to count (with and w/o charged specification, also simply tracks and neutrals) More...

std::vector< int >	fCombFsPdg
	particle codes for the lists used for combinatorics; !!! the codes are not selected from MC truth, just the corresponding mass hypos are applied to charged tracks More...

std::unordered_set< int >	fSetFsPdg
	set to identify particles from MC truth list which can be combined More...

std::map< TString, int >	fNameCodeMap
	mapes names to (pdg) codes More...

std::map< int, TString >	fCodeNameMap
	mapes (pdg) codes to names More...

Private Member Functions
	PndFilteredPrimaryGenerator (const PndFilteredPrimaryGenerator &)

PndFilteredPrimaryGenerator &	operator= (const PndFilteredPrimaryGenerator &)

	ClassDef (PndFilteredPrimaryGenerator, 1)

Detailed Description

Primary generator with added event filtering capabilities.

Author: Klaus Goetzen <k [dot] goetzen (at) gsi [dot] de>

This class adds event filtering capabilities to FairPrimaryGenerator which is used internally for handling the event generators. The event filtering is performed after the event generation and before the particle transport through the detector model.

The filter is configured with a string of the form "Filt_1.1 && Filt_1.2 || Filt_2.1 || Filt_3.1 && Filt_3.2 && Filt_3.2", where the logical && (AND) has higher priority than the logical || (OR). Each filter can be negated with a leading '!'.

There are two types of filters possible: particle type counters and invariant mass filters. The syntax is

count filter : C(<particle> ; <mult> ; [kinematic requirements]) ==> requires <mult> particles of type <particle> (counted in MC list) to meet certain kinematic conditions (p, pt, pz, theta, phi)

mass filter : M(<particle 1> <particle 2> [<particle 3> ... <particle 5>] ; <mult> ; <mass requirement>=""> ; [kinematic requirements]) ==> takes neutral particles and charged tracks, applies the mass hypotheses from <particle 1..5>, combines them and requires <mult> candidates with invariant mass to be in certain window

<particle> : e+, e-, e+-, mu+, mu-, mu+-, pi+, pi-, pi+-, K+, K-, K+-, p+, p-, p+-, n0, n0b, n00b (=n0 or n0b), gam, nt (neutral = gamma), t+ (positive tracks), t- (neg. tracks), t+- (chg. tracks), any (total multiplicity) <mult> : a..b, a (= a..a), ..a (= 0..a), a.. = (a..10000), not given (>=1 = 1..10000) <mass req.> : m[center,window] : center - window/2 < inv. mass < center + window/2 [kin. req.] : p[min,max] ; pt[min,max] ; pz[min,max] ; tht[min,max] ; phi[min,max]. Can be shortened to [,max] = [0.0,max] and [min,] = [min,1e8]

Example: C(any ; 4..) && C(tr+- ; 2 ; tht[5,20] ; p[1.0,]) && M(pi+ pi- gam gam ; m[0.782,0.1] ; pt[0.5,2.0])

at least >=4 particles in event
AND exactly two charged tracks with 5� < theta < 20� and p > 1.0 GeV/c
AND at least one possible pi+ pi- gamma gamma combination with |m - 0.782| < 0.1 GeV/c2 and 0.5 < pt < 2.0 GeV/c

Concerning FairPrimaryGenerator:

The PndFilteredPrimaryGenerator is responsible for the handling of the MC input. Several input generators can be registered to it; these have to be derived from the FairGenerator class. The PndFilteredPrimaryGenerator defines position and (optionally) smearing of the primary vertex. This class should be instantiated only once.

Definition at line 74 of file PndFilteredPrimaryGenerator.h.

Constructor & Destructor Documentation

PndFilteredPrimaryGenerator::PndFilteredPrimaryGenerator ( )

Default constructor.

Definition at line 15 of file PndFilteredPrimaryGenerator.cxx.

References fCodeNameMap, fCombFsPdg, fdbPdg, fNameCodeMap, fNamePdg, fPartNames, fSetFsPdg, and i.

 : FairPrimaryGenerator(),
 
   fEvtFilterStat(FairEvtFilterParams()),
   fVerbose(3),
   fEventNrFiltered(0),
   fEventPrintFreq(100)
 {
         fdbPdg = TDatabasePDG::Instance();
         
         std::vector<TString>    tmpfPartNames = {"e+", "e-", "mu+", "mu-", "pi+", "pi-", "k+", "k-", "p+",  "p-", "n0", "n0b", "e+-",   "mu+-",   "pi+-",   "k+-",    "p+-",   "n00b"  , "t+", "t-", "t+-", "nt", "any" };
         std::vector<int>        tmpfNamePdg   = {-11 , 11  , -13  , 13   , 211  , -211 , 321 , -321, 2212, -2212, 2112, -2112, 11000  , 13000   , 211000  , 321000 ,  2212000, 2112000,   1   , -1  ,  2    , 0   ,  3   };
         std::vector<int>        tmpfCombFsPdg = {-11 , 11  , -13  , 13   , 211  , -211 , 321 , -321, 2212, -2212 , 22};
         
         // initialize sets and vectors with names and codes
         fPartNames  = tmpfPartNames;
         fNamePdg    = tmpfNamePdg;
         fCombFsPdg  = tmpfCombFsPdg;
 
         std::unordered_set<int> tmpfSetFsPdg(fCombFsPdg.begin(), fCombFsPdg.end());
         fSetFsPdg   = tmpfSetFsPdg;
         
         // initialize maps between names and codes
         for (int i=0;i<(int)fPartNames.size();++i) fNameCodeMap[fPartNames[i]] = fNamePdg[i];
         for (int i=0;i<(int)fPartNames.size();++i) fCodeNameMap[fNamePdg[i]]   = fPartNames[i];
 }

PndFilteredPrimaryGenerator::PndFilteredPrimaryGenerator ( TString filterset )

Constructor with filter-set string.

Definition at line 44 of file PndFilteredPrimaryGenerator.cxx.

References AddFilter(), fCodeNameMap, fCombFsPdg, fdbPdg, fNameCodeMap, fNamePdg, fPartNames, fSetFsPdg, and i.

 : FairPrimaryGenerator(),
 
   fEvtFilterStat(FairEvtFilterParams()),
   fVerbose(3),
   fEventNrFiltered(0),
   fEventPrintFreq(100)
 {
         fdbPdg = TDatabasePDG::Instance();
         
         std::vector<TString>    tmpfPartNames = {"e+", "e-", "mu+", "mu-", "pi+", "pi-", "k+", "k-", "p+",  "p-", "n0", "n0b", "e+-",   "mu+-",   "pi+-",   "k+-",    "p+-",   "n00b"  , "t+", "t-", "t+-", "nt", "any" };
         std::vector<int>        tmpfNamePdg   = {-11 , 11  , -13  , 13   , 211  , -211 , 321 , -321, 2212, -2212, 2112, -2112, 11000  , 13000   , 211000  , 321000 ,  2212000, 2112000,   1   , -1  ,  2    , 0   ,  3   };
         std::vector<int>        tmpfCombFsPdg = {-11 , 11  , -13  , 13   , 211  , -211 , 321 , -321, 2212, -2212 , 22};
         
         // initialize sets and vectors with names and codes
         fPartNames  = tmpfPartNames;
         fNamePdg    = tmpfNamePdg;
         fCombFsPdg  = tmpfCombFsPdg;
         
         std::unordered_set<int> tmpfSetFsPdg(fCombFsPdg.begin(), fCombFsPdg.end());
         fSetFsPdg   = tmpfSetFsPdg;
 
         // initialize maps between names and codes
         for (int i=0;i<(int)fPartNames.size();++i) fNameCodeMap[fPartNames[i]] = fNamePdg[i];
         for (int i=0;i<(int)fPartNames.size();++i) fCodeNameMap[fNamePdg[i]]   = fPartNames[i];
         
         AddFilter(filterset);
 }

virtual PndFilteredPrimaryGenerator::~PndFilteredPrimaryGenerator ( )

inlinevirtual

Destructor.

Definition at line 86 of file PndFilteredPrimaryGenerator.h.

86 {};

PndFilteredPrimaryGenerator::PndFilteredPrimaryGenerator ( const PndFilteredPrimaryGenerator & )

private

Member Function Documentation

void PndFilteredPrimaryGenerator::AddFilter ( TString filterStr )

Registers a filter as a string to be parsed Each filter set consist of one filter definition or a number of filters connected with logical && (AND) (= a certain event type). All filter sets are combined with a logical || (OR) (all event types to be accepted). Syntax: "Filt_1.1 && Filt_1.2 || Filt_2.1 || Filt_3.1 && Filt_3.2 && Filt_3.2".

Definition at line 157 of file PndFilteredPrimaryGenerator.cxx.

References PndSmpFilt::compo, cuts, exit(), f, fFilterSets, fNameCodeMap, for(), fVerbose, GetRangeDouble(), GetRangeInt(), i, if(), jj, PndSmpFilt::mcntr, PndSmpFilt::mwin, PndSmpFilt::name, PndSmpFilt::ndau, PndSmpFilt::nmax, PndSmpFilt::nmin, PndSmpFilt::nocc, PndSmpFilt::pdg, PndSmpFilt::phimax, PndSmpFilt::phimin, PndSmpFilt::pmax, PndSmpFilt::pmin, PndSmpFilt::ptmax, PndSmpFilt::ptmin, PndSmpFilt::pzmax, PndSmpFilt::pzmin, SplitString(), PndSmpFilt::thtmax, PndSmpFilt::thtmin, TString, and PndSmpFilt::veto.

Referenced by PndFilteredPrimaryGenerator(), and prod_fsim().

 {
         // some replacements
         filterStr.ToLower();
         filterStr.ReplaceAll(":","&&");
         filterStr.ReplaceAll("gamma","nt");
         filterStr.ReplaceAll("gam","nt");
         filterStr.ReplaceAll("ch0","nt");
         filterStr.ReplaceAll("tr","t");
         filterStr.ReplaceAll("ch","t");
         filterStr.ReplaceAll("mass","m");
 
 
         // split filter string
         StrVec filtersets_str = SplitString(filterStr,"||");
         
         // *** loop over event type filters (the filter sets)
         for (int ifs=0; ifs < (int) filtersets_str.size(); ++ ifs)
         {
                 StrVec filters_str = SplitString(filtersets_str[ifs],"&&");
                 
                 // --------------------
                 // parse filter strings
                 // --------------------
                 
                 PndSmpFilterSet filters;
                 
                 for (int i=0;i<(int) filters_str.size();++i)
                 {
                         // some filter checks
                         TString fs = filters_str[i];
                         if ( (!fs.EndsWith(")")) || (fs.Contains(",") && !fs.Contains("[")) ) 
                         {
                                 cout <<" --> Please check filter string '"<<fs<<"'"<<endl;
                                 exit(0);
                         }
                         
                         // a copy without the X() stuff
                         TString tmpfs = fs(fs.Index("(")+1,fs.Length()-fs.Index("(")-2);
                         
                         StrVec cuts = SplitString(tmpfs,";");
 
                         if (fVerbose>3) cout<<endl<<fs<<endl<<tmpfs<<endl;
                         
                         PndSmpFilt f;
                         f.name = fs;
 
                         // negate filter?
                         if (fs.BeginsWith("!")) {f.veto = true; fs.ReplaceAll("!","");}
                         
                         // if the filter string has the form "M(...)" it's a mass filter, else a multiplicity filter ("(...)") 
                         if (fs.BeginsWith("m(")) f.compo=true;
                         
                         for (int j=0;j<(int)cuts.size();++j)
                         {
                                 TString ct = cuts[j];
                                 //cout <<"  "<<ct<<endl;
 
                                 // if specifier (p, tht, pt, ...) store it and remove the '[]'
                                 TString key = "";
                                 if (ct.Contains("[")) 
                                 {
                                         key = ct(0,ct.Index("["));
                                         ct = ct(ct.Index("[")+1,ct.Length()-ct.Index("[")-2);
                                 }
 
                                 // either multiplicity or particle names
                                 if (key=="")
                                 {
                                         // mult, either range delimiter '..' found
                                         if      (ct.Contains("..")) GetRangeInt(ct, f.nmin, f.nmax);
                                         // or single number, explicitely allowed to be 0 (kind of veto filtering)
                                         else if (ct.Atoi()!=0 || ct=="0")      f.nmin = f.nmax = ct.Atoi();
                                         // particle name of the particle to be counted or a list of particles for inv. mass filter (in case f.compo = true)
                                         else
                                         {
                                                 StrVec parts = SplitString(ct);
                                                 for (int jj=0;jj<(int)parts.size();++jj)
                                                 {
                                                         if (fNameCodeMap.find(parts[jj])!=fNameCodeMap.end() && f.ndau<5) f.pdg[f.ndau++] = fNameCodeMap[parts[jj]];
                                                         
                                                         // the key word 'nocc' prevents automatic charged conjugates for 
                                                         if (parts[jj]=="nocc") f.nocc = true;
                                                 }
                                                 for (int jj=0;jj<(int)parts.size();++jj) if (f.pdg[jj]==0) f.pdg[jj] = 22; 
                                         }
                                 }
                                 
                                 // other pdg code
                                 else if (key=="pdg")  {f.pdg[0] = ct.Atoi(); f.ndau=1;}
                                                                 
                                 // kinematic ranges
                                 else if (key=="p")    GetRangeDouble(ct, f.pmin,   f.pmax);
                                 else if (key=="pt")   GetRangeDouble(ct, f.ptmin,  f.ptmax);
                                 else if (key=="pz")   GetRangeDouble(ct, f.pzmin,  f.pzmax);
                                 else if (key=="tht")  GetRangeDouble(ct, f.thtmin, f.thtmax);
                                 else if (key=="phi")  GetRangeDouble(ct, f.phimin, f.phimax, ",", 1);
                                 else if (key=="m")    GetRangeDouble(ct, f.mcntr,  f.mwin);
                         }
                         
                         if ( (f.compo && f.ndau<2) || (!f.compo && f.ndau!=1))
                         {
                                 cout <<" --> Please check filter string '"<<fs<<"'"<<endl;
                                 exit(0);                                
                         }
                         else filters.push_back(f);
                 }
                 
                 fFilterSets.push_back(filters); 
         }
 }

int PndFilteredPrimaryGenerator::AntiPdgCode ( int pdg )

protected

Gets anti-pdg code, if exists. If not returns the code itself (particle is its anti-particle)

Definition at line 122 of file PndFilteredPrimaryGenerator.cxx.

References fdbPdg.

Referenced by GenerateEvent().

 {
         if (fdbPdg->GetParticle(pdg)!=0x0 && fdbPdg->GetParticle(pdg)->AntiParticle()!=0x0) 
                 return fdbPdg->GetParticle(pdg)->AntiParticle()->PdgCode(); 
         
         return pdg;
 }

bool PndFilteredPrimaryGenerator::CheckKinematic	(	const PndSmpFilt &	f,
		const TLorentzVector &	p4
	)

protected

Checks whether P4 kinematics match the criteria of a PndSmpFilt.

Definition at line 107 of file PndFilteredPrimaryGenerator.cxx.

References PndSmpFilt::phimax, PndSmpFilt::phimin, PndSmpFilt::pmax, PndSmpFilt::pmin, PndSmpFilt::ptmax, PndSmpFilt::ptmin, PndSmpFilt::pzmax, PndSmpFilt::pzmin, PndSmpFilt::thtmax, and PndSmpFilt::thtmin.

Referenced by GenerateEvent().

 {
         bool acc_trk = true;
         
         double curr_p   = p4.P(); acc_trk = (curr_p   >= f.pmin   && curr_p   <= f.pmax);       
         if (acc_trk)  { double curr_pt  = p4.Pt();            acc_trk = (curr_pt  >= f.ptmin  && curr_pt  <= f.ptmax);  } 
         if (acc_trk)  { double curr_pz  = p4.Pz();            acc_trk = (curr_pz  >= f.pzmin  && curr_pz  <= f.pzmax);  }
         if (acc_trk)  { double curr_tht = p4.Theta()*57.2958; acc_trk = (curr_tht >= f.thtmin && curr_tht <= f.thtmax); }
         if (acc_trk)  { double curr_phi = p4.Phi()*57.2958;   acc_trk = (curr_phi >= f.phimin && curr_phi <= f.phimax); }       
         
         return acc_trk;
 } 

PndFilteredPrimaryGenerator::ClassDef	(	PndFilteredPrimaryGenerator	,
		1
	)

private

PndSmpCandList PndFilteredPrimaryGenerator::CombineList	(	int	pdg,
		PndSmpCandList *	l0,
		PndSmpCandList *	l1,
		PndSmpCandList *	l2 = `0`,
		PndSmpCandList *	l3 = `0`,
		PndSmpCandList *	l4 = `0`
	)

protected

Combines upt to five particle lists of PndSmpCand with overlap and double counting prevention.

Definition at line 272 of file PndFilteredPrimaryGenerator.cxx.

References c1, c2, c3, c4, PndSmpCand::Overlap(), and res.

Referenced by GenerateEvent().

 {
         PndSmpCandList res;
 
         //-------------
         for (int i0=0; i0<(int)l0->size(); ++i0)
         {
                 PndSmpCand &c0 = (*l0)[i0];
                 
                 // if l0 == l1, start at i1 = i0+1 to avoid double counting
                 int st1 = 0;
                 if ((*l1) == (*l0)) st1 = i0+1;         
                 
                 // -------------
                 for (int i1=st1; i1<(int)l1->size(); ++i1)
                 {
                         PndSmpCand &c1 = (*l1)[i1];
                         if (c1.Overlap(c0)) continue;
                         
                         // only two lists? done!
                         if (l2==0) res.push_back(PndSmpCand(pdg, c0, c1));
                         else
                         {
                                 // if l2 == l0 or l1, start at i2 = i0/i1+1 to avoid double counting
                                 int st2 = 0;
                                 if (*l2==*l1) st2 = i1+1;
                                 else if (*l2==*l0) st2 = i0+1;
                                 
                                 // -------------
                                 for (int i2=st2; i2<(int)l2->size(); ++i2)
                                 {
                                         PndSmpCand &c2 = (*l2)[i2];
                                         if (c2.Overlap(c1) || c2.Overlap(c0)) continue;
                                         
                                         // only three lists? done!
                                         if (l3==0) res.push_back(PndSmpCand(pdg, c0, c1, c2));
                                         else
                                         {
                                                 // if l3 == l0, l1 or l2, start at i3 = i0/i1/i2+1 to avoid double counting
                                                 int st3 = 0;
                                                 if (*l3==*l2) st3 = i2+1; 
                                                 else if (*l3==*l1) st3 = i1+1; 
                                                 else if (*l3==*l0) st3 = i0+1;
                                                 
                                                 // -------------
                                                 for (int i3=st3; i3<(int)l3->size(); ++i3)
                                                 {
                                                         PndSmpCand &c3 = (*l3)[i3];
                                                         if (c3.Overlap(c2) || c3.Overlap(c1) || c3.Overlap(c0)) continue;
                                                         
                                                         // only four lists? done!
                                                         if (l4==0) res.push_back(PndSmpCand(pdg, c0, c1, c2, c3));
                                                         else
                                                         {
                                                                 // if l3 == l0, l1 or l2, start at i3 = i0/i1/i2+1 to avoid double counting
                                                                 int st4 = 0;
                                                                 if (*l4==*l3) st4 = i3+1; 
                                                                 else if (*l4==*l2) st3 = i2+1; 
                                                                 else if (*l4==*l1) st3 = i1+1;
                                                                 else if (*l4==*l0) st3 = i0+1;
                                                                 
                                                                 // -------------
                                                                 for (int i4=st4; i4<(int)l4->size(); ++i4)
                                                                 {
                                                                         PndSmpCand &c4 = (*l4)[i4];
                                                                         if (c4.Overlap(c3) || c4.Overlap(c2) || c4.Overlap(c1) || c4.Overlap(c0)) continue;
                                                                         
                                                                         res.push_back(PndSmpCand(pdg, c0, c1, c2, c3, c4));
                                                                 } // for i4
                                                         }  // else i3
                                                 }  // for i3
                                         }  // else i2
                                 }  // for i2            
                         }  // else i1
                 }  // for i1    
         }  // for i0
         
         return res;
 } 

Bool_t PndFilteredPrimaryGenerator::GenerateEvent ( FairGenericStack * pStack )

virtual

Calls event generators and the event filters.

To be called at the beginning of each event from FairMCApplication. Generates an event vertex and calls the ReadEvent methods from the registered generators. Calls defined event (veto) filters to decide whether to process the event or to call the event generators again.

Parameters

pStack The particle stack

Returns: kTRUE if successful, kFALSE if not

Definition at line 364 of file PndFilteredPrimaryGenerator.cxx.

References all, AntiPdgCode(), CheckKinematic(), cnt, CombineList(), PndSmpFilt::compo, f, fabs(), fCodeNameMap, fCombFsPdg, fdbPdg, fEventNrFiltered, fEventPrintFreq, fEvtFilterStat, FairEvtFilterParams::fFailedFilterEvents, FairEvtFilterParams::fFilterMaxTries, fFilterSets, FairEvtFilterParams::fGeneratedEvents, fSetFsPdg, fVerbose, i, PndSmpFilt::mcntr, PndSmpFilt::mwin, nc, PndSmpFilt::ndau, PndSmpFilt::nmax, PndSmpFilt::nocc, particle, PndSmpFilt::pdg, PrintSmpCandList(), PndSmpCand::SetP4(), and PndSmpFilt::veto.

 {
         Int_t iTry = 0; // number of attempts to find the next event that suits your filter
         
         bool acc_glob = false; // is kTRUE if the event is finally accepted
 
         bool active_filters = (fFilterSets.size()>0);
 
         while( !acc_glob && iTry < fEvtFilterStat.fFilterMaxTries)
         {
                 ++iTry; // count how often we run the generators before we accept an event
                 ++fEvtFilterStat.fGeneratedEvents; // total number of generated events
 
                 pStack->Reset(); // Clean the stack
                 FairPrimaryGenerator::GenerateEvent(pStack); // fill the stack
                 
                 // no filtering --> accept event
                 if (!active_filters) {acc_glob = true; continue;}
                 
                 TClonesArray* fParticleList = pStack->GetListOfParticles();
                 
                 
                 fEventNr = fEventNrFiltered; // Fix event numbering in FairPrimaryGenerator (otherwise fRun will stop too early...)
 
                 PndSmpCandList all;
                 
                 // init all lists to empty list
                 std::map<int,PndSmpCandList> pdgList;
                 for (unsigned int i=0;i<fCombFsPdg.size();++i)   pdgList[fCombFsPdg[i]] = all;
 
                 int cnt = 0;
                 
                 // --------------
                 // Fill the lists
                 // --------------
                 
                 //for (int i=0;i<*nTrk;++i)
                 //{
                 for (Int_t iPart=0; iPart<fParticleList->GetEntries(); ++iPart) 
                 {
                         TParticle *particle = (TParticle*)fParticleList->At(iPart);
                         
                         // the current MCT pdg code
                         Int_t pdg = particle->GetPdgCode();
 
                         // do we want to count multiplicity for this code?
                         //if (set_mult_pdg.find(pdg)==set_mult_pdg.end()) continue;
                         
                         // is yes, create SmpCand
                         TLorentzVector l(particle->Px(), particle->Py(), particle->Pz(), particle->Energy());
                         Float_t ch = fdbPdg->GetParticle(pdg) ? fdbPdg->GetParticle(pdg)->Charge()/3. : 0;
                         
                         // and push to to all-list
                         if (fSetFsPdg.find(pdg)==fSetFsPdg.end())
                         {
                                 // push to all list, but with uid = -1 for non-final states, since not needed for combinatorics
                                 all.push_back(PndSmpCand(l, ch, pdg, -1));
                                 continue;
                         }
 
                         // create candidate with marker (uid >= 0) only for final states used for combinatorics
                         PndSmpCand sc(l, ch, pdg, cnt++);
                         
                         // and push to to all-list
                         all.push_back(sc);
                         
                         // add particle to all lists with the same charge
                         for (unsigned int j=0;j<fCombFsPdg.size();++j)
                         {
                                 int    pdg_ch = fdbPdg->GetParticle(fCombFsPdg[j])->Charge()/3.;
                                 double pdg_m  = fdbPdg->GetParticle(fCombFsPdg[j])->Mass();
                                 
                                 if (pdg_ch == (int)ch) 
                                 {
                                         l.SetVectM(l.Vect(), pdg_m);
                                         sc.SetP4(l);
                                         pdgList[fCombFsPdg[j]].push_back(sc);
                                 }
                         }               
                 }
 
                 // --------------
                 // print all lists
                 // --------------
                 
                 if (fVerbose>3)
                 {
                         PrintSmpCandList(all,"All particles:");
                         
                         if (fVerbose>4)
                         for (int i=0;i<(int)fCombFsPdg.size();++i) 
                         {
                                 int pdg = fCombFsPdg[i];
                                 int nc = pdgList[pdg].size();
                                 if (nc>0) PrintSmpCandList(pdgList[pdg],Form("%s",fCodeNameMap[pdg].Data()));
                         }
                         cout <<endl<<endl;
                 }
                                 
                 
                 // -------------------
                 // Loop over filters
                 // -------------------
 
                 acc_glob=false;
                 
                 for (int ifs=0; ifs < (int) fFilterSets.size(); ++ifs)
                 {
                         // did another filterset accept? then we don't need to check anymore
                         if (acc_glob) continue;
                         
                         // initialize filter set accept with true
                         bool acc_fset = true;
 
                         // -------------------
                         // *** loop over count filters in filter set
                         // -------------------
                         for (int i=0; i<(int)fFilterSets[ifs].size(); ++i)
                         {
                                 // if the event cannot be accepted anymore, skip the next filter test
                                 if (!acc_fset) continue;
 
                                 // event accept flag
                                 bool acc_evt = true;
                                 
                                 PndSmpFilt &f = fFilterSets[ifs][i];
 
                                 // if mass filter, skip here (first apply count filters, which are probably much faster)
                                 if (f.compo) continue;
                                 
                                 int cnt_matches = 0;
                                 
                                 // *** loop over particles
                                 for (int j=0 ; j<(int)all.size()  ; ++j)
                                 {
                                         // *** check particle type first
                                         //vector<TString> names = {"e+", "e-", "mu+", "mu-", "pi+", "pi-", "k+", "k-", "p+",  "p-", "n0", "n0b", "e+-", "mu+-", "pi+-", "k+-",  "p+-", "n"  , "t+", "t-", "t+-", "nt", "any" };
                                         //vector<int> name_pdg  = {-11 , 11  , -13  , 13   , 211  , -211 , 321 , -321, 2212, -2212, 2112, -2112, 11000  , 13000   , 211000  , 321000 ,  2212000, 2112000, 1    , -1   ,  2    ,  0   ,  3   };
                                         
                                         // tracks pdg and charge
                                         int trpdg = all[j].Pdg();       
                                         int trchg = (int) all[j].Charge(); 
                                         
                                         // filter pdg (or special) code
                                         int flpdg = f.pdg[0];
 
                                         // the counters for any, t+-, t+, t- and nt/gam only for final states
                                         if (all[j].Marker()>0)
                                         {
                                                 // any particle: code = 3
                                                 if (flpdg == 3)   trpdg = flpdg;
 
                                                 // all charged tracks t+- : code = 2
                                                 if (flpdg == 2)   trpdg = abs(trchg)*2;
 
                                                 // simple track+- and neutral counting (t+,t-,nt) : code = charge (-1, +1, 0); 0 also used for photons with pdg = 22 
                                                 if (abs(flpdg)<2) trpdg = trchg;
                                         }
                                         
                                         // the pdg code ignoring charge : code = 1000*pdg, e.g. K+- = 321000, pi+- = 211000
                                         if (flpdg%1000==0) trpdg = abs(trpdg)*1000;
                                         
                                         // *** check code and kinematic variables and count track
                                         if ( (trpdg==flpdg)  &&  CheckKinematic(f,  all[j].P4()) )
                                         {
                                                 cnt_matches++;
                                                 if (fVerbose>2) cout  << "\033[1;34m ** ";
                                         }
                                         else if (fVerbose>2) cout << "    ";
                                         
                                         // *** (and print if verbose)
                                         if (fVerbose>2) {all[j].Print(); cout<<"\033[0;0m";}
                                 }
                                 
                                 // does multiplicity match?
                                 if (cnt_matches<f.nmin || cnt_matches>f.nmax) acc_evt = false;
                                 
                                 // is this a veto filter -> negate result
                                 acc_evt = acc_evt^f.veto;
                                 
                                 if (fVerbose>2) cout <<"Filterset "<<ifs<<" Filter "<<i<<" : N_matched = "<<cnt_matches<<" : event "<<(!acc_evt?"not ":"")<<"accepted"<<endl<<endl<<endl; 
 
                                 // global accept (of all filters)
                                 acc_fset = (acc_fset && acc_evt);               
                                 
                         } // **** count filters
                         
                         
                         // -------------------
                         // Loop over mass filters
                         // -------------------
                         for (int i=0; i<(int)fFilterSets[ifs].size(); ++i)
                         {
                                 // if the event cannot be accepted anymore, skip the next filter test
                                 if (!acc_fset) continue;
 
                                 // event accept flag
                                 bool acc_evt = true;
                                 
                                 PndSmpFilt &f = fFilterSets[ifs][i];
 
                                 // now only mass filters
                                 if (!f.compo) continue;
                                                 
                                 // ----------------
                                 // do combinatorics
                                 // ----------------
 
                                 // prepare lists of particles and the charged conjugation
                                 PndSmpCandList *l[5] = {0}, *al[5] = {0};
                                 std::vector<int> vpdg, vapdg;
                                 
                                 for (int j=0; j<f.ndau; ++j)
                                 {
                                         vpdg.push_back(f.pdg[j]); 
                                         vapdg.push_back(AntiPdgCode(f.pdg[j]));
                                         
                                         l[j]  = &(pdgList[f.pdg[j]]);
                                         al[j] = &(pdgList[AntiPdgCode(f.pdg[j])]);
                                 }
                                                                 
                                 // *** do combinatorics
                                 PndSmpCandList comblist = CombineList(0, l[0], l[1], l[2], l[3], l[4]);
 
                                 // do we want c.c. and is the cc'd list not the same as the original one (= a permutation)?
                                 if (!f.nocc && !is_permutation(vpdg.begin(), vpdg.end(), vapdg.begin())) 
                                 {
                                         PndSmpCandList anticomblist = CombineList(0, al[0], al[1], al[2], al[3], al[4]);
                                         comblist.insert(comblist.end(), anticomblist.begin(), anticomblist.end());
                                 }
                                 
                                 int cnt_matches = 0;
                                 
                                 // *** loop over composite particles
                                 for (int j=0 ; j<(int)comblist.size()  ; ++j)
                                 {
                                         TLorentzVector p4 = comblist[j].P4();
                                         
                                         // *** check mass window and kinematics of composite
                                         bool acc_cand = (fabs(p4.M() - f.mcntr) < f.mwin/2.)  &&  CheckKinematic(f, p4);   
 
                                         // *** count track
                                         if (acc_cand)
                                         {
                                                 cnt_matches++;
                                                 if (fVerbose>2) cout  << "\033[1;35m ** ";
                                         }
                                         else if (fVerbose>2) cout << "    ";
                                         
                                         // *** (and print if verbose)
                                         if (fVerbose>2) { comblist[j].Print(); cout<<"\033[0;0m";}
                                 }               
                                 
                                 // does multiplicity match?
                                 if (cnt_matches<f.nmin || cnt_matches>f.nmax) acc_evt = false;
                                 
                                 // is this a veto filter -> negate result
                                 acc_evt = acc_evt^f.veto;
                                 
                                 if (fVerbose>2) cout <<"Filterset "<<ifs<<" Filter "<<i<<" : N_matched = "<<cnt_matches<<" : event "<<(!acc_evt?"not ":"")<<"accepted"<<endl<<endl<<endl; 
 
                                 // global accept (of all filters)
                                 acc_fset = (acc_fset && acc_evt);               
                                 
                         } // **** mass filters
                         
                         
                         // the global event accept, being an OR connection of all filter sets ('event types' to be accepted)
                         
                         acc_glob = (acc_glob || acc_fset);
                 } // **** filter sets
                 
                 
                 //if (fVerbose>1) cout <<" --> EVENT "<<entry<<" "<<(!acc_glob?"NOT ":"")<<"ACCEPTED"<<endl<<endl<<endl; 
                 
                 //if (acc_glob) {cnt_acc++;}            
         } 
 
 
         // Set the event number ALWAYS when filtering
         ++fEventNrFiltered;
         fEvent->SetEventID(fEventNrFiltered);
 
         if ( !acc_glob ){
                 ++fEvtFilterStat.fFailedFilterEvents;
                 cout << "\n -E PndFilteredPrimaryGenerator: No event was found within " << iTry << " tries which satisfies your event filter.\n ";
                 cout << "I accept a random event as evtNr " << fEventNrFiltered << " to avoid infinite loops. \n";
                 cout <<  "Try increasing the max. number of tries or change your filter\n\n";
                 if (fVerbose > 3 ){
                         cout << iTry << " events simulated until I found a good one.\n";
                         cout << fEvtFilterStat.fGeneratedEvents << " events generated for finding " << fEventNr << " accepted events.\n";
                         cout << fEvtFilterStat.fFailedFilterEvents << " unsuccessful attempts in total to find an event that suits your filters\n\n";
                 }
         }
 
         if (0 < fVerbose && active_filters && (fEventNrFiltered%fEventPrintFreq)==0) cout <<"[PndFilteredPrimaryGenerator] " << fEventNrFiltered << " / " << fEvtFilterStat.fGeneratedEvents << " generated events accepted.\n";
 
         return kTRUE;
 }

Int_t PndFilteredPrimaryGenerator::GetNumberOfFilterFailedEvents ( )

inline

Returns the number of cases in which no matching event was found within the set max. tries.

This method returns 0 if everything works fine. If it returns a value >0 it means that you should set a higher limit in SetFilterMaxTries. If it returns a value which is equal to the number of events that you requested, it means that either the max. number of tries is set way too low or that the generator does not create such events that you are interested in or that your event filters cannot be satisfied at all (logical error).

Definition at line 145 of file PndFilteredPrimaryGenerator.h.

References fEvtFilterStat, and FairEvtFilterParams::fFailedFilterEvents.

Referenced by WriteEvtFilterStatsToRootFile().

                                              {
                 return fEvtFilterStat.fFailedFilterEvents;
         }

Int_t PndFilteredPrimaryGenerator::GetNumberOfFilterMaxTries ( )

inline

returns the maximum number of times that this object should try to find an event which suits all event filters.

Definition at line 122 of file PndFilteredPrimaryGenerator.h.

References fEvtFilterStat, and FairEvtFilterParams::fFilterMaxTries.

                                          {
                 return fEvtFilterStat.fFilterMaxTries;
         }

Int_t PndFilteredPrimaryGenerator::GetNumberOfGeneratedEvents ( )

inline

returns the total (accepted + rejected) number of events generated by the event generators. If no event filters are used this number is equal to the number of simulated events.

Definition at line 127 of file PndFilteredPrimaryGenerator.h.

References fEvtFilterStat, and FairEvtFilterParams::fGeneratedEvents.

Referenced by prod_fsim(), and WriteEvtFilterStatsToRootFile().

                                           {
                 return fEvtFilterStat.fGeneratedEvents;
         }

void PndFilteredPrimaryGenerator::GetRangeDouble	(	TString	s,
		double &	a,
		double &	b,
		TString	delim = `","`,
		bool	forceset = `false`
	)

protected

Turns a string of the form <float><delim><float> (e.g. "124.2,178.3") to a range a,b.

Definition at line 132 of file PndFilteredPrimaryGenerator.cxx.

References s, and TString.

Referenced by AddFilter().

 {
         int len = delim.Length();
         a = (TString(s(0,s.Index(delim)))).Atof();
         double tmpb = (TString(s(s.Index(delim)+len,10000))).Atof();
         if (tmpb>0. || forceset) b = tmpb;
 }

void PndFilteredPrimaryGenerator::GetRangeInt	(	TString	s,
		int &	a,
		int &	b,
		TString	delim = `".."`
	)

protected

Turns a string of the form <int><delim><int> (e.g. "3..8") to a range a,b.

Definition at line 142 of file PndFilteredPrimaryGenerator.cxx.

References s, and TString.

Referenced by AddFilter().

 {
         int len = delim.Length();
         a = (TString(s(0,s.Index(delim)))).Atoi();
         int tmpb = (TString(s(s.Index(delim)+len,10000))).Atoi();
         if (tmpb>0) b=tmpb;
 }

Bool_t PndFilteredPrimaryGenerator::Init ( )

virtual

Initialize the event generator(s) and the event (veto) filter(s).

Definition at line 75 of file PndFilteredPrimaryGenerator.cxx.

References fFilterSets, fVerbose, i, Init(), and Print().

 {
         FairPrimaryGenerator::Init();
 
         if (fVerbose) 
                 for (int ifs = 0; ifs < (int)fFilterSets.size(); ++ifs)
                 {
                         cout <<"FILTER SET "<< ifs<<endl<<endl;
                         for (int i = 0; i < (int)fFilterSets[ifs].size(); ++i) fFilterSets[ifs][i].Print();
                         cout <<endl;
                 }
 
         return kTRUE;
 }

PndFilteredPrimaryGenerator& PndFilteredPrimaryGenerator::operator= ( const PndFilteredPrimaryGenerator & )

private

void PndFilteredPrimaryGenerator::PrintSmpCandList	(	PndSmpCandList	l,
		TString	name = `""`
	)

protected

Prints a candidate lits.

Definition at line 354 of file PndFilteredPrimaryGenerator.cxx.

References i, and printf().

Referenced by GenerateEvent().

 {
         printf("LIST %s with %d candidates\n",name.Data(), (int)l.size());
         for (int i=0;i<(int)l.size(); ++i) {printf(" %2d : ",i); l[i].Print();}
 }

void PndFilteredPrimaryGenerator::SetEventPrintFrequency ( int freq )

inline

Sets the frequency (accepted events) for printout (verbose>0) of accepted and generated events.

Definition at line 132 of file PndFilteredPrimaryGenerator.h.

References fEventPrintFreq.

Referenced by prod_fsim().

         {
                 fEventPrintFreq = freq;
         }

void PndFilteredPrimaryGenerator::SetFilterMaxTries ( Int_t maxTries = 99999 )

inline

Define the maximum number of times that this object should try to find an event which suits all event filters.

Definition at line 111 of file PndFilteredPrimaryGenerator.h.

References fEvtFilterStat, and FairEvtFilterParams::fFilterMaxTries.

         {
                 if ( maxTries > 0 ){
                         fEvtFilterStat.fFilterMaxTries = maxTries;
                         std::cout << "PndFilteredPrimaryGenerator: maxTries is now set to " << fEvtFilterStat.fFilterMaxTries << "\n";
                 } else {
                         std::cout << "\n\n\n  -WARNING from PndFilteredPrimaryGenerator: maxTries must be a positive number! Check your SetFilterMaxTries call!\n\n\n";
                 }
         }

void PndFilteredPrimaryGenerator::SetVerbose ( Int_t verbose = 12 )

inline

Set the level of commenting output.

Parameters

verbose Level of commenting output, 0 means no output, higher gives more output.

Definition at line 168 of file PndFilteredPrimaryGenerator.h.

References fVerbose, and verbose.

Referenced by prod_fsim().

         {
                 if ( verbose >= 0 )
                 {
                         fVerbose = verbose;
                         std::cout << "PndFilteredPrimaryGenerator: fVerbose is now set to " << verbose << "\n";
                 }
                 else 
                 {
                         std::cout << "\n\n\n  -WARNING from PndFilteredPrimaryGenerator: verbose must be a positive number! Check your SetVerbose call!\n\n\n";
                 }
         }

StrVec PndFilteredPrimaryGenerator::SplitString	(	TString	s,
		TString	delim = `" "`
	)

protected

Splits a TString to substrings.

Definition at line 92 of file PndFilteredPrimaryGenerator.cxx.

References i, and TString.

Referenced by AddFilter().

 {
         StrVec toks;
         TObjArray *tok = s.Tokenize(delim);
         int N = tok->GetEntries();      
         for (int i=0;i<N;++i) 
         {
                 TString token = (((TObjString*)tok->At(i))->String()).Strip(TString::kBoth);
                 toks.push_back(token);
         }
         return toks;
 }

void PndFilteredPrimaryGenerator::WriteEvtFilterStatsToRootFile ( TFile * outputFile = NULL )

inline

Writes all relevant event filter information to the output root file.

Definition at line 150 of file PndFilteredPrimaryGenerator.h.

References fEvtFilterStat, GetNumberOfFilterFailedEvents(), and GetNumberOfGeneratedEvents().

Referenced by prod_fsim().

                                                                   {
                 std::cout << "\n\nGenerated Events = " << GetNumberOfGeneratedEvents() << "\n";
                 if (0 < GetNumberOfFilterFailedEvents() ) {
                         std::cout << "WARNING: Number of events where the event filter FAILED " << GetNumberOfFilterFailedEvents() << "\n\n\n";
                         std::cout << "Random events were accepted to avoid infinite loops. \n";
                         std::cout <<  "Try increasing the max. number of tries or change your filter (maybe the generators do not produce such events as you want).\n\n";
                 }
                 if(outputFile==NULL) outputFile = FairRunSim::Instance()->GetOutputFile();
                 outputFile->cd();
                 outputFile->mkdir("FairEvtFilter");
                 outputFile->cd("FairEvtFilter");
                 fEvtFilterStat.Write();
                 outputFile->cd();
         }

Member Data Documentation

std::map<int, TString> PndFilteredPrimaryGenerator::fCodeNameMap

protected

mapes (pdg) codes to names

Definition at line 244 of file PndFilteredPrimaryGenerator.h.

Referenced by GenerateEvent(), and PndFilteredPrimaryGenerator().

std::vector<int> PndFilteredPrimaryGenerator::fCombFsPdg

protected

particle codes for the lists used for combinatorics; !!! the codes are not selected from MC truth, just the corresponding mass hypos are applied to charged tracks

Definition at line 235 of file PndFilteredPrimaryGenerator.h.

Referenced by GenerateEvent(), and PndFilteredPrimaryGenerator().

TDatabasePDG* PndFilteredPrimaryGenerator::fdbPdg

protected

Shortcut to TDatabasePDG.

Definition at line 226 of file PndFilteredPrimaryGenerator.h.

Referenced by AntiPdgCode(), GenerateEvent(), and PndFilteredPrimaryGenerator().

Int_t PndFilteredPrimaryGenerator::fEventNrFiltered

protected

Event number (Set by the filtered primary generator.

Definition at line 220 of file PndFilteredPrimaryGenerator.h.

Referenced by GenerateEvent().

Int_t PndFilteredPrimaryGenerator::fEventPrintFreq

protected

Print frequency for filtered events.

Definition at line 223 of file PndFilteredPrimaryGenerator.h.

Referenced by GenerateEvent(), and SetEventPrintFrequency().

FairEvtFilterParams PndFilteredPrimaryGenerator::fEvtFilterStat

protected

Contains the statistics of the event filtering process.

Definition at line 214 of file PndFilteredPrimaryGenerator.h.

Referenced by GenerateEvent(), GetNumberOfFilterFailedEvents(), GetNumberOfFilterMaxTries(), GetNumberOfGeneratedEvents(), SetFilterMaxTries(), and WriteEvtFilterStatsToRootFile().

std::vector<PndSmpFilterSet> PndFilteredPrimaryGenerator::fFilterSets

protected

Contains the filter-sets. Each filter set consist of one filter definition or a number of filters connected with logical AND (= a certain event type). All filter sets are combined with a logical OR (all event types to be accepted). Structure is like: (Filt_1.1 AND Filt_1.2) OR (Filt_2.1) OR (Filt_3.1 AND Filt_3.2 AND Filt_3.2)

Definition at line 211 of file PndFilteredPrimaryGenerator.h.

Referenced by AddFilter(), GenerateEvent(), and Init().