// @(#)root/quadp:$Name:  $:$Id: TQpVar.cxx,v 1.5 2005/09/04 10:02:30 brun Exp $
// Author: Eddy Offermann   May 2004

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

/*************************************************************************
 * Parts of this file are copied from the OOQP distribution and          *
 * are subject to the following license:                                 *
 *                                                                       *
 * COPYRIGHT 2001 UNIVERSITY OF CHICAGO                                  *
 *                                                                       *
 * The copyright holder hereby grants you royalty-free rights to use,    *
 * reproduce, prepare derivative works, and to redistribute this software*
 * to others, provided that any changes are clearly documented. This     *
 * software was authored by:                                             *
 *                                                                       *
 *   E. MICHAEL GERTZ      gertz@mcs.anl.gov                             *
 *   Mathematics and Computer Science Division                           *
 *   Argonne National Laboratory                                         *
 *   9700 S. Cass Avenue                                                 *
 *   Argonne, IL 60439-4844                                              *
 *                                                                       *
 *   STEPHEN J. WRIGHT     swright@cs.wisc.edu                           *
 *   Computer Sciences Department                                        *
 *   University of Wisconsin                                             *
 *   1210 West Dayton Street                                             *
 *   Madison, WI 53706   FAX: (608)262-9777                              *
 *                                                                       *
 * Any questions or comments may be directed to one of the authors.      *
 *                                                                       *
 * ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES OF   *
 * ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT, AND     *
 * OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A CONTRACT   *
 * WITH THE DEPARTMENT OF ENERGY.                                        *
 *************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Class containing the variables for the general QP formulation         //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include "Riostream.h"
#include "TQpVar.h"
#include "TMatrixD.h"

ClassImp(TQpVar)

//______________________________________________________________________________
 TQpVar::TQpVar()
{
  fNx   = 0;
  fMy   = 0;
  fMz   = 0;
  fNxup = 0;
  fNxlo = 0;
  fMcup = 0;
  fMclo = 0;
}

//______________________________________________________________________________
 TQpVar::TQpVar(TVectorD &x_in,TVectorD &s_in,TVectorD &y_in,TVectorD &z_in,
               TVectorD &v_in,TVectorD &gamma_in,TVectorD &w_in,TVectorD &phi_in,
               TVectorD &t_in,TVectorD &lambda_in,TVectorD &u_in,TVectorD &pi_in,
               TVectorD &ixlow_in,TVectorD &ixupp_in,TVectorD &iclow_in,TVectorD &icupp_in)
{
  if (x_in     .GetNrows() > 0) fX.       Use(x_in     .GetNrows(),x_in     .GetMatrixArray());
  if (s_in     .GetNrows() > 0) fS.       Use(s_in     .GetNrows(),s_in     .GetMatrixArray());
  if (y_in     .GetNrows() > 0) fY.       Use(y_in     .GetNrows(),y_in     .GetMatrixArray());
  if (z_in     .GetNrows() > 0) fZ.       Use(z_in     .GetNrows(),z_in     .GetMatrixArray());
  if (v_in     .GetNrows() > 0) fV.       Use(v_in     .GetNrows(),v_in     .GetMatrixArray());
  if (phi_in   .GetNrows() > 0) fPhi.     Use(phi_in   .GetNrows(),phi_in   .GetMatrixArray());
  if (w_in     .GetNrows() > 0) fW.       Use(w_in     .GetNrows(),w_in     .GetMatrixArray());
  if (gamma_in .GetNrows() > 0) fGamma.   Use(gamma_in .GetNrows(),gamma_in .GetMatrixArray());
  if (t_in     .GetNrows() > 0) fT.       Use(t_in     .GetNrows(),t_in     .GetMatrixArray());
  if (lambda_in.GetNrows() > 0) fLambda.  Use(lambda_in.GetNrows(),lambda_in.GetMatrixArray());
  if (u_in     .GetNrows() > 0) fU.       Use(u_in     .GetNrows(),u_in     .GetMatrixArray());
  if (pi_in    .GetNrows() > 0) fPi.      Use(pi_in    .GetNrows(),pi_in    .GetMatrixArray());
  if (ixlow_in .GetNrows() > 0) fXloIndex.Use(ixlow_in .GetNrows(),ixlow_in .GetMatrixArray());
  if (ixupp_in .GetNrows() > 0) fXupIndex.Use(ixupp_in .GetNrows(),ixupp_in .GetMatrixArray());
  if (iclow_in .GetNrows() > 0) fCloIndex.Use(iclow_in .GetNrows(),iclow_in .GetMatrixArray());
  if (icupp_in .GetNrows() > 0) fCupIndex.Use(icupp_in .GetNrows(),icupp_in .GetMatrixArray());

  fNx = fX.GetNrows();
  fMy = fY.GetNrows();
  fMz = fZ.GetNrows();

  Assert(fNx == fXloIndex.GetNrows() || 0 == fXloIndex.GetNrows());
  Assert(fNx == fXloIndex.GetNrows() || 0 == fXloIndex.GetNrows());
  Assert(fMz == fCloIndex.GetNrows() || 0 == fCloIndex.GetNrows());
  Assert(fMz == fCupIndex.GetNrows() || 0 == fCupIndex.GetNrows());

  fNxlo = fXloIndex.NonZeros();
  fNxup = fXupIndex.NonZeros();
  fMclo = fCloIndex.NonZeros();
  fMcup = fCupIndex.NonZeros();
  fNComplementaryVariables = fMclo+fMcup+fNxlo+fNxup;

  Assert(fMz == fS.GetNrows());
  Assert(fNx == fV     .GetNrows() || (0 == fV     .GetNrows() && fNxlo == 0));
  Assert(fNx == fGamma .GetNrows() || (0 == fGamma .GetNrows() && fNxlo == 0));

  Assert(fNx == fW     .GetNrows() || (0 == fW     .GetNrows() && fNxup == 0));
  Assert(fNx == fPhi   .GetNrows() || (0 == fPhi   .GetNrows() && fNxup == 0));

  Assert(fMz == fT     .GetNrows() || (0 == fT     .GetNrows() && fMclo == 0));
  Assert(fMz == fLambda.GetNrows() || (0 == fLambda.GetNrows() && fMclo == 0));

  Assert(fMz == fU     .GetNrows() || (0 == fU     .GetNrows() && fMcup == 0));
  Assert(fMz == fPi    .GetNrows() || (0 == fPi    .GetNrows() && fMcup == 0));
}

//______________________________________________________________________________
 TQpVar::TQpVar(Int_t nx,Int_t my,Int_t mz,TVectorD &ixlow,TVectorD &ixupp,
                                          TVectorD &iclow,TVectorD &icupp)
{
  Assert(nx == ixlow.GetNrows() || 0 == ixlow.GetNrows());
  Assert(nx == ixlow.GetNrows() || 0 == ixlow.GetNrows());
  Assert(mz == iclow.GetNrows() || 0 == iclow.GetNrows());
  Assert(mz == icupp.GetNrows() || 0 == icupp.GetNrows());

  fNxlo = ixlow.NonZeros();
  fNxup = ixupp.NonZeros();
  fMclo = iclow.NonZeros();
  fMcup = icupp.NonZeros();

  if (ixlow.GetNrows() > 0) fXloIndex.Use(ixlow.GetNrows(),ixlow.GetMatrixArray()); 
  if (ixupp.GetNrows() > 0) fXupIndex.Use(ixupp.GetNrows(),ixupp.GetMatrixArray());  
  if (iclow.GetNrows() > 0) fCloIndex.Use(iclow.GetNrows(),iclow.GetMatrixArray()); 
  if (icupp.GetNrows() > 0) fCupIndex.Use(icupp.GetNrows(),icupp.GetMatrixArray());
  
  fNx = nx;
  fMy = my;
  fMz = mz;
  
  if (fMclo > 0) {
    fT.ResizeTo(fMz);
    fLambda.ResizeTo(fMz);
  }
  if (fMcup > 0) {
    fU.ResizeTo(fMz);
    fPi.ResizeTo(fMz);
  }
  if (fNxlo > 0) {
    fV.ResizeTo(fNx);
    fGamma.ResizeTo(fNx);
  }

  if (fNxup > 0) {
    fW.ResizeTo(fNx);
    fPhi.ResizeTo(fNx);
  }

  fS.ResizeTo(fMz);
  fX.ResizeTo(fNx);
  fY.ResizeTo(fMy);
  fZ.ResizeTo(fMz);
  fNComplementaryVariables = fMclo+fMcup+fNxlo+fNxup;
}

//______________________________________________________________________________
 TQpVar::TQpVar(const TQpVar &another) : TObject(another)
{
  *this = another;
}

//______________________________________________________________________________
 Double_t TQpVar::GetMu()
{   
  Double_t mu = 0.0;
  if (fNComplementaryVariables > 0 ) {
    if (fMclo > 0) mu += fT*fLambda;
    if (fMcup > 0) mu += fU*fPi;
    if (fNxlo > 0) mu += fV*fGamma;
    if (fNxup > 0) mu += fW*fPhi;
    
    mu /= fNComplementaryVariables;

  }
  return mu;
}

//______________________________________________________________________________
 Double_t TQpVar::MuStep(TQpVar *step,Double_t alpha)
{
  Double_t mu = 0.0;
  if (fNComplementaryVariables > 0) {
    if (fMclo > 0)
      mu += (fT+alpha*step->fT)*(fLambda+alpha*step->fLambda);
    if (fMcup > 0)
      mu += (fU+alpha*step->fU)*(fPi+alpha*step->fPi);
    if (fNxlo > 0)
      mu += (fV+alpha*step->fV)*(fGamma+alpha*step->fGamma);
    if (fNxup > 0)
      mu += (fW+alpha*step->fW)*(fPhi+alpha*step->fPhi);
    mu /= fNComplementaryVariables;
  }
  return mu;
}

//______________________________________________________________________________
 void TQpVar::Saxpy(TQpVar *b,Double_t alpha)
{
  Add(fX,alpha,b->fX);
  Add(fY,alpha,b->fY);
  Add(fZ,alpha,b->fZ);
  Add(fS,alpha,b->fS);
  if (fMclo > 0) {
    Assert((b->fT)     .MatchesNonZeroPattern(fCloIndex) &&
           (b->fLambda).MatchesNonZeroPattern(fCloIndex));

    Add(fT     ,alpha,b->fT);
    Add(fLambda,alpha,b->fLambda);
  }
  if (fMcup > 0) {
    Assert((b->fU) .MatchesNonZeroPattern(fCupIndex) &&
           (b->fPi).MatchesNonZeroPattern(fCupIndex));

    Add(fU ,alpha,b->fU);
    Add(fPi,alpha,b->fPi);
  }
  if (fNxlo > 0) {
    Assert((b->fV)    .MatchesNonZeroPattern(fXloIndex) &&
           (b->fGamma).MatchesNonZeroPattern(fXloIndex));

    Add(fV    ,alpha,b->fV);
    Add(fGamma,alpha,b->fGamma);
  }
  if (fNxup > 0) {
    Assert((b->fW)  .MatchesNonZeroPattern(fXupIndex) &&
           (b->fPhi).MatchesNonZeroPattern(fXupIndex));

    Add(fW  ,alpha,b->fW);
    Add(fPhi,alpha,b->fPhi);
  }
}

//______________________________________________________________________________
 void TQpVar::Negate()
{
  fS *= -1.;
  fX *= -1.;
  fY *= -1.;
  fZ *= -1.;
  if (fMclo > 0) {
    fT      *= -1.;
    fLambda *= -1.;
  }
  if (fMcup > 0) {
    fU  *= -1.;
    fPi *= -1.;
  }
  if (fNxlo > 0) {
    fV     *= -1.;
    fGamma *= -1.;
  }
  if (fNxup > 0) {
    fW   *= -1.;
    fPhi *= -1.;
  }
}

//______________________________________________________________________________
 Double_t TQpVar::StepBound(TQpVar *b)
{
  Double_t maxStep = 1.0;

  if (fMclo > 0 ) {
    Assert(fT     .SomePositive(fCloIndex));
    Assert(fLambda.SomePositive(fCloIndex));

    maxStep = this->StepBound(fT,     b->fT,     maxStep);
    maxStep = this->StepBound(fLambda,b->fLambda,maxStep);
  }

  if (fMcup > 0 ) {
    Assert(fU .SomePositive(fCupIndex));
    Assert(fPi.SomePositive(fCupIndex));

    maxStep = this->StepBound(fU, b->fU, maxStep);
    maxStep = this->StepBound(fPi,b->fPi,maxStep);
  }

  if (fNxlo > 0 ) {
    Assert(fV    .SomePositive(fXloIndex));
    Assert(fGamma.SomePositive(fXloIndex));

    maxStep = this->StepBound(fV,    b->fV,    maxStep);
    maxStep = this->StepBound(fGamma,b->fGamma,maxStep);
  }

  if (fNxup > 0 ) {
    Assert(fW  .SomePositive(fXupIndex));
    Assert(fPhi.SomePositive(fXupIndex));

    maxStep = this->StepBound(fW,  b->fW,  maxStep);
    maxStep = this->StepBound(fPhi,b->fPhi,maxStep);
  }

  return maxStep;
}

//______________________________________________________________________________
 Double_t TQpVar::StepBound(TVectorD &v,TVectorD &dir,Double_t maxStep)
{ 
  if (!AreCompatible(v,dir)) {
    ::Error("StepBound(TVectorD &,TVectorD &,Double_t)","vector's not compatible");
    return kFALSE;
  }

  const Int_t n = v.GetNrows();
  const Double_t * const pD = dir.GetMatrixArray();
  const Double_t * const pV = v.GetMatrixArray();

  Double_t bound = maxStep;
  for (Int_t i = 0; i < n; i++) {
    Double_t tmp = pD[i];
    if( pV[i] >= 0 && tmp < 0 ) {
      tmp = -pV[i]/tmp;
      if (tmp < bound)
        bound = tmp;
    }
  }
  return bound;
}

//______________________________________________________________________________
 Bool_t TQpVar::IsInteriorPoint()
{
  Bool_t interior = kTRUE;
  if (fMclo > 0)
    interior = interior &&
                   fT.SomePositive(fCloIndex) && fLambda.SomePositive(fCloIndex);

  if (fMcup > 0)
    interior = interior &&
                   fU.SomePositive(fCupIndex) && fPi.SomePositive(fCupIndex);

  if (fNxlo > 0)
    interior = interior &&
                   fV.SomePositive(fXloIndex) && fGamma.SomePositive(fXloIndex);

  if (fNxup > 0)
    interior = interior &&
                   fW.SomePositive(fXupIndex) && fPhi.SomePositive(fXupIndex);

  return interior;
}

//______________________________________________________________________________
 Double_t TQpVar::FindBlocking(TQpVar   *step, 
                              Double_t &primalValue,
                              Double_t &primalStep,
                              Double_t &dualValue,
                              Double_t &dualStep, 
                              Int_t    &fIrstOrSecond)
{ 
  fIrstOrSecond = 0;
  Double_t alpha = 1.0;
  if (fMclo > 0)
    alpha = FindBlocking(fT,step->fT,fLambda,step->fLambda,alpha,
                         primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
  
  if (fMcup > 0)
    alpha = FindBlocking(fU,step->fU,fPi,step->fPi,alpha,
                         primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
  
  if (fNxlo > 0)
    alpha = FindBlocking(fV,step->fV,fGamma,step->fGamma,alpha,
                         primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);

  if (fNxup > 0)
    alpha = FindBlocking(fW,step->fW,fPhi,step->fPhi,alpha,
                         primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);

  return alpha;
}

//______________________________________________________________________________
 Double_t TQpVar::FindBlocking(TVectorD &w,TVectorD &wstep,TVectorD &u,TVectorD &ustep,
                              Double_t maxStep,Double_t &w_elt,Double_t &wstep_elt,Double_t &u_elt,
                              Double_t &ustep_elt,int& fIrst_or_second)
{
  return FindBlockingSub(w.GetNrows(),
                         w.GetMatrixArray(),    1,
                         wstep.GetMatrixArray(),1,
                         u.GetMatrixArray(),    1,
                         ustep.GetMatrixArray(),1,
                         maxStep,
                         w_elt,wstep_elt,
                         u_elt,ustep_elt,
                         fIrst_or_second);
}

//______________________________________________________________________________
 Double_t TQpVar::FindBlockingSub(Int_t n,
                                 Double_t *w,    Int_t incw,
                                 Double_t *wstep,Int_t incwstep,
                                 Double_t *u,    Int_t incu,
                                 Double_t *ustep,Int_t incustep,
                                 Double_t maxStep,
                                 Double_t &w_elt,Double_t &wstep_elt,
                                 Double_t &u_elt,Double_t &ustep_elt,
                                 Int_t &fIrst_or_second)
{
  Double_t bound = maxStep;

  Int_t i = n-1;
  Int_t lastBlocking = -1;

  // Search backward so that we fInd the blocking constraint of lowest
  // index. We do this to make things consistent with MPI's MPI_MINLOC,
  // which returns the processor with smallest rank where a min occurs.
  //
  // Still, going backward is ugly!
  Double_t *pw     = w    +(n-1)*incw;
  Double_t *pwstep = wstep+(n-1)*incwstep;
  Double_t *pu     = u    +(n-1)*incu;
  Double_t *pustep = ustep+(n-1)*incustep;

  while (i >= 0) {
    Double_t tmp = *pwstep;
    if (*pw > 0 && tmp < 0) {
      tmp = -*pw/tmp;
      if( tmp <= bound ) {
        bound = tmp;
        lastBlocking = i;
        fIrst_or_second = 1;
      }
    }
    tmp = *pustep;
    if (*pu > 0 && tmp < 0) {
      tmp = -*pu/tmp;
      if( tmp <= bound ) {
        bound = tmp;
        lastBlocking = i;
        fIrst_or_second = 2;
      }
    }

    i--;
    if (i >= 0) {
      // It is safe to decrement the pointers
      pw     -= incw;
      pwstep -= incwstep;
      pu     -= incu;
      pustep -= incustep;
    }
  }

  if (lastBlocking > -1) {
    // fIll out the elements
    w_elt     = w[lastBlocking];
    wstep_elt = wstep[lastBlocking];
    u_elt     = u[lastBlocking];
    ustep_elt = ustep[lastBlocking];
  }
  return bound;
}

//______________________________________________________________________________
 void TQpVar::InteriorPoint(Double_t alpha,Double_t beta)
{
  fS.Zero();
  fX.Zero();
  fY.Zero();
  fZ.Zero();

  if (fNxlo > 0) {
    fV = alpha;
    fV.SelectNonZeros(fXloIndex);
    fGamma = beta;
    fGamma.SelectNonZeros(fXloIndex);
  }

  if (fNxup > 0) {
    fW = alpha;
    fW.SelectNonZeros(fXupIndex);
    fPhi = beta;
    fPhi.SelectNonZeros(fXupIndex);
  }

  if (fMclo > 0 ) {
    fT = alpha;
    fT.SelectNonZeros(fCloIndex);
    fLambda = beta;
    fLambda.SelectNonZeros(fCloIndex);
  }

  if (fMcup > 0) {
    fU = alpha;
    fU.SelectNonZeros(fCupIndex);
    fPi = beta;
    fPi.SelectNonZeros(fCupIndex);
  }
}

//______________________________________________________________________________
 Double_t TQpVar::Violation()
{
  Double_t viol = 0.0;
  Double_t cmin;

  if (fNxlo > 0) {
    cmin = fV.Min();
    if (cmin < viol) viol = cmin;
   
    cmin = fGamma.Min();
    if (cmin < viol) viol = cmin;
  }
  if (fNxup > 0) {
    cmin = fW.Min();
    if (cmin < viol) viol = cmin;
   
    cmin = fPhi.Min();
    if (cmin < viol) viol = cmin;
  }
  if (fMclo > 0) {
    cmin = fT.Min();
    if (cmin < viol) viol = cmin;
   
    cmin = fLambda.Min();
    if (cmin < viol) viol = cmin;
  }
  if (fMcup > 0) {
    cmin = fU.Min();
    if (cmin < viol) viol = cmin;
   
    cmin = fPi.Min();
    if (cmin < viol) viol = cmin;
  }

  return -viol;
}

//______________________________________________________________________________
 void TQpVar::ShiftBoundVariables(Double_t alpha,Double_t beta)
{
  if (fNxlo > 0) {
    fV    .AddSomeConstant(alpha,fXloIndex);
    fGamma.AddSomeConstant(beta, fXloIndex);
  }
  if (fNxup > 0) {
    fW  .AddSomeConstant(alpha,fXupIndex);
    fPhi.AddSomeConstant(beta, fXupIndex);
  }
  if (fMclo > 0) {
    fT     .AddSomeConstant(alpha,fCloIndex);
    fLambda.AddSomeConstant(beta, fCloIndex);
  }
  if (fMcup > 0) {
    fU .AddSomeConstant(alpha,fCupIndex);
    fPi.AddSomeConstant(beta, fCupIndex);
  }
}


//______________________________________________________________________________
 void TQpVar::Print(Option_t * /*option*/) const
{
  cout << "fNx  : " << fNx   << endl;
  cout << "fMy  : " << fMy   << endl;
  cout << "fMz  : " << fMz   << endl;
  cout << "fNxup: " << fNxup << endl;
  cout << "fNxlo: " << fNxlo << endl;
  cout << "fMcup: " << fMcup << endl;
  cout << "fMclo: " << fMclo << endl;

  fXloIndex.Print("fXloIndex");
  fXupIndex.Print("fXupIndex");
  fCupIndex.Print("fCupIndex");
  fCloIndex.Print("fCloIndex");

  fX.Print("fX");
  fS.Print("fS");
  fY.Print("fY");
  fZ.Print("fZ");

  fV.Print("fV");
  fPhi.Print("fPhi");

  fW.Print("fW");
  fGamma.Print("fGamma");

  fT.Print("fT");
  fLambda.Print("fLambda");

  fU.Print("fU");
  fPi.Print("fPi");
}

