// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Tools/ParticleIdUtils.hh" namespace Rivet { class ALICE_2012_I1181770 : public Analysis { public: ALICE_2012_I1181770() : Analysis("ALICE_2012_I1181770") { } public: void init() { // this will return all final state particles addProjection(FinalState(),"FS"); // create histograms 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); } } // Event loop void analyze(const Event& event) { const double weight = event.weight(); // now get the final state collection const FinalState& fs = applyProjection(event, "FS"); // 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 LRG(0.0), etapre(0.0), gapbwd(0.0), gapfwd(0.0); //gap(0.0); double Eslowest(0.0), pslowest(0.0), etaslowest(0.0); double Efastest(0.0), pfastest(0.0), etafastest(0.0); int pidfirst(0); unsigned int num_p(0); FourMomentum leadP(0.,0.,0.,0.); double Elead(0.0), plead(0.0); bool foundslowest = false; bool foundfastest = false; // Particle loop foreach(const Particle& p, fs.particlesByEta()) { // sorted from minus to plus const PdgId pid = p.pdgId(); double eta = p.momentum().eta(); num_p += 1; //SD case if (num_p==1 && pid==2212) { // if it is the first particle(closest to the minus edge) and if it is proton Eslowest = p.momentum().E(); pslowest = p.momentum().vector3().mod(); etaslowest = p.momentum().eta(); foundslowest = true; } if (num_p==fs.size() && pid==2212) { // if it is the last particle(closest to the plus edge) and if it is proton Efastest = p.momentum().E(); pfastest = p.momentum().vector3().mod(); etafastest = p.momentum().eta(); foundfastest = true; } //DD case //DD is only for charged particles. we simply remove neutral particles in the final state projection. //if (PID::threeCharge(p.pdgId()) == 0) continue; // gap calculation for DD events. For each event we find the largest gap (LRG). Also, if the first (last) particle is proton , we calculate the // backward gap (forward gap). if (num_p==1) { etapre = p.momentum().eta(); pidfirst = pid; } else if (num_p > 1) { if (num_p==2 && pidfirst==2212) gapbwd = fabs(eta-etapre); double gap = fabs(eta-etapre); LRG = (gap > LRG ? gap : LRG); // largest gap if (num_p==fs.size() && pid==2212) gapfwd = fabs(eta-etapre); etapre = eta; } } // end of particle loop // Now we are going to calculate Mx event by event if ( foundslowest && !foundfastest ) { // only slowest proton Elead = Eslowest; plead = pslowest; } else if ( !foundslowest && foundfastest ) { // only fastest proton Elead = Efastest; plead = pfastest; } else if ( foundslowest && foundfastest ) { // both slowest and fastest proton if ( abs(etaslowest) > abs(etafastest) ) { Elead = Eslowest; plead = pslowest; } else if (abs(etaslowest) < abs(etafastest) ) { Elead = Efastest; plead = pfastest; } else { // if they are equal generate random numbers.... Elead = Eslowest; plead = pslowest; if ( (double)rand() / (double)RAND_MAX > 0.5) { // generate random numbers in [0.,1.] range and make decision randomly Elead = Efastest; plead = pfastest; } } } //double Mx = -999.; // we don't need this line, please see the explanation below //if (Elead > 0) // we don't need this line, please see the explanation below // The initlation values of Elead and plead is 0. If we don't find any slowest / fastest proton in the event, Mx will be equal to 7000. and // SD plots will not be filled. double Mx = sqrt((sqrtS()/GeV-Elead-plead)*(sqrtS()/GeV-Elead+plead)); bool singleDiff = false; //if (Mx < 200. && Mx > 0) { // Fill SD plot 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); } } if ( singleDiff ) vetoEvent; // DD events are defined as NSD with large gap. vetoEvent means "return". // also remove SD-like events in NSD events if ( std::abs(gapbwd-LRG) < std::numeric_limits::epsilon() || std::abs(gapfwd-LRG) < std::numeric_limits::epsilon() ) 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); } } } void finalize() { // get the ratio plots AIDA::IDataPointSet* h = 0; const string dir = histoDir(); if (fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) { h = histogramFactory().divide( dir + "/d01-x01-y01", *_h_xsec_sd , *_h_xsec_inel ); h = histogramFactory().divide( dir + "/d02-x01-y01", *_h_xsec_dd , *_h_xsec_inel ); } else if (fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) { h = histogramFactory().divide( dir + "/d01-x01-y02", *_h_xsec_sd , *_h_xsec_inel ); h = histogramFactory().divide( dir + "/d02-x01-y02", *_h_xsec_dd , *_h_xsec_inel ); } else if (fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) { h = histogramFactory().divide( dir + "/d01-x01-y03", *_h_xsec_sd , *_h_xsec_inel ); h = histogramFactory().divide( dir + "/d02-x01-y03", *_h_xsec_dd , *_h_xsec_inel ); } // These are original plots 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; }; DECLARE_RIVET_PLUGIN(ALICE_2012_I1181770); }