EventShape.h

Go to the documentation of this file.

#include "TLorentzVector.h"
#include "TVector3.h"
#include "SimpleCand.C"
#include "TMatrixD.h"
#include "TMatrixDEigen.h"

#define FWMAX 6 // maximum Fox Wolfram moment

class EventShape {
public:
        EventShape(const CandList &l, TLorentzVector cms);      
        ~EventShape(){};

        // ******* multiplicities
        int NParticles() const {return fN;}                     // number of particle candidates
        int NCharged() const {return fnChrg;}           // number of charged candidates
        int NNeutral() const  {return fnNeut;}          // number of neutral candidates
        
        // ******* maxima of momenta
        double PmaxLab() const  {return fpmaxlab;}              // max momentum in lab system
        double PmaxCms() const {return fpmaxcms;}               // max momentum im cms system
        double PminLab() const  {return fpminlab;}              // min momentum in lab system
        double PminCms() const {return fpmincms;}               // min momentum im cms system
        double Ptmax() const {return fptmax;}                   // max pt (same for lab and cms)
        double Ptmin() const {return fptmin;}                   // min pt (same for lab and cms)
        double PRapmax() const {return fprapmax;}               // max pseudorapidity (lab)
        
        // ******* sum of energies/momenta (lab system)
        double PtSumLab() const {return fptsumlab;}                             // sum of pt in (lab)
        double NeutEtSumLab() const {return fneutetsumlab;}             // sum of transvers energys of neutrals (lab)
        double NeutESumLab()  const {return fneutesumlab;}              // sum of energys of neutrals (lab)
        double ChrgPtSumLab() const {return fchrgptsumlab;}             // sum of pt of charged (lab)
        double ChrgPSumLab() const {return fchrgpsumlab;}               // sum of momenta of charged (lab)
        
        // ******* sum of energies/momenta (cms system)
        double PtSumCms() const {return fptsumcms;}                             // sum of pt in (cms)
        double NeutEtSumCms() const {return fneutetsumcms;}             // sum of transvers energys of neutrals (cms)
        double NeutESumCms()  const {return fneutesumcms;}              // sum of energys of neutrals (cms)
        double ChrgPtSumCms() const {return fchrgptsumcms;}             // sum of pt of charged (cms)
        double ChrgPSumCms() const {return fchrgpsumcms;}               // sum of momenta of charged (cms)
        
        // ******* multiplicities with threshold
        int MultPminLab(double pmin);   // number of particles with p>pmin (lab frame)
        int MultPmaxLab(double pmax);   // number of particles with p<pmax (lab frame)
        int MultPminCms(double pmin);   // number of particles with p>pmin (cms frame)
        int MultPmaxCms(double pmax);   // number of particles with p<pmax (cms frame)
        
        // ******* multiplicities with threshold
        int MultPtminLab(double ptmin); // number of particles with pt>pmin (lab frame)
        int MultPtmaxLab(double ptmax); // number of particles with pt<pmax (lab frame)
        int MultPtminCms(double ptmin); // number of particles with pt>pmin (cms frame)
        int MultPtmaxCms(double ptmax); // number of particles with pt<pmax (cms frame)
        
        // ******* neutrals multiplicities with threshold
        int MultNeutEminLab(double emin);       // number of neutrals with E>emin (lab frame)
        int MultNeutEmaxLab(double emax);       // number of neutrals with E<emax (lab frame)
        int MultNeutEminCms(double emin);       // number of neutrals with E>emin (cms frame)
        int MultNeutEmaxCms(double emax);       // number of neutrals with E<emax (cms frame)
        
        // ******* charged multiplicities with threshold
        int MultChrgPminLab(double pmin);       // number of charged with p>pmin (lab frame)
        int MultChrgPmaxLab(double pmax);       // number of charged with p<pmax (lab frame)
        int MultChrgPminCms(double pmin);       // number of charged with p>pmin (cms frame)
        int MultChrgPmaxCms(double pmax);       // number of charged with p<pmax (cms frame)

