matching.hh

#ifndef Cone_H
#define Cone_H

#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <utility>



#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"

using namespace std;

namespace TemplateOverlap{

enum  {CONE, GAUSSIAN};
enum  {HADRONIC, SEMILEPTONIC};

class Settings {
private:
   double _subConeRadius;     // subCone radius
   double _coneRadius;    // characteristic jet radius
   double _sigma;  // Gaussian energy resolution relative to parton pT
   double _coneRadius2;  // subjet radius
   bool _variableCone; // Dynamically change cone radius based on parton  pT
   double _minPtParton; //For infrared safety, partons are not too soft

public:
   Settings(const double subConeRadiusIn=0.20, 
        const double sigmaIn=0.33, const double coneRadiusIn = 1.0,  const double coneRadius2In=0.40, const bool variableConeIn=false, 
         const double minPtPartonIn = 10.) :
   _subConeRadius(subConeRadiusIn), _coneRadius(coneRadiusIn), _sigma(sigmaIn),
   _coneRadius2(coneRadius2In),_variableCone(variableConeIn), _minPtParton(minPtPartonIn) {}
   double r() const {return _subConeRadius;}
   double R0() const {return _coneRadius;}
   double R1() const {return _coneRadius2;} 
   double sigma() const {return _sigma;}
   bool variableCone() const {return _variableCone;}
   double minPtParton() const {return _minPtParton;}
};



class SingleParticle {

public:

  // Constructors.
  SingleParticle(  double pTIn = 0., double yIn = 0., double phiIn = 0.)
  :  pT(pTIn), y(yIn), phi(phiIn), mult(1), isUsed(false), radius(0.), sigma(0.) { }
  SingleParticle(const SingleParticle& ssj) :  pT(ssj.pT),
    y(ssj.y), phi(ssj.phi), mult(ssj.mult), isUsed(ssj.isUsed), radius(ssj.radius), sigma(ssj.sigma) { }
  SingleParticle& operator=(const SingleParticle& ssj) { if (this != &ssj)
    {  pT = ssj.pT; y = ssj.y; phi = ssj.phi; 
    mult = ssj.mult; isUsed = ssj.isUsed; radius =ssj.radius; sigma=ssj.sigma;} return *this; }

  // Properties of particle.
  double pT, y, phi;
  int    mult; 
  bool   isUsed;
  double radius; //For templates
  double sigma;

  double deltaR2(const SingleParticle & other) const {
    double  dPhi = abs(phi - other.phi );
    if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
    double  dEta =  y - other.y;
    return (dEta * dEta + dPhi * dPhi );
  }
  
};

typedef std::vector<SingleParticle> jet_t;
typedef std::pair<double, std::vector<fastjet::PseudoJet> > templ_t;
bool checkTemplate (const jet_t & probe_template);
void matchTemplate( jet_t & jet, const  vector<jet_t> & templates, vector<double>  & result, int match_method=CONE, int normalization=HADRONIC);
void maximize(const vector<double>  & result, double & maxVal,  int & maxLoc );

std::vector<fastjet::PseudoJet> ConvertToPseudoJet(const jet_t& particles) ;
std::vector<fastjet::PseudoJet> ConvertToPseudoJet(const jet_t& particles, const fastjet::PseudoJet & tilt_axis) ;

jet_t ConvertToJet(const fastjet::PseudoJet & jet) ;
double DeltaPhi(double phi1, double phi2);



class MatchingMethod {
private:

   ifstream _mystream;
   vector<jet_t> _templates; 
   double _subConeRadius;
   double _sigma;
   double _minPtParton;
   double _maxVal;
   int _maxLoc; 
   
