/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Wolfgang Wittek, 04/2003 <mailto:wittek@mppmu.mpg.de>
!               David Paneque,   02/2004 <mailto:dpaneque@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MSupercutsCalcONOFF                                                     //
//                                                                         //
//   this class calculates the hadronness for the supercuts                //
//   the parameters of the supercuts are taken                             //
//                  from the container MSupercuts                       //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MSupercutsCalcONOFF.h"

#include <math.h>
#include <fstream>

#include "TFile.h"
#include "TArrayD.h"

#include "MParList.h"
#include "MHillasExt.h"
#include "MHillasSrc.h"
#include "MNewImagePar.h"
#include "MPointingPos.h"
#include "MCerPhotEvt.h"
#include "MGeomCam.h"
#include "MHadronness.h"
#include "MTSupercutsApplied.h"
#include "MHMatrix.h"
#include "MSupercuts.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MSupercutsCalcONOFF);

using namespace std;


// --------------------------------------------------------------------------
//
// constructor
//

MSupercutsCalcONOFF::MSupercutsCalcONOFF(const char *hilname, 
                                     const char *hilsrcname, 
                                     const char *name, const char *title)
  : fHadronnessName("MHadronness"), fSupercutsAppliedName("MSupercutsApplied"), 
    fHilName(hilname), fHilSrcName(hilsrcname),
    fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), 
    fSuperName("MSupercuts") 
{
    fName  = name  ? name  : "MSupercutsCalcONOFF";
    fTitle = title ? title : "Class to evaluate the Supercuts";

    fMatrix = NULL;

    fStoreAppliedSupercuts = kFALSE; // by default, applied SC parameters are not stored

    fNotUseTheta = kFALSE; // by default, theta info is used in the computation of the cuts

    fUseStaticCuts = kFALSE; // by default, dynamical cuts are used


    // Usage of DIST parameter in the parameterization of the cuts.
    // For the time being is in the constructor. If finally it comes out 
    // that it is important to disable the DIST parameter from the 
    // cut parameterization I will make a function that access this 
    // data variable from outside the class.
    fUseDist = kTRUE;



    // OFFSETS FOR THE DYNAMICAL CUTS
    // Values of Size (photons), Dist (degrees) and Theta (degrees) 
    // for which the value of the dynamical cut is equal to the 
    // non dependent parameter; i.e. to the static cut. 
    // By adjusting these offsets the user can set the value in size, 
    // dist and theta for which the dynamical cuts will be given by 
    // the term that DOES not depend on size, dist and theta.

    // For the time being, these quantities are set in the constructor. 
    // In future, if they show to be useful, I will make them available 
    // as external variables.

    fSizeOffset = 3000; // photons
    fDistOffset = 0.9; // degrees  
    fThetaOffset = 20; // degrees   NOT USED FOR THE TIME BEING

    




    // Variables defining upper limits for some of the hillas params.
    // The default values are set to conservative 
    // values so that no gamma showers are removed

    // If a cut value computed by function MSupercutsCalcONOFF::CtsMCut, 
    // (i.e. widthlow) exceeds (larger or smaller 
    // depending on wether is cutUP or cutLOW) one of these values
    // (i.e. fWidthLowerLimit), such cut value is replaced by the 
    // the corresponding Limit (i.e. widthlow = fWidthLowerLimit)



    fDistUpperLimit = 1.4; // in deg
    fLengthUpperLimit = 0.6;
    fWidthUpperLimit = 0.4;

    fLeakage1UpperLimit = 0.25; 

    fLengthOverWidthUpperLimit = 1.0;

    fDistLowerLimit = 0.1;
    fLengthLowerLimit = 0.09;
    fWidthLowerLimit = 0.05;

    




}

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


void MSupercutsCalcONOFF::SetStoreAppliedSupercuts(Bool_t b)
{
    fStoreAppliedSupercuts = b;
    if (fStoreAppliedSupercuts)
    {
	*fLog << "Supercuts applied to all the individual events will be stored " 
	      << "in a container of the class MSupercutsApplied." << endl;
    }

}