        // ******* sums with threshold
        double SumPminLab(double pmin);         // sum of momenta with p>pmin (lab frame)
        double SumPmaxLab(double pmax);         // sum of momenta with p<pmax (lab frame)
        double SumPminCms(double pmin);         // sum of momenta with p>pmin (cms frame)
        double SumPmaxCms(double pmax);         // sum of momenta with p<pmax (cms frame)
        
        // ******* sums with threshold
        double SumPtminLab(double ptmin);               // sum of pt with pt>pmin (lab frame)
        double SumPtmaxLab(double ptmax);               // sum of pt with pt<pmax (lab frame)
        double SumPtminCms(double ptmin);               // sum of pt with pt>pmin (cms frame)
        double SumPtmaxCms(double ptmax);               // sum of pt with pt<pmax (cms frame)
        
        // ******* neutrals sums with threshold
        double SumNeutEminLab(double emin);     // sum of energies of neutrals with E>emin (lab frame)
        double SumNeutEmaxLab(double emax);     // sum of energies of neutrals with E<emax (lab frame)
        double SumNeutEminCms(double emin);     // sum of energies of neutrals with E>emin (cms frame)
        double SumNeutEmaxCms(double emax);     // sum of energies of neutrals with E<emax (cms frame)
        
        // ******* charged sums with threshold
        double SumChrgPminLab(double pmin);     // sum of momenta of charged with p>pmin (lab frame)
        double SumChrgPmaxLab(double pmax);     // sum of momenta of charged with p<pmax (lab frame)
        double SumChrgPminCms(double pmin);     // sum of momenta of charged with p>pmin (cms frame)
        double SumChrgPmaxCms(double pmax);     // sum of momenta of charged with p<pmax (cms frame)

        // ******* shape variables
        double Sphericity();                    // Sphericity
        double Aplanarity();                    // Aplanarity
        double Planarity();                             // Planarity
        double Circularity();                   // Cirularity
        
        double FoxWolfMomH(int order);  // Fox Wolfram moment absolute H_i
        double FoxWolfMomR(int order);  // Fox Wolfram moment relative R_i = H_i/H_0
        
        double Thrust();
        TVector3 ThrustVector();
        
private:
        
        void ComputeSphericity();                               // compute sph, apl, pla
        double Eps(const TVector3 v1, const TVector3 v2) {return (v1*v2)>0. ? 1. : -1.;}  // aux for Thrust
        double Legendre( int l, double x );             // Legendre function; auxilliary for Fox Wolfram moments
        
        std::vector<TLorentzVector> fLabList;   // List of 4-vectors in lab frame
        std::vector<TLorentzVector> fCmsList;   // List of 4-vectors in cms frame
        std::vector<int>                        fCharge;        // List of charges of particles
        
        int    fnChrg;                          // number of charged particles
        int    fnNeut;                          // number of neutral particles
        int    fN;                                      // number of particles
        
        double fpmaxlab;                        // maximum momentum lab frame
        double fpmaxcms;                        // maximum momentum cms frame
        double fpminlab;                        // minimum momentum lab frame
        double fpmincms;                        // minimum momentum cms frame
        double fptmax;                          // maximum transvers momentum
        double fptmin;                          // minimum transvers momentum
        double fprapmax;                        // maximum pseudorapidity

        double fptsumlab;                       // sum of pt in (lab)
        double fneutetsumlab;           // sum of transvers energys of neutrals (lab)
        double fneutesumlab;            // sum of energys of neutrals (lab)
        double fchrgptsumlab;           // sum of pt of charged (lab)
        double fchrgpsumlab;            // sum of momenta of charged (lab)
        
        double fptsumcms;                       // sum of pt in (lab)
        double fneutetsumcms;           // sum of transvers energys of neutrals (lab)
        double fneutesumcms;            // sum of energys of neutrals (lab)
        double fchrgptsumcms;           // sum of pt of charged (lab)
        double fchrgpsumcms;            // sum of momenta of charged (lab)
        
        double fsph;                            // Sphericity
        double fapl;                            // Aplanarity
        double fpla;                            // Planarity
        double fcir;                            // Circularity
  
        double fFWmom[FWMAX];           // Fox Wolfram moments up to FWMAX
        bool   fFWready;                        // did we compute FW moments?
  