//______________________________________________________________________________
 Double_t TQpVar::Norm1()
{
  Double_t norm = 0.0;
  norm += fX.Norm1();
  norm += fS.Norm1();
  norm += fY.Norm1();
  norm += fZ.Norm1();

  norm += fV.Norm1();
  norm += fPhi.Norm1();
  norm += fW.Norm1();
  norm += fGamma.Norm1();
  norm += fT.Norm1();
  norm += fLambda.Norm1();
  norm += fU.Norm1();
  norm += fPi.Norm1();

  return norm;
}

//______________________________________________________________________________
 Double_t TQpVar::NormInf()
{
  Double_t norm = 0.0;

  Double_t tmp = fX.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fS.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fY.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fZ.NormInf();
  if (tmp > norm) norm = tmp;

  tmp = fV.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fPhi.NormInf();
  if (tmp > norm) norm = tmp;

  tmp = fW.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fGamma.NormInf();
  if (tmp > norm) norm = tmp;

  tmp = fT.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fLambda.NormInf();
  if (tmp > norm) norm = tmp;

  tmp = fU.NormInf();
  if (tmp > norm) norm = tmp;
  tmp = fPi.NormInf();
  if (tmp > norm) norm = tmp;

  return norm;
}

//______________________________________________________________________________
 Bool_t TQpVar::ValidNonZeroPattern()
{
  if (fNxlo > 0 && 
      ( !fV    .MatchesNonZeroPattern(fXloIndex) ||
        !fGamma.MatchesNonZeroPattern(fXloIndex) ) ) {
    return kFALSE;
  }

  if (fNxup > 0 &&
      ( !fW  .MatchesNonZeroPattern(fXupIndex) ||
        !fPhi.MatchesNonZeroPattern(fXupIndex) ) ) {
    return kFALSE;
  }
  if (fMclo > 0 &&
      ( !fT     .MatchesNonZeroPattern(fCloIndex) ||
        !fLambda.MatchesNonZeroPattern(fCloIndex) ) ) {
    return kFALSE;
  }

  if (fMcup > 0 &&
      ( !fU .MatchesNonZeroPattern(fCupIndex) ||
        !fPi.MatchesNonZeroPattern(fCupIndex) ) ) {
    return kFALSE;
  }
  
  return kTRUE;
}

