|
[Rivet-svn] r4069 - in branches/2012-06-aidarivet: . data/anainfo data/plotinfo src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgFri Dec 7 13:21:06 GMT 2012
Author: buckley Date: Fri Dec 7 13:21:06 2012 New Revision: 4069 Log: A rewrite of ALICE_2012_I1181770 -- needs to be checked by original authors to make sure logic is intact Modified: branches/2012-06-aidarivet/ChangeLog branches/2012-06-aidarivet/data/anainfo/ALICE_2012_I1181770.info branches/2012-06-aidarivet/data/plotinfo/ALICE_2012_I1181770.plot branches/2012-06-aidarivet/src/Analyses/ALICE_2012_I1181770.cc Modified: branches/2012-06-aidarivet/ChangeLog ============================================================================== --- branches/2012-06-aidarivet/ChangeLog Fri Dec 7 13:19:31 2012 (r4068) +++ branches/2012-06-aidarivet/ChangeLog Fri Dec 7 13:21:06 2012 (r4069) @@ -1,5 +1,9 @@ 2012-12-07 Andy Buckley <andy.buckley at cern.ch> + * Version number bump to 1.8.2 -- release approaching. + + * Rewrite of ALICE_2012_I1181770 analysis to make it a bit more sane and acceptable. + * Adding a note on FourVector and FourMomentum that operator- and operator-= invert both the space and time components: use of -= can result in a vector with negative energy. Modified: branches/2012-06-aidarivet/data/anainfo/ALICE_2012_I1181770.info ============================================================================== --- branches/2012-06-aidarivet/data/anainfo/ALICE_2012_I1181770.info Fri Dec 7 13:19:31 2012 (r4068) +++ branches/2012-06-aidarivet/data/anainfo/ALICE_2012_I1181770.info Fri Dec 7 13:21:06 2012 (r4069) @@ -10,18 +10,21 @@ - Martin Poghosyan <Martin.Poghosyan at cern.ch> - Sercan Sen <Sercan.Sen at cern.ch> - Burak Bilki <bbilki at gmail.com> + - Andy Buckley <andy.buckley at cern.ch> References: - arXiv:1208.4968 [hep-ex] RunInfo: - Inelastic events (non-diffractive and inelastic diffractive). + Inelastic events (non-diffractive and inelastic diffractive). NumEvents: 100K Beams: [p+, p+] Energies: [900, 2760, 7000] Description: - 'Measurements of cross sections of inelastic and diffractive processes in proton-proton collisions at $\sqrt{s}$=900, 2760 and 7000 GeV. - The fractions of diffractive processes in inelastic collisions were determined from a study of gaps in charged particle pseudorapidity distributions. - Single-diffractive events are selected with $M_{X} < 200$ GeV $/c^{2}$ and double-diffractive events defined as NSD events with $\Delta\eta >3$. - To measure the inelastic cross section, beam properties were determined with van der Meer scans using a simulation of diffraction adjusted to data.' + 'Measurements of cross-sections of inelastic and diffractive processes in proton-proton collisions at $\sqrt{s} = 900$, 2760 and 7000 GeV. + The fractions of diffractive processes in inelastic collisions were determined from a study of gaps in charged particle pseudorapidity distributions. + Single-diffractive events are selected with $M_{X} < 200\;\GeV/c^2$ and double-diffractive events defined as NSD events with $\Delta\eta > 3$. + To measure the inelastic cross-section, beam properties were determined with van der Meer scans using a simulation of diffraction adjusted to data. + + Note that these are experimental approximations to theoretical concepts -- it's not totally clear whether the data point values are model-independent.' BibKey: :2012sja BibTeX: '@article{:2012sja, author = "Abelev, Betty and others", Modified: branches/2012-06-aidarivet/data/plotinfo/ALICE_2012_I1181770.plot ============================================================================== --- branches/2012-06-aidarivet/data/plotinfo/ALICE_2012_I1181770.plot Fri Dec 7 13:19:31 2012 (r4068) +++ branches/2012-06-aidarivet/data/plotinfo/ALICE_2012_I1181770.plot Fri Dec 7 13:21:06 2012 (r4069) @@ -3,7 +3,7 @@ # END PLOT # BEGIN PLOT /ALICE_2012_I1181770/d01-x01-y0 -Title=Production ratios of SD with $M_{X} < 200$ GeV$/c^{2}$ to INEL +Title=Production ratios of SD with $M_{X} < 200\;\GeV/c^2$ to INEL YLabel=$\sigma_\text{SD} / \sigma_\text{inel}$ # END PLOT @@ -13,16 +13,16 @@ # END PLOT # BEGIN PLOT /ALICE_2012_I1181770/d03-x01-y0 -Title=Single diffraction cross-section for $M_{X} < 200$ GeV$/c^{2}$ +Title=Single diffraction cross-section for $M_{X} < 200\;\GeV/c^2$ YLabel=$\sigma_\text{SD}$ [mb] # END PLOT # BEGIN PLOT /ALICE_2012_I1181770/d04-x01-y0 -Title=Double diffraction cross-section for $\Delta\eta >3$ +Title=Double diffraction cross-section for $\Delta\eta > 3$ YLabel=$\sigma_\text{DD}$ [mb] # END PLOT # BEGIN PLOT /ALICE_2012_I1181770/d05-x01-y0 -Title=Inleastic cross-section +Title=Inelastic cross-section YLabel=$\sigma_\text{inel}$ [mb] # END PLOT Modified: branches/2012-06-aidarivet/src/Analyses/ALICE_2012_I1181770.cc ============================================================================== --- branches/2012-06-aidarivet/src/Analyses/ALICE_2012_I1181770.cc Fri Dec 7 13:19:31 2012 (r4068) +++ branches/2012-06-aidarivet/src/Analyses/ALICE_2012_I1181770.cc Fri Dec 7 13:21:06 2012 (r4069) @@ -6,6 +6,7 @@ namespace Rivet { + class ALICE_2012_I1181770 : public Analysis { public: @@ -13,170 +14,104 @@ : Analysis("ALICE_2012_I1181770") { } - public: void init() { - addProjection(ChargedFinalState(),"CFS"); + // Projection setup + addProjection(ChargedFinalState(), "CFS"); - if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { - _h_xsec_sd = bookHistogram1D(3, 1, 1); - _h_xsec_dd = bookHistogram1D(4, 1, 1); - _h_xsec_inel = bookHistogram1D(5, 1, 1); - } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { - _h_xsec_sd = bookHistogram1D(3, 1, 2); - _h_xsec_dd = bookHistogram1D(4, 1, 2); - _h_xsec_inel = bookHistogram1D(5, 1, 2); - } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { - _h_xsec_sd = bookHistogram1D(3, 1, 3); - _h_xsec_dd = bookHistogram1D(4, 1, 3); - _h_xsec_inel = bookHistogram1D(5, 1, 3); - } + // Book (energy-specific) histograms + _isqrts = -1; + if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) _isqrts = 1; + else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) _isqrts = 2; + else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) _isqrts = 3; + assert(_isqrts > 0); + _h_xsec_sd = bookHistogram1D(3, 1, _isqrts); + _h_xsec_dd = bookHistogram1D(4, 1, _isqrts); + _h_xsec_inel = bookHistogram1D(5, 1, _isqrts); } - void analyze(const Event& event) { - const double weight = event.weight(); + void analyze(const Event& event) { const ChargedFinalState& cfs = applyProjection<ChargedFinalState>(event, "CFS"); + if (cfs.size() < 2) vetoEvent; // need at least two particles to calculate gaps - // fill INEL plots for each event - if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { - _h_xsec_inel->fill( 900/GeV, weight ); - } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { - _h_xsec_inel->fill( 2760/GeV, weight ); - } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { - _h_xsec_inel->fill( 7000/GeV, weight ); - } - - double Eslowest(0.0), pslowest(0.0), yslowest(999.); - double Efastest(0.0), pfastest(0.0), yfastest(-999.); - int pidslowest(0), pidfastest(0); - - double LRG(0.0), etapre(0.0), gapbwd(0.0), gapfwd(0.0); - unsigned int num_p(0); - - FourMomentum leadP(0.,0.,0.,0.); - double Elead(0.0), plead(0.0); - - foreach(const Particle& p, cfs.particlesByEta()) { //sorted from minus to plus - const PdgId pid = p.pdgId(); - double y = p.momentum().rapidity(); - double eta = p.momentum().eta(); - //SD case - if (y < yslowest) { - Eslowest = p.momentum().E(); - pslowest = p.momentum().vector3().mod(); - yslowest = p.momentum().rapidity(); - pidslowest = pid; - } - if (y > yfastest) { - Efastest = p.momentum().E(); - pfastest = p.momentum().vector3().mod(); - yfastest = p.momentum().rapidity(); - pidfastest = pid; - } - - num_p += 1; - // DD case - if (num_p==1) { - etapre = p.momentum().eta(); - } else if (num_p > 1) { - if (num_p==2) gapbwd = fabs(eta-etapre); - - double gap = fabs(eta-etapre); - LRG = (gap > LRG ? gap : LRG); // largest gap - - if (num_p==cfs.size()) gapfwd = fabs(eta-etapre); - etapre = eta; - } - } + const double weight = event.weight(); - // Mx calculation - if (pidslowest==2212 && pidfastest==2212) { - if (fabs(yslowest) > fabs(yfastest)) { - Elead = Eslowest; - plead = pslowest; - } else if (fabs(yslowest) < fabs(yfastest)) { - Elead = Efastest; - plead = pfastest; - } else { - Elead = Eslowest; - plead = pslowest; - if ( (double)rand() / (double)RAND_MAX > 0.5) { // generate random number in [0.,1.] range and make decision randomly - Elead = Efastest; - plead = pfastest; - } - } - } else if (pidslowest==2212) { - Elead = Eslowest; - plead = pslowest; - } else if (pidfastest==2212) { - Elead = Efastest; - plead = pfastest; - } + // Fill INEL plots for each event + _h_xsec_inel->fill(sqrtS()/GeV, weight); - double Mx = sqrt((sqrtS()/GeV-Elead-plead)*(sqrtS()/GeV-Elead+plead)); - bool singleDiff = false; + // Identify particles with most positive/most negative rapidities + const ParticleVector particlesByRap = cfs.particlesByRapidity(); + const Particle pslowest = particlesByRap.front(); + const Particle pfastest = particlesByRap.back(); + + // Find gap sizes + const ParticleVector particlesByEta = cfs.particlesByEta(); // sorted from minus to plus + const size_t num_particles = particlesByEta.size(); + vector<double> gaps; + for (size_t ip = 1; ip < num_particles; ++ip) { + const Particle& p1 = particlesByEta[ip-1]; + const Particle& p2 = particlesByEta[ip]; + const double gap = p2.momentum().eta() - p1.momentum().eta(); + assert(gap >= 0); + gaps.push_back(gap); + } + + // First, last, and largest gaps + const double gapmax = *max_element(gaps.begin(), gaps.end()); + const double gapbwd = gaps.front(); + const double gapfwd = gaps.back(); - // Fill SD - if (Mx < 200.) { - singleDiff = true; - if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { - _h_xsec_sd->fill( 900/GeV, weight); - } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { - _h_xsec_sd->fill( 2760/GeV, weight); - } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { - _h_xsec_sd->fill( 7000/GeV, weight); - } + // Mx calculation + FourMomentum p4lead; + if (pslowest.pdgId() == PROTON && pfastest.pdgId() == PROTON) { + p4lead = (fabs(pslowest.momentum().rapidity()) > fabs(pfastest.momentum().rapidity())) ? pslowest.momentum() : pfastest.momentum(); + } else if (pslowest.pdgId() == PROTON) { + p4lead = pslowest.momentum(); + } else if (pfastest.pdgId() == PROTON) { + p4lead = pfastest.momentum(); + } + const double Mx = sqrt( (sqrtS()-p4lead.E()-p4lead.vector3().mod()) * (sqrtS()-p4lead.E()+p4lead.vector3().mod()) ); + + // Fill SD (and escape) if Mx is sufficiently low + if (Mx < 200*GeV) { + _h_xsec_sd->fill(sqrtS()/GeV, weight); + return; } - if ( singleDiff ) vetoEvent; // DD events are defined as NSD with large gap. - - // also remove SD-like events in NSD events - if ( std::abs(gapbwd-LRG) < std::numeric_limits<double>::epsilon() || std::abs(gapfwd-LRG) < std::numeric_limits<double>::epsilon() ) vetoEvent; + // Also remove SD-like events in NSD events + if (fuzzyEquals(gapbwd, gapmax) || fuzzyEquals(gapfwd, gapmax)) vetoEvent; // Fill DD plots - if (LRG > 3.) { - if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { - _h_xsec_dd->fill( 900/GeV, weight); - } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { - _h_xsec_dd->fill( 2760/GeV, weight); - } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { - _h_xsec_dd->fill( 7000/GeV, weight); - } - } + if (gapmax > 3) _h_xsec_dd->fill(sqrtS()/GeV, weight); } - void finalize() { - // get the ratio plots: SD/inel, DD/inel + void finalize() { + // Make the ratio plots: SD/inel, DD/inel const string dir = histoDir(); + const string ypart = "y0" + lexical_cast<string>(_isqrts); + histogramFactory().divide( dir + "/d01-x01-" + ypart, *_h_xsec_sd, *_h_xsec_inel); + histogramFactory().divide( dir + "/d02-x01-" + ypart, *_h_xsec_dd, *_h_xsec_inel); - if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { - histogramFactory().divide( dir + "/d01-x01-y01", *_h_xsec_sd , *_h_xsec_inel ); - histogramFactory().divide( dir + "/d02-x01-y01", *_h_xsec_dd , *_h_xsec_inel ); - } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { - histogramFactory().divide( dir + "/d01-x01-y02", *_h_xsec_sd , *_h_xsec_inel ); - histogramFactory().divide( dir + "/d02-x01-y02", *_h_xsec_dd , *_h_xsec_inel ); - } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { - histogramFactory().divide( dir + "/d01-x01-y03", *_h_xsec_sd , *_h_xsec_inel ); - histogramFactory().divide( dir + "/d02-x01-y03", *_h_xsec_dd , *_h_xsec_inel ); - } - - scale(_h_xsec_sd, crossSection()/millibarn/sumOfWeights()); - scale(_h_xsec_dd, crossSection()/millibarn/sumOfWeights()); + scale(_h_xsec_sd, crossSection()/millibarn/sumOfWeights()); + scale(_h_xsec_dd, crossSection()/millibarn/sumOfWeights()); scale(_h_xsec_inel, crossSection()/millibarn/sumOfWeights()); - } + private: AIDA::IHistogram1D *_h_xsec_sd; AIDA::IHistogram1D *_h_xsec_dd; AIDA::IHistogram1D *_h_xsec_inel; + /// Energy index for histogram IDs + int _isqrts; + }; - // Hook for the plugin system + DECLARE_RIVET_PLUGIN(ALICE_2012_I1181770); }
More information about the Rivet-svn mailing list |