        double fthr;                            // thrust
        TVector3 fThrVect;                      // direction of thrust
        TVector3 fBoost;                        // boost vector to go to requested frame

};

// ************************************ IMPLEMENTATION

EventShape::EventShape(const CandList &l, TLorentzVector cms) : 
  fnChrg(0), fnNeut(0), fN(0), fpmaxlab(0.), fpmaxcms(0.),fpminlab(0.), fpmincms(0.), fptmax(0.), fptmin(0.),
  fptsumlab(0.),fneutetsumlab(0.),fneutesumlab(0.),fchrgptsumlab(0.),fchrgpsumlab(0.),
  fptsumcms(0.),fneutetsumcms(0.),fneutesumcms(0.),fchrgptsumcms(0.),fchrgpsumcms(0.),
  fsph(-1.), fapl(-1.), fpla(-1.), fcir(-1.), fFWready(false), fthr(-1.)
{
  int i;

  // initialize more complex things
  fThrVect.SetXYZ(0.,0.,0.);
  fBoost=cms.BoostVector();  
  for (i=0;i<=FWMAX;++i) fFWmom[i]=0.;
  
  
  double pmax=0., ptmax=0., pmaxcms=0.;
  double pmin=1000., ptmin=1000., pmincms=1000.;
  double prapmax=0.;
  
  fLabList.clear();
  fCmsList.clear();
  
  for (i=0;i<l.size();++i)
  {
        TLorentzVector lv(l[i].P4());
        int chrg(l[i].Charge());

        fN++;
        
        // cache multiplicities
        if (chrg==0) fnNeut++;
        else fnChrg++;
        
        // cache unboosted 4-vectors
        fLabList.push_back(lv);
        // cache charges
        fCharge.push_back(chrg);
        
        // sum momentum variables (lab)
        fptsumlab += lv.Pt();
        if (chrg==0)
        {
          fneutetsumlab += lv.Pt();  
          fneutesumlab  += lv.E();
        }
        else
        {
          fchrgptsumlab += lv.Pt();
          fchrgpsumlab  += lv.P();
        }
        
        // cache maximum momenta in lab
        if (lv.P()>pmax) pmax=lv.P();
        if (lv.Pt()>ptmax) ptmax=lv.Pt();
        
        // cache minimum momenta in lab
        if (lv.P()<pmin) pmin=lv.P();
        if (lv.Pt()<ptmin) ptmin=lv.Pt();
        
        // cache pseudorapidity
        if (lv.PseudoRapidity()>prapmax) prapmax = lv.PseudoRapidity();
        
        // cache boosted vectors
        lv.Boost(-fBoost);
        fCmsList.push_back(lv);

        // sum momentum variables (cms)
        fptsumcms += lv.Pt();
        if (chrg==0)
        {
          fneutetsumcms += lv.Pt();  
          fneutesumcms  += lv.E();
        }
        else
        {
          fchrgptsumcms += lv.Pt();
          fchrgpsumcms  += lv.P();
        }
        
        // cache maximum momenta in cms
        if (lv.P()>pmaxcms) pmaxcms=lv.P();
        // cache minimum momenta in cms
        if (lv.P()<pmincms) pmincms=lv.P();
  }
  
  fpmaxlab = pmax;
  fptmax   = ptmax;
  fpmaxcms = pmaxcms;
  fpminlab = pmin;
  fptmin   = ptmin;
  fpmincms = pmincms;
  fprapmax = prapmax;
}

// ----------------------------

void EventShape::ComputeSphericity()
{
  if( fN==0 ) return;
  
  double stot=0, sxx=0, sxy=0, sxz=0, syy=0, syz=0, szz=0;      
  int i;
  
  for (i=0;i<fN;++i)
  {
        TVector3 v(fCmsList[i].Vect());
        
        sxx += v.X()*v.X(); sxy += v.X()*v.Y(); sxz += v.X()*v.Z();
        syy += v.Y()*v.Y(); syz += v.Y()*v.Z();
        szz += v.Z()*v.Z();
        
        stot += v.Mag2();
  }
  
  TMatrixD sm(3,3);
  sm(0,0) = sxx/stot; sm(0,1) = sxy/stot; sm(0,2) = sxz/stot; 
  sm(1,0) = sxy/stot; sm(1,1) = syy/stot; sm(1,2) = syz/stot; 
  sm(2,0) = sxz/stot; sm(2,1) = syz/stot; sm(2,2) = szz/stot; 
  
  TMatrixDEigen ei(sm);
  TMatrixD eiv=ei.GetEigenValues();
  
  fsph = 1.5 * (eiv(1,1) + eiv(2,2));           // 3/2 (lam_2+lam_3)
  fapl = 1.5 * eiv(2,2);                                        // 3/2 lam_3
  fpla = eiv(1,1) - eiv(2,2);                           // lam_2-lam_3
  fcir = 2.*eiv(1,1)/(eiv(1,1)+eiv(0,0));       // 2* min(lam_1, lam_2)/(lam_1+lam_2)
}

