FillHistograms.h

//__________________________________________________________________________________________________________________________
void FillSimHistograms(const evt::sim::VertexTrack& simTrack,
                              const std::vector<int>& vParticleETypes,
                              const std::vector<double>& vParticleMass,
                              const TVector3& boost_cms,
                              std::vector<TH2D*> vHist2D,
                              std::vector<double>& vParticleCounters,
                              const double& weight = 1)
{
    const double px = simTrack.GetMomentum().GetX();
    const double py = simTrack.GetMomentum().GetY();
    const double pz = simTrack.GetMomentum().GetZ();
    const double p  = simTrack.GetMomentum().GetMag();
    
    for (int ip = 0; ip < int(vParticleETypes.size()); ip++)
    {
        if (simTrack.GetParticleId() == vParticleETypes[ip])                        //eDeuteron, eProton etc; from utl::ParticleConst 
        {
            TLorentzVector tv;   //four momentum of particle in LAB frame of colliding p+p
            tv.SetPxPyPzE(px, py, pz, sqrt(p*p + vParticleMass[ip]*vParticleMass[ip]));
            
            tv.Boost(-boost_cms);
            vHist2D[ip]->Fill(tv.Rapidity(), tv.Pt(), weight);
            
            
            if (simTrack.GetParticleId() == ParticleConst::eDeuteron)
            {
                if (tv.Rapidity() < -0.4 && tv.Rapidity() > -1.2 && tv.Pt() < 0.8)
                    vParticleCounters[ip] = vParticleCounters[ip] + weight;
            }
            else if (simTrack.GetParticleId() == ParticleConst::eProton)
            {
                if(tv.Rapidity() < 0.2 && tv.Rapidity() > -0.6 && tv.Pt() < 1.0)
                    vParticleCounters[ip] = vParticleCounters[ip] + weight;
            }
            else
                vParticleCounters[ip] = vParticleCounters[ip] + weight;
            
            
            break;
        }
    }
}




//__________________________________________________________________________________________________________________________
void FillRecHistograms(const evt::sim::VertexTrack& simTrack,
                       const evt::rec::VertexTrack& recTrack,
                       const std::vector<int>& vParticleETypes,
                       const std::vector<double>& vParticleMass,
                       const TVector3& boost_cms,
                       std::vector<TH2D*> vHist2D,
                       std::vector<double>& vParticleCounters,
                       const double& weight = 1)
{
    const double px = recTrack.GetMomentum().GetX();
    const double py = recTrack.GetMomentum().GetY();
    const double pz = recTrack.GetMomentum().GetZ();
    const double p  = recTrack.GetMomentum().GetMag();
    
    for (int ip = 0; ip < int(vParticleETypes.size()); ip++)
    {
        if (simTrack.GetParticleId() == vParticleETypes[ip])                        //eDeuteron, eProton etc; from utl::ParticleConst 
        {
            TLorentzVector tv;   //four momentum of particle in LAB frame of colliding p+p
            tv.SetPxPyPzE(px, py, pz, sqrt(p*p + vParticleMass[ip]*vParticleMass[ip]));
            
            tv.Boost(-boost_cms);
            vHist2D[ip]->Fill(tv.Rapidity(), tv.Pt(), weight);
            
            
            if (simTrack.GetParticleId() == ParticleConst::eDeuteron)
            {
                if (tv.Rapidity() < -0.4 && tv.Rapidity() > -1.2 && tv.Pt() < 0.8)
                    vParticleCounters[ip] = vParticleCounters[ip] + weight;
            }
            
            else if (simTrack.GetParticleId() == ParticleConst::eProton)
            {
                if(tv.Rapidity() < 0.2 && tv.Rapidity() > -0.6 && tv.Pt() < 1.0)
                    vParticleCounters[ip] = vParticleCounters[ip] + weight;
            }
            
            else
                vParticleCounters[ip] = vParticleCounters[ip] + weight;
            
            
            break;
        }
    }
}