ATLAS_2012_I1094061.cc

Go to the documentation of this file.

// -*- C++ -*-

#include "Rivet/Analysis.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Tools/Logging.hh"

#include "Rivet/Projections/ChargedFinalState.hh"

#include <boost/assign/std/vector.hpp>

namespace Rivet{
    
  class ATLAS_2012_I1094061;
  
  AnalysisBuilder<ATLAS_2012_I1094061> plugin_ATLAS_2012_I1094061;
  
  class ATLAS_2012_I1094061: public Analysis{
    
  private: 
    
    //////////////////////////////////////////////////////////////////////////// 
    /**
     * Little container to hold a pair of foreground and background histos and then 
     * divide them at the end of the analysis
     */
    class HistoPair{
      
    public:
      
      enum HistoType{FOREGROUND, BACKGROUND};
      
      HistoPair(): _analysis(0),
      _h_foreground(0), _h_background(0), _d_final(0){}
      
      void init(int ds, int xaxis, int yaxis, ATLAS_2012_I1094061 *analysis){
        
        _ds = ds;
        _xaxis = xaxis;
        _yaxis = yaxis;
        
        _analysis = analysis;
        
        ++HistoPair::_s_counter;
        
        const BinEdges &edges = _analysis->binEdges(_ds, _xaxis, _yaxis);
        
        string sCount = boost::lexical_cast<string>(HistoPair::_s_counter);
        
        _h_foreground = analysis->bookHistogram1D("tmpForeground" + sCount, edges);
        _h_background = analysis->bookHistogram1D("tmpBackground" + sCount, edges);
      }
      
      void fillForeground(double value, double weight){
        _h_foreground->fill(value, weight);
        _h_foreground->fill(-value, weight);        
        return;
      }
      
      void fillBackground(double value, double weight){
        _h_background->fill(value, weight);
        _h_background->fill(-value, weight);
        return;
      }
      
      void fill(double value, double weight, HistoType type){
        
        switch(type){
          case FOREGROUND:
            fillForeground(value, weight);
            break;
          case BACKGROUND:
            fillBackground(value, weight);
            break;
        }
        
        return;
      }
      
      void finalize(double wgtSum, double bgWeight, double avNTracks){
        
        _h_foreground->scale(1. / wgtSum);
        _h_background->scale(1. / bgWeight);
        
        string histoPath = _analysis->histoPath(_ds, _xaxis, _yaxis);
        
        AIDA::IDataPointSet *final = _analysis->histogramFactory().divide(histoPath, *_h_foreground, *_h_background);
        
        for(int ii=0; ii!= final->size(); ++ii){
          AIDA::IDataPoint *pt = final->point(ii);
          double y=pt->coordinate(1)->value();
          pt->coordinate(1)->setValue(y-(avNTracks - 1));
        }
        
        _analysis->histogramFactory().destroy(_h_foreground);
        _analysis->histogramFactory().destroy(_h_background);
        
        return;
      }
      
    private:
      
      int _ds, _xaxis, _yaxis;
      
      ATLAS_2012_I1094061 *_analysis;
      
      AIDA::IHistogram1D* _h_foreground;
      AIDA::IHistogram1D* _h_background;
      AIDA::IDataPointSet* _d_final;
      
      static short _s_counter;
      
    };
    
    ////////////////////////////////////////////////////////////////////////////
    
  public:
    
    ATLAS_2012_I1094061(): Analysis("ATLAS_2012_I1094061"), 
    _minpT(100.*MeV), _etaMax(2.5), _nVersions(5), _version(0),
    _etaCut(2.), _phiCut(0.5*M_PI),
    _historyInclusive(_nVersions, ParticleVector()), _historyN20(_nVersions, ParticleVector()),
    _historyInclusiveWgts(_nVersions, 0.), _historyN20Wgts(_nVersions, 0.),
    _particleCountInclusive(0.), _particleCountN20(0.),
    _weightInclusive(0.), _weightN20(0.),
    _bgWeightInclusive(0.), _bgWeightN20(0.){
      
    }
    
