PhotosBranch.cxx

#include <vector>
#include <list>
#include "PH_HEPEVT_Interface.h"
#include "PhotosParticle.h"
#include "PhotosBranch.h"
#include "Photos.h"
#include "Log.h"
using std::vector;
using std::list;
using std::endl;

namespace Photospp
{

PhotosBranch::PhotosBranch(PhotosParticle* p)
{
        daughters = p->getDaughters();

        //Suppress if somehow got stable particle
        if(daughters.size()==0)
        {
                Log::Debug(1)<<"Stable particle."<<endl;
                suppression = 1;
                forcing     = 0;
                particle    = NULL;
                return;
        }
        else if(daughters.at(0)->getMothers().size()==1)
        {
                // Regular case - one mother
                Log::Debug(1)<<"Regular case."<<endl;
                particle  = p;
                mothers   = p->findProductionMothers();
        }
        else
        {
                // Advanced case - branch with multiple mothers - no mid-particle
                Log::Debug(1)<<"Advanced case."<<endl;
                particle  = NULL;
                mothers   = daughters.at(0)->getMothers();
        }
        forcing = checkForcingLevel();
        if(!forcing) suppression = checkSuppressionLevel();
        else         suppression = 0;

        // Even if forced or passed suppression check, we still have to check few things
        if(!suppression)
        {
                // Check momentum conservation
                suppression=!checkMomentumConservation();
                if(suppression) Log::Warning()<<"Branching ignored due to 4-momentum non conservation"<<endl;

                // Check if advanced case has only one daughter
                if(!particle && daughters.size()==1) suppression=-1;

                // If any of special cases is true, we're not forcing this branch
                if(suppression) forcing=0;
        }
}

void PhotosBranch::process()
{
        Log::Debug(703)<<"   Processing barcode: "<<( (particle) ? particle->getBarcode() : ( (mothers.size()) ? mothers.at(0)->getBarcode() : -1) )<<endl;
        /*
        cout<<"Particles send to photos (with barcodes in brackets):"<<endl;
        vector<PhotosParticle *> get = getParticles();
        for(int i=0;i<(int)get.size();i++) cout<<"ID: "<<get.at(i)->getPdgID()<<" ("<<get.at(i)->getBarcode()<<"), "; cout<<endl;
        */
        int index = PH_HEPEVT_Interface::set(this);
        PH_HEPEVT_Interface::prepare();
        photos_make_c_(&index);
        PH_HEPEVT_Interface::complete();
        PH_HEPEVT_Interface::get();
        checkMomentumConservation();
}

vector<PhotosParticle *> PhotosBranch::getParticles()
{
        vector<PhotosParticle *> ret = mothers;
        if(particle) ret.push_back(particle);
        ret.insert(ret.end(),daughters.begin(),daughters.end());
        return ret;
}

bool PhotosBranch::checkMomentumConservation()
{
        if(particle)           return particle->checkMomentumConservation();
        if(mothers.size()>0)   return mothers.at(0)->checkMomentumConservation();
        return true;
}

vector<PhotosBranch *> PhotosBranch::createBranches(vector<PhotosParticle *> particles)
{
        Log::Debug(700)<<"PhotosBranch::createBranches - filtering started"<<endl;
        list<PhotosParticle *> list(particles.begin(),particles.end());
        vector<PhotosBranch *> branches;

        // First - add all forced decays
        if(Photos::forceBremList)
        {
                std::list<PhotosParticle *>::iterator it;
                for(it=list.begin();it!=list.end();it++)
                {
                        PhotosBranch *branch = new PhotosBranch(*it);
                        int forcing = branch->getForcingStatus();
                        if(forcing)
                        {
                                Log::Debug(701)<<" Forced: "<<(*it)->getPdgID()<<" (barcode: "<<(*it)->getBarcode()<<") with forcing status= "<<forcing<<endl;
                                branches.push_back(branch);
                                it = list.erase(it);
                                --it;
                                // If forcing consecutive decays
                                if(forcing==2)
                                {
                                        PhotosParticle *p = branch->getDecayingParticle();
                                        if(!p)
                                                if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
                                                else continue;
                                        vector<PhotosParticle *> tree = p->getDecayTree();
                                        //Add branches for all particles from the list - max O(n*m)
                                        std::list<PhotosParticle *>::iterator it2;
                                        for(it2=list.begin();it2!=list.end();it2++)
                                        {
                                                for(int i=0;i<(int)tree.size();i++)
                                                {
                                                        if(tree.at(i)->getBarcode()==(*it2)->getBarcode())
                                                        {
                                                                PhotosBranch *b = new PhotosBranch(*it2);
                                                                branches.push_back(b);
                                                                // If we were to delete our next particle in line
                                                                if(it==it2) --it;
                                                                it2 = list.erase(it2);
                                                                --it2;
                                                                break;
                                                        }
                                                }
                                        }
                                }
                        }
                        else delete branch;
                }
        }
        // Quit if we're suppressing everything
        if(Photos::isSuppressed) return branches;
        // Now - check if remaining decays are suppressed
        while(!list.empty())
        {
                PhotosParticle *particle = list.front();
                list.pop_front();
                if(!particle) continue;

                PhotosBranch *branch = new PhotosBranch(particle);
                int suppression = branch->getSuppressionStatus();
                if(!suppression) branches.push_back(branch);
                else
                {
                        Log::Debug(702)<<"  Suppressed: "<<particle->getPdgID()<<" (barcode: "<<particle->getBarcode()<<") with suppression status= "<<suppression<<endl;
                        //If suppressing consecutive decays
                        if(suppression==2)
                        {
                                PhotosParticle *p = branch->getDecayingParticle();
                                if(!p)
                                        if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
                                        else continue;
                                vector<PhotosParticle *> tree = p->getDecayTree();
                                //Remove all particles from the list - max O(n*m)
                                std::list<PhotosParticle *>::iterator it;
                                for(it=list.begin();it!=list.end();it++)
                                {
                                        for(int i=0;i<(int)tree.size();i++)
                                        {
                                                if(tree.at(i)->getBarcode()==(*it)->getBarcode())
                                                {
                                                        it = list.erase(it);
                                                        --it;
                                                        break;
                                                }
                                        }
                                }
                        }
                        delete branch;
                        continue;
                }

                //In case we don't have mid-particle erase rest of the mothers from list
                if(!branch->getDecayingParticle())
                {
                        vector<PhotosParticle *> mothers = branch->getMothers();
                        for(int i=0;i<(int)mothers.size();i++)
                        {
                                PhotosParticle *m = mothers.at(i);
                                if(m->getBarcode()==particle->getBarcode()) continue;
                                std::list<PhotosParticle *>::iterator it;
                                for(it=list.begin();it!=list.end();it++)
                                        if(m->getBarcode()==(*it)->getBarcode())
                                        {
                                                it = list.erase(it);
                                                break;
                                        }
                        }
                }
        }
        return branches;
}

int PhotosBranch::checkList(bool forceOrSuppress)
{
        vector< vector<int>* > *list = (forceOrSuppress) ? Photos::forceBremList : Photos::supBremList;
        if(!list) return 0;

        // Can't check without pdgid
        int motherID;
        if(particle) motherID = particle->getPdgID();
        else
        {
                if(mothers.size()==0) return 0;
                motherID = mothers.at(0)->getPdgID();
        }

        // Create list of daughters
        vector<int> dID;
        for(int j=0;j<(int)daughters.size();j++) dID.push_back(daughters[j]->getPdgID());

        vector< vector<int> *> &patternList = *list;

        // Check if the mother and list of daughters matches any of the declared patterns
        for(int j=0; j<(int)patternList.size();j++)
        {
                // Skip patterns that don't have our mother
                if(motherID!=(*patternList[j])[0]) continue;

                // Compare decay daughters with pattern - max O(n*m)
                vector<int> &pattern = *patternList[j];
                bool fullMatch=true;
                for(int k = 1; k<(int)pattern.size()-1; k++)
                {
                        bool oneMatch=false;
                        for(int l=0;l<(int)dID.size(); l++)
                                if(pattern[k]==dID[l]) { oneMatch=true; break; }
                        if(!oneMatch) { fullMatch=false; break; }
                }
                // Check if the matching pattern is set for consecutive suppression
                /*
                  Currently minimal matching is used.
                  If e.g. 25 -> -15 is suppressed, then 25 -> 15,-15 and 25 -> 15,-15,22 etc. is suppressed too
                  For exact matching (then suppress 25 -> 15,-15 ; 25 -> 15,-15,22 etc. must be done separately) uncoment line ...:
                */
                // if(pattern.size()<=2 || (fullMatch && dID.size()==pattern.size()-2) )
                // ...and comment out line:
                if(pattern.size()<=2 || fullMatch)
                        return (pattern.back()==1) ? 2 : 1;
        }
        return 0;
}

} // namespace Photospp