void MSupercutsCalcONOFF::SetVariableNotUseTheta(Bool_t b)
{
    fNotUseTheta = b;
    if (fNotUseTheta)
    {
	*fLog << "Theta variable is NOT used in the computation of the the  " 
	      << "dynamical cuts" << endl;
    }
    else
    {
	*fLog << "Theta variable is used in the computation of the the  " 
	      << "dynamical cuts" << endl;
    }

}


void MSupercutsCalcONOFF::SetVariableUseStaticCuts(Bool_t b)
{
    fUseStaticCuts = b;
    if (fUseStaticCuts)
    {
	*fLog << "Supercuts DO NOT take into account Theta, Size and Dist dependence; "
	      << "i.e, they are STATIC CUTS." << endl;
    }
    else
    {
	*fLog << "Supercuts are computed taking into account Theta, size and Dist dependence."
	      << "i.e., they are DYNAMICAL CUTS" << endl;

    }



}


void MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits(Double_t distup, 
							      Double_t lengthup, 
							      Double_t widthup, 
							      Double_t distlow, 
							      Double_t lengthlow, 
							      Double_t widthlow)
{
  
  fDistUpperLimit = distup;
  fLengthUpperLimit = lengthup;
  fWidthUpperLimit = widthup;


  fDistLowerLimit = distlow;
  fLengthLowerLimit = lengthlow;
  fWidthLowerLimit = widthlow;


  *fLog << "MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits" << endl
	<< "Upper limits to Hillas parameters Dist, Length and Width " 
	<< "are set respectively to : " 
	<<  fDistUpperLimit << ", " 
	<< fLengthUpperLimit << ", "
	<< fWidthUpperLimit << " (in degrees!!)" << endl

	<< "Lower limits to Hillas parameters Dist, Length and Width " 
	<< "are set respectively to : " 
	<<  fDistLowerLimit << ", " 
	<< fLengthLowerLimit << ", "
	<< fWidthLowerLimit << " (in degrees!!)" << endl

	<< "The Hadronnes of those images exceeding one of these limits will " 
	<< "be set to 0.75." << endl << endl;


} 



Int_t MSupercutsCalcONOFF::PreProcess(MParList *pList)
{
    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
        *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
        return kFALSE;
    }

    fMm2Deg = cam->GetConvMm2Deg();

    fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
    if (!fHadronness)
    {
        *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
        return kFALSE;
    }

    fSuper = (MSupercuts*)pList->FindObject(fSuperName, "MSupercuts");
    if (!fSuper)
    {
        *fLog << err << fSuperName << " [MSupercuts] not found... aborting." << endl;
        return kFALSE;
    }

    if (fStoreAppliedSupercuts)
    {
	
	fSupercutsApplied = (MTSupercutsApplied*)pList->FindObject(fSupercutsAppliedName,
								  "MTSupercutsApplied");
	if(!fSupercutsApplied)
	{
	    *fLog << err << fSupercutsAppliedName 
		  << " [MTSupercutsApplied] not found... aborting." << endl;
	    return kFALSE;
	}
	

    }

    if (fMatrix)
        return kTRUE;

    //-----------------------------------------------------------
    fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
    if (!fHil)
    {
        *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
        return kFALSE;
    }

    fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
    if (!fHilExt)
    {
        *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
        return kFALSE;
    }

    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
    if (!fHilSrc)
    {
        *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
        return kFALSE;
    }

    fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
    if (!fNewPar)
    {
        *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
        return kFALSE;
    }

   fPointPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
   if (!fPointPos)
   {
       *fLog << err << "MSupercutsCalcONOFF::PreProcess; MPointingPos not found... aborting." << endl;
       return kFALSE;
   }



    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculation of upper and lower limits