// ---------------------------------

double EventShape::Sphericity()
{
  if (fsph<0.) ComputeSphericity();
  
  return fsph;
}

// ---------------------------------

double EventShape::Planarity()
{
  if (fpla<0.) ComputeSphericity();
  
  return fpla;
}

// ---------------------------------

double EventShape::Aplanarity()
{
  if (fapl<0.) ComputeSphericity();
  
  return fapl;
}

// ---------------------------------

double EventShape::Circularity()
{
  if (fcir<0.) ComputeSphericity();
  
  return fcir;
}

// ---------------------------------------
// use iterative formula for thrust vector
//
//             sum eps( n(j)*p_i) * p_i
//   n(j+1) = --------------------------
//            |        ''             | 
//
// with eps(x)=-1 for x<0 and =+1 for x>0

double EventShape::Thrust()
{
  // did we already compute?
  if (fthr>-1.) return fthr; 
  
  TVector3 n0(0,0,0);
  if( fN==0 ) return -1.;
  
  int i,j;
  double pmax=0;
  
  //find starting vector as maximum momentum vector
  for (i=0;i<fN;++i)
  {
        if (fCmsList[i].Vect().Mag()>pmax)
        {
          n0=fCmsList[i].Vect();
          pmax=fCmsList[i].Vect().Mag();
        }
  }
  
  TVector3 nnew(0,0,0);

  // find thrust axis (10 iterations)
  for (i=0;i<10;++i)
  {
        for (j=0;j<fN;++j) 
          nnew += Eps(n0,fCmsList[j].Vect()) * fCmsList[j].Vect();
        
        n0=nnew.Unit();
  }
  
  fThrVect = n0;
  
  double thr=0, sum=0;
  for (i=0;i<fN;++i)
  {
        thr += fabs(fThrVect.Dot(fCmsList[i].Vect()));
        sum += fCmsList[i].Vect().Mag();
  }
  
  fthr = thr/sum;
  
  return fthr;
}

// ---------------------------------------

TVector3 EventShape::ThrustVector()
{
  if (fthr<0.) Thrust();
  
  return fThrVect;
}

// ---------------------------------------

double EventShape::Legendre( int l, double x )
{
    if (fabs(x)>1.) return -999.;

        if (l==0) return 1.;
        
    double pmm = 1.;
        double pmmp1 = x;

        if (l>1)
        {
                for(int ll=2; ll<=l; ll++)
                {
                  double pll = (x * (2 * ll - 1) * pmmp1 - (ll - 1) * pmm) / (ll);
                  pmm = pmmp1;
                  pmmp1 = pll;
                }
        }
   return pmmp1;
}

// ---------------------------------------

