|
[Rivet-svn] r3605 - in trunk: . src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Mar 15 11:14:16 GMT 2012
Author: hoeth Date: Thu Mar 15 11:14:16 2012 New Revision: 3605 Log: Whitespace consistency in ATLAS_2012_I1084540 Modified: trunk/ChangeLog trunk/src/Analyses/ATLAS_2012_I1084540.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Thu Mar 15 11:03:40 2012 (r3604) +++ trunk/ChangeLog Thu Mar 15 11:14:16 2012 (r3605) @@ -1,6 +1,7 @@ 2012-03-15 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * Make ATLAS_2012_I1084540 less verbose. + * Make ATLAS_2012_I1084540 less verbose, and general code + cleanup of that analysis. * New-style plugin hook in ATLAS_2011_I926145, ATLAS_2011_I944826 and ATLAS_2012_I1084540 Modified: trunk/src/Analyses/ATLAS_2012_I1084540.cc ============================================================================== --- trunk/src/Analyses/ATLAS_2012_I1084540.cc Thu Mar 15 11:03:40 2012 (r3604) +++ trunk/src/Analyses/ATLAS_2012_I1084540.cc Thu Mar 15 11:14:16 2012 (r3605) @@ -1,7 +1,7 @@ // -*- C++ -*- -/** +/** * @name ATLAS Diffractive Gaps Rivet Analysis - * @author Tim Martin, tim.martin at cern.ch + * @author Tim Martin, tim.martin at cern.ch * @version 1.0 * @date 16/01/2012 * @see http://arxiv.org/abs/1201.2808 @@ -15,7 +15,7 @@ * detector at eta = |4.9| which contains zero particles above * pT Cut. Soft diffractive topologies are isolated at large * gap sizes. - * + * */ #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" @@ -25,214 +25,212 @@ namespace Rivet { -class ATLAS_2012_I1084540 : public Analysis { + class ATLAS_2012_I1084540 : public Analysis { public: - /// @name Constructor - //@{ - ATLAS_2012_I1084540() : Analysis("ATLAS_2012_I1084540"), _etaBins(49), _etaMax(4.9) { - setNeedsCrossSection(true); - } - //@} - - /// @name Analysis methods - //@{ - /// Book histograms and initialise projections before the run - void init() { - //All final states. Rapidity range = ATLAS calorimetry. Lowest pT cut = 200 MeV. - const FinalState cnfs2(-_etaMax, _etaMax, 0.2 * GeV); - const FinalState cnfs4(-_etaMax, _etaMax, 0.4 * GeV); - const FinalState cnfs6(-_etaMax, _etaMax, 0.6 * GeV); - const FinalState cnfs8(-_etaMax, _etaMax, 0.8 * GeV); - addProjection(cnfs2, "CNFS2"); - addProjection(cnfs4, "CNFS4"); - addProjection(cnfs6, "CNFS6"); - addProjection(cnfs8, "CNFS8"); - - _etaBinSize = (2. * _etaMax)/(double)_etaBins; - - //Book histogram - _h_DeltaEtaF_200 = bookHistogram1D(1, 1, 1); - _h_DeltaEtaF_400 = bookHistogram1D(2, 1, 1); - _h_DeltaEtaF_600 = bookHistogram1D(3, 1, 1); - _h_DeltaEtaF_800 = bookHistogram1D(4, 1, 1); - } - - void fillMap(const FinalState& fs, bool* energyMap, double pTcut) { - // Fill true/false array by iterating over particles and compare their - // pT with pTcut - foreach (const Particle& p, fs.particlesByEta()) { - int checkBin = -1; - double checkEta = -_etaMax; - while (1) { - checkEta += _etaBinSize; - ++checkBin; - if (p.momentum().eta() < checkEta) { - energyMap[checkBin] = (p.momentum().pT() > pTcut * GeV); - break; - } - } - } - } - - /// Perform the per-event analysis - void analyze(const Event& event) { - static unsigned int event_count = 0; - ++event_count; - const double weight = event.weight(); - const FinalState& fs2 = applyProjection<FinalState>(event, "CNFS2"); - const FinalState& fs4 = applyProjection<FinalState>(event, "CNFS4"); - const FinalState& fs6 = applyProjection<FinalState>(event, "CNFS6"); - const FinalState& fs8 = applyProjection<FinalState>(event, "CNFS8"); - - // Set up Yes/No arrays for energy in each eta bin at each pT cut - bool* energyMap_200 = new bool[_etaBins]; - bool* energyMap_400 = new bool[_etaBins]; - bool* energyMap_600 = new bool[_etaBins]; - bool* energyMap_800 = new bool[_etaBins]; - for (int i = 0; i < _etaBins; ++i) { - energyMap_200[i] = false; - energyMap_400[i] = false; - energyMap_600[i] = false; - energyMap_800[i] = false; + ATLAS_2012_I1084540() : Analysis("ATLAS_2012_I1084540"), _etaBins(49), _etaMax(4.9) { + setNeedsCrossSection(true); } - // Veto bins based on final state particles > Cut (Where Cut = 200 - 800 MeV pT) - fillMap(fs2, energyMap_200, 0.2); - fillMap(fs4, energyMap_400, 0.4); - fillMap(fs6, energyMap_600, 0.6); - fillMap(fs8, energyMap_800, 0.8); - - // Apply gap finding algorithm - // Detector layout follows... - // -Eta [Proton --- DetectorCSide --- DetectorBarrel --- DetectorASide --- Proton] +Eta - bool gapDirectionAt200 = false; //False is gap on C size, True is gap on A side. - double largestEdgeGap_200 = 0.; - double largestEdgeGap_400 = 0.; - double largestEdgeGap_600 = 0.; - double largestEdgeGap_800 = 0.; - - for (int E = 200; E <= 800; E += 200) { - double EdgeGapSizeA = -1, EdgeGapSizeC = -1; - bool* energyMap = 0; - switch (E) { - case 200: energyMap = energyMap_200; break; - case 400: energyMap = energyMap_400; break; - case 600: energyMap = energyMap_600; break; - case 800: energyMap = energyMap_800; break; - } - // Look left to right - for (int a = 0; a < _etaBins; ++a) { - if (energyMap[a] == true) { - EdgeGapSizeA = (_etaBinSize * a); - break; - } - } + /// @name Analysis methods + //@{ + /// Book histograms and initialise projections before the run + void init() { + //All final states. Rapidity range = ATLAS calorimetry. Lowest pT cut = 200 MeV. + const FinalState cnfs2(-_etaMax, _etaMax, 0.2 * GeV); + const FinalState cnfs4(-_etaMax, _etaMax, 0.4 * GeV); + const FinalState cnfs6(-_etaMax, _etaMax, 0.6 * GeV); + const FinalState cnfs8(-_etaMax, _etaMax, 0.8 * GeV); + addProjection(cnfs2, "CNFS2"); + addProjection(cnfs4, "CNFS4"); + addProjection(cnfs6, "CNFS6"); + addProjection(cnfs8, "CNFS8"); + + _etaBinSize = (2. * _etaMax)/(double)_etaBins; + + //Book histogram + _h_DeltaEtaF_200 = bookHistogram1D(1, 1, 1); + _h_DeltaEtaF_400 = bookHistogram1D(2, 1, 1); + _h_DeltaEtaF_600 = bookHistogram1D(3, 1, 1); + _h_DeltaEtaF_800 = bookHistogram1D(4, 1, 1); + } - // And look right to left - for (int c = _etaBins-1; c >= 0; --c) { - if (energyMap[c] == true) { - EdgeGapSizeC = (2 * _etaMax) - (_etaBinSize * (c+1)); - if (fuzzyEquals(EdgeGapSizeC, 4.47035e-08)) EdgeGapSizeC = 0.0; - break; - } - } - // Put your hands on your hips - - // Find the largest gap - double largestEdgeGap = 0.; - if (E == 200) { - // If the 200 MeV pass, take the biggest of the two gaps. Make note of which side for higher pT cuts. - largestEdgeGap = std::max(EdgeGapSizeA,EdgeGapSizeC); - gapDirectionAt200 = (EdgeGapSizeA > EdgeGapSizeC); - } else { - // Use the direction from 200 MeV pass, most accurate measure of which side gap is on. - if (gapDirectionAt200) { - largestEdgeGap = EdgeGapSizeA; + void fillMap(const FinalState& fs, bool* energyMap, double pTcut) { + // Fill true/false array by iterating over particles and compare their + // pT with pTcut + foreach (const Particle& p, fs.particlesByEta()) { + int checkBin = -1; + double checkEta = -_etaMax; + while (1) { + checkEta += _etaBinSize; + ++checkBin; + if (p.momentum().eta() < checkEta) { + energyMap[checkBin] = (p.momentum().pT() > pTcut * GeV); + break; + } } - else largestEdgeGap = EdgeGapSizeC; } - - // Check case of empty detector - if (largestEdgeGap < 0.0) largestEdgeGap = 2.0 * _etaMax; - - // Fill bin centre - switch (E) { - case 200: _h_DeltaEtaF_200->fill(largestEdgeGap + _etaBinSize/2., weight); break; - case 400: _h_DeltaEtaF_400->fill(largestEdgeGap + _etaBinSize/2., weight); break; - case 600: _h_DeltaEtaF_600->fill(largestEdgeGap + _etaBinSize/2., weight); break; - case 800: _h_DeltaEtaF_800->fill(largestEdgeGap + _etaBinSize/2., weight); break; - } - - if (E == 200) largestEdgeGap_200 = largestEdgeGap; - if (E == 400) largestEdgeGap_400 = largestEdgeGap; - if (E == 600) largestEdgeGap_600 = largestEdgeGap; - if (E == 800) largestEdgeGap_800 = largestEdgeGap; } - // Algorithm result every 1000 events - if (event_count % 1000 == 0) { + /// Perform the per-event analysis + void analyze(const Event& event) { + static unsigned int event_count = 0; + ++event_count; + const double weight = event.weight(); + const FinalState& fs2 = applyProjection<FinalState>(event, "CNFS2"); + const FinalState& fs4 = applyProjection<FinalState>(event, "CNFS4"); + const FinalState& fs6 = applyProjection<FinalState>(event, "CNFS6"); + const FinalState& fs8 = applyProjection<FinalState>(event, "CNFS8"); + + // Set up Yes/No arrays for energy in each eta bin at each pT cut + bool* energyMap_200 = new bool[_etaBins]; + bool* energyMap_400 = new bool[_etaBins]; + bool* energyMap_600 = new bool[_etaBins]; + bool* energyMap_800 = new bool[_etaBins]; + for (int i = 0; i < _etaBins; ++i) { + energyMap_200[i] = false; + energyMap_400[i] = false; + energyMap_600[i] = false; + energyMap_800[i] = false; + } + + // Veto bins based on final state particles > Cut (Where Cut = 200 - 800 MeV pT) + fillMap(fs2, energyMap_200, 0.2); + fillMap(fs4, energyMap_400, 0.4); + fillMap(fs6, energyMap_600, 0.6); + fillMap(fs8, energyMap_800, 0.8); + + // Apply gap finding algorithm + // Detector layout follows... + // -Eta [Proton --- DetectorCSide --- DetectorBarrel --- DetectorASide --- Proton] +Eta + bool gapDirectionAt200 = false; //False is gap on C size, True is gap on A side. + double largestEdgeGap_200 = 0.; + double largestEdgeGap_400 = 0.; + double largestEdgeGap_600 = 0.; + double largestEdgeGap_800 = 0.; + for (int E = 200; E <= 800; E += 200) { + double EdgeGapSizeA = -1, EdgeGapSizeC = -1; bool* energyMap = 0; - double largestEdgeGap = 0; switch (E) { - case 200: energyMap = energyMap_200; largestEdgeGap = largestEdgeGap_200; break; - case 400: energyMap = energyMap_400; largestEdgeGap = largestEdgeGap_400; break; - case 600: energyMap = energyMap_600; largestEdgeGap = largestEdgeGap_600; break; - case 800: energyMap = energyMap_800; largestEdgeGap = largestEdgeGap_800; break; - } - MSG_DEBUG("Largest Forward Gap at pT Cut " << E << " MeV=" << largestEdgeGap - << " eta, NFinalState pT > 200 in ATLAS acceptance:" << fs2.particles().size()); - std::string hitPattern = "Detector HitPattern=-4.9["; + case 200: energyMap = energyMap_200; break; + case 400: energyMap = energyMap_400; break; + case 600: energyMap = energyMap_600; break; + case 800: energyMap = energyMap_800; break; + } + + // Look left to right for (int a = 0; a < _etaBins; ++a) { - if (energyMap[a] == true) hitPattern += "X"; - else hitPattern += "_"; + if (energyMap[a] == true) { + EdgeGapSizeA = (_etaBinSize * a); + break; + } + } + + // And look right to left + for (int c = _etaBins-1; c >= 0; --c) { + if (energyMap[c] == true) { + EdgeGapSizeC = (2 * _etaMax) - (_etaBinSize * (c+1)); + if (fuzzyEquals(EdgeGapSizeC, 4.47035e-08)) EdgeGapSizeC = 0.0; + break; + } + } + // Put your hands on your hips + + // Find the largest gap + double largestEdgeGap = 0.; + if (E == 200) { + // If the 200 MeV pass, take the biggest of the two gaps. Make note of which side for higher pT cuts. + largestEdgeGap = std::max(EdgeGapSizeA,EdgeGapSizeC); + gapDirectionAt200 = (EdgeGapSizeA > EdgeGapSizeC); + } else { + // Use the direction from 200 MeV pass, most accurate measure of which side gap is on. + if (gapDirectionAt200) { + largestEdgeGap = EdgeGapSizeA; + } + else largestEdgeGap = EdgeGapSizeC; + } + + // Check case of empty detector + if (largestEdgeGap < 0.0) largestEdgeGap = 2.0 * _etaMax; + + // Fill bin centre + switch (E) { + case 200: _h_DeltaEtaF_200->fill(largestEdgeGap + _etaBinSize/2., weight); break; + case 400: _h_DeltaEtaF_400->fill(largestEdgeGap + _etaBinSize/2., weight); break; + case 600: _h_DeltaEtaF_600->fill(largestEdgeGap + _etaBinSize/2., weight); break; + case 800: _h_DeltaEtaF_800->fill(largestEdgeGap + _etaBinSize/2., weight); break; + } + + if (E == 200) largestEdgeGap_200 = largestEdgeGap; + if (E == 400) largestEdgeGap_400 = largestEdgeGap; + if (E == 600) largestEdgeGap_600 = largestEdgeGap; + if (E == 800) largestEdgeGap_800 = largestEdgeGap; + } + + // Algorithm result every 1000 events + if (event_count % 1000 == 0) { + for (int E = 200; E <= 800; E += 200) { + bool* energyMap = 0; + double largestEdgeGap = 0; + switch (E) { + case 200: energyMap = energyMap_200; largestEdgeGap = largestEdgeGap_200; break; + case 400: energyMap = energyMap_400; largestEdgeGap = largestEdgeGap_400; break; + case 600: energyMap = energyMap_600; largestEdgeGap = largestEdgeGap_600; break; + case 800: energyMap = energyMap_800; largestEdgeGap = largestEdgeGap_800; break; + } + MSG_DEBUG("Largest Forward Gap at pT Cut " << E << " MeV=" << largestEdgeGap + << " eta, NFinalState pT > 200 in ATLAS acceptance:" << fs2.particles().size()); + std::string hitPattern = "Detector HitPattern=-4.9["; + for (int a = 0; a < _etaBins; ++a) { + if (energyMap[a] == true) hitPattern += "X"; + else hitPattern += "_"; + } + hitPattern += "]4.9"; + MSG_DEBUG(hitPattern); + std::string gapArrow = " "; + if (!gapDirectionAt200) { + int drawSpaces = (int)(_etaBins - (largestEdgeGap/_etaBinSize) + 0.5); + for (int i = 0; i < drawSpaces; ++i) gapArrow += " "; + } + int drawArrows = (int)((largestEdgeGap/_etaBinSize) + 0.5); + for (int i = 0; i < drawArrows; ++i) gapArrow += "^"; + MSG_DEBUG(gapArrow); } - hitPattern += "]4.9"; - MSG_DEBUG(hitPattern); - std::string gapArrow = " "; - if (!gapDirectionAt200) { - int drawSpaces = (int)(_etaBins - (largestEdgeGap/_etaBinSize) + 0.5); - for (int i = 0; i < drawSpaces; ++i) gapArrow += " "; - } - int drawArrows = (int)((largestEdgeGap/_etaBinSize) + 0.5); - for (int i = 0; i < drawArrows; ++i) gapArrow += "^"; - MSG_DEBUG(gapArrow); } + + delete[] energyMap_200; + delete[] energyMap_400; + delete[] energyMap_600; + delete[] energyMap_800; } - delete[] energyMap_200; - delete[] energyMap_400; - delete[] energyMap_600; - delete[] energyMap_800; - } - - /// Normalise histograms after the run, Scale to cross section - void finalize() { - MSG_DEBUG("Cross Section=" << crossSection() / millibarn << "mb, SumOfWeights=" << sumOfWeights()); - scale(_h_DeltaEtaF_200, (crossSection() / millibarn)/sumOfWeights()); - scale(_h_DeltaEtaF_400, (crossSection() / millibarn)/sumOfWeights()); - scale(_h_DeltaEtaF_600, (crossSection() / millibarn)/sumOfWeights()); - scale(_h_DeltaEtaF_800, (crossSection() / millibarn)/sumOfWeights()); - } - //@} + /// Normalise histograms after the run, Scale to cross section + void finalize() { + MSG_DEBUG("Cross Section=" << crossSection() / millibarn << "mb, SumOfWeights=" << sumOfWeights()); + scale(_h_DeltaEtaF_200, (crossSection() / millibarn)/sumOfWeights()); + scale(_h_DeltaEtaF_400, (crossSection() / millibarn)/sumOfWeights()); + scale(_h_DeltaEtaF_600, (crossSection() / millibarn)/sumOfWeights()); + scale(_h_DeltaEtaF_800, (crossSection() / millibarn)/sumOfWeights()); + } + //@} private: - /// @name Histograms - //@{ - AIDA::IHistogram1D* _h_DeltaEtaF_200; - AIDA::IHistogram1D* _h_DeltaEtaF_400; - AIDA::IHistogram1D* _h_DeltaEtaF_600; - AIDA::IHistogram1D* _h_DeltaEtaF_800; - //@} - /// @name Private variables - //@{ - int _etaBins; - double _etaMax; - double _etaBinSize; - //@} -}; + /// @name Histograms + //@{ + AIDA::IHistogram1D* _h_DeltaEtaF_200; + AIDA::IHistogram1D* _h_DeltaEtaF_400; + AIDA::IHistogram1D* _h_DeltaEtaF_600; + AIDA::IHistogram1D* _h_DeltaEtaF_800; + //@} + /// @name Private variables + //@{ + int _etaBins; + double _etaMax; + double _etaBinSize; + //@} + }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2012_I1084540);
More information about the Rivet-svn mailing list |