NtupleReader/src/MyTruthParticleManager.cxx

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// 07.12.2006, AUTHOR: MANFRED GROH
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS MyTruthParticleManager ::
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//*****************************************************************************
#include "NtupleReader/MyTruthParticleManager.h"

//:::::::::::::::::::::::::
//:: NAME SPACE SETTINGS ::
//:::::::::::::::::::::::::

using namespace std;


bool MyTruthParticleManager::init(vector<MyParticle*> in_truth_particles){
    // temporary variable
//    unsigned int tmp_barcode(0);
    // reset data
    n_of_decay_chains = 0;
    primary_intpart_0 = NULL;
    v_intpart_0_children_gen1.clear();
    vv_intpart_0_children_gen2.clear();
    primary_intpart_1 = NULL;
    v_intpart_1_children_gen1.clear();
    vv_intpart_1_children_gen2.clear();
        genstat_e=9999;
        genstat_mu=9999;
        genstat_tau=9999;
     
    // loop over given vector
    for (unsigned int i=0; i<in_truth_particles.size(); i++){
         // check if a barcode is stored twice
/*
         tmp_barcode=(unsigned int)( (MyTruthParticle*)(in_truth_particles[i]))->barcode();
         if (getParticle_with_barcode(tmp_barcode)!=NULL){
             cerr<<"Error in MyTruthParticleManager::init"<<endl
                 <<"vector of TruthParticles contains barcode "
                 <<tmp_barcode<<" twice"<<endl;
             return false;
         }
//*/
         // store the particle
         v_truth_particles.push_back( (MyTruthParticle*)(in_truth_particles[i]));
    }
    // return true if all particles were stored
    if (in_truth_particles.size()==v_truth_particles.size()){
        return true;
    }
    // return false otherwise
    else{
        cerr<<"Error in MyTruthParticleManager::init"<<endl;
        return false;
    }
}
//*****************************************************************************

//:::::::::::::::::::::::::::::::::
//:: METHOD set_generator ::
//:::::::::::::::::::::::::::::::::

void MyTruthParticleManager::set_generator(string gen){
        generator = gen;
    if (generator=="Alpgen"){
                genstat_e = 1;
                genstat_mu = 1;
                genstat_tau = 195;
        }
        else if (generator=="Herwig"){
                genstat_e = 1;
                genstat_mu = 1;
                genstat_tau = 195;
        }
        else if (generator=="Pythia"){
                genstat_e = 1;
                genstat_mu = 1;
                genstat_tau = 2;
        }
        else if (generator=="McAtNlo"){
                genstat_e = 1;
                genstat_mu = 1;
                genstat_tau = 195;
        }
        else if (generator=="Sherpa"){
                genstat_e = 1;
                genstat_mu = 1;
                genstat_tau = 2; //has to be checked !!!
        }
        else {
                cerr<<"ERROR: MyTruthParticleManager::set_generator:"<<endl
                <<"Unknown Generator: "<<generator<<endl;
                exit(1);
        }
//* debug info */       cout<<"Set Genstats to:"<<endl
//* debug info */           <<"e:   "<<genstat_e<<endl
//* debug info */           <<"mu:  "<<genstat_mu<<endl
//* debug info */           <<"tau: "<<genstat_tau<<endl;
                
        
}

//*****************************************************************************

//:::::::::::::::::::::::::::::::::
//:: METHOD getParticles_with_ID ::
//:::::::::::::::::::::::::::::::::

vector<MyParticle*> MyTruthParticleManager::getParticles_with_ID(unsigned int pdgId) const {
        vector<MyParticle*> tmp_v;
        tmp_v.clear();
        
        // if given ID==0, return all particles
        if (pdgId==0){
                for (unsigned int i=0; i<v_truth_particles.size(); i++){
                        tmp_v.push_back(v_truth_particles[i]);
                }
                return tmp_v;
        }
        
        // scan full truth info for the given ID
        for (unsigned int i=0; i<v_truth_particles.size(); i++){
                if ( (unsigned int)TMath::Abs((v_truth_particles[i])->pdgId())==pdgId){
                        tmp_v.push_back(v_truth_particles[i]);
                }
        }
        return tmp_v;
}