//______________________________________________________________________________
TQpVar &TQpVar::operator=(const TQpVar &source)
{
  if (this != &source) {
    TObject::operator=(source);
    fNx       = source.fNx;
    fMy       = source.fMy;
    fMz       = source.fMz;
    fNxup     = source.fNxup;
    fNxlo     = source.fNxlo;
    fMcup     = source.fMcup;
    fMclo     = source.fMclo;

    fXloIndex = source.fXloIndex;
    fXupIndex = source.fXupIndex;
    fCupIndex = source.fCupIndex;
    fCloIndex = source.fCloIndex;

    fX     .ResizeTo(source.fX);      fX      = source.fX;
    fS     .ResizeTo(source.fS);      fS      = source.fS;
    fY     .ResizeTo(source.fY);      fY      = source.fY;
    fZ     .ResizeTo(source.fZ);      fZ      = source.fZ;
  
    fV     .ResizeTo(source.fV);      fV      = source.fV;
    fPhi   .ResizeTo(source.fPhi);    fPhi    = source.fPhi;
  
    fW     .ResizeTo(source.fW);      fW      = source.fW;
    fGamma .ResizeTo(source.fGamma) ; fGamma  = source.fGamma;

    fT     .ResizeTo(source.fT);      fT      = source.fT;
    fLambda.ResizeTo(source.fLambda); fLambda = source.fLambda;
  
    fU     .ResizeTo(source.fU);      fU      = source.fU;
    fPi    .ResizeTo(source.fPi);     fPi     = source.fPi;
  }
  return *this;
}