double EventShape::FoxWolfMomH(int order)
{
  if (order>FWMAX) return -1.;
 
  // already computed FW moments
  if (fFWready) return fFWmom[order];
  
  if( fN==0 ) return -1.;

  double s = 0.;
  int i,j,l;

  for (i=0; i<fN-1; ++i)
  {
     // this candidate's 3-momentum
          TVector3 p1(fCmsList[i].Vect());
      double pmag1 = p1.Mag();

      // loop over other candidates, starting at the next one in the list
          
          for (j=i+1; j<fN; ++j)
          {
                // this candidate's 3-momentum
                TVector3 p2(fCmsList[j].Vect());
                double pmag2 = p2.Mag();

                // the cosine of the angle between the two candidates
                double cosPhi =  cos ( p1.Angle(p2) );

                // the contribution of this pair of track 
                // (note the factor 2 : the pair enters the sum twice)
                for( l=0; l<=FWMAX; l++ )
                  fFWmom[l] += 2 * pmag1 * pmag2 * Legendre( l, cosPhi );
                  
          }
      // contribution of this track
      for( l=0; l<=FWMAX; l++ )
                fFWmom[l] += pmag1 * pmag1 * Legendre( l, 1. );
      
      // total energy
      s += fCmsList[i].Energy();
   }
  
  // well ...
  if( s<=0. ) return -1.;
  double s2=s*s;
  
  // normalize Fox Wolfram Moments
  for( i=0; i<=FWMAX; i++) fFWmom[i]/=s2 ;  
  
  fFWready = true;
  
  return fFWmom[order];
}


// ---------------------------------------

double EventShape::FoxWolfMomR(int order)
{
  return FoxWolfMomH(order)/FoxWolfMomH(0);
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

int EventShape::MultPminLab(double pmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].P()>pmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPmaxLab(double pmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].P()<pmax) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPminCms(double pmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].P()>pmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPmaxCms(double pmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].P()<pmax) cnt++;
 
  return cnt;
}
// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

int EventShape::MultPtminLab(double ptmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].Pt()>ptmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPtmaxLab(double ptmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].Pt()<ptmax) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPtminCms(double ptmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].Pt()>ptmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultPtmaxCms(double ptmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].Pt()<ptmax) cnt++;
 
  return cnt;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

int EventShape::MultNeutEminLab(double emin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fLabList[i].E()>emin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultNeutEmaxLab(double emax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fLabList[i].E()<emax) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultNeutEminCms(double emin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fCmsList[i].E()>emin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultNeutEmaxCms(double emax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fCmsList[i].E()<emax) cnt++;
 
  return cnt;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

int EventShape::MultChrgPminLab(double pmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fLabList[i].P()>pmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultChrgPmaxLab(double pmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fLabList[i].P()<pmax) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultChrgPminCms(double pmin)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fCmsList[i].P()>pmin) cnt++;
 
  return cnt;
}

// ---------------------------------------

int EventShape::MultChrgPmaxCms(double pmax)
{
  int cnt=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fCmsList[i].P()<pmax) cnt++;
 
  return cnt;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

double EventShape::SumPminLab(double pmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].P()>pmin) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPmaxLab(double pmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].P()<pmax) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPminCms(double pmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].P()>pmin) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPmaxCms(double pmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].P()<pmax) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

double EventShape::SumPtminLab(double ptmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].Pt()>ptmin) sum+=fLabList[i].Pt();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPtmaxLab(double ptmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fLabList[i].Pt()<ptmax) sum+=fLabList[i].Pt();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPtminCms(double ptmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].Pt()>ptmin) sum+=fLabList[i].Pt();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumPtmaxCms(double ptmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCmsList[i].Pt()<ptmax) sum+=fLabList[i].Pt();
 
  return sum;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

double EventShape::SumNeutEminLab(double emin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fLabList[i].E()>emin) sum+=fLabList[i].E();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumNeutEmaxLab(double emax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fLabList[i].E()<emax) sum+=fLabList[i].E();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumNeutEminCms(double emin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fCmsList[i].E()>emin) sum+=fLabList[i].E();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumNeutEmaxCms(double emax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]==0 && fCmsList[i].E()<emax) sum+=fLabList[i].E();
 
  return sum;
}

// ---------------------------------------
// ---------------------------------------
// ---------------------------------------

double EventShape::SumChrgPminLab(double pmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fLabList[i].P()>pmin) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumChrgPmaxLab(double pmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fLabList[i].P()<pmax) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumChrgPminCms(double pmin)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fCmsList[i].P()>pmin) sum+=fLabList[i].P();
 
  return sum;
}

// ---------------------------------------

double EventShape::SumChrgPmaxCms(double pmax)
{
  double sum=0;
  for (int i=0;i<fN;++i) 
        if (fCharge[i]!=0 && fCmsList[i].P()<pmax) sum+=fLabList[i].P();
 
  return sum;
}