1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7#include "Rivet/Projections/Sphericity.hh"
  8#include "Rivet/Projections/Hemispheres.hh"
  9#include "Rivet/Projections/ParisiTensor.hh"
 10#include "Rivet/Tools/BinnedHistogram.hh"
 11#include "fastjet/EECambridgePlugin.hh"
 12
 13namespace Rivet {
 14
 15
 16  /// @brief event shapes vs thrust direction
 17  class DELPHI_2000_I522656 : public Analysis {
 18  public:
 19
 20    /// Constructor
 21    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2000_I522656);
 22
 23
 24    /// @name Analysis methods
 25    ///@{
 26
 27    /// Book histograms and initialise projections before the run
 28    void init() {
 29      // Initialise and register projections.
 30      declare(Beam(), "Beams");
 31      const FinalState fs;
 32      declare(fs, "FS");
 33      const Thrust thrust(fs);
 34      declare(thrust, "Thrust");
 35      declare(Sphericity(fs), "Sphericity");
 36      declare(ParisiTensor(fs), "Parisi");
 37      declare(Hemispheres(thrust), "Hemispheres");
 38      declare(FastJets(fs, FastJets::DURHAM, 0.7), "DurhamJets");
 39      declare(FastJets(fs, FastJets::JADE  , 0.7), "JadeJets"  );
 40
 41      // book histograms
 42      vector<double> bins={0.00,0.12,0.24,0.36,0.48,0.60,0.72,0.84,0.96};
 43      // thrust angle binned
 44      unsigned int iy=1,ioff=0;
 45      for(unsigned int ix=0;ix<8;++ix) {
 46	{Histo1DPtr tmp; _h_EEC.add(bins[ix],bins[ix+1],book(tmp,21+ioff,1,iy));}
 47	{Histo1DPtr tmp;_h_AEEC.add(bins[ix],bins[ix+1],book(tmp,25+ioff,1,iy));}
 48	{Histo1DPtr tmp;_h_cone.add(bins[ix],bins[ix+1],book(tmp,29+ioff,1,iy));}
 49	if(ioff==0) {
 50	  Histo1DPtr tmp;_h_thrust.add(bins[iy],bins[iy+1],book(tmp,33,1,iy));
 51	}
 52	++iy;
 53	if(iy==3) {
 54	  ++ioff;
 55	  iy=1;
 56	}
 57      }
 58      // total values
 59      book(_h_EEC_all   , 3,1,1);
 60      book(_h_AEEC_all  , 4,1,1);
 61      book(_h_cone_all  , 5,1,1);
 62      book(_h_thrust_all, 6,1,1);
 63      book(_h_Oblateness, 7,1,1);
 64      book(_h_C         , 8,1,1);
 65      book(_h_heavy     , 9,1,1);
 66      book(_h_sum       ,10,1,1);
 67      book(_h_diff      ,11,1,1);
 68      book(_h_wide      ,12,1,1);
 69      book(_h_total     ,13,1,1);
 70      book(_h_jade      ,17,1,1);
 71      book(_h_dur       ,18,1,1);
 72      book(_h_cam       ,20,1,1);
 73      book(_h_bin,"/TMP/hbin",bins);
 74    }
 75
 76
 77    /// Perform the per-event analysis
 78    void analyze(const Event& event) {
 79      const FinalState& fs = apply<FinalState>(event, "FS");
 80      if ( fs.particles().size() < 2) vetoEvent;
 81      // Get beams and average beam momentum
 82      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 83      // thrust
 84      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 85      // angle bettwen thrust and beam
 86      double cosThrust = abs(beams.first.p3().unit().dot(thrust.thrustAxis()));
 87      _h_bin->fill(cosThrust);
 88      // thrust and related
 89      _h_thrust_all->fill(           1.-thrust.thrust());
 90      _h_thrust     .fill(cosThrust, 1.-thrust.thrust());
 91      _h_Oblateness->fill(thrust.oblateness() );
 92      
 93      // visible energy and make pseudojets
 94      double Evis = 0.0;
 95      PseudoJets pjs;
 96      for (const Particle& p : fs.particles()) {
 97        Evis += p.E();
 98	fastjet::PseudoJet pj = p;
 99	pjs.push_back(pj);
100      }
101      double Evis2 = sqr(Evis);
102      // (A)EEC
103      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
104      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
105        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
106	  if(p_i == p_j) continue; 
107          const Vector3 mom3_i = p_i->momentum().p3();
108          const Vector3 mom3_j = p_j->momentum().p3();
109          const double energy_i = p_i->momentum().E();
110          const double energy_j = p_j->momentum().E();
111          const double thetaij = 180.*mom3_i.unit().angle(mom3_j.unit())/M_PI;
112          double eec = (energy_i*energy_j) / Evis2;
113	  eec *= 2.;
114          _h_EEC_all->fill(           thetaij, eec);
115          _h_EEC     .fill(cosThrust, thetaij, eec);
116          if (thetaij <90.){
117	    _h_AEEC_all->fill(           thetaij, -eec);
118	    _h_AEEC     .fill(cosThrust, thetaij, -eec);
119	  }
120          else {
121	    _h_AEEC_all->fill(          180.-thetaij, eec);
122	    _h_AEEC     .fill(cosThrust,180.-thetaij, eec);
123	  }
124        }
125      }
126      // hemisphere related
127      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
128      _h_heavy->fill(hemi.scaledM2high());
129      _h_diff ->fill(hemi.scaledM2diff());
130      _h_sum  ->fill(hemi.scaledM2low()+hemi.scaledM2high());
131      _h_wide ->fill(hemi.Bmax() );
132      _h_total->fill(hemi.Bsum() );
133      // C-parameter
134      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
135      _h_C->fill(parisi.C());
136      // jets
137      const FastJets&  durjet = apply<FastJets>(event, "DurhamJets");
138      const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
139      if (durjet .clusterSeq()) _h_dur ->fill( durjet.clusterSeq()->exclusive_ymerge_max(2));
140      if (jadejet.clusterSeq()) _h_jade->fill(jadejet.clusterSeq()->exclusive_ymerge_max(2));
141      // Cambridge is more complicated, inclusive defn
142      for (size_t i = 0; i < _h_cam->numBins(); ++i) {
143	double ycut = _h_cam->bin(i).xMax();
144	// 	  double width = _h_y_2_Cambridge->bin(i).width();
145	fastjet::EECambridgePlugin plugin(ycut);
146	fastjet::JetDefinition jdef(&plugin);
147	fastjet::ClusterSequence cseq(pjs, jdef);
148	unsigned int njet = cseq.inclusive_jets().size();
149	if(njet==2) {
150	  _h_cam->fill( _h_cam->bin(i).xMid());
151	  break;
152	}
153      }
154      // jet cone
155      Vector3 jetAxis=thrust.thrustAxis();
156      if(hemi.highMassDirection()) jetAxis *=-1.;
157      for(const Particle & p : fs.particles()) {
158	const double thetaij = 180.*jetAxis.angle(p.p3().unit())/M_PI;
159	double jcef = p.E()/ Evis;
160	_h_cone_all->fill(          thetaij,jcef);
161	_h_cone     .fill(cosThrust,thetaij,jcef);
162      }
163    }
164
165
166    /// Normalise histograms etc., after the run
167    void finalize() {
168      for(unsigned int ix=0;ix<8;++ix) {
169	if(ix<2) scale(_h_thrust.histos()[ix],1./_h_bin->bins()[ix].area());
170	scale(_h_EEC.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
171	scale(_h_AEEC.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
172	scale(_h_cone.histos()[ix],180./M_PI/_h_bin->bins()[ix].area());
173      }
174
175      
176      // _h_thrust.scale(1./sumOfWeights(),this);
177      // _h_EEC.scale(180./M_PI/sumOfWeights(),this);
178      // _h_AEEC.scale(180./M_PI/sumOfWeights(),this);
179      // _h_cone.scale(180./M_PI/sumOfWeights(),this);
180
181
182
183      
184      scale(_h_thrust_all, 1./sumOfWeights());
185      scale(_h_EEC_all, 180./M_PI/sumOfWeights());
186      scale(_h_AEEC_all, 180./M_PI/sumOfWeights());
187      scale(_h_cone_all, 180./M_PI/sumOfWeights());
188      scale(_h_Oblateness, 1./sumOfWeights());
189      scale(_h_C         , 1./sumOfWeights());
190      scale(_h_heavy     , 1./sumOfWeights());
191      scale(_h_sum       , 1./sumOfWeights());
192      scale(_h_diff      , 1./sumOfWeights());
193      scale(_h_wide      , 1./sumOfWeights());
194      scale(_h_total     , 1./sumOfWeights());
195      scale(_h_dur       , 1./sumOfWeights());
196      scale(_h_jade      , 1./sumOfWeights());
197      scale(_h_cam       , 1./sumOfWeights());
198    }
199
200    ///@}
201
202
203    /// @name Histograms
204    ///@{
205    Histo1DPtr _h_thrust_all,_h_EEC_all,_h_AEEC_all,_h_cone_all;
206    Histo1DPtr _h_Oblateness,_h_C,_h_heavy,_h_sum,_h_diff,_h_wide,_h_total;
207    Histo1DPtr _h_jade,_h_dur,_h_cam;
208    BinnedHistogram _h_thrust,_h_EEC,_h_AEEC,_h_cone;
209    Histo1DPtr _h_bin;
210    ///@}
211
212
213  };
214
215
216  RIVET_DECLARE_PLUGIN(DELPHI_2000_I522656);
217
218}