    //////////////////////////////////////////////////////////////////////////// 
    void init(){
      
      const ChargedFinalState cfs(-2.5, 2.5, _minpT);
      addProjection(cfs, "ChargedParticles");
      
      // Only do the multiplicity > 20 plots for 7 TeV collisions
      _doN20 = (fabs(sqrtS() - 7000.*GeV) < 0.1*GeV);
      
      int yaxis = (_doN20)? 2: 1;

      _hp_DEta_0_pi.init(1, 1, yaxis, this);
      _hp_DEta_0_pi2.init(2, 1, yaxis, this);
      _hp_DEta_pi2_pi.init(3, 1, yaxis, this);
      
      _hp_DPhi_0_2.init(4, 1, yaxis, this);
      _hp_DPhi_2_5.init(5, 1, yaxis, this);
      
      if(_doN20){
        
        yaxis = 3;
        
        _hp_N20_DEta_0_pi.init(1, 1, yaxis, this);
        _hp_N20_DEta_0_pi2.init(2, 1, yaxis, this);
        _hp_N20_DEta_pi2_pi.init(3, 1, yaxis, this);
        
        _hp_N20_DPhi_0_2.init(4, 1, yaxis, this);
        _hp_N20_DPhi_2_5.init(5, 1, yaxis, this);
              
      }
      return;
    }
    
    //////////////////////////////////////////////////////////////////////////// 
    void analyze(const Event &evt){
           
      const ChargedFinalState &cfsProj = applyProjection<ChargedFinalState>(evt, "ChargedParticles");
      
      ParticleVector chargedParticles = cfsProj.particles();
      
      if(chargedParticles.size() < 2) vetoEvent;
      
      bool hasN20 = (_doN20 && chargedParticles.size() >= 20);
      
      double dMultiplicity = (double)chargedParticles.size();
      
      double multiplicityWeightIncr = dMultiplicity * evt.weight();
      
      _weightInclusive += evt.weight();
      _particleCountInclusive += multiplicityWeightIncr;
      
      if(hasN20){
        _weightN20 += evt.weight();
        _particleCountN20 += multiplicityWeightIncr;
      }
      
      double fgWeight = 2.*evt.weight() / dMultiplicity;
            
      for(ParticleVector::const_iterator p1 = chargedParticles.begin();
          p1 != chargedParticles.end(); ++p1){

        ParticleVector::const_iterator p2 = p1;
        ++p2;
        
        // fill the foreground distributions
        while(p2 != chargedParticles.end()){
          fillHistosInclusive(*p1, *p2, fgWeight, HistoPair::FOREGROUND);
          if(hasN20) fillHistosN20(*p1, *p2, fgWeight, HistoPair::FOREGROUND);
          ++p2;
        }// end filling the foreground distributions
        
        // loop over the history of particles from previous events and fill the background 
        // by correlating those particles with the current event
                
        for(size_t version = 0; version != _nVersions; ++version){
          
          const ParticleVector &bgParticles = _historyInclusive[version];
          double bgWeight = evt.weight() * _historyInclusiveWgts[version];
          
          for(ParticleVector::const_iterator p2 = bgParticles.begin();
              p2 != bgParticles.end(); ++p2){
            fillHistosInclusive(*p1, *p2, bgWeight, HistoPair::BACKGROUND);
            _bgWeightInclusive += bgWeight;
          }
          
          if(!hasN20) continue;
          
          const ParticleVector &bgParticlesN20 = _historyN20[version];
          bgWeight = evt.weight() * _historyN20Wgts[version];
          
          for(ParticleVector::const_iterator p2 = bgParticlesN20.begin();
              p2 != bgParticlesN20.end(); ++p2){
            fillHistosN20(*p1, *p2, bgWeight, HistoPair::BACKGROUND);
            _bgWeightN20 += bgWeight;
          }
          
        }//end loop over particle history for background fill
        
      }// end particle loop
        
      // Overwrite the history for the version count number
      _historyInclusive[_version] = chargedParticles;
      _historyInclusiveWgts[_version] = evt.weight();
      
      if(hasN20){
        _historyN20[_version] = chargedParticles;
        _historyN20Wgts[_version] = evt.weight();
      }
      
      ++_version;
      if(_version == _nVersions) _version = 0;
      
      return;
    }
    
