/***************************************************************************
 * Author: Dominik Flierl, flierl@ikf.uni-frankfurt.de
 ***************************************************************************/

// Revision 2004/09/10 skniege 
// Add Ceres methods  T49HbtFit_BP_LS_CHI2::TheoreticalValue()


using namespace std ;

#include <iostream>
#include <stdlib.h>
#include "T49HbtFitDivers.h"

ClassImp(T49HbtFit_BP_LS_CHI2)
ClassImp(T49HbtFit_BP_LS_MML)
ClassImp(T49HbtFit_BP_ULS_CHI2)
ClassImp(T49HbtFit_YKP_LS_CHI2)


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtFit_BP_LS_CHI2::mfcn(Int_t &nParameters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag) 
{
  // calculate chi2 of the given paramset 'parameterSet'
            
  // loop over Internal Array and caluclate chi2
  double chi2 = 0.0 ;
  for (Int_t index = 0 ;  index < mInternalArraySize ; index++)
    {
      // calculate the expected52 value at this point (variablePosition) for this set of paramters, 
      // attention the return value might be coulomb corrected !
      
     
      
      Double_t theoreticalValue  = TheoreticalValue( &mVariablePositionArray[index], parameterSet) ; 
      //cout <<  &mVariablePositionArray[index] << endl ;

      // now we have to consider the normalization
      theoreticalValue = theoreticalValue *  parameterSet[0] ;
      
      // get the difference between measured and predicted value
      double delta = (mNumeratorInternalArray[index]/mDenominatorInternalArray[index]-theoreticalValue)/mErrorInternalArray[index];
     
      // sum the chi2
      chi2 += delta*delta;
    }
    
  // chi2 is calculated
  finalChi2 = chi2 ;
};
//____________________________
Double_t T49HbtFit_BP_LS_CHI2::TheoreticalValue(TVector3* position, double* parameterSet )
{

  //if(position->X()<0.0) { cout << "unmoeglich !" << endl ; } 

  // return C2 at 'position' with 'parameterSet' for Bertsch-Pratt parametrization, if set, with coulomb correction

  // position->x = qside
  // position->y = qout 
  // position->z = qlong  
  // parameterSet[0] = normalization
  // parameterSet[1] = lambda        // corresponds to the measured reduction of the correlationfunction for CERES fit procedure 
  // parameterSet[2] = Rside 
  // parameterSet[3] = Rout  
  // parameterSet[4] = Rlong
  // parameterSet[5] = Routlong
  
  Double_t valueFromCorrelationFunction = 1.0;
  Double_t coulombWeight = 1.0;

  switch (fCoulombFit) 
    {
      
    case 0 :   //Sinyukov coulomb weight (Na49 fit function)
	
      valueFromCorrelationFunction  = 1.0 + parameterSet[1] * exp(
								  (
								   (-1.0) * pow(position->x(),2) * pow(parameterSet[2],2) 
								   - pow(position->y(),2) * pow(parameterSet[3],2) 
								   - pow(position->z(),2) * pow(parameterSet[4],2) 
								   - 2* position->y() * position->z() * pow(parameterSet[5],2) 
								   ) 
								  / mhc2 );
      
      if(unCorrectedBackground == 1 && mCoulombCalculator)
	{
	  coulombWeight =  
	    mCoulombCalculator->Weight(mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
								       mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
								       mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z()))) ;
	}
      
      // see equation 7 : C++ = p A Cqs + (1-p) 
      //valueFromCorrelationFunction =  mPurity * coulombWeight * valueFromCorrelationFunction +  (1.0-mPurity)  ;
      // added mom resolution to compare to ceres procedure  
      valueFromCorrelationFunction =  mPurity * valueFromCorrelationFunction*(mMomentumRes*coulombWeight  +(1-mMomentumRes)) +  (1.0-mPurity)  ;
     
      
      break;
      
    case 1 : // PBM coulomb weight (Na49 fit function)
     
      valueFromCorrelationFunction  = 1.0 + parameterSet[1] * exp(
								  (
								   (-1.0) * pow(position->x(),2) * pow(parameterSet[2],2) 
								   - pow(position->y(),2) * pow(parameterSet[3],2) 
								   - pow(position->z(),2) * pow(parameterSet[4],2) 
								   - 2* position->y() * position->z() * pow(parameterSet[5],2) 
								   ) 
								  / mhc2 );
      if(unCorrectedBackground == 1 && mCoulombCalculator)
	{
	  coulombWeight = 
	    mCoulombCalculator->ClassicalWeight(mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
										mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
										mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z()))) ;
	}
      valueFromCorrelationFunction =  mPurity * coulombWeight * valueFromCorrelationFunction +  (1.0-mPurity)  ;
      
    
      
      break;
      
    case 2 :  // PBM coulomb weight ( CERES fit function )
      
      valueFromCorrelationFunction  = 1.0 +  exp(
						 (
						   (-1.0) * pow(position->x(),2) * pow(parameterSet[2],2) 
						   - pow(position->y(),2) * pow(parameterSet[3],2) 
						   - pow(position->z(),2) * pow(parameterSet[4],2) 
						   - 2* position->y() * position->z() * pow(parameterSet[5],2) 
						   ) 
						  / mhc2 );
      if(unCorrectedBackground == 1 && mCoulombCalculator)
	{
	  coulombWeight = 
	    mCoulombCalculator->ClassicalWeight(mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
										mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
										mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z()))) ;
	}

      valueFromCorrelationFunction =  parameterSet[1]*(valueFromCorrelationFunction*(mMomentumRes*coulombWeight  +(1-mMomentumRes))) +  (1.0-parameterSet[1])  ;
      
      break;
      
    }
  
  if(valueFromCorrelationFunction==0)
    {       cout << " coulombWeight " << coulombWeight << endl ;
    cout << " valueFromCorrelationFunction    " <<  valueFromCorrelationFunction << endl;
    }
  // return result
  return valueFromCorrelationFunction ;
}