//*****************************************************************************

//:::::::::::::::::::::::::::::::::::::
//:: METHOD getParticle_with_barcode ::
//:::::::::::::::::::::::::::::::::::::

MyTruthParticle* MyTruthParticleManager::getParticle_with_barcode(unsigned int tmp_barcode) const {
    for (unsigned int i=0; i<v_truth_particles.size(); i++){
                if ( (unsigned int)((v_truth_particles[i])->barcode())==tmp_barcode){
            return (v_truth_particles[i]);
        }
    }
    // if no particle with the given barcode found, return NULL
    return NULL;
}

//*****************************************************************************

//:::::::::::::::::::::::
//:: METHOD get_mother ::
//:::::::::::::::::::::::

MyTruthParticle* MyTruthParticleManager::get_mother(MyTruthParticle* truth_particle) const {
    int tmp_nParents(0);
    int tmp_barcode(0);
    MyTruthParticle* tmp_truth_part;

    // get number of parents
    tmp_nParents=truth_particle->nParents();

    // if number == 0, or >1 return NULL
    if (tmp_nParents==0 || tmp_nParents>1){
        return NULL;
    }

    // if number == 1, search for the barcode of the mother
    tmp_barcode=truth_particle->mother0_barcode();
    tmp_truth_part = getParticle_with_barcode(tmp_barcode);
    if (tmp_truth_part==NULL){
        //cerr<<"ERROR: Mother not found!"<<endl;
    }
    return tmp_truth_part;
}

//*****************************************************************************

//:::::::::::::::::::::::::::::::
//:: METHOD get_primary_mother ::
//:::::::::::::::::::::::::::::::

MyTruthParticle* MyTruthParticleManager::get_primary_mother(MyTruthParticle* truth_particle) const {
    // check number of parents
    // return particle itself if no mother, or more than 1 mother
    if (truth_particle->nParents()==0 || truth_particle->nParents()>1){
        return truth_particle;
    }
    // stop at a Higgs
    if (truth_particle->pdgId()==25){
        return truth_particle;
    }

//* debug info */ cerr<<"truth_particle->pdgId()"<<truth_particle->pdgId()<<endl;
//* debug info */ cerr<<""<<endl;
//* debug info */ cerr<<""<<endl;

    // if number of parents equals 1, get mother particle
    MyTruthParticle* truth_particle2;
//* debug info */ cerr<<"in MyTruthParticleManager::get_primary_mother: before get_mother"<<endl;
    truth_particle2 = get_mother(truth_particle);
    if (truth_particle2==NULL){
        //cerr<<"ERROR: Mother not found!"<<endl;
        return NULL;
    }
//     // If mother is a quark: 
//     if (TMath::Abs(truth_particle2->pdgId())<9){
//         return truth_particle;
//     }



    // call this routine recoursively with this particle
    // return output
    return get_primary_mother(truth_particle2);

}

//*****************************************************************************

//:::::::::::::::::::::::::
//:: METHOD  get_children::
//:::::::::::::::::::::::::

