H1_2000_S4129130.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/Math/Constants.hh"
00005 #include "Rivet/Tools/ParticleIdUtils.hh"
00006 #include "Rivet/Projections/FinalState.hh"
00007 #include "Rivet/Projections/DISKinematics.hh"
00008 
00009 namespace Rivet {
00010 
00011 
00012   /// @brief H1 energy flow and charged particle spectra
00013   /// @author Peter Richardson
00014   /// Based on the HZtool analysis hz99091
00015   class H1_2000_S4129130 : public Analysis {
00016   public:
00017 
00018     /// Constructor
00019     H1_2000_S4129130() : Analysis("H1_2000_S4129130")
00020     {
00021       setBeams(POSITRON, PROTON);
00022     }
00023  
00024  
00025     /// @name Analysis methods
00026     //@{
00027 
00028     void analyze(const Event& event) {
00029       // Get the projections
00030       const FinalState & fs = applyProjection<FinalState>(event, "FS");
00031       const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
00032       const DISLepton    & dl = applyProjection<DISLepton>(event,"Lepton");
00033 
00034       // Get the DIS kinematics
00035       double q2  = dk.Q2();
00036       double x   = dk.x();
00037       double y   = dk.y();
00038       double w2  = dk.W2();
00039    
00040       // Momentum of the scattered lepton
00041       FourMomentum leptonMom = dl.out().momentum();
00042       // pT energy and angle
00043       double enel = leptonMom.E();
00044       double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
00045 
00046       // Extract the particles other than the lepton
00047       ParticleVector particles;
00048       particles.reserve(fs.particles().size());
00049       const GenParticle& dislepGP = dl.out().genParticle();
00050       for (ParticleVector::const_iterator p = fs.particles().begin();
00051            p != fs.particles().end(); ++p) {
00052         const GenParticle& loopGP = p->genParticle();
00053         if (&loopGP == &dislepGP) continue;
00054         particles.push_back(*p);
00055       }
00056    
00057       // Cut on the forward energy
00058       double efwd = 0.;
00059       foreach (const Particle& p, particles) {
00060         double th = 180.-p.momentum().angle(dl.in().momentum())/degree;
00061         if (th > 4.4 && th < 15.0) efwd += p.momentum().E();
00062       }
00063       // There are four possible selections for events
00064       bool evcut[4];
00065       // Low  Q2 selection a
00066       evcut[0] = enel/GeV > 12. && w2 >= 4400.*GeV2 && efwd/GeV > 0.5 &&
00067         inRange(thel,157.,176.);
00068       // Low  Q2 selection b
00069       evcut[1] = enel/GeV > 12. && inRange(y,0.3,0.5);
00070       // High Q2 selection a
00071       evcut[2] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) &&
00072         w2 >= 4400.*GeV2 && efwd > 0.5;
00073       // High Q2 selection b
00074       evcut[3] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) &&
00075     inRange(w2,27110.*GeV2,45182.*GeV2);
00076   
00077       // Veto if fails all cuts
00078       if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3]) ) {
00079         vetoEvent;
00080       }
00081    
00082       // Find the bins
00083       int bin[4] = {-1,-1,-1,-1};
00084       // For the low Q2 selection a)
00085       if (q2 > 2.5*GeV && q2 <= 5.*GeV) {
00086         if (x > 0.00005 && x <= 0.0001 ) bin[0] = 0;
00087         if (x > 0.0001  && x <= 0.0002 ) bin[0] = 1;
00088         if (x > 0.0002  && x <= 0.00035) bin[0] = 2;
00089         if (x > 0.00035 && x <= 0.0010 ) bin[0] = 3;
00090       }
00091       else if(q2 > 5.*GeV && q2 <= 10.*GeV) {
00092         if (x > 0.0001  && x <= 0.0002 ) bin[0] = 4;
00093         if (x > 0.0002  && x <= 0.00035) bin[0] = 5;
00094         if (x > 0.00035 && x <= 0.0007 ) bin[0] = 6;
00095         if (x > 0.0007  && x <= 0.0020 ) bin[0] = 7;
00096       }
00097       else if(q2 > 10.*GeV && q2 <= 20.*GeV) {
00098         if (x > 0.0002 && x <= 0.0005) bin[0] = 8;
00099         if (x > 0.0005 && x <= 0.0008) bin[0] = 9;
00100         if (x > 0.0008 && x <= 0.0015) bin[0] = 10;
00101         if (x > 0.0015 && x <= 0.040 ) bin[0] = 11;
00102       }
00103       else if(q2 > 20.*GeV && q2 <= 50.*GeV) {
00104         if (x > 0.0005 && x <= 0.0014) bin[0] = 12;
00105         if (x > 0.0014 && x <= 0.0030) bin[0] = 13;
00106         if (x > 0.0030 && x <= 0.0100) bin[0] = 14;
00107       }
00108       else if (q2 > 50.*GeV && q2 <= 100.*GeV) {
00109         if (x >0.0008 && x <= 0.0030) bin[0] = 15;
00110         if (x >0.0030 && x <= 0.0200) bin[0] = 16;
00111       }
00112       // check in one of the bins
00113       evcut[0] &= bin[0] >= 0;
00114       // For the low Q2 selection b)
00115       if (q2 > 2.5*GeV && q2 <= 5.  *GeV) bin[1] = 0;
00116       if (q2 > 5. *GeV && q2 <= 10. *GeV) bin[1] = 1;
00117       if (q2 > 10.*GeV && q2 <= 20. *GeV) bin[1] = 2;
00118       if (q2 > 20.*GeV && q2 <= 50. *GeV) bin[1] = 3;
00119       if (q2 > 50.*GeV && q2 <= 100.*GeV) bin[1] = 4;
00120       // check in one of the bins
00121       evcut[1] &= bin[1] >= 0;
00122       // for the high Q2 selection a)
00123       if (q2 > 100.*GeV && q2 <= 400.*GeV) {
00124         if (x > 0.00251 && x <= 0.00631) bin[2] = 0;
00125         if (x > 0.00631 && x <= 0.0158 ) bin[2] = 1;
00126         if (x > 0.0158  && x <= 0.0398 ) bin[2] = 2;
00127       }
00128       else if (q2 > 400.*GeV && q2 <= 1100.*GeV) {
00129         if (x > 0.00631 && x <= 0.0158 ) bin[2] = 3;
00130         if (x > 0.0158  && x <= 0.0398 ) bin[2] = 4;
00131         if (x > 0.0398  && x <= 1.     ) bin[2] = 5;
00132       }
00133       else if (q2 > 1100.*GeV && q2 <= 100000.*GeV) {
00134         if (x > 0. && x <= 1.) bin[2] = 6;
00135       }
00136       // check in one of the bins
00137       evcut[2] &= bin[2] >= 0;
00138       // for the high Q2 selection b)
00139       if      (q2 > 100.*GeV && q2 <= 220.*GeV) bin[3] = 0;
00140       else if (q2 > 220.*GeV && q2 <= 400.*GeV) bin[3] = 1;
00141       else if (q2 > 400.              ) bin[3] = 2;
00142       // check in one of*GeV the bins
00143       evcut[3] &= bin[3] >= 0;
00144    
00145       // Veto if fails all cuts after bin selection
00146       if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3])) {
00147         vetoEvent;
00148       }
00149    
00150       // Increment the count for normalisation
00151       const double weight = event.weight();
00152       if (evcut[0]) _weightETLowQa [bin[0]] += weight;
00153       if (evcut[1]) _weightETLowQb [bin[1]] += weight;
00154       if (evcut[2]) _weightETHighQa[bin[2]] += weight;
00155       if (evcut[3]) _weightETHighQb[bin[3]] += weight;
00156    
00157       // Boost to hadronicCM
00158       const LorentzTransform hcmboost = dk.boostHCM();
00159    
00160       // Loop over the particles
00161       double etcent = 0;
00162       double etfrag = 0;
00163       foreach (const Particle& p, particles) {
00164         // Boost momentum to CMS
00165         const FourMomentum hcmMom = hcmboost.transform(p.momentum());
00166         double et = fabs(Et(hcmMom));
00167         double eta = hcmMom.pseudorapidity();
00168         // Averages in central and forward region
00169         if (fabs(eta) < .5 ) etcent += et;
00170         if (eta > 2 && eta <= 3.) etfrag += et;
00171         // Histograms of Et flow
00172         if (evcut[0]) _histETLowQa [bin[0]]->fill(eta, et*weight);
00173         if (evcut[1]) _histETLowQb [bin[1]]->fill(eta, et*weight);
00174         if (evcut[2]) _histETHighQa[bin[2]]->fill(eta, et*weight);
00175         if (evcut[3]) _histETHighQb[bin[3]]->fill(eta, et*weight);
00176       }
00177       // Fill histograms for the average quantities
00178       if (evcut[1] || evcut[3]) {
00179         _histAverETCentral->fill(q2, etcent*weight,weight);
00180         _histAverETFrag   ->fill(q2, etfrag*weight,weight);
00181       }
00182     }
00183  
00184  
00185     void init() {
00186       // Projections
00187       addProjection(DISLepton(), "Lepton");
00188       addProjection(DISKinematics(), "Kinematics");
00189       addProjection(FinalState(), "FS");
00190    
00191       // Histos
00192       IHistogram1D* h = 0;
00193    
00194       // Histograms and weight vectors for low Q^2 a
00195       _histETLowQa.reserve(17);
00196       _weightETLowQa.reserve(17);
00197       for (size_t ix = 0; ix < 17; ++ix) {
00198         h = bookHistogram1D(ix+1, 1, 1);
00199         _histETLowQa.push_back(h);
00200         _weightETLowQa.push_back(0.);
00201       }
00202    
00203       // Histograms and weight vectors for high Q^2 a
00204       _histETHighQa.reserve(7);
00205       _weightETHighQa.reserve(7);
00206       for (size_t ix = 0; ix < 7; ++ix) {
00207         h = bookHistogram1D(ix+18, 1, 1);
00208         _histETHighQa.push_back(h);
00209         _weightETHighQa.push_back(0.);
00210       }
00211    
00212       // Histograms and weight vectors for low Q^2 b
00213       _histETLowQb.reserve(5);
00214       _weightETLowQb.reserve(5);
00215       for (size_t ix = 0; ix < 5; ++ix) {
00216         h = bookHistogram1D(ix+25, 1, 1);
00217         _histETLowQb.push_back(h);
00218         _weightETLowQb.push_back(0.);
00219       }
00220    
00221       // Histograms and weight vectors for high Q^2 b
00222       _histETHighQb.reserve(3);
00223       _weightETHighQb.reserve(3);
00224       for (size_t ix = 0; ix < 3; ++ix) {
00225         h = bookHistogram1D(30+ix, 1, 1);
00226         _histETHighQb.push_back(h);
00227         _weightETHighQb.push_back(0.0);
00228       }
00229    
00230       // Histograms for the averages
00231       _histAverETCentral = bookProfile1D(33,  1, 1);
00232       _histAverETFrag = bookProfile1D(34,  1, 1);
00233     }
00234  
00235  
00236     // Finalize
00237     void finalize() {
00238       // Normalization of the Et distributions
00239       for (size_t ix=0; ix<17; ++ix) {
00240         scale(_histETLowQa[ix], 1./_weightETLowQa[ix]);
00241       }
00242       for(size_t ix=0; ix<7; ++ix) {
00243         scale(_histETHighQa[ix], 1./_weightETHighQa[ix]);
00244       }
00245       for(size_t ix=0; ix<5; ++ix) {
00246         scale(_histETLowQb[ix], 1./_weightETLowQb[ix]);
00247       }
00248       for(size_t ix=0; ix<3; ++ix) {
00249         scale(_histETHighQb[ix], 1./_weightETHighQb[ix]);
00250       }
00251     }
00252  
00253 
00254     //@}
00255 
00256 
00257   private:
00258  
00259     /// @name Histograms
00260     //@{
00261     vector<AIDA::IHistogram1D *> _histETLowQa;
00262     vector<AIDA::IHistogram1D *> _histETHighQa;
00263     vector<AIDA::IHistogram1D *> _histETLowQb;
00264     vector<AIDA::IHistogram1D *> _histETHighQb;
00265     AIDA::IProfile1D * _histAverETCentral;
00266     AIDA::IProfile1D * _histAverETFrag;
00267     //@}
00268 
00269     /// @name storage of weights for normalisation
00270     //@{
00271     vector<double> _weightETLowQa;
00272     vector<double> _weightETHighQa;
00273     vector<double> _weightETLowQb;
00274     vector<double> _weightETHighQb;
00275     //@}
00276   };
00277 
00278 
00279 
00280   // This global object acts as a hook for the plugin system
00281   AnalysisBuilder<H1_2000_S4129130> plugin_H1_2000_S4129130;
00282 
00283 }