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,324 +1,355 @@ // -*- C++ -*- #include "pluginALICE/HeavyIonAnalysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Tools/Cuts.hh" #include #define _USE_MATH_DEFINES #include namespace Rivet { class ALICE_2012_I930312 : public HeavyIonAnalysis { public: ALICE_2012_I930312() : HeavyIonAnalysis("ALICE_2012_I930312") { } void init() { HeavyIonAnalysis::init(); - + + // Set centrality method 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; - + // Set limit values of pT bins 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]; + // Chaged final states with |eta| < 1.0 and different pT bins for associated particles 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)); } - // 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)"); - + // Create event strings + event_string[0] = "pp"; + event_string[1] = "central"; + event_string[2] = "peripheral"; + event_string[3] = "other"; + + // For each event type + for (int itype = 0; itype < EVENT_TYPES; itype++) { + // For each pT range + for (int ipt = 0; ipt < PT_BINS; ipt++) { + // Initialize yield histograms + _histYield[itype][ipt] = bookHisto1D("Yield_" + event_string[itype] + "_" + 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[itype][ipt] = bookHisto1D("Yield_" + event_string[itype] + "_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)"); + } } + + // Histogram for counting trigger particles for each event type _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(2, 1, 1); - _histIAA[2] = bookHisto1D(5, 1, 1); + // Initialize IAA and ICP histograms + _histIAA[0] = bookScatter2D(1, 1, 1); + _histIAA[1] = bookScatter2D(2, 1, 1); + _histIAA[2] = bookScatter2D(5, 1, 1); - _histIAA[3] = bookHisto1D(3, 1, 1); - _histIAA[4] = bookHisto1D(4, 1, 1); - _histIAA[5] = bookHisto1D(6, 1, 1); + _histIAA[3] = bookScatter2D(3, 1, 1); + _histIAA[4] = bookScatter2D(4, 1, 1); + _histIAA[5] = bookScatter2D(6, 1, 1); + } void analyze(const Event& event) { + // Check if heavy-ion event 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 "; + std::cout << "HI EVENT: "; } else { - std::cout << "PP EVENT "; + std::cout << "PP EVENT: "; } std::cout << event.genEvent()->heavy_ion(); + // Create charged final state for trigger particle const ChargedFinalState& triggerFinalState = applyProjection(event, "CFSTrigger"); Particles triggerParticles = triggerFinalState.particlesByPt(); std::cout << "trig: " << triggerParticles.size(); - + + // Create charged final state for associated particle ChargedFinalState associatedFinalState[PT_BINS]; Particles associatedParticles[PT_BINS]; 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 - string event_string = ""; if (is_heavy_ion) { double centr = centrality(event); std::cout << "centrality: " << centr << std::endl; - if (centr > 0.0 && centr < 5.0) { + if (centr > 0.0 && centr < 5.0) event_type = 1; // PbPb, central - event_string = "PbPb central"; - } - else if (centr > 60.0 && centr < 90.0) { + else if (centr > 60.0 && centr < 90.0) event_type = 2; // PbPb, peripherial - event_string = "PbPb peripheral"; - } - else { + else event_type = 3; // PbPb, other - event_string = "PbPb other"; - } } else { event_type = 0; // pp - event_string = "pp"; } - std::cout << "ev type: " << event_string << " (" << event_type << ")" << std::endl; + std::cout << "ev type: " << event_string[event_type] << " (" << event_type << ")" << std::endl; + // Veto event if not valid event type if (event_type == 3) vetoEvent; + // Fill trigger histogram for a proper event type _histTriggerCounter->fill(event_type, triggerParticles.size()); // Loop over trigger particles foreach (const Particle& triggerParticle, triggerParticles) { std::cout << "Trigger particle" << std::endl; + // For each pt bin 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]) { - if (associatedParticle != triggerParticle) { - if (triggerParticle.pt() > associatedParticle.pt()) { - _histTriggerCounter[event_type]->fill(ipt, 1); - std::cout << "Adding trigger particle" << std::endl; - break; - } - } - }*/ + // Loop over associated particles foreach (const Particle& associatedParticle, associatedParticles[ipt]) { std::cout << "Associated particle" << std::endl; + // If associated and trigger particle are not the same particles if (associatedParticle != triggerParticle) { + // If pt of trigger particle is greater than pt of associated particle if (triggerParticle.pt() > associatedParticle.pt()) { + // Calculate delta phi in range (-0.5*PI, 1.5*PI) 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; } + // Fill yield histogram for calculated delta phi 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 { - /*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])); - } - }*/ + std::cout << "Finalizing..." << std::endl; } void post() { + // Variables for near and away side peak calculation double nearSide[EVENT_TYPES][PT_BINS] = { {0.0} }; double awaySide[EVENT_TYPES][PT_BINS] = { {0.0} }; + // Variables for background error calculation + double background[EVENT_TYPES][PT_BINS] = { {0.0} }; + double backgroundError[EVENT_TYPES][PT_BINS] = { {0.0} }; + + // Variables for integration error calculation + double scalingFactor[EVENT_TYPES] = {0.0}; + double numberOfEntries[EVENT_TYPES][PT_BINS][2] = { { {0.0} } }; + int numberOfBins[EVENT_TYPES][PT_BINS][2] = { { {0} } }; + 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->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->bin(itype).sumW() != 0) { + if ((_histTriggerCounter->bin(itype).sumW() != 0)) { std::cout << "Scaling histogram type " << itype << " pt bin " << ipt << "..." << std::endl; std::cout << "Scaling factor: " << _histTriggerCounter->bin(itype).sumW() << std::endl; + scalingFactor[itype] = 1. / _histTriggerCounter->bin(itype).sumW(); 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; + background[itype][ipt] = sum / nbins; + std::cout << "background: " << background[itype][ipt] << std::endl; + + // Calculate background error + sum = 0.0; + nbins = 0; + for (unsigned int ibin = 0; ibin < _histYield[itype][ipt]->numBins(); ibin++) { + 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)) { + sum += (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]) * + (_histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); + nbins++; + } + } + backgroundError[itype][ipt] = sqrt(sum / (nbins - 1)); + std::cout << "backgroundError: " << backgroundError[itype][ipt] << std::endl; // 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 << "Bin " << ibin << ", value " << _histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt] << std::endl; + _histYieldBkgRemoved[itype][ipt]->fillBin(ibin, _histYield[itype][ipt]->bin(ibin).sumW() - background[itype][ipt]); } 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; + unsigned int lowerBin = _histYield[itype][ipt]->binIndexAt(-0.7); + unsigned int upperBin = _histYield[itype][ipt]->binIndexAt(0.7); + nbins = upperBin - lowerBin; + numberOfBins[itype][ipt][0] = nbins; + nearSide[itype][ipt] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) - nbins * background[itype][ipt]; + numberOfEntries[itype][ipt][0] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) * _histTriggerCounter->bin(itype).sumW(); + std::cout << "_histYield[" << itype << "][" << ipt << "]->integralRange(" << lowerBin << ", " << upperBin << ") - nbins * bkg = " << + _histYield[itype][ipt]->integralRange(lowerBin, upperBin) << " - " << nbins << " * " << background[itype][ipt] << " = " << + 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; + lowerBin = _histYield[itype][ipt]->binIndexAt(M_PI - 0.7); + upperBin = _histYield[itype][ipt]->binIndexAt(M_PI + 0.7); + nbins = upperBin - lowerBin; + numberOfBins[itype][ipt][1] = nbins; + awaySide[itype][ipt] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) - nbins * background[itype][ipt]; + numberOfEntries[itype][ipt][1] = _histYield[itype][ipt]->integralRange(lowerBin, upperBin) * _histTriggerCounter->bin(itype).sumW(); + std::cout << "_histYield[" << itype << "][" << ipt << "]->integralRange(" << lowerBin << ", " << upperBin << ") - nbins * bkg = " << + _histYield[itype][ipt]->integralRange(lowerBin, upperBin) << " - " << nbins << " * " << background[itype][ipt] << " = " << + awaySide[itype][ipt] << std::endl; } } - // Print I_CP results + // Variables for IAA/ICP plots + double dI = 0.0; + int near = 0; + int away = 1; + double xval[PT_BINS] = { 3.5, 5.0, 7.0, 9.0 }; + double xerr[PT_BINS] = { 0.5, 1.0, 1.0, 1.0 }; + + int types1[3] = {1, 2, 1}; + int types2[3] = {0, 0, 2}; + string type_string[3] = {"I_AA 0-5%/pp", "I_AA 60-90%/pp", "I_CP 0-5%/60-90%"}; + + // Fill IAA/ICP plots for near side peak std::cout << "Near side" << std::endl; - for (int ipt = 0; ipt < PT_BINS; 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]); + for (int ihist = 0; ihist < 3; ihist++) { + int type1 = types1[ihist]; + int type2 = types2[ihist]; + for (int ipt = 0; ipt < PT_BINS; ipt++) { + std::cout << type_string[ihist] << ", pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << nearSide[type1][ipt] / nearSide[type2][ipt]; + dI = scalingFactor[type1] * sqrt(numberOfEntries[type1][ipt][near]) / nearSide[type2][ipt] + + scalingFactor[type2] * sqrt(numberOfEntries[type2][ipt][near]) * nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]) + + numberOfBins[type1][ipt][near] * backgroundError[type1][ipt] / nearSide[type2][ipt] + + numberOfBins[type2][ipt][near] * backgroundError[type2][ipt] * nearSide[type1][ipt] / (nearSide[type2][ipt] * nearSide[type2][ipt]); + std::cout << " +- " << dI << std::endl; + _histIAA[ihist]->addPoint(xval[ipt], nearSide[type1][ipt] / nearSide[type2][ipt], xerr[ipt], dI); + } } + + // Fill IAA/ICP plots for away side peak std::cout << "Away side" << std::endl; - for (int ipt = 0; ipt < PT_BINS; 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]); + for (int ihist = 0; ihist < 3; ihist++) { + int type1 = types1[ihist]; + int type2 = types2[ihist]; + for (int ipt = 0; ipt < PT_BINS; ipt++) { + std::cout << type_string[ihist] << ", pt=(" << pt_limits[ipt] << ", " << pt_limits[ipt+1] << "): " << awaySide[type1][ipt] / awaySide[type2][ipt]; + dI = scalingFactor[type1] * sqrt(numberOfEntries[type1][ipt][away]) / awaySide[type2][ipt] + + scalingFactor[type2] * sqrt(numberOfEntries[type2][ipt][away]) * awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]) + + numberOfBins[type1][ipt][away] * backgroundError[type1][ipt] / awaySide[type2][ipt] + + numberOfBins[type2][ipt][away] * backgroundError[type2][ipt] * awaySide[type1][ipt] / (awaySide[type2][ipt] * awaySide[type2][ipt]); + std::cout << " +- " << dI << std::endl; + _histIAA[ihist + 3]->addPoint(xval[ipt], awaySide[type1][ipt] / awaySide[type2][ipt], xerr[ipt], dI); + } } } private: 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; - Histo1DPtr _histIAA[6]; - int num_trig[PT_BINS]; + Scatter2DPtr _histIAA[6]; int pt_limits[5]; int event_type; + string event_string[EVENT_TYPES + 1]; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I930312); }