/Users/jmonk/Physics/ForIA/src/AnalysisTools/SimpleTruthSelection.cxx

Go to the documentation of this file.

#include "ForIA/AnalysisTools/SimpleTruthSelection.hh"
#include "ForIA/Event.hh"
#include <cmath>

namespace ForIA{
 
  SimpleTruthSelection::SimpleTruthSelection(): FourVectorSelection(), m_charge(ANY), m_haveSetPT(false), 
                                                m_haveSetE(false),m_haveSetET(false),m_haveSetEta(false), 
                                                m_haveSetCharged(false),m_haveSetChargedET(false), m_haveSetNeutralET(false),
                                                m_haveSetChargedP(false), m_haveSetNeutralP(false),
                                                m_haveSetSpeciesDepSel(false), m_haveSetBarcode(false), 
                                                m_haveSetIsPrimary(false), m_haveSetID(false){}
    
  
  const TruthParticleVector &SimpleTruthSelection::particles()const{
    if(m_freshEvent) fillVectors();
    return m_passedParticles;
  }
  
  const TruthParticleVector &SimpleTruthSelection::rejectedParticles()const{
    if(m_freshEvent) fillVectors();
    return m_rejectedParticles;
  }
  
  bool SimpleTruthSelection::setEMin(double e){
    if(!m_canSetCuts) return false;
    m_haveSetE = true;
    m_eMin = e;
    return true;
  }

 bool SimpleTruthSelection::setETMin(double et){
    if(!m_canSetCuts) return false;
    m_haveSetET = true;
    m_etMin = et;
    return true; 
  }

  bool SimpleTruthSelection::setPTMin(double pt){
    if(!m_canSetCuts) return false;
    m_haveSetPT = true;
    m_ptMin = pt;
    return true; 
  }

  bool SimpleTruthSelection::setEtaMin(double eta){
    if(!m_canSetCuts) return false;
    m_haveSetEta = true;
    m_etaMin = eta;
    return true;
  }
  
  bool SimpleTruthSelection::setEtaMax(double eta){
    if(!m_canSetCuts) return false;
    m_haveSetEta = true;
    m_etaMax = eta;
    return true;    
  }
  
  bool SimpleTruthSelection::setCharged(charge chrg){
    if(!m_canSetCuts) return false;
    m_haveSetCharged = true;
    m_charge = chrg;
    return true;
  }

  bool SimpleTruthSelection::setChargedETMin(double et){
    if(!m_canSetCuts) return false;
    m_haveSetChargedET = true;
    m_chargedEtMin = et;
    return true;
  }
 
  bool SimpleTruthSelection::setNeutralETMin(double et){
    if(!m_canSetCuts) return false;
    m_haveSetNeutralET = true;
    m_neutralEtMin = et;
    return true;
  }

  bool SimpleTruthSelection::setChargedPMin(double p){
    if(!m_canSetCuts) return false;
    m_haveSetChargedP = true;
    m_chargedPMin = p;
    return true;
  }

  bool SimpleTruthSelection::setNeutralPMin(double p){
    if(!m_canSetCuts) return false;
    m_haveSetNeutralP = true;
    m_neutralPMin = p;
    return true;
  }

  bool SimpleTruthSelection::setSpeciesDepETMin(double etaMin, double etaMax, double etCharged, double etNeutral){
    if(!m_canSetCuts) return false;
    if(etaMax != 0.0F) m_haveSetSpeciesDepSel = true;
    m_speciesDepEtaMin = etaMin;
    m_speciesDepEtaMax = etaMax;
    m_chargedEtMin = etCharged;
    m_neutralEtMin = etNeutral;
    return true;
  }

  bool SimpleTruthSelection::setSpeciesDepPMin(double etaMin, double etaMax, double pCharged, double pNeutral){
    if(!m_canSetCuts) return false;
    m_haveSetSpeciesDepSel = true;
    m_haveSetChargedP = true;
    m_haveSetNeutralP = true;
    m_speciesDepEtaMin = etaMin;
    m_speciesDepEtaMax = etaMax;
    m_chargedPMin = pCharged;
    m_neutralPMin = pNeutral;
    return true;
  }
  
  bool SimpleTruthSelection::setBarcodeMax(int bcode){
    if(!m_canSetCuts) return false;
    m_haveSetBarcode = true;
    m_barcodeMax = bcode;
    return true;
  }
  
  bool SimpleTruthSelection::setIsPrimary(bool isPrimary){
    if(!m_canSetCuts) return false;
    m_haveSetIsPrimary = true;
    m_isPrimary = isPrimary;
    return true;
  }

  bool SimpleTruthSelection::setPdgId(int pdgID){
    if(!m_canSetCuts) return false;
    m_haveSetID = true;
    m_pdgID = pdgID;
    return true;
  }
  