//____________________________
void T49HbtFit_BP_LS_CHI2::InitMinuit() 
{
  // initialization with default start values

  // Errorflag
  Int_t ierflg = 0 ;
  mMinuit->mnparm(0, "Normalization", 0.1, 0.01, 0,0,ierflg);
  mMinuit->mnparm(1, "lambda",        0.5, 0.1, 0,0,ierflg);
  mMinuit->mnparm(2, "Rside",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(3, "Rout",          5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(4, "Rlong",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(5, "Routlong",      5.0,  0.1, 0,0,ierflg);

}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtFit_BP_LS_MML::mfcn(Int_t &nParameters, Double_t *gin, Double_t &finalMLH, Double_t *parameterSet, Int_t iflag) 
{
  // calculate maximum likelihood asuming an underlying poisson distribution, for details see for example
  // http://nedwww.ipac.caltech.edu/level5/Leo/Stats_contents.html
  /*
    poison distribution : p(r) = y^r * exp(-y) / r!      y = mean
    likelihood distribution to be maximized L, n number of measurements (x1....xn) :
    L = -n*y + Sum(xi * ln(y)) - Sum(ln(xi!)) = Sum( -y + xi * ln(y) - ln(xi!) )
  */
      
  // some definitions
  Double_t MLH  = 0.0 ;
             
  // loop over Internal Array and caluclate chi2 
  for (Int_t index = 0 ;  index < mInternalArraySize ; index++)
	{
	  // calculate the expected value at this point (variablePosition) for this set of paramters
	  Double_t expectedValue  = TheoreticalValue( &mVariablePositionArray[index] , parameterSet) * mDenominatorInternalArray[index] ; 

	  // now we have to consider the normalization
	  expectedValue = expectedValue  *  parameterSet[0] ;
	  
	  // measured value
	  Double_t measueredValue = mNumeratorInternalArray[index] ;
	  // sum up
	  MLH += -expectedValue + measueredValue * log(expectedValue) - Lnfactorial(measueredValue) ;
	}
    // likelihood for this set of paramters
    finalMLH = - MLH ;
}
//____________________________
Double_t T49HbtFit_BP_LS_MML::TheoreticalValue(TVector3* position, double* parameterSet )
{
  // return C2 at 'position' with 'parameterSet' for Bertsch-Pratt parametrization, if set, with coulomb correction

  // position.x = qside
  // position.y = qout 
  // position.z = qlong  
  // parameterSet[0] = normalization
  // parameterSet[1] = lambda 
  // parameterSet[2] = Rside 
  // parameterSet[3] = Rout  
  // parameterSet[4] = Rlong
  // parameterSet[5] = Routlong
  Double_t valueFromCorrelationFunction  = 1.0 + parameterSet[1] * exp(
						       (
							(-1.0) * pow(position->x(),2) * pow(parameterSet[2],2) 
							- pow(position->y(),2) * pow(parameterSet[3],2) 
							- pow(position->z(),2) * pow(parameterSet[4],2) 
							- 2* position->y() * position->z() * pow(parameterSet[5],2) 
							) 
						       / mhc2 );
  
  // now check coulomb must be included
  if(unCorrectedBackground == 1 && mCoulombCalculator)
    {
      Double_t coulombWeight = 
	mCoulombCalculator->Weight(mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
								   mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
								   mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z()))) ;

      // see equation 7 : C++ = p A Cqs + (1-p) 
      valueFromCorrelationFunction =  mPurity * coulombWeight * valueFromCorrelationFunction +  (1.0-mPurity)  ;
    
    }

  // return result
  return valueFromCorrelationFunction ;
}
//____________________________
void T49HbtFit_BP_LS_MML::InitMinuit() 
{  
  // initialization with default start values

  // Errorflag
  Int_t ierflg = 0 ;
  mMinuit->mnparm(0, "Normalization", 0.05, 0.01, 0,0,ierflg);
  mMinuit->mnparm(1, "lambda",        0.5, 0.1, 0,0,ierflg);
  mMinuit->mnparm(2, "Rside",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(3, "Rout",          5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(4, "Rlong",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(5, "Routlong",      5.0,  0.1, 0,0,ierflg);
  
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void T49HbtFit_BP_ULS_CHI2::mfcn(Int_t &nParameters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag) 
{
  // calculate chi2 of the given paramset 'parameterSet'
  Double_t chi2 = 0.0 ;

  // set coulomb radius, according to parameterSet (the second parameter (normalization) is varied below)
  mCoulombCalculator->SetMeanSeparation(parameterSet[1]) ;

  // loop over Internal Array and caluclate chi2 
  for (Int_t index = 0 ;  index < mInternalArraySize ; index++)
    {
      // get the coulomb weight for the mean qinv in this (qout,qside,qlong)-bin for a given coulomb radius and vary normalization
      Double_t theoreticalValue  =  parameterSet[0] * TheoreticalValue(&mVariablePositionArray[index], parameterSet) ;
     	  
      // get the difference between measured and predicted value
      double delta = (mNumeratorInternalArray[index]/mDenominatorInternalArray[index]-theoreticalValue)/mErrorInternalArray[index];
      
      // sum the chi2
      chi2 += delta*delta;
    }
   
  // chi2 has been calculated:
  finalChi2 = chi2 ;
};
//____________________________
 Double_t T49HbtFit_BP_ULS_CHI2::TheoreticalValue(TVector3* position, double* parameterSet )
{
  // calculate the expected value in this (qout,qside,qlong) cell
  // in the coulomb case (unlikesign pairs) the expected value only depends on qinv, so first we go and get the mean qinv
  Double_t qinv = mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
						  mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
						  mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z())) ;

  // get the coulomb weight for this qinv, where the coulomb radius has been set before !
  Double_t theoreticalValue  = mCoulombCalculator->WeightUnlikeSignPairs(qinv) ;
  
  // if you want it the classical way, use this one !
  // Double_t theoreticalValue  = mCoulombCalculator->ClassicalWeightUnlikeSignPairs(qinv) ;   

  // take the purity into account
  theoreticalValue =  mPurity * theoreticalValue + (1.0-mPurity)  ;
  
  return theoreticalValue ;
}
//____________________________
 void T49HbtFit_BP_ULS_CHI2::InitMinuit() 
{
  // initialization with default start values

  // Errorflag
  Int_t ierflg = 0 ;
  mMinuit->mnparm(0, "Normalization",  0.11, 0.01,  0,0,ierflg) ;
  mMinuit->mnparm(1, "MeanSeparation", 9.0,  1.0, 0,0,ierflg) ;
  
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
void T49HbtFit_YKP_LS_CHI2::mfcn(Int_t &nParameters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag) 
{
  // calculate chi2 of the given paramset 'parameterSet'
            
  // loop over Internal Array and caluclate chi2
  double chi2 = 0.0 ;
  for (Int_t index = 0 ;  index < mInternalArraySize ; index++)
    {
      // calculate the expected value at this point (variablePosition) for this set of paramters, 
      // attention the return value might be coulomb corrected !
      Double_t theoreticalValue  = TheoreticalValue( &mVariablePositionArray[index], parameterSet) ; 

      // now we have to consider the normalization
      theoreticalValue = theoreticalValue *  parameterSet[0] ;
      
      // get the difference between measured and predicted value
      double delta = (mNumeratorInternalArray[index]/mDenominatorInternalArray[index]-theoreticalValue)/mErrorInternalArray[index];
     
      // sum the chi2
      chi2 += delta*delta;
    }
    
  // chi2 is calculated
  finalChi2 = chi2 ;
};
//____________________________
Double_t T49HbtFit_YKP_LS_CHI2::TheoreticalValue(TVector3* position, double* parameterSet )
{
  // return C2 at 'position' with Parameterset 'parameterSet' for YKP parametrization
    
  // position.x = qlong 
  // position.y = qperp 
  // position.z = qnull    
  // parameterSet[0] = constant 
  // parameterSet[1] = lambda 
  // parameterSet[2] = Rlong  
  // parameterSet[3] = Rperp  
  // parameterSet[4] = Rnull
  // parameterSet[5] = beta
  
  Double_t gamma2 = 1.0/fabs(1.0-pow(parameterSet[5],2)) ;
  Double_t valueFromCorrelationFunction =  parameterSet[0] +  parameterSet[1] * exp(
						    (
						     (-1.0) * pow(position->y(),2) * pow(parameterSet[3],2) 
						     - gamma2 * pow((position->x() - parameterSet[5]*position->z()),2) * pow(parameterSet[2],2)
						     - gamma2 * pow((position->z() - parameterSet[5]*position->x()),2) * pow(parameterSet[4],2)
						     ) 
						    / double(mhc2) );
   
  // now check coulomb must be included
  if(unCorrectedBackground == 1 && mCoulombCalculator)
    {
      Double_t coulombWeight = 
	mCoulombCalculator->Weight(mMeanQinv3dHisto->GetBinContent(mMeanQinv3dHisto->GetXaxis()->FindBin(position->X()),
								   mMeanQinv3dHisto->GetYaxis()->FindBin(position->Y()),
								   mMeanQinv3dHisto->GetZaxis()->FindBin(position->Z()))) ;

      // see equation 7 : C++ = p A Cqs + (1-p) 
      valueFromCorrelationFunction =  mPurity * coulombWeight * valueFromCorrelationFunction +  (1.0-mPurity)  ;
    
    }

  // return result
  return valueFromCorrelationFunction ;
}
//____________________________
void T49HbtFit_YKP_LS_CHI2::InitMinuit() 
{
  // initialization with default start values
  
  // Errorflag
  Int_t ierflg = 0 ;
  mMinuit->mnparm(0, "Normalization", 0.1, 0.01, 0,0,ierflg);
  mMinuit->mnparm(1, "lambda",        0.05, 0.1, 0,0,ierflg);
  mMinuit->mnparm(2, "Rlong",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(3, "Rperp",          5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(4, "Rnull",         5.0,  0.1, 0,0,ierflg);
  mMinuit->mnparm(5, "beta",      5.0,  0.1, 0,0,ierflg);
  
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////


ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.