//
Double_t MSupercutsCalcONOFF::CtsMCut(const Double_t* a,  Double_t ls, Double_t ct,
                                    Double_t ls2, Double_t dd2) const
{
    // define cut-function
    //
    //   
    //    dNOMCOSZA   = 1.0
    //
    //      a: array of cut parameters
    //     ls: log(SIZE) - log(fSizeOffset)
    //    ls2: ls^2
    //     ct: Cos(ZA.) - Cos(fThetaOffset)
    //    dd2: (DIST- fDistOffset)^2

    Double_t limit;

    if(fUseStaticCuts)
    { // static cuts are used
	limit = a[0];
    }
    else
    {
	if (fNotUseTheta)
	{ // Theta info is NOT used in the computation of the dynamical cuts
	  // For the time being dist info will not be used

	  if(fUseDist)
	    {

	     
	      limit = a[0] + a[1] * dd2 + 
		      ls  * (a[3] + a[4] * dd2) + 
		      ls2 * (a[6] + a[7] * dd2);
	      

	    }
	  else
	    {
	      limit =  a[0] + ls  * (a[3])+   ls2 * (a[6]);
	    }


	}
	else
	{// Theta info IS used in the computation of the dynamical cuts
	    
	    limit =
		a[0] + a[1] * dd2 + a[2] * ct  +
		ls  * (a[3] + a[4] * dd2 + a[5] * ct) +
		ls2 * (a[6] + a[7] * dd2);
	}

    }



    //*fLog << "MSupercutsCalcONOFF::CtsMCut; *a = "
    //      << *a     << ",  " << *(a+1) << ",  " << *(a+2) << ",  "
    //      << *(a+3) << ",  " << *(a+4) << ",  " << *(a+5) << ",  "
    //      << *(a+6) << ",  " << *(a+7) << endl;

    //*fLog << "MSupercutsCalcONOFF::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
    //      << ",  " << ls2 << ",  " << ct << ",  " << dd2 << ",  "
    //      << limit << endl;

    return limit;
}

// --------------------------------------------------------------------------
//
// Returns the mapped value from the Matrix
//
Double_t MSupercutsCalcONOFF::GetVal(Int_t i) const
{

    Double_t val = (*fMatrix)[fMap[i]];


    //*fLog << "MSupercutsCalcONOFF::GetVal; i, fMatrix, fMap, val = "
    //    << i << ",  " << fMatrix << ",  " << fMap[i] << ",  " << val << endl;


    return val;
}

// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of the
// given containers. This function adds all necessary columns to the
// given matrix. Afterward you should fill the matrix with the corresponding
// data (eg from a file by using MHMatrix::Fill). If you now loop
// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
// will take the values from the matrix instead of the containers.
//


void MSupercutsCalcONOFF::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
      return;

    fMatrix = mat;

    fMap[0] = fMatrix->AddColumn("MPointingPos.fZd");  //deg
    fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
    fMap[2] = fMatrix->AddColumn("MHillas.fLength");
    fMap[3] = fMatrix->AddColumn("MHillas.fSize");
    fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
    fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
    fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
    fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
    fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
    fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
    fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
}