  void SimpleTruthSelection::fillVectors()const{

    m_freshEvent = false;
    
    m_passedParticles.clear();
    m_passedFourVectors.clear();
    m_rejectedParticles.clear();
    m_rejectedFourVectors.clear();

    if(m_event == 0) return;

    TruthParticleVector input = m_event->truthParticles();
    //    sort( input.begin(), input.end(), SortEt );
    //    std::cout<<"NEw ECENT " << std::endl;
    for(TruthParticleVector::const_iterator particle = input.begin();
        particle != input.end(); ++particle){

      bool isGenStable = true;
      if(m_haveSetIsPrimary ){

        isGenStable = (  ((*particle)->status()%1000 == 1) ||
                        // need vertex bar code ((*particle)->status()==2 && vertex_barcode <-200000) ||
                        ((*particle)->status()%1000 == 2 &&  (*particle)->status() > 1000)
                         );

        isGenStable = isGenStable  &&  (  (*particle)->barcode()<200000 );
        isGenStable = isGenStable && ! (abs( (*particle)->pdgId() ) == 21 && (*particle)->E()==0);
      }

    

      bool passedCuts = isGenStable;
      //first check if Eta/Species dependent cuts are requested
      bool isTagged = false;
      if(m_haveSetSpeciesDepSel && passedCuts){
        if ( fabs((*particle)->eta()) > m_speciesDepEtaMin && fabs((*particle)->eta()) < m_speciesDepEtaMax ){
          double p = 0;
          if( sin((*particle)->theta()) != 0 || fabs(sin((*particle)->theta())) > 0.000001 ) 
            p = (*particle)->PT() / sin((*particle)->theta());
          //check if particles is charged
          //if ( (*particle)->charge() != 0 && (*particle)->ET() < m_chargedEtMin ){
          if( (*particle)->charge() != 0 ){
            if( m_haveSetChargedET && (*particle)->ET() < m_chargedEtMin ){
              passedCuts = false; 
              isTagged = true;
            }
            if( m_haveSetChargedP && fabs(p) < m_chargedPMin ){
              passedCuts = false;
              isTagged = true;
            }
          }
          //check if particle is neutral
          //if ( (*particle)->charge() == 0 && (*particle)->ET() < m_neutralEtMin ){
          if ( (*particle)->charge() == 0 ){ 
            if( m_haveSetNeutralET && (*particle)->ET() < m_neutralEtMin ){
              passedCuts = false;
              isTagged = true;
            }
            if( m_haveSetNeutralP && fabs(p) < m_neutralPMin ){
              passedCuts = false;
              isTagged = true;
            }
          }
        }//end species dependent cuts
        //else? what about forward particles
        /*
        else if( (*particle)->ET() < m_neutralEtMin ){
          passedCuts = false;
          isTagged = true;
        }
        */
      }
      //apply generic cuts otherwise
      if(!isTagged){
        if(m_haveSetPT  && passedCuts) passedCuts = ((*particle)->PT() > m_ptMin); //std::cout<<"here 6a " << passedCuts <<std::endl; 
        if(m_haveSetET  && passedCuts) passedCuts = ((*particle)->ET() > m_etMin);  //std::cout<<"here 6b "<< m_etMin<<" " << passedCuts <<std::endl;
        if(m_haveSetE   && passedCuts) passedCuts = ((*particle)->E() > m_eMin);  //std::cout<<"here 6c " << passedCuts <<std::endl;
        if(m_haveSetEta && passedCuts) passedCuts = (fabs((*particle)->eta()) < m_etaMax);  //std::cout<<"here 6d " << passedCuts <<std::endl;
      }
      if(m_haveSetBarcode && passedCuts){ passedCuts = ((*particle)->barcode() < m_barcodeMax);
        std::cout << "ZXCV barcodeMax = " << m_barcodeMax << std::endl;}
      // if(m_haveSetIsPrimary && passedCuts) passedCuts = (m_isPrimary == ((*particle)->barcode() < 200000 && (*particle)->barcode() != 0));
      
      
      if(m_haveSetCharged && passedCuts){
        
        switch(m_charge){
        case NEUTRAL:
          if((*particle)->charge() != 0) passedCuts = false;
          break;
        case CHARGED:
          if((*particle)->charge() == 0) passedCuts = false;
          break;
        case POSITIVE:
          if((*particle)->charge() <= 0) passedCuts = false;
          break;
        case NEGATIVE:
          if((*particle)->charge() >= 0) passedCuts = false;
          break;
          
        case ANY:
          
        default:
          break;
        }
      }
      if(passedCuts){
        
        m_passedParticles.push_back(*particle);
        m_passedFourVectors.push_back(*particle);
      }else if(isGenStable){
        m_rejectedParticles.push_back(*particle);
        m_rejectedFourVectors.push_back(*particle);
      }
    }
    
    return;
  }
  
}