// ************************************************************************ // // PANDA Simple Particle Combiner Class // // Offers simple combinatorics // // ************************************************************************ // // Parameters: // - fAna : PndAnalysis instance // - decay : decay specification, e.g. "phi -> K+ K-; D_s+ -> phi pi+ cc" (cc indicates charged conjugate; particle used have to be defined beforehand) // - params : configuration parameters, e.g. "mwin=0.5:mwin(phi)=0.2:emin=0.1:pid=Loose:pide=Tight:algo=PidChargedProbability" // - mwin : mass window for all composites // - mwin(X) : mass window for composite 'X'; X (or cc of X) has to appear in decay // - emin : minimum energy threshold for neutrals // - pmin : minimum momentum threshold for charged // - pid : PID criterion ("All","VeryLoose","Loose","Tight","VeryTight","Best") for all species // - pide, pidmu, pidpi, pidk, pidp : PID criterion for specific species // - algo : PID algorithm (e.g. "PidAlgoEmcBayes;PidAlgoDrc" or "PidChargedProbability") for all species // - algoe, ... algop : PID algorithm for specific species // // // K.Goetzen 01/2015 // // ************************************************************************ // The header file #include "PndSimpleCombiner.h" // C++ headers #include #include // ROOT headers #include "TParticlePDG.h" #include "TObjArray.h" #include "TObjString.h" // RHO headers #include "RhoMassParticleSelector.h" #include "RhoEnergyParticleSelector.h" #include "RhoMomentumParticleSelector.h" #include "RhoTuple.h" // PANDA headers #include "PndAnalysis.h" using std::cout; using std::endl; // ------------------------------------------------------------------------- // constructor // ------------------------------------------------------------------------- PndSimpleCombiner::PndSimpleCombiner(PndAnalysis *fAna, TString decay, TString params, double Ecm) : fAnalysis(fAna), fDecay(decay), fGlobParams(params), fNLists(11), fVerbose(0), fEmin(0.), fPmin(0.), fEcm(Ecm), fESel(0), fPSel(0) { // initialize mapping pdg -> list index and list name int pdgcodes[] = { -11, 11, -13, 13, 211, -211, 321, -321, 2212, -2212, 22}; TString pdgnames[] = {"ElectronPlus", "ElectronMinus", "MuonPlus", "MuonMinus", "PionPlus", "PionMinus", "KaonPlus", "KaonMinus", "ProtonPlus", "ProtonMinus", "Neutral"}; fPdgIdxMap.clear(); fIdxPdgMap.clear(); fIdxListNameMap.clear(); for (int i=0;i<11;++i) { fPdgIdxMap[pdgcodes[i]] = i; // maps pdg code -> list index fIdxPdgMap[i] = pdgcodes[i]; // maps list index -> pdg code fIdxListNameMap[i] = pdgnames[i]; // maps list index -> generic list name (ElectronPlus, PionMinus, ...; see above) } // set initial pid configuration SetPid("All", "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts;PidAlgoSciT;PidAlgoRich;PidAlgoFtof"); if (fEcm==0) fEcm=-999; assert(ParseDecay(decay)); ParseParams(params); } // ------------------------------------------------------------------------- // destructor PndSimpleCombiner::~PndSimpleCombiner() { if (fESel) delete fESel; if (fPSel) delete fPSel; for (size_t i=0;iGetEntries(); for (int i=0;iAt(i))->String(); st.ReplaceAll("\t",""); st = st.Strip(TString::kBoth); if (st != "") toks.push_back(st); } return toks.size(); } // ------------------------------------------------------------------------- void PndSimpleCombiner::InitDecayInfo(SCDecayInfo &info, int pdg, int idx) { info.mpdg = pdg; // mother(composite) pdg code info.midx = idx; // add another list for the new composite info.daucc = false; // in case only the daughters have a cc, but not the mother, we have to care (e.g. etac -> Ks K+ pi-) info.mwin = 0; // default: no mass selection info.mwinlo = 0; // default: no mass selection info.mwinhi = 0; // default: no mass selection info.msel = 0; // default: no mass selection info.dpdg.clear(); // pdgs of daughters info.didx.clear(); // list index of daughters } // ------------------------------------------------------------------------- bool PndSimpleCombiner::ParseDecay(TString decay) { StringList subdec; int ndec = SplitString(decay, ";", subdec); // loop over subdecays for (int i=0;i' in the middle of the string? subdec[i].ReplaceAll("->",">"); SplitString(subdec[i],">",dectoks); if ( dectoks.size()!=2 || !TDatabasePDG::Instance()->GetParticle(dectoks[0]) ) {cout <<"[PndSimpleCombiner] **** ERROR : Invalid decay pattern '"<10) if (dectoks.size()>11) {cout <<"[PndSimpleCombiner] **** ERROR : Exceeding max number of daughters (10): '"<GetParticle(dectoks[0])->PdgCode(), fNLists++); fPdgIdxMap[info.mpdg] = info.midx; // add new list to the pdg <-> list index maps fIdxPdgMap[info.midx] = info.mpdg; bool cc=true; // loop over daughters for (size_t j=1;jGetParticle(dectoks[j])) {cout <<"[PndSimpleCombiner] **** ERROR : Unknown particle '"<GetParticle(dectoks[j])->PdgCode(); // add corresponding index if exists if ( fPdgIdxMap.find(dpdg) == fPdgIdxMap.end() ) {cout <<"[PndSimpleCombiner] **** ERROR : Undefined list '"< K- pi+ has cc -> D0bar -> K+ pi- (different FS and mother) // 2) D0 -> KS pi+ pi- has cc -> D0bar -> KS pi+ pi- (different mother) // 3) etac -> KS K- pi+ has cc -> etac -> KS K+ pi- (different FS) // 4) phi -> K+ K- has _no_ cc! (identical FS and mother) // does the mother have a cc? (case 1 + 2) // then we add a new list if (!CCInvariant(info.mpdg)) { SCDecayInfo ainfo; InitDecayInfo(ainfo, AntiPdg(info.mpdg), fNLists++); fPdgIdxMap[ainfo.mpdg] = ainfo.midx; fIdxPdgMap[ainfo.midx] = ainfo.mpdg; for (size_t j=0;j cut +- 0.25 around nominal mass of part // or : mwin(part)=0.25|0.75 -> cut 0.25 < mpart < 0.75 // **** if (pair[0]=="mwin") { double window = pair[1].Atof(); for (size_t j=0;jGetParticle(info.mpdg)->Mass(); if (info.mpdg/10 == 8888) { mean = fEcm;} info.mwin = window; info.mwinlo = mean - window/2; info.mwinhi = mean + window/2; // check if Ecm was not set and particle is pbarpSystem if (mean>-999) info.msel = new RhoMassParticleSelector("msel",mean,window); else { info.mwin = -1; // indicates that we skip the mass window cut for pbarpSystem if Ecm was not set } } } // **** // mass window for one composite // **** else if (pair[0].BeginsWith("mwin")) { // extract particle name from strings like 'mwin(D0)' pair[0] = pair[0](5,pair[0].Length()-6); if (!TDatabasePDG::Instance()->GetParticle(pair[0])) {cout <<"[PndSimpleCombiner] **** WARNING : Unknown particle type '"<GetParticle(pair[0])->PdgCode(); if (fPdgIdxMap.find(pdg) == fPdgIdxMap.end()) {cout <<"[PndSimpleCombiner] **** WARNING : Unknown particle list '"<GetParticle(pdg)->Mass(), window = 0; if (pdg/10 == 8888) { mean = fEcm; } // range setting like 1.2|1.6 ? if (pair[1].Contains("|")) { double lower = TString((pair[1])(0,pair[1].Index("|"))).Atof(); double upper = TString((pair[1])(pair[1].Index("|")+1,1000)).Atof(); mean = (upper+lower)/2.; window = upper-lower; if (mean<0 || window <=0) {cout <<"[PndSimpleCombiner] **** WARNING : Invalid mass window defintion '"<-999) info.msel = new RhoMassParticleSelector("msel",mean,window); else { info.mwin = -1; // indicates that we skip the mass window cut for pbarpSystem if Ecm was not set } } } } // **** // check for pid setting // **** if (pair[0] == "pid") SetPid(pair[1]); if (pair[0] == "pide") SetPidElectron(pair[1]); if (pair[0] == "pidmu") SetPidMuon(pair[1]); if (pair[0] == "pidpi") SetPidPion(pair[1]); if (pair[0] == "pidk") SetPidKaon(pair[1]); if (pair[0] == "pidp") SetPidProton(pair[1]); if (pair[0] == "algo") SetPid("",pair[1]); if (pair[0] == "algoe") SetPidElectron("",pair[1]); if (pair[0] == "algomu") SetPidMuon("",pair[1]); if (pair[0] == "algopi") SetPidPion("",pair[1]); if (pair[0] == "algok") SetPidKaon("",pair[1]); if (pair[0] == "algop") SetPidProton("",pair[1]); // **** // E_min or p_min // **** if (pair[0] == "emin") {fEmin = pair[1].Atof(); fESel = new RhoEnergyParticleSelector("eSel",50.+fEmin,100.);} if (pair[0] == "pmin") {fPmin = pair[1].Atof(); fPSel = new RhoMomentumParticleSelector("pSel",50.+fPmin,100.);} } return kTRUE; } // ------------------------------------------------------------------------- void PndSimpleCombiner::SetPid(TString crit, TString algo) { if (crit!="") fIdxPidCritMap.clear(); if (algo!="") fIdxPidAlgoMap.clear(); SetPidElectron(crit, algo); SetPidMuon(crit, algo); SetPidPion(crit, algo); SetPidKaon(crit, algo); SetPidProton(crit, algo); } // ------------------------------------------------------------------------- void PndSimpleCombiner::SetPidElectron(TString crit, TString algo) { if (crit!="") { fIdxPidCritMap[0] = crit; fIdxPidCritMap[1] = crit;} if (algo!="") { fIdxPidAlgoMap[0] = algo; fIdxPidAlgoMap[1] = algo;} } void PndSimpleCombiner::SetPidMuon(TString crit, TString algo) { if (crit!="") {fIdxPidCritMap[2] = crit; fIdxPidCritMap[3] = crit;} if (algo!="") {fIdxPidAlgoMap[2] = algo; fIdxPidAlgoMap[3] = algo;} } void PndSimpleCombiner::SetPidPion(TString crit, TString algo) { if (crit!="") {fIdxPidCritMap[5] = crit; fIdxPidCritMap[4] = crit;} if (algo!="") {fIdxPidAlgoMap[5] = algo; fIdxPidAlgoMap[4] = algo;} } void PndSimpleCombiner::SetPidKaon(TString crit, TString algo) { if (crit!="") {fIdxPidCritMap[7] = crit; fIdxPidCritMap[6] = crit;} if (algo!="") {fIdxPidAlgoMap[7] = algo; fIdxPidAlgoMap[6] = algo;} } void PndSimpleCombiner::SetPidProton(TString crit, TString algo) { if (crit!="") {fIdxPidCritMap[9] = crit; fIdxPidCritMap[8] = crit;} if (algo!="") {fIdxPidAlgoMap[9] = algo; fIdxPidAlgoMap[8] = algo;} } // ------------------------------------------------------------------------- // fill all generic lists, which are used by any of the composites void PndSimpleCombiner::FillGenericLists() { // fill all generic lists, which are used by any of the composites // loop through composites bool filled[11] = {0}; // marker, whether list is already filled int n = fDecayInfoArray.size(); for (int i=0;i10 || filled[lidx]) continue; // this list is not generic or already filled filled[lidx] = true; // mark list as filled int laidx = fPdgIdxMap[AntiPdg(info.dpdg[j])]; // index of antiparticle list (only relevant for charged) if (lidx<10) // charged lists { fAnalysis->FillList(fList[lidx],fIdxListNameMap[lidx]+fIdxPidCritMap[lidx],fIdxPidAlgoMap[lidx]); if (fPSel) fList[lidx].Select(fPSel); // also fill anti pdg generic list for charged generic ... fAnalysis->FillList(fList[laidx],fIdxListNameMap[laidx]+fIdxPidCritMap[laidx],fIdxPidAlgoMap[laidx]); if (fPSel) fList[laidx].Select(fPSel); filled[laidx] = true; // ... and mark list as filled } else if (lidx==10) //neutrals { fAnalysis->FillList(fList[lidx],fIdxListNameMap[lidx]); if (fESel) fList[lidx].Select(fESel); } if (lidx<11 && fVerbose) { cout <<"idx :"< &vpdg) { int sum=0, nminus=0, nplus=0; for (size_t i=0;i0) nplus++; else nminus++; } } return (sum==0 && nplus==nminus); } // ------------------------------------------------------------------------- // Returns pdg code of anti-particle int PndSimpleCombiner::AntiPdg(int pdg) { if (TDatabasePDG::Instance()->GetParticle(pdg)) { if (!TDatabasePDG::Instance()->GetParticle(pdg)->AntiParticle()) return pdg; else return TDatabasePDG::Instance()->GetParticle(pdg)->AntiParticle()->PdgCode(); } return -999999; } // ------------------------------------------------------------------------- // Prints the configuration void PndSimpleCombiner::Print() { cout <<"\n------------------------------------------"<GetParticle(info.mpdg)->GetName()<<"("< "; for (size_t j=0;jGetParticle(info.dpdg[j])->GetName()<<"("<0) cout <<" ["< aidx; for (size_t j=0;j &idx) { l.Cleanup(); int nd = idx.size(); switch (nd) { case 2: l.Combine(fList[idx[0]], fList[idx[1]], mpdg); break; case 3: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], mpdg); break; case 4: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], mpdg); break; case 5: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], mpdg); break; case 6: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], fList[idx[5]], mpdg); break; case 7: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], fList[idx[5]], fList[idx[6]], mpdg); break; case 8: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], fList[idx[5]], fList[idx[6]], fList[idx[7]], mpdg); break; case 9: l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], fList[idx[5]], fList[idx[6]], fList[idx[7]], fList[idx[8]], mpdg); break; case 10:l.Combine(fList[idx[0]], fList[idx[1]], fList[idx[2]], fList[idx[3]], fList[idx[4]], fList[idx[5]], fList[idx[6]], fList[idx[7]], fList[idx[8]], fList[idx[9]], mpdg); break; } return l.GetLength(); } // ------------------------------------------------------------------------- // access lists by name, pdg or index bool PndSimpleCombiner::GetList(RhoCandList &l, TString comp) { if (!TDatabasePDG::Instance()->GetParticle(comp)) return false; return GetList(l, TDatabasePDG::Instance()->GetParticle(comp)->PdgCode()); } // ------------------------------------------------------------------------- bool PndSimpleCombiner::GetList(RhoCandList &l, int pdg) { l.Cleanup(); if (fPdgIdxMap.find(pdg) == fPdgIdxMap.end()) return false; l = fList[fPdgIdxMap[pdg]]; return true; } // ------------------------------------------------------------------------- bool PndSimpleCombiner::GetListN(RhoCandList &l, int idx) { l.Cleanup(); if (idx>=0 && idx