diff --git a/analyses/pluginALICE/ALICE_2012_I930312.cc b/analyses/pluginALICE/ALICE_2012_I930312.cc --- a/analyses/pluginALICE/ALICE_2012_I930312.cc +++ b/analyses/pluginALICE/ALICE_2012_I930312.cc @@ -1,321 +1,324 @@ // -*- C++ -*- -#include "Rivet/Analysis.hh" +#include "pluginALICE/HeavyIonAnalysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Tools/Cuts.hh" #include #define _USE_MATH_DEFINES #include -#define PT_BINS 4 -#define EVENT_TYPES 3 - namespace Rivet { - double quantile(double value, const Histo1DPtr hist); - - class ALICE_2012_I930312 : public Analysis { + class ALICE_2012_I930312 : public HeavyIonAnalysis { public: ALICE_2012_I930312() : - Analysis("ALICE_2012_I930312") + HeavyIonAnalysis("ALICE_2012_I930312") { } void init() { + HeavyIonAnalysis::init(); + + setCentralityMethod(HeavyIonAnalysis::ImpactParameter, 10000); // Charged final states with |eta| < 1.0 and 8 < pT < 15 GeV/c for trigger particles const Cut& cutTrigger = Cuts::abseta < 1.0 && Cuts::pT > 8*GeV && Cuts::pT < 15*GeV; const ChargedFinalState cfsTrigger(cutTrigger); addProjection(cfsTrigger,"CFSTrigger"); // Chaged final states with |eta| < 1.0 and 3 < pT < 15 GeV/c for associated particles for (int ipt = 0; ipt < PT_BINS; ipt++) num_trig[ipt] = 0; pt_limits[0] = 3; pt_limits[1] = 4; pt_limits[2] = 6; pt_limits[3] = 8; pt_limits[4] = 10; // Cut cutAssoc[4]; // ChargedFinalState cfsAssoc[4]; for (int ipt = 0; ipt < PT_BINS; ipt++) { Cut mycut = Cuts::abseta < 1.0 && Cuts::pT > pt_limits[ipt]*GeV && Cuts::pT < pt_limits[ipt + 1]*GeV; addProjection(ChargedFinalState(mycut), "CFSAssoc" + std::to_string(ipt)); //ChargedFinalState mycfs = ChargedFinalState(mycut); // const ChargedFinalState &mycfs(mycut); //addProjection(mycfs, "CFSAssoc" + std::to_string(ipt)); } - // Impact parameter distribution histogram for centrality calculation. If the histogram with the same name - // is loaded in readData() function, it will bind such preloaded histogram to the following histogram. - // Otherwise, empty histogram will be created - // @note the name was not discussed yet, this is only an example: calib_mult... - std::cout << "Initialization..." << std::endl; - _histCalibration = bookHisto1D("calib_impactpar_(eta<0.5)&&(pT>50*MeV)", 200, 0, 20.0, "Calibration histogram", "xlabel", "ylabel"); - _histControl = bookHisto1D("control_impactpar_(eta<0.5)&&(pT>50*MeV)", 200, 0, 20.0, "Control histogram", "xlabel", "ylabel"); - // Histograms and variables initialization. Do it for each centrality range for (int ipt = 0; ipt < PT_BINS; ipt++) { _histYield[0][ipt] = bookHisto1D("Yield_pp_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); _histYield[1][ipt] = bookHisto1D("Yield_central_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); _histYield[2][ipt] = bookHisto1D("Yield_peripheral_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); _histYieldBkgRemoved[0][ipt] = bookHisto1D("Yield_pp_nobkg_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield no bkg", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); _histYieldBkgRemoved[1][ipt] = bookHisto1D("Yield_central_nobkg_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield no bkg", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); _histYieldBkgRemoved[2][ipt] = bookHisto1D("Yield_peripheral_nobkg_" + std::to_string(ipt), 36, -0.5 * M_PI, 1.5 * M_PI, "Associated particle per trigger particle yield no bkg", "#Delta#eta (rad)", "1 / N_trig dN_assoc / d#Delta#eta (rad^-1)"); } - _histTriggerCounter[0] = bookHisto1D("Trigger_pp", PT_BINS, 0.0, PT_BINS, "Trigger counter pp", "pt bin", "N"); - _histTriggerCounter[1] = bookHisto1D("Trigger_central", PT_BINS, 0.0, PT_BINS, "Trigger counter central", "pt bin", "N"); - _histTriggerCounter[2] = bookHisto1D("Trigger_peripheral", PT_BINS, 0.0, PT_BINS, "Trigger counter peripheral", "pt bin", "N"); + _histTriggerCounter = bookHisto1D("Trigger", EVENT_TYPES, 0.0, EVENT_TYPES, "Trigger counter", "event type", "N"); + //_histTriggerCounter[0] = bookHisto1D("Trigger_pp", PT_BINS, 0.0, PT_BINS, "Trigger counter pp", "pt bin", "N"); + //_histTriggerCounter[1] = bookHisto1D("Trigger_central", PT_BINS, 0.0, PT_BINS, "Trigger counter central", "pt bin", "N"); + //_histTriggerCounter[2] = bookHisto1D("Trigger_peripheral", PT_BINS, 0.0, PT_BINS, "Trigger counter peripheral", "pt bin", "N"); _histIAA[0] = bookHisto1D(1, 1, 1); - _histIAA[1] = bookHisto1D(4, 1, 1); + _histIAA[1] = bookHisto1D(2, 1, 1); _histIAA[2] = bookHisto1D(5, 1, 1); - _histIAA[3] = bookHisto1D(6, 1, 1); - + + _histIAA[3] = bookHisto1D(3, 1, 1); + _histIAA[4] = bookHisto1D(4, 1, 1); + _histIAA[5] = bookHisto1D(6, 1, 1); } void analyze(const Event& event) { - if (event.genEvent()->heavy_ion()) { - std::cout << "HI EVENT" << std::endl; + std::pair beam_pair = event.genEvent()->beam_particles(); + std::cout << "PDG ID 1: " << (std::get<0>(beam_pair))->pdg_id() << std::endl; + std::cout << "PDG ID 2: " << (std::get<1>(beam_pair))->pdg_id() << std::endl; + bool is_heavy_ion = ((std::get<0>(beam_pair))->pdg_id() == 1000822080) && ((std::get<1>(beam_pair))->pdg_id() == 1000822080); + if (is_heavy_ion) { + std::cout << "HI EVENT "; } else { - std::cout << "PP EVENT" << std::endl; + std::cout << "PP EVENT "; } - std::cout << event.genEvent()->heavy_ion() << std::endl; + std::cout << event.genEvent()->heavy_ion(); const ChargedFinalState& triggerFinalState = applyProjection(event, "CFSTrigger"); Particles triggerParticles = triggerFinalState.particlesByPt(); - std::cout << "Trigger particles: " << triggerParticles.size() << std::endl; + std::cout << "trig: " << triggerParticles.size(); ChargedFinalState associatedFinalState[PT_BINS]; Particles associatedParticles[PT_BINS]; - std::cout << "Associated particles: "; + std::cout << ", assoc: "; for (int ipt = 0; ipt < PT_BINS; ipt++) { associatedFinalState[ipt] = applyProjection(event, "CFSAssoc" + std::to_string(ipt)); associatedParticles[ipt] = associatedFinalState[ipt].particlesByPt(); std::cout << associatedParticles[ipt].size() << " "; } std::cout << std::endl; // Check event type - if (event.genEvent()->heavy_ion()) { - // For PbPb events fill calibration and control histogram - if (_histCalibration->numEntries() >= 10000) { - _histControl->fill(event.genEvent()->heavy_ion() ? event.genEvent()->heavy_ion()->impact_parameter() : -1.0, 1); - } else { - _histCalibration->fill(event.genEvent()->heavy_ion() ? event.genEvent()->heavy_ion()->impact_parameter() : -1.0, 1); - // Veto event in case there are not enough events in calibration histogram - vetoEvent; + string event_string = ""; + if (is_heavy_ion) { + double centr = centrality(event); + std::cout << "centrality: " << centr << std::endl; + if (centr > 0.0 && centr < 5.0) { + event_type = 1; // PbPb, central + event_string = "PbPb central"; } - double centr = quantile(event.genEvent()->heavy_ion()->impact_parameter(), _histCalibration) * 100.0; - if (centr > 0.0 && centr < 5.0) - event_type = 1; // PbPb, central - else if (centr > 60.0 && centr < 90.0) + else if (centr > 60.0 && centr < 90.0) { event_type = 2; // PbPb, peripherial - else + event_string = "PbPb peripheral"; + } + else { event_type = 3; // PbPb, other + event_string = "PbPb other"; + } } else { event_type = 0; // pp + event_string = "pp"; } - std::cout << "Event type: " << event_type << std::endl; + std::cout << "ev type: " << event_string << " (" << event_type << ")" << std::endl; if (event_type == 3) vetoEvent; + _histTriggerCounter->fill(event_type, triggerParticles.size()); + // Loop over trigger particles foreach (const Particle& triggerParticle, triggerParticles) { - - for (int ipt = 0; ipt < PT_BINS; ipt++) { - if (associatedParticles[ipt].size() > 1 || + std::cout << "Trigger particle" << std::endl; + for (int ipt = 0; ipt < PT_BINS; ipt++) { + std::cout << "pt bin " << ipt << std::endl; + /*if (associatedParticles[ipt].size() > 1 || (associatedParticles[ipt].size() == 1 && associatedParticles[ipt].at(0) != triggerParticle)) { num_trig[ipt] += 1; _histTriggerCounter[event_type]->fill(ipt, 1); - std::cout << "Found some associated particles" << std::endl; - foreach (const Particle& associatedParticle, associatedParticles[ipt]) { + std::cout << "Found some associated particles" << std::endl;*/ + /*foreach (const Particle& associatedParticle, associatedParticles[ipt]) { + if (associatedParticle != triggerParticle) { + if (triggerParticle.pt() > associatedParticle.pt()) { + _histTriggerCounter[event_type]->fill(ipt, 1); + std::cout << "Adding trigger particle" << std::endl; + break; + } + } + }*/ + foreach (const Particle& associatedParticle, associatedParticles[ipt]) { + std::cout << "Associated particle" << std::endl; + if (associatedParticle != triggerParticle) { if (triggerParticle.pt() > associatedParticle.pt()) { double dPhi = triggerParticle.phi() - associatedParticle.phi(); while (dPhi > 1.5 * M_PI) { dPhi -= 2 * M_PI; } while (dPhi < -0.5 * M_PI) { dPhi += 2 * M_PI; } std::cout << "Filling yield histogram for pt = " << pt_limits[ipt] << "-" << pt_limits[ipt + 1] << " GeV/c for delta phi = " << dPhi << std::endl; _histYield[event_type][ipt]->fill(dPhi, 1); } } + //} } - else { + //else { + /*if (trigger_added == false) { std::cout << "No associated particles with pt = " << pt_limits[ipt] << "-" << pt_limits[ipt + 1] << " GeV/c found for that trigger particle" << std::endl; - } + }*/ } } + std::cout << std::endl; } void finalize() { std::cout << "Finalising..." << std::endl; // Normalize all histograms by number of trigger particles used to create each one of them /*for (int ipt = 0; ipt < PT_BINS; ipt++) { std::cout << "finalize(): number of trigger particles: " << num_trig[ipt] << std::endl; if (num_trig[ipt] != 0) { scale(_histYield[event_type][ipt], (1. / num_trig[ipt])); } }*/ } void post() { double nearSide[EVENT_TYPES][PT_BINS] = { {0.0} }; double awaySide[EVENT_TYPES][PT_BINS] = { {0.0} }; - std::cout << "Calibration entries: " << _histCalibration->numEntries() << std::endl; - std::cout << "Control entries: " << _histControl->numEntries() << std::endl; - std::cout << "pp trigger counter histogram entries: " <<_histTriggerCounter[0]->numEntries() << std::endl; - std::cout << "Central trigger counter histogram entries: " << _histTriggerCounter[1]->numEntries() << std::endl; - std::cout << "Peripheral trigger counter histogram entries: " << _histTriggerCounter[2]->numEntries() << std::endl; + std::cout << "Trigger counter histogram entries: " << + _histTriggerCounter->bin(0).sumW() << " " << + _histTriggerCounter->bin(1).sumW() << " " << + _histTriggerCounter->bin(2).sumW() << std::endl; + //std::cout << "pp trigger counter histogram entries: " << _histTriggerCounter[0]->numEntries() << std::endl; + //std::cout << "Central trigger counter histogram entries: " << _histTriggerCounter[1]->numEntries() << std::endl; + //std::cout << "Peripheral trigger counter histogram entries: " << _histTriggerCounter[2]->numEntries() << std::endl; std::cout << "pp yield pt bin 0, entries: " << _histYield[0][0]->numEntries() << std::endl; std::cout << "pp yield pt bin 1, entries: " << _histYield[0][1]->numEntries() << std::endl; std::cout << "pp yield pt bin 2, entries: " << _histYield[0][2]->numEntries() << std::endl; std::cout << "pp yield pt bin 3, entries: " << _histYield[0][3]->numEntries() << std::endl; std::cout << "Central yield pt bin 0, entries: " << _histYield[1][0]->numEntries() << std::endl; std::cout << "Central yield pt bin 1, entries: " << _histYield[1][1]->numEntries() << std::endl; std::cout << "Central yield pt bin 2, entries: " << _histYield[1][2]->numEntries() << std::endl; std::cout << "Central yield pt bin 3, entries: " << _histYield[1][3]->numEntries() << std::endl; std::cout << "Peripheral yield pt bin 0, entries: " << _histYield[2][0]->numEntries() << std::endl; std::cout << "Peripheral yield pt bin 1, entries: " << _histYield[2][1]->numEntries() << std::endl; std::cout << "Peripheral yield pt bin 2, entries: " << _histYield[2][2]->numEntries() << std::endl; std::cout << "Peripheral yield pt bin 3, entries: " << _histYield[2][3]->numEntries() << std::endl; // For each event type for (int itype = 0; itype < EVENT_TYPES; itype++) { // For each pT range for (int ipt = 0; ipt < PT_BINS; ipt++) { // Check if histograms are fine - if (_histTriggerCounter[itype]->numEntries() == 0 || _histYield[itype][ipt]->numEntries() == 0) { - std::cout << _histTriggerCounter[itype]->numEntries() << std::endl; + if (_histTriggerCounter->numEntries() == 0 || _histYield[itype][ipt]->numEntries() == 0) { + std::cout << _histTriggerCounter->numEntries() << std::endl; std::cout << _histYield[itype][ipt]->numEntries() << std::endl; std::cout << "There are no entries in one of the histograms" << std::endl; continue; } // Scale yield histogram std::cout << "Scaling yield histograms..." << std::endl; - if (_histTriggerCounter[itype]->bin(ipt).sumW() != 0) { + if (_histTriggerCounter->bin(itype).sumW() != 0) { std::cout << "Scaling histogram type " << itype << " pt bin " << ipt << "..." << std::endl; - scale(_histYield[itype][ipt], (1. / _histTriggerCounter[itype]->bin(ipt).sumW())); + std::cout << "Scaling factor: " << _histTriggerCounter->bin(itype).sumW() << std::endl; + scale(_histYield[itype][ipt], (1. / _histTriggerCounter->bin(itype).sumW())); } // Calculate background std::cout << "Calculating background" << std::endl; double sum = 0.0; int nbins = 0; for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { std::cout << "Bin " << ibin << std::endl; if (_histYield[itype][ipt]->bin(ibin).xMid() > (0.5 * M_PI - 0.4) && _histYield[itype][ipt]->bin(ibin).xMid() < (0.5 * M_PI + 0.4)) { std::cout << "Adding " << _histYield[itype][ipt]->bin(ibin).sumW() << std::endl; sum += _histYield[itype][ipt]->bin(ibin).sumW(); nbins++; } } if (nbins == 0) { std::cout << "Failed to estimate background!" << std::endl; continue; } std::cout << "Dividing " << sum << " / " << nbins << std::endl; double bkg = sum / nbins; // Fill histograms with removed background std::cout << "Filling histogram with removed background..." << std::endl; for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { std::cout << "Bin " << ibin << ", value " << _histYield[itype][ipt]->bin(ibin).sumW() - bkg << std::endl; _histYieldBkgRemoved[itype][ipt]->fillBin(ibin, _histYield[itype][ipt]->bin(ibin).sumW() - bkg); } std::cout << "Integrating the whole histogram: " << _histYieldBkgRemoved[itype][ipt]->integral() << std::endl; // Integrate near-side yield std::cout << "Integrating near-side yield..." << std::endl; unsigned int lowerBin = _histYieldBkgRemoved[itype][ipt]->binIndexAt(-0.7); unsigned int upperBin = _histYieldBkgRemoved[itype][ipt]->binIndexAt(0.7); nearSide[itype][ipt] = _histYieldBkgRemoved[itype][ipt]->integralRange(lowerBin, upperBin); std::cout << "Integrating bins (" << lowerBin << ", " << upperBin << "): " << nearSide[itype][ipt] << std::endl; // Integrate away-side yield std::cout << "Integrating away-side yield..." << std::endl; lowerBin = _histYieldBkgRemoved[itype][ipt]->binIndexAt(M_PI - 0.7); upperBin = _histYieldBkgRemoved[itype][ipt]->binIndexAt(M_PI + 0.7); awaySide[itype][ipt] = _histYieldBkgRemoved[itype][ipt]->integralRange(lowerBin, upperBin); std::cout << "Integrating bins (" << lowerBin << ", " << upperBin << "): " << awaySide[itype][ipt] << std::endl; } } // Print I_CP results std::cout << "Near side" << std::endl; for (int ipt = 0; ipt < PT_BINS; ipt++) { - std::cout << "I_AA, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[0][ipt] / nearSide[1][ipt] << std::endl; - _histIAA[0]->fillBin(ipt, nearSide[0][ipt] / nearSide[1][ipt]); - std::cout << "I_CP, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[1][ipt] / nearSide[2][ipt] << std::endl; - _histIAA[1]->fillBin(ipt, nearSide[1][ipt] / nearSide[2][ipt]); + std::cout << "I_AA 0-5%/pp, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[1][ipt] / nearSide[0][ipt] << std::endl; + _histIAA[0]->fillBin(ipt, nearSide[1][ipt] / nearSide[0][ipt]); + std::cout << "I_AA 60-90%/pp, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[2][ipt] / nearSide[0][ipt] << std::endl; + _histIAA[1]->fillBin(ipt, nearSide[2][ipt] / nearSide[0][ipt]); + std::cout << "I_CP 0-5%/60-90%, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[1][ipt] / nearSide[2][ipt] << std::endl; + _histIAA[2]->fillBin(ipt, nearSide[1][ipt] / nearSide[2][ipt]); } std::cout << "Away side" << std::endl; for (int ipt = 0; ipt < PT_BINS; ipt++) { - std::cout << "I_AA, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[0][ipt] / awaySide[1][ipt] << std::endl; - _histIAA[2]->fillBin(ipt, awaySide[0][ipt] / awaySide[1][ipt]); - std::cout << "I_CP, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[1][ipt] / awaySide[2][ipt] << std::endl; - _histIAA[3]->fillBin(ipt, awaySide[1][ipt] / awaySide[2][ipt]); + std::cout << "I_AA 0-5%/pp, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[1][ipt] / awaySide[0][ipt] << std::endl; + _histIAA[3]->fillBin(ipt, awaySide[1][ipt] / awaySide[0][ipt]); + std::cout << "I_AA 60-90%/pp, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[2][ipt] / awaySide[0][ipt] << std::endl; + _histIAA[4]->fillBin(ipt, awaySide[2][ipt] / awaySide[0][ipt]); + std::cout << "I_CP 0-5%/60-90%, pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[1][ipt] / awaySide[2][ipt] << std::endl; + _histIAA[5]->fillBin(ipt, awaySide[1][ipt] / awaySide[2][ipt]); } } private: - Histo1DPtr _histCalibration; - Histo1DPtr _histControl; + static const int PT_BINS = 4; + static const int EVENT_TYPES = 3; + Histo1DPtr _histYield[EVENT_TYPES][PT_BINS]; Histo1DPtr _histYieldBkgRemoved[EVENT_TYPES][PT_BINS]; - Histo1DPtr _histTriggerCounter[EVENT_TYPES]; - Histo1DPtr _histIAA[4]; + Histo1DPtr _histTriggerCounter; + Histo1DPtr _histIAA[6]; int num_trig[PT_BINS]; int pt_limits[5]; int event_type; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I930312); - - // Function for calculating the quantile from histogram at selected value - // (below this value, including overflow) - // @note this function is here temporarily, it will be moved and modified - // to create a general function for calculating the quantile of a histogram - double quantile(double value, const Histo1DPtr hist) { - - // Check if there are entries in the histogram - if (hist->numEntries() == 0) { - throw WeightError("There are no entires in the histogram!"); - } - - // Integration ranges - size_t upperBin = hist->binIndexAt(value); - - // Calculate centrality as percentile - return hist->integralTo(upperBin, true) / hist->integral(true); - - } - } diff --git a/analyses/pluginALICE/HeavyIonAnalysis.cc b/analyses/pluginALICE/HeavyIonAnalysis.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/HeavyIonAnalysis.cc @@ -0,0 +1,149 @@ +// -*- C++ -*- +#include "pluginALICE/HeavyIonAnalysis.hh" + +namespace Rivet { + HeavyIonAnalysis::HeavyIonAnalysis(const std::string &name) : + Analysis(name), + _centmethod(Undefined) + { } + + void HeavyIonAnalysis::setCentralityMethod(const CentralityMethod method, + const size_t numEventsRequired, + const FinalState *fs) { + switch (method) { + + case ImpactParameter: + std::cout << "Method: impact parameter" << std::endl; + _methodName = "ImpactParameter"; + _histCentralityCalibration = bookHisto1D("calib_impactpar", 200, 0, 20.0, "Calibration histogram", "xlabel", "ylabel"); + std::cout << "Calibration histogram entries: " << _histCentralityCalibration->numEntries() << std::endl; + _histCentralityControl = bookHisto1D("control_impactpar", 200, 0, 20.0, "Control histogram", "xlabel", "ylabel"); + break; + + case Multiplicity: + if (!fs) + throw Exception("something bad happened ..."); + addProjection(*fs, "FS"); + _methodName = "Multiplicity"; + _histCentralityCalibration = bookHisto1D("calib_multiplicity_" + trimString(fs->getCuts()->description()), 10000, 0, 10000.0, "Calibration histogram", "xlabel", "ylabel"); + _histCentralityControl = bookHisto1D("control_multiplicity_" + trimString(fs->getCuts()->description()), 10000, 0, 10000.0, "Control histogram", "xlabel", "ylabel"); + break; + + default: + centmethod_not_found(); + } + + _numEventsRequired = numEventsRequired; + _centmethod = method; + } + + void HeavyIonAnalysis::centmethod_not_found() const { + throw Exception("Unimplemented method of centrality calibration."); + } + + double HeavyIonAnalysis::centrality(const Event& event) const { + // Initialize value of the event observable to be filled to the histograms + if (_centmethod == Undefined) + throw Exception("no centrality method defined"); + + const float observable = calculateObservable(event); + float centrality = -1.; + + // Check if there are enough events in the calibration histogram + if (_histCentralityCalibration->numEntries() >= _numEventsRequired) { + // Calculate centrality as percentile using observable calculated with the selected method + centrality = quantile(observable, _histCentralityCalibration) * 100.; + if (_centmethod == Multiplicity) + centrality = 100. - centrality; + // Fill the control histogram with the impact parameter value and event weight + MSG_INFO("Adding control point nr " << _histCentralityControl->numEntries() << ": (" << observable << ", " << event.weight() << ")"); + _histCentralityControl->fill(observable, event.weight()); + } + // Otherwise, if there are not enough events in the calibration histogram + else { + // Fill the calibration histogram with the impact parameter value and event weight + MSG_INFO("Adding calibration point nr " << _histCentralityCalibration->numEntries() << ": (" << observable << ", " << event.weight() << ")"); + _histCentralityCalibration->fill(observable, event.weight()); + } + return centrality; + } + + float HeavyIonAnalysis::calculateObservable(const Event& e) const { + + // Initialize observable + float observable = -1.; + + // Calculate its value according to the selected method + switch (_centmethod) { + + case ImpactParameter: + // Get impact parameter of the event + observable = e.genEvent()->heavy_ion() ? e.genEvent()->heavy_ion()->impact_parameter() : -1.; + break; + + case Multiplicity: + { + // Get multiplicity of the event + const FinalState& fs = applyProjection(e, "FS"); + observable = e.genEvent()->heavy_ion() ? fs.particles().size() : -1.; + break; + } + + default: + centmethod_not_found(); + } + + // Check if observable is correct. If not, skip this event + checkObservable(observable); + + return observable; + } + + + void HeavyIonAnalysis::checkObservable(const float observable) const { + + try { + // Check if the observable is greater than 0 + if (observable < 0.) + throw UserError("Calculated observable is lower than 0!"); + // Check if the observable is inside histogram ranges + if ((observable < _histCentralityCalibration->xMin()) || + (observable > _histCentralityCalibration->xMax()) || + (_histCentralityCalibration->binIndexAt(observable) < 0)) + throw RangeError("Calculated observable is out of bounds!"); + } catch (...) { + MSG_WARNING("Skipping this event..."); + } + + } + + string HeavyIonAnalysis::trimString(const string& str) const { + // @note Which characters should be removed? + string trimmedString = str; + + const char characters[] = "\\ "; + for (const auto character : characters) { + trimmedString.erase (std::remove (trimmedString.begin(), trimmedString.end(), character), trimmedString.end()); + } + return trimmedString; + } + + // Function for calculating the quantile from histogram at selected value + // (below this value, including overflow) + // @note this function is here temporarily, it will be moved and modified + // to create a general function for calculating the quantile of a histogram + double HeavyIonAnalysis::quantile(double value, const Histo1DPtr hist) const { + + // Check if there are entries in the histogram + if (hist->numEntries() == 0) { + throw WeightError("There are no entires in the histogram!"); + } + + // Integration ranges + size_t upperBin = hist->binIndexAt(value); + + // Calculate centrality as percentile + return hist->integralTo(upperBin, true) / hist->integral(true); + + } +} diff --git a/analyses/pluginALICE/HeavyIonAnalysis.hh b/analyses/pluginALICE/HeavyIonAnalysis.hh new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/HeavyIonAnalysis.hh @@ -0,0 +1,54 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" + +namespace Rivet { + + class HeavyIonAnalysis : public Analysis { + + public: + HeavyIonAnalysis(const std::string &name = ""); + + // methods specific to heavy-ion analyses + enum CentralityMethod { + ImpactParameter, + Multiplicity, + Undefined + }; + + double quantile(double value, const Histo1DPtr hist) const; + + /// centrality implementation + void setCentralityMethod(const CentralityMethod method, + const size_t numEventsRequired, + const FinalState *fs = 0x0); + + double centrality(const Event& event) const; + + void centmethod_not_found() const; + + float calculateObservable(const Event &event) const; + + void checkObservable(const float observable) const; + + private: + string trimString(const string& str) const; + + /// Histogram for centrality calibration + Histo1DPtr _histCentralityCalibration; + /// Histogram for centrality control. It may be used to compare distribution + /// in the current run to the one provided in calibration histogram. + Histo1DPtr _histCentralityControl; + + /// String with the cuts + string _cutString; + + /// Method of the centrality calibration + CentralityMethod _centmethod; + + /// Number of events required for a selected method of centrality calibration + size_t _numEventsRequired; + + /// Name of the centrality calibration method + string _methodName; + }; +}