// ---------------------------------------------------------------------------
//
// Evaluate dynamical supercuts 
// 
//          set hadronness to 0.25 if cuts are fullfilled
//                            0.75 otherwise
//
Int_t MSupercutsCalcONOFF::Process()
{
  //    const Double_t kNomLogSize = 4.1;
    const Double_t kNomCosZA   = 1.0;

    //const Double_t theta   = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
    // For the time being I do it in a way that theta must be 
    // the value used in the row 0 of fMatrix. That means that if there is 
    // no matrix, the program will complain (and crash). And tehn I will know 
    // that value 0 (supposed to be fThetaOrig.Val) could not be taken.
    const Double_t theta   = GetVal(0);
    const Double_t width0  = fMatrix ? GetVal(1) : fHil->GetWidth();
    const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
    const Double_t size    = fMatrix ? GetVal(3) : fHil->GetSize();
    const Double_t meanx   = fMatrix ? GetVal(4) : fHil->GetMeanX();
    const Double_t meany   = fMatrix ? GetVal(5) : fHil->GetMeanY();
    const Double_t dist0   = fMatrix ? GetVal(6) : fHilSrc->GetDist();

    const Double_t alpha   = fMatrix ? GetVal(7) : fHilSrc->GetAlpha();

    Double_t help;
    if (!fMatrix)
      help = TMath::Sign(fHilExt->GetM3Long(), 
	      		 fHilSrc->GetCosDeltaAlpha());
    const Double_t asym0   = fMatrix ? GetVal(8) : help;
    const Double_t conc    = fMatrix ? GetVal(9) : fNewPar->GetConc();
    const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();

    const Double_t newdist = dist0 * fMm2Deg;

    const Double_t dist2   = meanx*meanx + meany*meany;
    const Double_t dist    = sqrt(dist2) * fMm2Deg;
    
    // const Double_t dd2     = dist*dist;


    // in the parameterization of the cuts 
    // the dist parameter used is the one computed from the source position, 
    // and not the one computed from the camera center.

    // Actually the value used in the parameterization is newdist^2, 
    // minus the offset_in_dist^2
    
    const Double_t dd2     = newdist*newdist - fDistOffset*fDistOffset;


    
    const Double_t dmls    = log(size) - log(fSizeOffset);
    const Double_t dmls2   = dmls * dmls;

    const Double_t dmcza   = cos(theta) - kNomCosZA;

    const Double_t length  = length0 * fMm2Deg;
    const Double_t width   = width0  * fMm2Deg;
    const Double_t asym    = asym0   * fMm2Deg;


    
    // computation of the cut limits

    Double_t lengthup  = CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2);    
    Double_t lengthlow = CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2);

    Double_t widthup   = CtsMCut (fSuper->GetWidthUp(),  dmls, dmcza, dmls2, dd2);     
    Double_t widthlow  =  CtsMCut (fSuper->GetWidthLo(),  dmls, dmcza, dmls2, dd2);

    Double_t distup   =  CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2);
    Double_t distlow  =  CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2);

    Double_t asymup   = CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2);
    Double_t asymlow  = CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2);

    Double_t concup   = CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2);
    Double_t conclow  = CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2);

    Double_t  leakageup = CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2);
    Double_t  leakagelow = CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2); 

    
    Double_t lengthoverwidth = length/width;

    Double_t hadronness = 0.0;


    // If the cut values computed before for Dist, Length and Width exceed the upper limits set 
    // by the user through function MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits,
    // (or the default values defined in constructor); such cut values are replaced by the 
    // the limit values

    

    if (distup > fDistUpperLimit)
      {
	distup = fDistUpperLimit;
      }

   
    if (lengthup > fLengthUpperLimit)
      {
	lengthup = fLengthUpperLimit;
      }

    if (widthup > fWidthUpperLimit)
      {
	widthup = fWidthUpperLimit;
      }

    

    if (distlow < fDistLowerLimit)
      {
	distlow = fDistLowerLimit;
      }
                   
    if (lengthlow < fLengthLowerLimit)
      {
	lengthlow = fLengthLowerLimit;
      }

    if (widthlow < fWidthLowerLimit)
      {            
	widthlow = fWidthLowerLimit;
      }

    
    if (leakageup > fLeakage1UpperLimit)
      {
	leakageup = fLeakage1UpperLimit;
      }
    

   


    // Upper cut in leakage 1 also implemented

    
    
   
    if (// 
	newdist > distup ||
	newdist < distlow ||	    
	length  > lengthup ||
	length  < lengthlow ||	    
	width   > widthup ||
	width   < widthlow ||	    
	leakage > leakageup ||
	leakage < leakagelow ||
	lengthoverwidth < fLengthOverWidthUpperLimit
	//
	//asym    < asymup &&
	//asym    > asymlow &&
	
	//dist    < distup &&
	//dist    > distlow &&
	
	//conc    < concup &&
	//conc    > conclow &&	
	//
	) 
	  
      {hadronness = 0.75;}
    else
      {
	hadronness = 0.25;
      }
      
    
    fHadronness->SetHadronness(hadronness);
    fHadronness->SetReadyToSave();
    
    
    if(fStoreAppliedSupercuts)
    {

	// TArrayD vector with the shower parameters (matching the ones 
	// specified in function MTSupercutsApplied::CreateTreeBranches)
	// is created and filled with this event features

	// Eventually this definition might be done from outside this 
	// class, allowing for adding removing parameters without 
	// recompiling mars. yet for the time being, let's keep it simple.
	
	TArrayD ShowerParams(10);
	ShowerParams[0] = length;
	ShowerParams[1] = width;
	ShowerParams[2] = dist;
	ShowerParams[3] = newdist;
	ShowerParams[4] = asym;
	ShowerParams[5] = conc;
	ShowerParams[6] = leakage;
	ShowerParams[7] = theta;
	ShowerParams[8] = size;
	ShowerParams[9] = alpha;
	
	
	// TArrayD vector with the cut parameters (matching the ones 
	// specified in function MTSupercutsApplied::CreateTreeBranches)
	// is created and filled with this event features

	// Eventually this definition might be done from outside this 
	// class, allowing for adding removing parameters without 
	// recompiling mars. yet for the time being, let's keep it simple.
	
	
	TArrayD SuperCutParams(13);
	SuperCutParams[0] = lengthup;
	SuperCutParams[1] = lengthlow;
	SuperCutParams[2] = widthup;
	SuperCutParams[3] = widthlow;
	SuperCutParams[4] = distup;
	SuperCutParams[5] = distlow;
	SuperCutParams[6] = asymup;
	SuperCutParams[7] = asymlow;
	SuperCutParams[8] = concup;
	SuperCutParams[9] = conclow;
	SuperCutParams[10] = leakageup;
	SuperCutParams[11] = leakagelow;
	SuperCutParams[12] = hadronness;

	// SC parameters applied to this event, as well as the 
	// shower parameters, are stored in the branches of the
	// TTree object of a MTSupercutsApplied object

   
	if (!fSupercutsApplied ->FillTreeBranches(SuperCutParams, ShowerParams))
	{
	    *fLog << "MSupercutsCalcONOFF::Process()" << endl
		  << "Supercuts applied could not be stored in tree..."
		  << endl;
	    
	    return kFALSE;
	}
	

    }
     


    return kTRUE;


    // OLD STUFF

    /*
    *fLog << "newdist, length, width, asym, dist, conc, leakage = " 
          << newdist << ",  " << length << ",  " << width << ",  "
          << asym    << ",  " << dist   << ",  " << conc  << ",  " << leakage
          << endl;
  
    *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " 
          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetLengthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetLengthLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetWidthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetWidthLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) << ",  "

          << CtsMCut (fSuper->GetLeakage1Up(),   dmls, dmcza, dmls2, dd2) << ",  "
          << CtsMCut (fSuper->GetLeakage1Lo(),   dmls, dmcza, dmls2, dd2) << ",  "
          << endl;
    */

    /*
    if (
        //dist    < 1.05                                                     &&
        //newdist < 1.05                                                     &&

        newdist < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
        newdist > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&

        length  < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
        length  > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&

        width   < CtsMCut (fSuper->GetWidthUp(),  dmls, dmcza, dmls2, dd2) &&
        width   > CtsMCut (fSuper->GetWidthLo(),  dmls, dmcza, dmls2, dd2) &&

        asym    < CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) &&
        asym    > CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) &&

        dist    < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
        dist    > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&

        conc    < CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) &&
        conc    > CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) &&

        leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
        leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2)  ) 

      fHadronness->SetHadronness(0.25);
    else
      fHadronness->SetHadronness(0.75);

   

    // *fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;

    fHadronness->SetReadyToSave();
    */

    // END OF OLD STUFF
}

/*
Bool_t MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent(TArrayD CutParams, 
							   TArrayD ShowerParams)
{
    // SC parameters applied to this event are stored in 
    // TTree object of MTSupercutsApplied object

   
    if (!fSupercutsApplied ->FillTreeBranches(CutParams, ShowerParams))
    {
	*fLog << "MSupercutsCalcONOFF::StoreSupercutsAppliedToThisEvent" << endl
	      << "Supercuts applied could not be stored in tree..."
	      << endl;

	return kFALSE;
    }

    
    return kTRUE;
}

*/