   int templateRead(ifstream& is);
   bool checkTemplate (const jet_t & probe_template);

public:
   MatchingMethod(const string & templateFile, Settings mySettings) : 
   _mystream(templateFile.c_str()), _subConeRadius(mySettings.r()), _sigma(mySettings.sigma()), _minPtParton(mySettings.minPtParton()) {templateRead(_mystream);}
  templ_t getOv (const fastjet::PseudoJet & jet);
};

templ_t  MatchingMethod::getOv(const fastjet::PseudoJet& jet)
{ 
  int match_method = CONE;
  int normalize = HADRONIC;
  
  jet_t jet_mat; 
  jet_mat = ConvertToJet(jet);
  
  vector<double>  result; 
  
  matchTemplate(jet_mat, _templates, result, match_method, normalize);
  
  maximize(result, _maxVal, _maxLoc );
  
  std::vector<fastjet::PseudoJet> peak_template = ConvertToPseudoJet(_templates[_maxLoc], jet);  
  return std::make_pair<double, std::vector<fastjet::PseudoJet> >(_maxVal, peak_template);
} 

int MatchingMethod::templateRead(ifstream& is) {

  string line, tag;
 
 //Event particlesSave;
 
jet_t particlesSave; 
_templates.clear();

  while (getline(is,line)){
  istringstream getpro(line);
  int nupSave;
  getpro >> nupSave;
  if (!getpro) return false;

  // Reset particlesSave vector.
  particlesSave.clear(); 

  // Read in particle info one by one, and store it.
  double pup1, pup2, pup3, pup4;
  
   if (!getline(is, line)) return false;
    istringstream getall(line);
    getall  >> pup1 >> pup2 >> pup3 >> pup4;
  
  for (int ip = 0; ip < nupSave; ++ip) { 
    if (!getline(is, line)) return false;
    istringstream getall(line);
    getall  >> pup1 >> pup2 >> pup3 >> pup4;
    if (!getall) return false;   
    
   fastjet::PseudoJet pTemp(pup1,pup2,pup3,pup4);
   SingleParticle pAux= SingleParticle( pTemp.perp(), pTemp.rapidity(), pTemp.phi());
   pAux.radius = _subConeRadius;
   pAux.sigma = _sigma;
   particlesSave.push_back(pAux);
  }
 
 if (checkTemplate (particlesSave) ) _templates.push_back(particlesSave);

        
 
  // Reading worked.
  
}
  return (_templates.size());

}

bool MatchingMethod::checkTemplate(const jet_t& probe_template)
{  
  bool isOK = true; 
  for (jet_t::const_iterator i = probe_template.begin() ; i != probe_template.end() && isOK ; ++i)  {  
    if ((*i).pT < _minPtParton) { isOK = false; break;}
    for (jet_t::const_iterator j = i + 1 ; j != probe_template.end(); ++j){ 
      if((*i).deltaR2((*j)) <   ((*i).radius+(*j).radius) *  ((*i).radius+(*j).radius))
    { isOK = false; break ;}
    }
  } 
  return (isOK);
}



void reset(jet_t & jet){
  for (jet_t::iterator j = jet.begin(); j != jet.end(); ++j) 
    if ((*j).isUsed) (*j).isUsed = false;   
        
}


double overlapDistance (jet_t & jet1, jet_t & jet2, double R0) 
{
  double sum = 0.;
    for(jet_t::const_iterator k = jet1.begin(); k != jet1.end(); ++k) {
            double eTjetNow = 0.;              
       //  Sum up unused particles within required distance of parton.
    for (jet_t::iterator j = jet2.begin(); j !=jet2.end(); ++j) {
      if ((*j).isUsed) continue;
    if ((*j).deltaR2(*k) > R0 * R0 ) continue;
         eTjetNow += (*j).pT;
        (*j).isUsed = true;
    }              
    double eTcenterNow = (*k).pT;
    //change sigma to E/3
    sum += 4.5 * (eTjetNow - eTcenterNow) * (eTjetNow - eTcenterNow)/(eTcenterNow * eTcenterNow);

}

reset(jet2);
return (exp(-sum));
    }


void matchTemplate(jet_t & jet, const  vector<jet_t> & templates, vector<double>  & result, int match_method, int normalization)
{ 
  for (int i = 0; i < int (templates.size()); ++i){
      double sum = 0.;

    for(jet_t::const_iterator t_particle = templates[i].begin(); t_particle != 
      templates[i].end(); ++t_particle) 
    {
      double eTjetNow = 0.;
      //  Sum up unused particles within required distance of parton.
      for ( jet_t::iterator particle = jet.begin(); particle != jet.end(); ++particle) 
      {
    if ((*particle).isUsed) continue;
    if ((*particle).deltaR2(*t_particle) > (*t_particle).radius * (*t_particle).radius ) continue;
        eTjetNow += (*particle).pT;
    (*particle).isUsed = true;
      }
      double eTcenterNow = (*t_particle).pT;

    // change sigma to E/3
      sum += 0.5*1/ ( ((*t_particle).sigma)*((*t_particle).sigma) ) * (eTjetNow - eTcenterNow) * (eTjetNow - eTcenterNow)/(eTcenterNow * eTcenterNow); 
    // if one overlap is too low, we skip the rest. 
      if (sum > 4.6 ) break;
      
    } 
// Clean flags. 
  reset(jet);
    
  result.push_back(exp(-sum));
  }
  

}

void maximize(const vector<double>  & result, double & maxVal, int & maxLoc )
{
   maxVal = 0;
   maxLoc = 0;

 for (int i = 0; i < int (result.size()); ++i){
                double thisVal = result[i];
        if (thisVal > maxVal){ maxLoc = i; maxVal = thisVal;}
 }
}


//--------------------------------------------------------------------------


std::vector<fastjet::PseudoJet> ConvertToPseudoJet(const jet_t& particles) 
{
   std::vector<fastjet::PseudoJet> Jets;
   for (jet_t::const_iterator part = particles.begin(); part != particles.end(); part++) {    
      fastjet::PseudoJet temp = fastjet::PtYPhiM((*part).pT, (*part).y, (*part).phi);
      Jets.push_back(temp);
   }
   return Jets;
}

std::vector<fastjet::PseudoJet> ConvertToPseudoJet(const jet_t& particles, const fastjet::PseudoJet & tilt_axis) 
{
   std::vector<fastjet::PseudoJet> Jets;
   double etaJet = tilt_axis.rapidity();
   double phiJet = tilt_axis.phi();
   for (jet_t::const_iterator part = particles.begin(); part != particles.end(); part++) {    
      fastjet::PseudoJet temp = fastjet::PtYPhiM((*part).pT, (*part).y + etaJet, (*part).phi + phiJet);
      Jets.push_back(temp);
   }
   return Jets;
}

jet_t ConvertToJet(const fastjet::PseudoJet& jet)
{
  std::vector<fastjet::PseudoJet> particles = jet.constituents();
  jet_t jet_mat; 
  
  double jetRapidity = jet.rapidity() ;   
  double jetPhi =      jet.phi()  ;

  for (unsigned int i = 0; i < particles.size(); ++i) { 
    double etaNow = particles[i].eta();
    double phiNow = particles[i].phi();
    double pTnow = particles[i].perp();
    
    jet_mat.push_back( SingleParticle( pTnow, etaNow - jetRapidity, phiNow - jetPhi) );
  }
  return jet_mat;
}

//--------------------------------------------------------------------------
double DeltaPhi(double phi1, double phi2)
{
  double dphi = phi1-phi2;
  
  if (dphi >  M_PI) dphi -= 2* M_PI;
  if (dphi < -M_PI) dphi += 2* M_PI;
 // if (abs (temp) > M_PI) temp = 2 * M_PI - abs(temp);
return dphi;  
}



} //end namespace


#endif // end Event_H