JetAnalyzer.hh

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetAnalyzer.hh

/*
  Author self-assessment statement: It gets the job done for A_J.
*/

#ifndef JETANALYZER_H
#define  JETANALYZER_H

// Includes and namespaces
//#include "FastJet3.h"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/ClusterSequencePassiveArea.hh"
#include "fastjet/ClusterSequenceActiveArea.hh"
#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
#include "fastjet/Selector.hh"

#include <TLorentzVector.h>
#include <TClonesArray.h>

#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
#include "fastjet/tools/Subtractor.hh"
#include "fastjet/tools/Filter.hh"
#include "fastjet/FunctionOfPseudoJet.hh"

#include <iostream>
#include <string>
#include <sstream>

// For pairs and sorting in the dijet finding
#include <utility>
#include <algorithm>


class JetAnalyzer : public fastjet::ClusterSequenceArea {

private :  
  std::vector<fastjet::PseudoJet>& OrigParticles;

  bool CanDoBackground;
  
  fastjet::JetMedianBackgroundEstimator* bkgd_estimator;
  fastjet::Subtractor* bkgd_subtractor;

  fastjet::Selector selector_bkgd;
  fastjet::JetDefinition* jet_def_bkgd;  
  fastjet::AreaDefinition* area_def_bkgd;

public :
  // ------------
  // Constructors
  // ------------
  JetAnalyzer ( std::vector<fastjet::PseudoJet>& InOrigParticles, fastjet::JetDefinition& JetDef, fastjet::AreaDefinition& AreaDef,
                fastjet::Selector selector_bkgd=fastjet::SelectorAbsRapMax( 0.6 ) * (!fastjet::SelectorNHardest(2)) );

  JetAnalyzer ( std::vector<fastjet::PseudoJet>& InOrigParticles, fastjet::JetDefinition& JetDef );
  
  ~JetAnalyzer (){
    if ( area_def_bkgd )  {
      delete area_def_bkgd;
      area_def_bkgd = 0 ;
    }
    if ( jet_def_bkgd )  {
      delete jet_def_bkgd;
      jet_def_bkgd = 0 ;
    }
    if ( bkgd_estimator ) {
      delete bkgd_estimator;
      bkgd_estimator = 0 ;
    }
    if ( bkgd_subtractor ) {
      delete bkgd_subtractor;
      bkgd_subtractor = 0 ;
    }
    
  };
  

  // ----------------
  // Analysis methods
  // ----------------

  // -------------
  // Other Methods
  // -------------
  fastjet::Subtractor* GetBackgroundSubtractor();
  fastjet::JetMedianBackgroundEstimator* GetBackgroundEstimator() { return bkgd_estimator; };
  // ---------
  // Operators 
  // ---------
  // May eventually be added.
  // bool operator==(JetAnalyzer& rhs);

  // -------
  // Helpers
  // -------
  static const double pi;

  static double phimod2pi ( double phi );  
  
};

// NOT part of the class!
// =============================================================================
class SelectorDijetWorker : public fastjet::SelectorWorker{
public:
  SelectorDijetWorker( double dPhi )
    : dPhi(dPhi){};
  
  std::string description() const{
    std::ostringstream oss;
    oss << "Searches for and returns dijet pairs within |phi1 - phi2 - pi| < " << dPhi;    
    return oss.str() ; 
  };

  bool applies_jet_by_jet() const { return false;};

  bool pass(const fastjet::PseudoJet& pj) const{
    if (!applies_jet_by_jet())
      throw ( std::string("Cannot apply this selector worker to an individual jet") );
    return false;
  };

  void terminator(std::vector<const fastjet::PseudoJet *> & jets) const;
 
private:
  const double dPhi;                   
};

struct sort_IntDoubleByDouble {
  bool operator()(const std::pair<int,double> &left, const std::pair<int,double> &right) {
    return left.second < right.second;
  }
};
fastjet::Selector SelectorDijets( const double dPhi=0.4);



// =============================================================================
bool IsMatched ( const std::vector<fastjet::PseudoJet>& jetset1, const std::vector<fastjet::PseudoJet>& jetset2, const double Rmax );

bool IsMatched ( const std::vector<fastjet::PseudoJet>& jetset1, const fastjet::PseudoJet& reference, const double Rmax );

bool IsMatched ( const fastjet::PseudoJet& jet1, const fastjet::PseudoJet& jet2, const double Rmax );



// =============================================================================
// ------------------------------------------------------------------------
TLorentzVector MakeTLorentzVector ( const fastjet::PseudoJet& pj );
// ------------------------------------------------------------------------

// // =============================================================================
// Experimental, not working.
// Should work if this class is _instantiated_ and then the instance is used.
// That's not really the usage scenario.
// class MakeTLorentzVector : public fastjet::FunctionOfPseudoJet<TLorentzVector>{
//   TLorentzVector result(const fastjet::PseudoJet &pj) const{
//     return TLorentzVector( pj.px(), pj.py(), pj.pz(), pj.E() );
//   };
//   std::string description() const {return "Translates Pseudojet into ROOT TLorentzVector";}
//   TLorentzVector operator()(const fastjet::PseudoJet & pj) const {return result(pj);};
// };


// =============================================================================
fastjet::PseudoJet MakePseudoJet ( const TLorentzVector* const lv );
// =============================================================================
// =============================================================================

class JetAnalysisUserInfo: public fastjet::PseudoJet::UserInfoBase {
public:
  JetAnalysisUserInfo(int quarkcharge=-999) :  quarkcharge( quarkcharge )  {};
  
  int GetQuarkCharge() const { return quarkcharge; };

private:
  const int quarkcharge;   
};

// =============================================================================
class SelectorChargeWorker : public fastjet::SelectorWorker{
public:
  SelectorChargeWorker( const int cmin, const int cmax)  : cmin(cmin), cmax(cmax) {};
  
  std::string description() const{
    std::ostringstream oss;
    oss.str("");
    oss << cmin << " <= quark charge <= " << cmax;    
    return oss.str() ; 
  };
  
  bool pass(const fastjet::PseudoJet &p) const{
    const int & quarkcharge = p.user_info<JetAnalysisUserInfo>().GetQuarkCharge();
    return (quarkcharge >= cmin) &&  (quarkcharge<= cmax);
  };

private: 
  const int cmin; 
  const int cmax; 
  
};
// =============================================================================
fastjet::Selector SelectorChargeRange( const int cmin=-999, const int cmax=999);

// =============================================================================
// =============================================================================
// =============================================================================
fastjet::JetAlgorithm AlgoFromString( std::string s);


#endif // JETANALYZER_H