sPhenixAj/src/JetAnalyzer.cxx

0001 /** @file JetAnalyzer.cxx
0002     @author Kauder:Kolja
0003     @version Revision 0.1
0004     @brief Light FastJet wrapper for Heavy Ion analysis
0005     program file and the functions in that file.
0006     @details
0007     @date Mar 04, 2015
0008 */
0009
0010 #include "JetAnalyzer.hh"
0011
0012 // ------------
0013 // Constructors
0014 // ------------
0015
0016 // ------------------------------------------------------------------------
0017 // Standard ctor
0018 JetAnalyzer::JetAnalyzer( std::vector<fastjet::PseudoJet>& InOrigParticles, fastjet::JetDefinition& JetDef,
0019               fastjet::AreaDefinition& AreaDef, fastjet::Selector selector_bkgd )
0020   : fastjet::ClusterSequenceArea ( InOrigParticles, JetDef, AreaDef ), OrigParticles ( InOrigParticles ),
0021     selector_bkgd (selector_bkgd)
0022 {
0023
0024   CanDoBackground=true;
0025
0026   // this should not be necessary :-/
0027   bkgd_subtractor=0;
0028   jet_def_bkgd = 0;
0029   bkgd_estimator = 0;
0030   bkgd_subtractor = 0;
0031   area_def_bkgd = new fastjet::AreaDefinition( AreaDef );
0032   // DEBUG
0033   // area_def_bkgd = new fastjet::AreaDefinition( fastjet::VoronoiAreaSpec( 1.0 ) );
0034
0035   // // std::cout << "Area Spec " << area_spec.description() << std::endl;
0036   // std::cout << "Area Def given: " << std::endl << AreaDef.description() << std::endl;
0037   // std::cout << "Area Def Bkgd:  " << std::endl << area_def_bkgd->description() << std::endl;
0038   // std::cout << "This:  " << std::endl << this->area_def().description() << std::endl;
0039
0040 }
0041 // ------------------------------------------------------------------------
0042
0043 // ------------------------------------------------------------------------
0044 // Standard ctor for use as ClusterSequence
0045 JetAnalyzer::JetAnalyzer( std::vector<fastjet::PseudoJet>& InOrigParticles, fastjet::JetDefinition& JetDef )
0046   : fastjet::ClusterSequenceArea ( InOrigParticles, JetDef, 0 ), OrigParticles ( InOrigParticles )
0047 {
0048   CanDoBackground=false;
0049
0050   // this should not be necessary :-/
0051   bkgd_subtractor=0;
0052   area_def_bkgd=0;
0053   jet_def_bkgd = 0;
0054   bkgd_estimator = 0;
0055   bkgd_subtractor = 0;
0056
0057 }
0058 // ------------------------------------------------------------------------
0059
0060
0061
0062 // ------------
0063 // Methods
0064 // ------------
0065
0066 // ------------------------------------------------------------------------
0067 // Background functionality
0068 fastjet::Subtractor* JetAnalyzer::GetBackgroundSubtractor(){
0069   if ( !CanDoBackground ) {
0070     // throw std::string("Should not be called unless we can actually do background subtraction.");
0071     return 0;
0072   }
0073
0074   // // Area
0075   // // ----
0076   // // Construct from jet definition and radius
0077   // double ymin, ymax;
0078   // selector_bkgd.get_rapidity_extent (ymin, ymax);
0079   // assert ( !std::isinf(ymin) && !std::isinf(ymax)   && "Can't handle unrestricted area." );
0080   // assert ( ymin == -ymax                            && "Can't construct ghosted area for asymmetric rapidity cuts.");
0081
0082   // double ghost_maxrap = ymax + 2.0*jet_def().R();
0083   // area_def_bkgd = new fastjet::AreaDefinition (fastjet::active_area_explicit_ghosts,
0084   //                           fastjet::GhostedAreaSpec( ghost_maxrap ));
0085
0086   // Background jet definition
0087   // -------------------------
0088   jet_def_bkgd = new fastjet::JetDefinition (fastjet::kt_algorithm, jet_def().R());
0089
0090   // Estimator and subtractor
0091   // ------------------------
0092   // std::cout << area_def_bkgd->description() << std::endl;
0093   // std::cout << selector_bkgd.description() << std::endl;
0094   // std::cout << "selector.has_finite_area() = "<< selector_bkgd.has_finite_area() << std::endl;
0095   bkgd_estimator  = new fastjet::JetMedianBackgroundEstimator(selector_bkgd, *jet_def_bkgd, *area_def_bkgd);
0096   bkgd_subtractor = new fastjet::Subtractor(bkgd_estimator);
0097
0098   // since FastJet 3.1.0, rho_m is supported natively in background
0099   // estimation (both JetMedianBackgroundEstimator and
0100   // GridMedianBackgroundEstimator).
0101   //
0102   // For backward-compatibility reasons its use is by default switched off
0103   // (as is the enforcement of m>0 for the subtracted jets). The
0104   // following 2 lines of code switch these on. They are strongly
0105   // recommended and should become the default in future versions of
0106   // FastJet.
0107 #if FASTJET_VERSION_NUMBER >= 30100
0108   bkgd_subtractor->set_use_rho_m(true);
0109   bkgd_subtractor->set_safe_mass(true);
0110 #endif
0111
0112   bkgd_estimator->set_particles( OrigParticles );
0113
0114   // std::cout << OrigParticles.size() <<  std::endl;
0115   // std::cout << "    rho     = " << bkgd_estimator->rho()   << std::endl;
0116   // std::cout << "    sigma   = " << bkgd_estimator->sigma() << std::endl;
0117
0118   return bkgd_subtractor;
0119 }
0120 // ------------------------------------------------------------------------
0121 // ----------------
0122 // Analysis methods
0123 // ----------------
0124 // Dijet finding as a Selector
0125 void SelectorDijetWorker::terminator(std::vector<const fastjet::PseudoJet *> & jets) const{
0126   // For each jet that does not pass the cuts, this routine sets the
0127   // pointer to 0.
0128
0129   // Save these ones from being zeroed
0130   int i0 = -1;
0131   int i1 = -1;
0132
0133   // For dPhi<0 or too few jets, we don't even need to try
0134   if ( dPhi>=0 && jets.size() >= 2 ) {
0135     // need to allow unsorted input
0136     // So find the two largest pT jets. Sigh.
0137
0138     // Pair negative pT^2 with the index
0139     // That takes care of pT sign and of stl::sorting order
0140     std::vector< std::pair<int,double> > IndexPt;
0141     for (unsigned int i=0; i<jets.size(); ++i ){
0142       // Note that jets.at(i) might be NULL
0143       IndexPt.push_back( std::pair<int,double>( i, jets.at(i) ? -jets.at(i)->perp2() : 0.0) );
0144     }
0145
0146     // and sort
0147     std::sort(IndexPt.begin(), IndexPt.end(), sort_IntDoubleByDouble());
0148
0149     // We only accept top two jets
0150     // Now Index 0 and 1 should point to the largest ones
0151     // Tentatively set them to be keepers
0152     i0 = IndexPt.at(0).first;
0153     i1 = IndexPt.at(1).first;
0154
0155     // clarified for readability
0156     double phi0 = jets.at(i0)->phi();
0157     double phi1 = jets.at(i1)->phi();
0158
0159     if ( ! (fabs ( JetAnalyzer::phimod2pi( phi0 - phi1 - JetAnalyzer::pi) ) < dPhi) ) {
0160       // Nope, don't save them after all
0161       i0=i1=-1;
0162     }
0163   } else {
0164     // std::cout << "  -------------------- dPhi<0 or Not enough jets" << std::endl;
0165   }
0166
0167   // Zero out everything else
0168   for ( int i=0; i< int( jets.size() ) ; ++i ){
0169     if ( i==i0 ) continue;
0170     if ( i==i1 ) continue;
0171     jets.at(i)= NULL;
0172   }
0173
0174   return;
0175 }
0176 fastjet::Selector SelectorDijets( const double dPhi ){
0177   return fastjet::Selector(new SelectorDijetWorker( dPhi ) );
0178 }
0179
0180 // ------------------------------------------------------------------------
0181
0182 // -------
0183 // Helpers
0184 // -------
0185 double JetAnalyzer::phimod2pi ( double phi ){
0186   while ( phi <-pi ) phi += 2.0*pi;
0187   while ( phi > pi ) phi -= 2.0*pi;
0188   return phi;
0189 }
0190
0191 // ------------------------------------------------------------------------
0192
0193 // -------------------
0194 // Getters and Setters
0195 // -------------------
0196
0197 // ------------------------------------------------------------------------
0198
0199
0200 // --------------
0201 // Other funtions
0202 // --------------
0203 // ------------------------------------------------------------------------
0204 bool IsMatched ( const std::vector<fastjet::PseudoJet>& jetset1, const std::vector<fastjet::PseudoJet>& jetset2, const double Rmax )
0205 {
0206   // could also return lists of matches or so
0207   // As it is, I'm using a primitive 1-to-1 matching
0208   // Enforcing 1 to 1 to avoid stupid pathologies
0209   // For dijet A_J, the following is already overkill, but I want to
0210   // lay the groundwork for more complex tasks
0211   if ( jetset1.size() != jetset2.size() ) return false;
0212
0213   std::vector<fastjet::PseudoJet> j2copy = jetset2; // cause we're going to destroy it
0214   for ( int j1 = 0; j1 < int( jetset1.size() ) ; ++j1 ){
0215     bool j1matched=false;
0216     for ( std::vector<fastjet::PseudoJet>::iterator pj2=j2copy.begin() ; pj2!=j2copy.end(); ++pj2 ){
0217       if ( jetset1.at(j1).delta_R( *pj2 ) < Rmax ){
0218     j1matched=true;
0219
0220     // std::cout << jetset1.at(j1).delta_R( *pj2 ) << std::endl; // DEBUG
0221     // could pop it into a result list, like so:
0222     // result.push_back (*pj2);
0223
0224     // remove from jetset2
0225     j2copy.erase( pj2 );
0226
0227     // done with this one
0228     break;
0229       } // if (R<max)
0230     } // for jet2
0231     // not matched? bad
0232     if ( !j1matched ) return false;
0233   } // for jet1
0234
0235   // Now we should have cleanly removed all jets from set 2;
0236   if ( j2copy.size() >0 ) return false;
0237
0238   return true;
0239 }
0240 // ------------------------------------------------------------------------
0241 bool IsMatched ( const std::vector<fastjet::PseudoJet>& jetset1, const fastjet::PseudoJet& reference, const double Rmax )
0242 {
0243   for ( int j1 = 0; j1 < int (jetset1.size()) ; ++ j1 ){
0244     if ( jetset1.at(j1).delta_R( reference ) < Rmax ){
0245       return true;
0246     } // if (R<max)
0247   } // for jet1
0248
0249   return false;
0250 }
0251
0252 // ------------------------------------------------------------------------
0253 bool IsMatched ( const fastjet::PseudoJet& jet1, const fastjet::PseudoJet& jet2, const double Rmax )
0254 {
0255   return jet1.delta_R( jet2 ) < Rmax;
0256 }
0257
0258 // ------------------------------------------------------------------------
0259 TLorentzVector MakeTLorentzVector ( const fastjet::PseudoJet& pj ){
0260   return TLorentzVector( pj.px(), pj.py(), pj.pz(), pj.E() );
0261   // Below could be slightly faster but circumvents encapsulation.
0262   // four_mom() is a valarray<double> which is NOT guaranteed to be contiguous
0263   // return TLorentzVector( (Double_t*) &pj.four_mom()[0] );
0264 }
0265 // ------------------------------------------------------------------------
0266
0267 // ------------------------------------------------------------------------
0268 fastjet::PseudoJet MakePseudoJet ( const TLorentzVector* const lv ){
0269   return fastjet::PseudoJet ( *lv );
0270 }
0271
0272
0273 // ------------------------------------------------------------------------
0274 // the function that allows to write simply
0275 // Selector sel = SelectorChargeRange( cmin, cmax );
0276 fastjet::Selector SelectorChargeRange( const int cmin, const int cmax){
0277   return fastjet::Selector(new SelectorChargeWorker(cmin, cmax));
0278 };
0279 // ------------------------------------------------------------------------
0280 // Helper to get an enum from a string
0281 fastjet::JetAlgorithm AlgoFromString( std::string s){
0282   std::transform(s.begin(), s.end(), s.begin(), ::tolower);
0283
0284   if ( s.substr(0,2) == "kt" )       return fastjet::kt_algorithm;
0285   if ( s.substr(0,6) == "antikt" )   return fastjet::antikt_algorithm;
0286   if ( s.substr(0,6) == "cambri" )   return fastjet::cambridge_algorithm;
0287
0288   return fastjet::undefined_jet_algorithm;
0289
0290 }
0291 // ------------------------------------------------------------------------
0292 const double JetAnalyzer::pi = 3.141592653589793238462643383279502884;