diff --git a/analyses/pluginALICE/ALICE_2010_I880049.cc b/analyses/pluginALICE/ALICE_2010_I880049.cc --- a/analyses/pluginALICE/ALICE_2010_I880049.cc +++ b/analyses/pluginALICE/ALICE_2010_I880049.cc @@ -1,77 +1,77 @@ // -*- 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_2010_I880049 : public HeavyIonAnalysis { public: ALICE_2010_I880049() : HeavyIonAnalysis("ALICE_2010_I880049") { } void init() { HeavyIonAnalysis::init(); // Set centrality method - addCentralityMethod(HeavyIonAnalysis::ImpactParameter, 10000, "method1"); + addCentralityMethod(HeavyIonAnalysis::ImpactParameter, 10000, "ImpactParameterMethod"); // Charged final states with |eta| < 0.5 and pT > 50 MeV const Cut& cut = Cuts::abseta < 0.5 && Cuts::pT > 50*MeV; const ChargedFinalState cfs(cut); addProjection(cfs,"CFS"); // Histograms and variables initialization. Do it for each centrality range _histNchVsCentr = bookHisto1D(1, 1, 1); } void analyze(const Event& event) { const ChargedFinalState& charged = applyProjection(event, "CFS"); Particles chargedParticles = charged.particlesByPt(); const double c = centrality(event, "ImpactParameterMethod"); cout << "Centrality: " << c << endl; if ((c < 0.) || (c > 80.)) vetoEvent; _histNchVsCentr->fill(c, charged.particles().size() * event.weight()); } void finalize() { // Finishing the centrality calculations //const Centrality& centralityProjection = getProjection("CENTR"); //centralityProjection.finalize(); // Right scaling of the histograms with their individual weights. for (size_t ii = 0; ii < _histNchVsCentr->numBins(); ii++) { double scale = _histNchVsCentr->bin(ii).xWidth() / _histNchVsCentr->bin(ii).numEntries(); _histNchVsCentr->bin(ii).scaleW( scale ); } } // @note First time a post() method is implemented and used by Rivet. This is only available on the // mcplots-alice-dev.cern.ch and mcplots-alice-dev2.cern.ch machine because the Rivet installation here // is the only one which implements Rivet::Analysis::post() and Rivet::AnalysisHandler::post() void post() { } private: Histo1DPtr _histNchVsCentr; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2010_I880049); } diff --git a/analyses/pluginALICE/ALICE_2012_I1127497.cc b/analyses/pluginALICE/ALICE_2012_I1127497.cc --- a/analyses/pluginALICE/ALICE_2012_I1127497.cc +++ b/analyses/pluginALICE/ALICE_2012_I1127497.cc @@ -1,176 +1,191 @@ // -*- C++ -*- #include "pluginALICE/HeavyIonAnalysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include #define NHISTOS 15 #define _USE_MATH_DEFINES #include // maybe use boost namespace Rivet { class ALICE_2012_I1127497 : public HeavyIonAnalysis { public: ALICE_2012_I1127497() : HeavyIonAnalysis("ALICE_2012_I1127497") { } void init() { - - // charged final states - const ChargedFinalState cfs(-0.8, 0.8, 150*MeV); + HeavyIonAnalysis::init(); + + // Set centrality method + addCentralityMethod(HeavyIonAnalysis::ImpactParameter, 10000, "method1"); + + // Charged final states with |eta| < 0.5 and pT > 150 MeV + const Cut& cut = Cuts::abseta < 0.5 && Cuts::pT > 150*MeV; + const ChargedFinalState cfs(cut); addProjection(cfs,"CFS"); - + // @debug Member "log" is declared below. //log.open("/data/mcplots/mcplots/scripts/mcprod/testarea/logs/log.log"); // @note and @todo In principle, only ONE pp histogram needed since centrality does not make a difference here. // However, in some cases it had to be re-binned since some of the binnings of histograms in this analysis // differ from each other. This is therefore the "brute-force" but easy-to-implement way. Making this more // efficient this could be part of the optimization procedure in the end. Now, initialize AA and pp dummy // histograms. for( size_t d = 0; d < NHISTOS; ++d ) { m_sumOfWeights[1][d] = 0; _histNch[1][d] = bookHisto1D(d + 1, 1, 1); _histRAA[d] = bookScatter2D(d+16, 1, 1); m_sumOfWeights[0][d] = 0; // @note Give a proper name for the dummy histograms. They have to incorporate at least the name of the // analysis ALICE_2012_I1127497 and must of course not have the same path as another analysis object used // in this analysis. std::string name = _histNch[1][d]->name(); //cout << "PbPb name: " << _histNch[1][d]->name() << endl; std::string namePP = name + "-pp"; _histNch[0][d] = bookHisto1D( namePP, refData(d + 1, 1, 1) ); // binning is also taken from ref data //cout << "pp name: " << _histNch[0][d]->name() << endl; } // Centrality regions, 2 following numbers are boundaries for a certain region. Note, that some // regions overlap with other regions. The current implementation of centrality bins may not be // the most efficient and flexible one. Take this into account during the optimization procedure. m_centrRegions.clear(); m_centrRegions += {{0., 0.05, 0.05, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5, 0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0., 0.1, 0., 0.2, 0.2, 0.4, 0.4, 0.6, 0.4, 0.8, 0.6, 0.8}}; } void analyze(const Event& event) { const double weight = event.weight(); // final state particles with at least pT = 50 MeV in eta range of |eta| < 0.8 const ChargedFinalState& charged = applyProjection(event, "CFS"); Particles chargedParticles = charged.particlesByPt(); // @note Choose this convention here since JEWEL return a centrality value of -1 for a vacuum run. - float centr = -1; + //float centr = -1; // @note It has to be checked explicitly, if the pointer is sane and is no null_ptr. The standard // constructor of HepMC does not initialize the pointer and returns a null_ptr in case no heavy-ion // was written to the HepMC file. This is e.g. the case for the current pp generators. Trying to // access this naivly may therefore lead to a segmentation fault. - if( event.genEvent()->heavy_ion() ) { + //if( event.genEvent()->heavy_ion() ) { // JEWEL saves the centrality as the impact parameter. This should be changed later in JEWEL directly. // Note, that changes of HepMC may be required as well. However, given the current implementation, // the centrality can be accessed via - centr = event.genEvent()->heavy_ion()->impact_parameter(); + //centr = event.genEvent()->heavy_ion()->impact_parameter(); + //} + + // Check event type + double centr = -1.; + if (HeavyIonAnalysis::is_heavy_ion(event)) { + centr = centrality(event, "method1") / 100.0; } + std::cout << "Centrality: " << centr << std::endl; // Veto event for too large centralities since those are not used in the analysis at all. if( centr > 0.8 ) vetoEvent; size_t pp_AA; float pT; vector< size_t > indices; indices.clear(); if( centr < 0. ) { // Use this as a flag to decide later which histograms should be filled. pp_AA = 0 corresponds // to pp beam (and in case of JEWEL to a vacuum run). + std::cout << "Filling pp indices" << std::endl; pp_AA = 0; for(size_t i = 0; i < NHISTOS; ++i) // Push all indices in case of pp dummy histograms. indices.push_back(i); } else { + std::cout << "Filling PbPb indices" << std::endl; pp_AA = 1; for( size_t i = 0; i < NHISTOS; ++i ) { // Check centrality bins and push corresponding indices of AA histograms. if( inRange( centr, m_centrRegions[ 2 * i ], m_centrRegions[ 2 * i + 1] ) ) { indices.push_back(i); } } } // Fill the right histograms and add weights based on the flag pp_AA and indices pushed to "indices". for( size_t i = 0; i < indices.size(); ++i ) { m_sumOfWeights[pp_AA][indices.at(i)] += weight; foreach (const Particle& p, chargedParticles) { pT = p.pT()/GeV; if(pT < 50.) { + std::cout << "Filling histogram with pT=" << pT << " and weight=" << weight * ( 1 / pT ) << std::endl; _histNch[pp_AA][indices.at(i)]->fill( pT , weight * ( 1 / pT ) ); //cout << "Weight: " << weight << endl; //cout << "pT: " << pT << endl; } } } } void finalize() { // Right scaling of the histograms with their individual weights. for( size_t pp_AA = 0; pp_AA < 2; ++pp_AA ) { for( size_t i = 0; i < NHISTOS; ++i ) { if( m_sumOfWeights[pp_AA][i] > 0 ) //@note is this necessary or can this be handeled by the scaling method? scale( _histNch[pp_AA][i], ( 1./m_sumOfWeights[pp_AA][i] / 2. / M_PI / 1.6 ) ); // TODO PI } } } // @note First time a post() method is implemented and used by Rivet. This is only available on the // mcplots-alice-dev.cern.ch machine because the Rivet installation here is the only one which implements // Rivet::Analysis::post() and Rivet::AnalysisHandler::post() void post() { // Divide AA and pp histograms to obtain R_AA histograms. for( size_t i = 0; i < NHISTOS; ++i) { //cout << "_histNch[1][" << i << "]->bin(0).numEntries(): " << _histNch[1][i]->bin(0).numEntries() << endl; //cout << "_histNch[0][" << i << "]->bin(0).numEntries(): " << _histNch[0][i]->bin(0).numEntries() << endl; //cout << "_histRAA[" << i << "]->bin(0): " << _histRAA[i]->bin(0) << endl << endl; divide( _histNch[1][i], _histNch[0][i], _histRAA[i] ); //cout << "_histRAA[" << i << "]->bin(0): " << _histRAA[i]->bin(0) << endl << endl; } // @debug //log.close(); } private: //std::ofstream validationFile; Histo1DPtr _histNch[2][NHISTOS]; double m_sumOfWeights[2][NHISTOS]; Scatter2DPtr _histRAA[NHISTOS]; std::vector m_centrRegions; // @debug //std::ofstream log; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I1127497); }