    //////////////////////////////////////////////////////////////////////////// 
    void finalize(){
      
      double avMultiplicity = _particleCountInclusive / _weightInclusive;
      
      _hp_DEta_0_pi.finalize(_weightInclusive,  _bgWeightInclusive, avMultiplicity);
      _hp_DEta_0_pi2.finalize(_weightInclusive, _bgWeightInclusive, avMultiplicity);
      _hp_DEta_pi2_pi.finalize(_weightInclusive,_bgWeightInclusive, avMultiplicity);
      
      _hp_DPhi_0_2.finalize(_weightInclusive, _bgWeightInclusive, avMultiplicity);
      _hp_DPhi_2_5.finalize(_weightInclusive, _bgWeightInclusive, avMultiplicity);
      
      if(_doN20){
        avMultiplicity = _particleCountN20 / _weightN20;
        _hp_N20_DEta_0_pi.finalize(_weightN20,   _bgWeightN20, avMultiplicity);
        _hp_N20_DEta_0_pi2.finalize(_weightN20,  _bgWeightN20, avMultiplicity);
        _hp_N20_DEta_pi2_pi.finalize(_weightN20, _bgWeightN20, avMultiplicity);
        
        _hp_N20_DPhi_0_2.finalize(_weightN20, _bgWeightN20, avMultiplicity);
        _hp_N20_DPhi_2_5.finalize(_weightN20, _bgWeightN20, avMultiplicity);
      }
      
      return;
    }
  
    //////////////////////////////////////////////////////////////////////////// 
    
    void fillHistos(const Particle &p1, const Particle &p2, double weight, 
                    HistoPair::HistoType type, bool inclusive){
            
      double dEta = fabs(p1.momentum().eta() - p2.momentum().eta());
      double dPhi = mapAngle0ToPi(p1.momentum().phi() - p2.momentum().phi());
      double dPhiShift = TWOPI - dPhi;
            
      HistoPair &dEta_0_pi   = (inclusive)? _hp_DEta_0_pi   :_hp_N20_DEta_0_pi;
      HistoPair &dPhi_0_2    = (inclusive)? _hp_DPhi_0_2    :_hp_N20_DPhi_0_2;
      HistoPair &dPhi_2_5    = (inclusive)? _hp_DPhi_2_5    :_hp_N20_DPhi_2_5;
      HistoPair &dEta_0_pi2  = (inclusive)? _hp_DEta_0_pi2  :_hp_N20_DEta_0_pi2;
      HistoPair &dEta_pi2_pi = (inclusive)? _hp_DEta_pi2_pi :_hp_N20_DEta_pi2_pi;
      
      dEta_0_pi.fill(dEta, weight, type);
      
      if(dEta < _etaCut){
        dPhi_0_2.fill(dPhi, weight, type);
        dPhi_0_2.fill(dPhiShift, weight, type);        
      }else{
        dPhi_2_5.fill(dPhi, weight, type);
        dPhi_2_5.fill(dPhiShift, weight, type);
      }
      
      if(dPhi < _phiCut){
        dEta_0_pi2.fill(dEta, weight, type);
      }else{
        dEta_pi2_pi.fill(dEta, weight, type);        
      }

      return;
      
    }
    //////////////////////////////////////////////////////////////////////////// 

    void fillHistosInclusive(const Particle &p1, const Particle &p2, double weight, 
                             HistoPair::HistoType type){
      
      fillHistos(p1, p2, weight, type, true);
      return;
    }
    
    void fillHistosN20(const Particle &p1, const Particle &p2, double weight, 
                       HistoPair::HistoType type){
      
      fillHistos(p1, p2, weight, type, false);
      return;
    }
    //////////////////////////////////////////////////////////////////////////// 
    
    
    double _minpT;
    double _etaMax;
    
    size_t _nVersions;
    size_t _version;
    
    double _etaCut;
    double _phiCut;
    
    /// The "history" vectors contain the history of particles from _nVersions previous events
    /// These are used to construct the background correlation.
    vector<ParticleVector> _historyInclusive;
    vector<ParticleVector> _historyN20;
    
    vector<double>         _historyInclusiveWgts;
    vector<double>         _historyN20Wgts;
    
    double _particleCountInclusive;
    double _particleCountN20;
    double _weightInclusive;
    double _weightN20;
    
    double _bgWeightInclusive;
    double _bgWeightN20;
    
    bool _doN20;
    
    HistoPair _hp_DEta_0_pi;
    HistoPair _hp_DEta_0_pi2;
    HistoPair _hp_DEta_pi2_pi;
    
    HistoPair _hp_DPhi_0_2;
    HistoPair _hp_DPhi_2_5;
    
    HistoPair _hp_N20_DEta_0_pi;
    HistoPair _hp_N20_DEta_0_pi2;
    HistoPair _hp_N20_DEta_pi2_pi;
    
    HistoPair _hp_N20_DPhi_0_2;
    HistoPair _hp_N20_DPhi_2_5;
    
  };
 
  short ATLAS_2012_I1094061::HistoPair::_s_counter = 0;
  
}