vector<MyTruthParticle*>
    MyTruthParticleManager::get_children(MyTruthParticle* truth_particle)const{
    unsigned int tmp_barcode(0);
    vector<MyTruthParticle*> v_tmp_truth_particles;
    v_tmp_truth_particles.clear();
    tmp_barcode = truth_particle->barcode();
    int mother_pdgId(0);
    mother_pdgId = truth_particle->pdgId();
    // loop over all stored truth particles
    // if mother0_barcode of a particle==barcode
    // of the given particle, store the particle 
    // in a vector
    for (unsigned int i=0; i<v_truth_particles.size(); i++){
        if ((unsigned int)((v_truth_particles[i])->mother0_barcode())
            ==tmp_barcode){
            v_tmp_truth_particles.push_back(v_truth_particles[i]);
        }
    }
/*      
    for (unsigned int i=0; i<v_tmp_truth_particles.size(); i++){
        cout<<"child "<<i<<" ID: "
            <<(v_tmp_truth_particles[i])->pdgId()<<endl;
    }
    //*/
    // return the vector
    return v_tmp_truth_particles;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::
//:: METHOD get_last_instance ::
//::::::::::::::::::::::::::::::

MyTruthParticle* MyTruthParticleManager::get_last_instance(MyTruthParticle* part_in){
    vector<MyTruthParticle*> v_tmp_children;
    int nb_children_with_same_id;
     nb_children_with_same_id = 0;
   
    // get the direct children of the given particle
    v_tmp_children = get_children(part_in);

    // if the given particle has no children
    // return given particle
    if (v_tmp_children.size() ==0 ){
        //cerr<<"No more children"<<endl;
        return part_in;
    }

    // loop over children
    for (unsigned int i=0; i<v_tmp_children.size(); i++){
        // if a child has the same ID like the given particle,
        // call this method recursive
        if ( ((v_tmp_children[i])->pdgId())==part_in->pdgId()){
            nb_children_with_same_id++;
            return get_last_instance(v_tmp_children[i]);
        }
    }
    // if there was no child with the same ID like the given particle
    // return the given particle
    if (nb_children_with_same_id==0){
        //cerr<<"No more particles of same type within children"<<endl;
        return part_in;
    }
    // one should not arrive here...
    cerr<<"shouldn't arrive here: MyTruthParticleManager::get_last_instance"<<endl;
    return NULL;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::::::::::::::::
//:: METHOD remove_multiples_wo_using_genstat ::
//::::::::::::::::::::::::::::::::::::::::::::::

vector<MyParticle*> MyTruthParticleManager::remove_multiples_wo_using_genstat(const vector<MyParticle*> v_input){
    // a temporary particle
    MyTruthParticle* tmp_part1;
    tmp_part1=NULL;

    // a new vector that contains the cleaned up sample of particles
    vector<MyParticle*> v_out;
    v_out.clear();
    
    // tmeporary variables
    int child_bc, tmp_part1_bc;
    child_bc = 0;
    tmp_part1_bc = 0;
    bool already_stored;
    already_stored = false;

    // iterate over vector from the beginning
    for (int j=0; j<((int)v_input.size()); j++){
        // get the last instance of this particle
        tmp_part1 = get_last_instance((MyTruthParticle*)v_input[j]);
        if (tmp_part1==NULL){
            // error
            cerr<<"ERROR in MyTruthParticleManager::remove_multiples"<<endl;
        }
        else {
            // check if particle was already stored in the output vector
            already_stored = false;
            for (unsigned int k=0; k<v_out.size(); k++){
                child_bc = ((MyTruthParticle*)v_out[k])->barcode();
                tmp_part1_bc = tmp_part1->barcode();
                if( child_bc==tmp_part1_bc ){
                    already_stored = true;
                }
            }
            
            // if not, store it
            if (!already_stored){
                v_out.push_back(tmp_part1);
            }
            // if yes, go on to the next particle
        }

    }

    // remove taus with no daughters:
    for (int k=v_out.size()-1; k>-1; k--){
//        if (TMath::Abs(((MyTruthParticle*)v_out[k])->pdgId())==15){
//            cerr<<"Id: "<<<<"nDaughters: "
//                <<((MyTruthParticle*)v_out[k])->nDaughters()<<endl;
//        }
        // check if it is a tau and has no daughters
        if (TMath::Abs(((MyTruthParticle*)v_out[k])->pdgId())==15 &&
                ((MyTruthParticle*)v_out[k])->nDaughters()==0){
            // remove this tau
            v_out.erase(v_out.begin()+k);
            }
        }
//        return v_out;

// check isolation of particles:
// inform user about content of vector

        
        // iterate over vector from the beginning
        for (int j=0; j<((int)v_out.size()); ){
        // j is not increased automaitcally, only if particle is not removed
        // if the particle is removed, the non increased index
        // automaically points to the next particle
            bool to_be_removed;
            to_be_removed=false;
/*
            cerr<<"part1-Id: "<<((MyTruthParticle*)(v_out[j]))->pdgId()
                <<", BC: "<<((MyTruthParticle*)(v_out[j]))->barcode()
                <<", Status: "<<((MyTruthParticle*)(v_out[j]))->status()
                <<", nP: "<<((MyTruthParticle*)(v_out[j]))->nParents()
                <<", PId: "<<((MyTruthParticle*)(v_out[j]))->mother0_pdgId()
                <<", PBC: "<<((MyTruthParticle*)(v_out[j]))->mother0_barcode()
                <<", nD: "<<((MyTruthParticle*)(v_out[j]))->nDaughters()<<endl;
//*/    
            
            for (int k=1; k<((int)v_out.size()); k++){
                double dr(100000000);
                double dE(100000000);
                double dpdgId(1);
                double dstatus(0);
                dr = (v_out[j])->tlv().DeltaR((v_out[k])->tlv());
                dE = TMath::Abs((v_out[j])->tlv().E()-(v_out[k])->tlv().E());
                dpdgId = ((MyTruthParticle*)(v_out[j]))->pdgId()-((MyTruthParticle*)(v_out[k]))->pdgId();
                dstatus = ((MyTruthParticle*)(v_out[j]))->status()-((MyTruthParticle*)(v_out[k]))->status();
/*              
                cerr<<"part2-Id: "<<((MyTruthParticle*)(v_out[k]))->pdgId()
                    <<", BC: "<<((MyTruthParticle*)(v_out[k]))->barcode()
                    <<", Status: "<<((MyTruthParticle*)(v_out[k]))->status()
                    <<", nP: "<<((MyTruthParticle*)(v_out[k]))->nParents()
                    <<", PId: "<<((MyTruthParticle*)(v_out[k]))->mother0_pdgId()
                    <<", PBC: "<<((MyTruthParticle*)(v_out[k]))->mother0_barcode()
                    <<", nD: "<<((MyTruthParticle*)(v_out[k]))->nDaughters()<<endl;
                    cerr<<"dR: "<<dr<<", dE: "<<dE<<endl;
//*/
                if (dr<2 && dE<500000 && dpdgId==0 && dstatus>0){
                    //*
//                  cerr<<"found double, pos 1,2: "<<j<<", "<<k<<endl;
                if (((MyTruthParticle*)(v_out[j]))->status()==1 || ((MyTruthParticle*)(v_out[j]))->status()==195){
                    cerr<<"Removing a truth particle with status 1 or 195"<<endl;
                }
/*
                    cerr<<"dR: "<<dr<<", dE: "<<dE<<", dstatus: "<<dstatus<<endl;
                    cerr<<"part1-Id: "<<((MyTruthParticle*)(v_out[j]))->pdgId()
                        <<", BC: "<<((MyTruthParticle*)(v_out[j]))->barcode()
                        <<", Status: "<<((MyTruthParticle*)(v_out[j]))->status()
                        <<", nP: "<<((MyTruthParticle*)(v_out[j]))->nParents()
                        <<", PId: "<<((MyTruthParticle*)(v_out[j]))->mother0_pdgId()
                        <<", PBC: "<<((MyTruthParticle*)(v_out[j]))->mother0_barcode()
                        <<", nD: "<<((MyTruthParticle*)(v_out[j]))->nDaughters()<<endl;
                    cerr<<"part2-Id: "<<((MyTruthParticle*)(v_out[k]))->pdgId()
                        <<", BC: "<<((MyTruthParticle*)(v_out[k]))->barcode()
                        <<", Status: "<<((MyTruthParticle*)(v_out[k]))->status()
                        <<", nP: "<<((MyTruthParticle*)(v_out[k]))->nParents()
                        <<", PId: "<<((MyTruthParticle*)(v_out[k]))->mother0_pdgId()
                        <<", PBC: "<<((MyTruthParticle*)(v_out[k]))->mother0_barcode()
                        <<", nD: "<<((MyTruthParticle*)(v_out[k]))->nDaughters()<<endl;
                    //*/
                    to_be_removed=true;
                }
            }
            if(to_be_removed){
                // remove this pair
                v_out.erase(v_out.begin()+j);
                }
            else{
                // check the next particle
                j++;
            }
        }
        return v_out;
}

//*****************************************************************************

//:::::::::::::::::::::::::::::
//:: METHOD remove_multiples ::
//:::::::::::::::::::::::::::::

vector<MyParticle*> MyTruthParticleManager::remove_multiples(const vector<MyParticle*> v_input){
    
    // check if generator was set
    if(genstat_e==9999 || genstat_mu==9999 || genstat_tau==9999){
        cerr<<"Generator not set!"<<endl
            <<"Set it before using remove_multiples!"<<endl;
        exit(1);
    }
    
    // a new vector that contains the cleaned up sample of particles
    vector<MyParticle*> v_out;
    v_out.clear();
    
    // tmeporary variables
    int tmp_abs_pdgId;
    tmp_abs_pdgId = 0;
    unsigned int tmp_genStat;
    tmp_genStat = 0;

    // iterate over vector from the beginning
    for (int j=0; j<((int)v_input.size()); j++){
        // get the abs(pdgId) and status of the particle
        tmp_abs_pdgId = TMath::Abs(((MyTruthParticle*)v_input[j])->pdgId());
        tmp_genStat = (unsigned int)TMath::Abs(((MyTruthParticle*)v_input[j])->status());
        
//      cout<<"genstat_e required: "<<genstat_e<<endl;
//      cout<<"genstat_mu required: "<<genstat_mu<<endl;
//      cout<<"genstat_tau required: "<<genstat_tau<<endl;
//      cout<<"genstat: "<<tmp_genStat<<endl;
//      cout<<"pdgId: "<<tmp_abs_pdgId<<endl;
        
        if((tmp_abs_pdgId == 11) && (tmp_genStat == genstat_e)){
            v_out.push_back(v_input[j]);
        }
        if((tmp_abs_pdgId == 13) && (tmp_genStat == genstat_mu)){
            v_out.push_back(v_input[j]);
        }
        if((tmp_abs_pdgId == 15) && (tmp_genStat == genstat_tau)){
            v_out.push_back(v_input[j]);
        }
    }
    return v_out;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::::
//:: METHOD get_n_of_decay_chains ::
//::::::::::::::::::::::::::::::::::

// get number of interesting decay chains
unsigned int MyTruthParticleManager::get_n_of_decay_chains(void) const {
    return n_of_decay_chains;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::
//:: METHOD get_primary_intpart ::
//::::::::::::::::::::::::::::::::

// get the primary interesting particle of decay chain n_decay_chain
MyTruthParticle* 
    MyTruthParticleManager::get_primary_intpart(int n_decay_chain) const {
    if (n_decay_chain>n_of_decay_chains){
        cerr<<"decay chain "<<n_decay_chain<<" does not exist"<<endl;
        return primary_intpart_0;
    }
    if (n_decay_chain==0){
        return primary_intpart_0;
    }
    if (n_decay_chain==1){
        return primary_intpart_1;
    }
    return primary_intpart_0;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::::::::
//:: METHOD get_gen1_children_intpart ::
//::::::::::::::::::::::::::::::::::::::

// get the 1st generation children of the
// primary interesting particle of decay chain n_decay_chain
vector<MyTruthParticle*> 
    MyTruthParticleManager::get_gen1_children_intpart(int n_decay_chain) const {
    if (n_decay_chain>n_of_decay_chains){
        cerr<<"decay chain "<<n_decay_chain<<" does not exist"<<endl;
        return v_intpart_0_children_gen1;
    }
    if (n_decay_chain==0){
        return v_intpart_0_children_gen1;
    }
    if (n_decay_chain==1){
        return v_intpart_1_children_gen1;
    }
    return v_intpart_0_children_gen1;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::::::::
//:: METHOD get_gen2_children_intpart ::
//::::::::::::::::::::::::::::::::::::::

// get the 2nd generation children of the
// primary interesting particle of decay chain n_decay_chain
vector< vector<MyTruthParticle*> > 
    MyTruthParticleManager::get_gen2_children_intpart(int n_decay_chain) const {
    if (n_decay_chain>n_of_decay_chains){
        cerr<<"decay chain "<<n_decay_chain<<" does not exist"<<endl;
        return vv_intpart_0_children_gen2;
    }
    if (n_decay_chain==0){
        return vv_intpart_0_children_gen2;
    }
    if (n_decay_chain==1){
        return vv_intpart_1_children_gen2;
    }
    return vv_intpart_0_children_gen2;
}

//*****************************************************************************

//::::::::::::::::::::::::::::::::
//:: METHOD find_decay_chain_of ::
//::::::::::::::::::::::::::::::::

int MyTruthParticleManager::find_decay_chain_of(vector<unsigned int> interesting_pdgIds){
    // reset the number of interesting decay chains
    n_of_decay_chains = 0;

    // temporary variables
    vector<MyTruthParticle*> v_int_part;
    MyTruthParticle* tmp_truth_part;
    int primary_mother_barcode(0);
    int int_part_0_barcode(0);  // barcode of the primary mother 
                                // of first decay chain
    int int_part_1_barcode(0);  // and of second chain


    // scan through all the truth particles for interesting PDG-IDs
    // store the pointers to all the particles with an interesting ID
    // in a temporary vector
    for (unsigned int i=0; i<interesting_pdgIds.size(); i++){
        for (unsigned int j=0; j<v_truth_particles.size(); j++){
            if ( (unsigned int)TMath::Abs((v_truth_particles[j])->pdgId())==interesting_pdgIds[i]
                && (v_truth_particles[j])->barcode()<50){
                v_int_part.push_back(v_truth_particles[j]);
            }
        }
    }

    // check wether there were interesting particles
    // by checking the size of the temporary vector
    // that was just built
    if (v_int_part.size()==0){
        //cout<<"No interesting particle found"<<endl;
        n_of_decay_chains = 0;
        return 0;
    }

//* debug info */       for (unsigned int i=0; i<v_int_part.size(); i++){
//* debug info */               cout<<"found PDG-ID: "<<(v_int_part[i])->pdgId()<<endl;
//* debug info */       }

    // loop over the temporary vector, assign the particles to decay chains
    // store the primary particle of each found deacy chain
    for(unsigned int i=0; i<v_int_part.size(); i++){
//* debug info */ cerr<<"in MyTruthParticleManager::find_decay_chain_of: tempvector loop"<<endl;
        
        // get the primary mother of the current particle
        tmp_truth_part = get_primary_mother((v_int_part[i]));
//* debug info */ cerr<<"in MyTruthParticleManager::find_decay_chain_of: after get_primary_mother"<<endl;

        tmp_truth_part = get_last_instance(tmp_truth_part);
        
        // throw a warning if no mother was found
        if (tmp_truth_part==NULL){
           //cerr<<"ERROR: Mother not found!"<<endl;
           n_of_decay_chains= -1;
           return n_of_decay_chains;
        }
    
        
        // check if generator was set
        if(genstat_e==9999 || genstat_mu==9999 || genstat_tau==9999){
            cerr<<"Generator not set!"<<endl
                <<"Set it before using find_decay_chain_of!"<<endl;
                exit(1);
        }


        // Fix to get only one Z in Alpgen samples:
        if(generator=="Alpgen"){
            if(tmp_truth_part->pdgId()==23 && tmp_truth_part->status()!=155){
//                cout<<"Generator is Alpgen and found Z with status!=155: skipping..."<<endl;
                continue;
            }
        }
        
        // get the barcode of the primary mother
        primary_mother_barcode = tmp_truth_part->barcode();
//* debug info */ cout<<"primary_mother_barcode: "<<primary_mother_barcode<<endl;
    
        // if there is no primary particle stored yet,
        // store this particle as first one
        if (int_part_0_barcode==0){
            int_part_0_barcode=primary_mother_barcode;
            primary_intpart_0=tmp_truth_part;
            continue;
        }
    
        // continue if this primary particle was already stored
        // as first interesting particle
        if (int_part_0_barcode==primary_mother_barcode){
            continue;
        }
    
        // from here on the current particle does NOT belong
        // to the first decay chain
    
        // if there is no second interesting primary particle
        // stored yet, store this particle as second one
        if (int_part_1_barcode==0){
            int_part_1_barcode=primary_mother_barcode;
            primary_intpart_1=tmp_truth_part;
            continue;
        }
    
        // continue if this primary particle was already stored
        // as second interesting particle
        if (int_part_1_barcode==primary_mother_barcode){
            continue;
        }
    
        // since from here on an interesting particle was found
        // that does not belong to a decay chain of the two found
        // interesting decay chains, throw a warning that more than
        // 2 interesting decay chains found
        n_of_decay_chains = 999;
        cout<<"More than 2 interesting decay chains found!!!"<<endl;
        // end method here
        return n_of_decay_chains;
    
    } // end of loop over all interesting truth particles

    // evaluate the number of decay chains
    if (int_part_0_barcode==0 && int_part_1_barcode==0){
        n_of_decay_chains=0;
    }       
    if (int_part_0_barcode>0 && int_part_1_barcode==0){
        n_of_decay_chains=1;
    }       
    if (int_part_0_barcode>0 && int_part_1_barcode>0){
        n_of_decay_chains=2;
    }
    return n_of_decay_chains;
}

//*****************************************************************************

//:::::::::::::::::::::::::::::::::
//:: METHOD process_decay_chains ::
//:::::::::::::::::::::::::::::::::

bool MyTruthParticleManager::process_decay_chains(void){
    // fills the first and second generation children
    // for the primary interesting particles found

    // first chain
    if (n_of_decay_chains>0){
    /*
        cout<<"primary particle 0 ID: "<<primary_intpart_0->pdgId()<<endl;
        //*/
        v_intpart_0_children_gen1.clear();
        MyTruthParticle* primary_intpart_0_li;
        primary_intpart_0_li = get_last_instance(primary_intpart_0);
        v_intpart_0_children_gen1=get_children(primary_intpart_0_li);
        // 2nd generation
        for (unsigned int i=0; i<v_intpart_0_children_gen1.size(); i++){
            vv_intpart_0_children_gen2.push_back(
            get_children(get_last_instance(v_intpart_0_children_gen1[i])) );
        }
    }

    // second chain
    if (n_of_decay_chains>1){
    /*
        cout<<"primary particle 1 ID: "<<primary_intpart_1->pdgId()<<endl;
        //*/
        v_intpart_1_children_gen1.clear();
        MyTruthParticle* primary_intpart_1_li;
        primary_intpart_1_li = get_last_instance(primary_intpart_1);
        v_intpart_1_children_gen1=get_children(primary_intpart_1_li);
        // 2nd generation
        for (unsigned int i=0; i<v_intpart_1_children_gen1.size(); i++){
            vv_intpart_1_children_gen2.push_back( get_children(get_last_instance(v_intpart_1_children_gen1[i])) );
        }
    }
    return true;
}

---