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

// Revision 2004/09/10 skniege 
// pass pointer to TheoreticalValue() instead of reference ( in get1dProjection(..) ( in get2dProjection(..)) ) 
// changed boarders for projections (errors) in get1dProjection(..) (for -loop)
// changed core to qube  
// changed protection against fluctuation 

using namespace std ;

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

ClassImp(T49HbtFitBase)

// this is the stupid way ROOT wants it: the fcn which is going to be minimzed must
// either be static or a global function. here it's global to make sur that we got
// the connection to this class we need the  TMinuit->SetObjectFit(this) statment in doFit()  
void gfcn(Int_t &nParamters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag)
{
  T49HbtFitBase* dummy = dynamic_cast<T49HbtFitBase*>(gMinuit->GetObjectFit()) ;
  dummy->mfcn(nParamters, gin, finalChi2 , parameterSet, iflag) ;
}

//____________________________________________________________________________________________________________________________________________
T49HbtFitBase::T49HbtFitBase(TMinuit* gMinuit,TH3D* numerator, TH3D* denominator)
{
  // Get TMinuit in
  mMinuit = gMinuit ;
       
  // Get numerator and denominator in
  mNumerator = numerator ; 
  mDenominator = denominator ;
    
  // Create ratio which will be filled later
  mRatio = (TH3D*) mNumerator->Clone("Correlationfunction") ;
  mRatio->Reset() ;
    
  // Set constants
  mhc  = 0.197326960277 ; // GeV * fm
  mhc2 = 0.038937929230 ; // GeV^2 * fm^2 ;
 
  // Set defaults for thresholds = 1 means accept every bin with at least one entry
  mThresholdNumerator = 1.0 ;
  mThresholdDenominator = 1.0 ;

  // Per default the normalization is a free fit parameter 
  mfixNormalization = 0 ;

  // Set defaults for fitting regions
  mSphereLimit = 0.5 ;
  mCore = 0.0 ;

  // Set default for coulomb 
  unCorrectedBackground = 0 ;
  mMeanQinv3dHisto = NULL ;
  mPurity = 1.0 ;
  
  // CERES like fitting needs these parameters
  mMomentumRes = 1.0;
  fCoulombFit  = 0  ;
  
  
  // Set arrays to null
  mNumeratorInternalArray   = NULL ;
  mDenominatorInternalArray = NULL ;
  mErrorInternalArray       = NULL ;
  mVariablePositionArray    = NULL ;
  mCenterOfGravityInQX      = NULL ;
  mCenterOfGravityInQY      = NULL ;
  mCenterOfGravityInQZ      = NULL ;

}
//_______________
 T49HbtFitBase::~T49HbtFitBase()
{
  // Delete arrays;
  if(mNumeratorInternalArray)   delete [] mNumeratorInternalArray ;
  if(mDenominatorInternalArray) delete [] mDenominatorInternalArray ;
  if(mErrorInternalArray)       delete [] mErrorInternalArray ;
  if(mVariablePositionArray)    delete [] mVariablePositionArray ;
}
//____________________________
 Int_t T49HbtFitBase::FillInternalArrays(Double_t previouslyUsedNormalization)
{
  //  Some output
  cout << " Filling internal arrays with numerator threshold : " << mThresholdNumerator 
       << "  denominator threshold : " << mThresholdDenominator << endl ;
    
    
  // Get the maximum number of bins
  Int_t maximumInternalArraySize = (mNumerator->GetNbinsX()+2) * (mNumerator->GetNbinsY()+2) * (mNumerator->GetNbinsZ()+2) ;
  if ( maximumInternalArraySize != (mDenominator->GetNbinsX()+2) * (mDenominator->GetNbinsY()+2) * (mDenominator->GetNbinsZ()+2) )
    cerr << "Warning : different histogram sizes in numerator and denominatorn " ;
    
  // Now create arrays
  mNumeratorInternalArray   = new Double_t[maximumInternalArraySize] ;
  mDenominatorInternalArray = new Double_t[maximumInternalArraySize] ;
  mErrorInternalArray       = new Double_t[maximumInternalArraySize] ;
  mVariablePositionArray    = new TVector3[maximumInternalArraySize] ;
        
  // Count number of actually used bins
  mInternalArraySize = 0 ;
    
  // Loop over histogram and fill values into vectors 
  for (Int_t index = 0 ; index < maximumInternalArraySize ; index++)
    {
      if ( mNumerator->GetBinContent(index) > mThresholdNumerator && mDenominator->GetBinContent(index) > mThresholdDenominator )
	{
	  // Take center of the BIN ( as corresponding position in qlong,qside,qnull )
	  Int_t binX = 0 ; 
	  Int_t binY = 0 ; 
	  Int_t binZ = 0 ;
	  Bin1ToBin3(mNumerator, index, binX, binY, binZ) ;
	  Double_t qx = mNumerator->GetXaxis()->GetBinCenter(binX) ;
	  Double_t qy = mNumerator->GetYaxis()->GetBinCenter(binY) ;
	  Double_t qz = mNumerator->GetZaxis()->GetBinCenter(binZ) ;

	  //  Check if bin is in fitting range
	  if ( TMath::Abs(qx) < mSphereLimit && TMath::Abs(qy)< mSphereLimit &&  TMath::Abs(qz)< mSphereLimit && 
	      (TMath::Abs(qx) > mCore        || TMath::Abs(qy)> mCore        ||  TMath::Abs(qz) > mCore) )
	    {
	      // Fill the 3d histograms into vectors
	      Double_t num   = mNumerator->GetBinContent(index) ;
	      Double_t denom = mDenominator->GetBinContent(index) ;
	      // Fill numerator
	      mNumeratorInternalArray[mInternalArraySize]    = num ;
	      // Normalize and fill background
	      mDenominatorInternalArray[mInternalArraySize]  = denom ;
			    
	      // The error is cacluated using Gauss' error propagation : e = num/denom * sqrt(1/num + 1/(denom))
	      Double_t error = num/(denom) * sqrt( (1.0/num) + (1.0/denom) ) ;
	      // In case we did the normalization before
	      if(previouslyUsedNormalization!=0.0)
		{
		  if(previouslyUsedNormalization>1.0) cerr << "I need the normalziation, you gave me 1.0/normalziation" << endl ;
		  mfixNormalization = 1 ;
		  // error = num/(denom*(1.0/previouslyUsedNormalization)) * sqrt( (1.0/num) + 1.0/(denom*(1.0/previouslyUsedNormalization)) ) ; // flierl
		  
		  // first calculate the error of the unscaled ratio
		  error = num/(denom*(1.0/previouslyUsedNormalization)) * sqrt( (1.0/num) + 1.0/(denom*(1.0/previouslyUsedNormalization)) ) ;
		  // the rescale the errors with the normalization constant
		  error *= 1.0/previouslyUsedNormalization;
		  
		}
	      mErrorInternalArray[mInternalArraySize]  = error ;

	      if(mCenterOfGravityInQX == NULL || mCenterOfGravityInQY == NULL || mCenterOfGravityInQZ == NULL )
		{
		  // The most simple case : take center of the BIN ( as according position in qlong,qside,qnull )
		  mVariablePositionArray[mInternalArraySize].SetX(qx) ;
		  mVariablePositionArray[mInternalArraySize].SetY(qy) ;
		  mVariablePositionArray[mInternalArraySize].SetZ(qz) ;
		}
	      else
		{
		  // More realisic: take center of gravity in qside,qlong,qout cell
		  mVariablePositionArray[mInternalArraySize].SetX(mCenterOfGravityInQX->GetBinContent(index)) ;
		  mVariablePositionArray[mInternalArraySize].SetY(mCenterOfGravityInQY->GetBinContent(index)) ;
		  mVariablePositionArray[mInternalArraySize].SetZ(mCenterOfGravityInQZ->GetBinContent(index)) ;
		}

	      // Some output for debugging only
	      if (num/(denom) <0.8 && num/(denom)>0.0 && 0) 
		{
		  cout << (Double_t) num << "t" << (Double_t) (denom) ; 
		  cout << "t" << (Double_t) (num/(denom)) ;
		  Double_t xx = mNumerator->GetXaxis()->GetBinCenter(binX) ;
		  Double_t yy = mNumerator->GetYaxis()->GetBinCenter(binY) ;
		  Double_t zz = mNumerator->GetZaxis()->GetBinCenter(binZ) ;
		  cout << "xa " << xx << "t" << binX << "t" ;
		  cout << "ya " << yy << "t" << binY << "t" ;
		  cout << "za " << zz << "t" << binZ << endl ;
		}
			  
	      // Fill 3d ratio for projection purposes
	      mRatio->SetBinContent(index,num/(denom)) ;
	      mRatio->SetBinError(index, error) ;
	      
	      // Count entries
	      mInternalArraySize++ ;
	    }
	}
    }
  // Done
  cout << "Internal array size :" << mInternalArraySize << endl ;
  return mInternalArraySize; // for chi^2/NDF
}
//________________________________
 void T49HbtFitBase::SetHistos(TH3D* numerator, TH3D* denominator)
{
  // Set the numerator and denominator
  mNumerator = numerator ;
  mDenominator = denominator ;

  // Calculate the ratio 
  mRatio = (TH3D*) numerator->Clone() ;
  mRatio->Reset() ;
  mRatio->Divide(numerator,denominator) ;
}
//________________________________
 void T49HbtFitBase::SetCentersOfGravity(TH3D* CenterOfGravityInQX, TH3D* CenterOfGravityInQY, TH3D* CenterOfGravityInQZ) 
{
  mCenterOfGravityInQX = CenterOfGravityInQX ;
  mCenterOfGravityInQY = CenterOfGravityInQY ;
  mCenterOfGravityInQZ = CenterOfGravityInQZ ;
}
//____________________________
 void T49HbtFitBase::Bin1ToBin3(TH3D* histo, Int_t bin, Int_t& binx, Int_t& biny, Int_t&  binz)
{
  // Actually this should be provided by ROOT ... grrrr
  // return bin number of every axis (binx,biny,binz) if the general bin number (bin) is known
  // general bin = binX + binY * (nbinX+2) + binZ * ((nbinx+2) * (nbiny+2) )   (+2 acounts for over and underflow!)
  Int_t nbbinX = (histo->GetNbinsX()) + 2 ;
  Int_t nbbinY = (histo->GetNbinsY()) + 2 ;
    
  binx = bin%nbbinX ;
  biny = (bin/nbbinX) % nbbinY ;
  binz = (bin/nbbinX) / nbbinY ;
    
}
//____________________________
 TH3D* T49HbtFitBase::Ratio()
{
  // If normalization was part of the fit we have to scale the ratio histogram for projection etc. ...
  Double_t currentValue ;
  Double_t currentError ;
  mMinuit->GetParameter(0, currentValue, currentError) ;
  mRatio->Scale(1.0/currentValue) ; 
  return mRatio ;
}
//____________________________
 void T49HbtFitBase::SetCoulombCalculator(TH3D* MeanQinv3dHisto, T49HbtCoulomb* CoulombCalculator, Double_t purity , Double_t momentumres) 
{
  // Switch coulomb correction on
  unCorrectedBackground = 1 ; 
  mMeanQinv3dHisto = MeanQinv3dHisto ; 
  mCoulombCalculator = CoulombCalculator ;  
  if(!mMeanQinv3dHisto || !mCoulombCalculator) cerr << "Error: Didn't get histos correctly ! " << endl ; 
  mPurity = purity ;
  mMomentumRes = momentumres;
}
//________________________________
 Double_t T49HbtFitBase::Lnfactorial(Double_t arg)
{
    // Return the factorial of arg
    Double_t fac = 1.0 ;
    Int_t iarg = (Int_t) arg ;
    if(iarg<50)
	{
	    for (int index = 1; index < iarg+1 ; index++)
		{ fac *=(Double_t)index; } ;
	    fac = log(fac) ;
	}
    else 
	{
	  // If the argument is too large : use Stirlings formula
	    fac = 0.91 + (iarg+0.5) * log((float)iarg) - iarg ;
	}

    return fac;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////

T49FitProjector::T49FitProjector(T49HbtFitBase* dFitter) 
{ 
  // Init my fitter
  mFitter = dFitter ;
  
  // Get the 3d hist
  mRatio = mFitter->Ratio() ;
  
  // If the ratio histogram does not exist 
  if (!mRatio)
    {
      cerr << " Attention no ratio histogram ! " << endl ;
    }

  // Get minuit to extract the fit result
  TMinuit* tMinuit = dFitter->getMinuit() ;

  // Fill fit result into array
  mFitParameters = new Double_t[6] ;
  mFitParameterErrors = new Double_t[6] ;
  cout << "Found fit parameters : " << "t" ;
  for( Int_t index = 0 ; index < 6 ; index++ )
    {
      Double_t value = 0 ;
      Double_t error = 0 ;
      tMinuit->GetParameter(index, value, error) ;
      mFitParameters[index] = value ;
      mFitParameterErrors[index]= error;
      cout << value << "t" ;
    }
  cout << endl ;
} ;
//___________________
T49FitProjector::~T49FitProjector() 
{ 
  delete [] mFitParameters ;
  delete [] mFitParameterErrors ;
} ;
//___________________
void T49FitProjector::Clear()
{
  if(m1dpro)      delete m1dpro ;
  if(m1dfit)      delete m1dfit ;
  if(m1dweight)   delete m1dweight ;
  if(m1dcounter)  delete m1dcounter ;
  if(m1derror)    delete m1derror ;
}
//___________________
TH1D* T49FitProjector::get1dProjection(TString axis, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax)
{
  // The 1d projection is done by summing over the two axis not equal to "axis" in the function call.
  // The summing is done by weighting each addend with its inverse error  S = sigma(a_i * (1.0/e_i) / E) with E = sigma(1.0/e_i)
  // The error is calculated based on regular error propagation :  totalerror = sqrt(n) / E  with n = number of addends

  // Get ranges of the ratio 3d histogram 
  Int_t xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
  Int_t xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;  Int_t yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
  Int_t ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
  Int_t zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
  Int_t zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;

  // Calculate size of the projection histo
  Int_t proNbins = 0 ;
  Double_t proMin = 0.0 ;
  Double_t proMax = 0.0 ;
  char* axistitle = new char[60] ;
  if (axis == "x")        
    {  
      proNbins = abs(xmaxBin-xminBin)+1 ;  
      proMin   = mRatio->GetXaxis()->GetBinLowEdge(xminBin) ; 
      proMax   = mRatio->GetXaxis()->GetBinUpEdge(xmaxBin) ; 
      strcpy(axistitle,mRatio->GetXaxis()->GetTitle()) ; 
    }
  else if (axis == "y")   
    {  
      proNbins = abs(ymaxBin-yminBin)+1 ;  
      proMin   = mRatio->GetYaxis()->GetBinLowEdge(yminBin) ; 
      proMax   = mRatio->GetYaxis()->GetBinUpEdge(ymaxBin) ; 
      strcpy(axistitle,mRatio->GetYaxis()->GetTitle()) ; 
    }
  else if (axis == "z")   
    {  
      proNbins = abs(zmaxBin-zminBin)+1 ;  
      proMin   =  mRatio->GetZaxis()->GetBinLowEdge(zminBin) ; 
      proMax   =  mRatio->GetZaxis()->GetBinUpEdge(zmaxBin) ; 
      strcpy(axistitle,mRatio->GetZaxis()->GetTitle()) ; }
  else 
    { 
      cerr << "Wrong axis in projector, this should not happen ! "<< endl ; 
      return 0; 
    } ;

  
  // Create 1d histos
  Int_t ran = rand() ; // to distinguish between objects give them a random id ...
 
  // Result histo
  m1dpro = new TH1D("my1dpro","my1dpro",proNbins,proMin,proMax) ; 
  m1dpro->SetXTitle(axistitle) ;
  m1dpro->SetTitle(axistitle) ;
  Char_t restitle[100] ; 
  sprintf(restitle,"fit%s%dn",axis.Data(),ran) ;
  m1dpro->SetName(restitle) ;

  // Fit histo
  m1dfit = new TH1D("my1dfit","my1dfit",proNbins,proMin,proMax) ;
  Char_t fittitle[100] ; 
  sprintf(fittitle,"fit%s%dn",axis.Data(),ran) ;
  m1dfit->SetName(fittitle);
  m1dfit->SetLineColor(2) ;
  m1dfit->SetLineStyle(1) ;
  m1dfit->SetLineWidth(1) ;

  // Weight histo
  m1dweight = new TH1D("my1dweight","my1dweight",proNbins,proMin,proMax) ;
  Char_t weightmtitle[100] ; 
  sprintf(weightmtitle,"weightm%s%dn",axis.Data(),ran) ;
  m1dweight->SetName(weightmtitle);

  // Counter histo
  m1dcounter = new TH1D("my1dcounter","my1dcounter",proNbins,proMin,proMax) ;
  Char_t countermtitle[100] ; 
  sprintf(countermtitle,"counterm%s%dn",axis.Data(),ran) ;
  m1dcounter->SetName(countermtitle);

  // Fill histo
  for(Int_t indexX = xminBin ; indexX <= xmaxBin ; indexX++)
    {
      for(Int_t indexY = yminBin ; indexY <= ymaxBin ; indexY++)
	{
	  for(Int_t indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
	    {
	      
	      // Get content of the 3d cell in ratio histo
	      Double_t cellContent = mRatio->GetBinContent(indexX,indexY,indexZ) ;
	      
	      
	      // Get error of this cell
	      Double_t cellContentError = mRatio->GetBinError(indexX,indexY,indexZ) ;
	      
	      // Get the content of this 3d cell using the fit function
	      TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
				mRatio->GetYaxis()->GetBinCenter(indexY),
				mRatio->GetZaxis()->GetBinCenter(indexZ));
	      Double_t cellContentFromFit = mFitter->TheoreticalValue(&position,mFitParameters) ;
	      
	      // Fill the 1d histos
	      Int_t proBin = -100 ;
	      if (axis == "x")        { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetXaxis()->GetBinCenter(indexX)); }  
	      else if (axis == "y")   { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetYaxis()->GetBinCenter(indexY)); }
	      else if (axis == "z")   { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetZaxis()->GetBinCenter(indexZ)); }
	      
	      // Accept only meaningfull cells
	      if( cellContent > 0.0 && cellContentError > 0.0 )
		{
		  // Fill projection 1d histo, weighted with 1/error
		  m1dpro->AddBinContent(proBin,cellContent/cellContentError ) ;
		  
		  // Fill fit 1d histo equally weighted with 1/error
		  m1dfit->AddBinContent(proBin,cellContentFromFit/cellContentError ) ;
		  
		  // Fill 1/error histogram to normalize
		  m1dweight->AddBinContent(proBin,1.0/cellContentError) ;
		  
		  // Count number of addends (German: Summanden :))
		  m1dcounter->AddBinContent(proBin,1.0) ;
		  
		  // Some ouput for debugging only
	   	  if ( proBin == 10 && 0) 
		    {
		       cout << proBin << " " ;
		      cout << "posx " << mRatio->GetXaxis()->GetBinCenter(indexX) << "  ";
		      cout << "posy " << mRatio->GetYaxis()->GetBinCenter(indexY) << "  ";
		      cout << "posz " << mRatio->GetZaxis()->GetBinCenter(indexZ) << "  ";
		      cout << "content " << cellContent << "  " ;
		      cout << "fit " <<  cellContentFromFit << "  " ;
		      cout << "error " << cellContentError  << "   " ;
		      cout << endl ;
		    }
		}
	    }
	}
    }
  
  // Normalize fit and projection histos 
  m1dpro->Divide(m1dweight) ;
  m1dfit->Divide(m1dweight) ;
 
  // Now attach an error bar to each data point, with error = number of addens / (sum over 1 over error) ;
  for(Int_t binIndex=1 ; binIndex <= m1dpro->GetNbinsX() ;  binIndex++)
    {
      if( m1dcounter->GetBinContent(binIndex) > 0.0 && m1dweight->GetBinContent(binIndex) != 0.0 )
	{
	  // Use sqrt since the errors are squard 
	  m1dpro->SetBinError(binIndex,sqrt(m1dcounter->GetBinContent(binIndex))/m1dweight->GetBinContent(binIndex)) ;
	 
	}
    }
  
  // Return projection
  return m1dpro ;
}
//___________________
TH2D* T49FitProjector::get2dProjection(TString axis, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax)
{
  // The projection is done by summing over entries in direction "axis" specified when the function is called
  // The entries are weighted by 1 over the statistical error

  // Get ranges of the ratio 3d histogram 
  Int_t xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
  Int_t xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;
  Int_t yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
  Int_t ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
  Int_t zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
  Int_t zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;
    
  // Get dimensions of the projection histos
  Int_t proNbinsX = 0 ;
  Double_t proMinX = 0.0 ;
  Double_t proMaxX = 0.0 ;
  Char_t* axistitleX = new char[100] ;
  Int_t proNbinsY = 0 ;
  Double_t proMinY = 0.0 ;
  Double_t proMaxY = 0.0 ;
  Char_t* axistitleY = new char[100] ;

  // Determine the dimension of the 2d projection histogram
  if (axis == "x")        	
    {  
      proNbinsX = abs(ymaxBin-yminBin)+1 ;  
      proMinX   = mRatio->GetYaxis()->GetBinLowEdge(yminBin) ; 
      proMaxX   = mRatio->GetYaxis()->GetBinUpEdge(ymaxBin) ;
      strcpy(axistitleX,mRatio->GetYaxis()->GetTitle()) ; 
      proNbinsY = abs(zmaxBin-zminBin)+1 ;  
      proMinY   = mRatio->GetZaxis()->GetBinLowEdge(zminBin) ; 
      proMaxY   = mRatio->GetZaxis()->GetBinUpEdge(zmaxBin) ; 
      strcpy(axistitleY,mRatio->GetZaxis()->GetTitle()) ; 
    }
  else if (axis == "y")   
    {  
      proNbinsX =  abs(zmaxBin-zminBin)+1 ;  
      proMinX   =  mRatio->GetZaxis()->GetBinLowEdge(zminBin) ;
      proMaxX   =  mRatio->GetZaxis()->GetBinUpEdge(zmaxBin) ; 
      strcpy(axistitleX,mRatio->GetZaxis()->GetTitle()) ; 
      proNbinsY =  abs(xmaxBin-xminBin)+1 ;  
      proMinY   =  mRatio->GetXaxis()->GetBinLowEdge(xminBin) ; 
      proMaxY   =  mRatio->GetXaxis()->GetBinUpEdge(xmaxBin) ; 
      strcpy(axistitleY,mRatio->GetXaxis()->GetTitle()) ;
    }
  else if (axis == "z")   
    {  
      proNbinsX = abs(xmaxBin-xminBin)+1 ;  
      proMinX   = mRatio->GetXaxis()->GetBinLowEdge(xminBin) ; 
      proMaxX   = mRatio->GetXaxis()->GetBinUpEdge(xmaxBin) ; 
      strcpy(axistitleX,mRatio->GetXaxis()->GetTitle()) ; 
      proNbinsY = abs(ymaxBin-yminBin)+1 ;  
      proMinY   = mRatio->GetYaxis()->GetBinLowEdge(yminBin) ; 
      proMaxY   = mRatio->GetYaxis()->GetBinUpEdge(ymaxBin) ;
      strcpy(axistitleY,mRatio->GetYaxis()->GetTitle()) ;
    }
  else 
    { 
      cerr << "Wrong axis in projector, this should not happen ! "<< endl ; 
      return 0; 
    } ;
    
  //
  // Create 2d histos
  //
  Int_t ran = (Int_t) rand() ; // to distinguish between objects give them a special random id ...

  // Projection
  m2dpro = new TH2D("my2dpro","my2dpro",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ; 
  Char_t resName[100] ; sprintf(resName,"res%s%dn",axis.Data(),ran) ;  m2dpro->SetName(resName) ;
  m2dpro->SetXTitle(axistitleX) ;
  m2dpro->SetYTitle(axistitleY) ;
  Char_t restitle[100] ; sprintf(restitle,"%sn",axis.Data()) ;m2dpro->SetTitle(restitle) ;

  // Fit
  m2dfit = new TH2D("my2dfit","my2dfit",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
  Char_t fittitle[100] ; sprintf(fittitle,"fit%s%dn",axis.Data(),ran) ; m2dfit->SetName(fittitle) ;

  // Norm
  m2dweight = new TH2D("my2dweight","my2dweight",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
  char weightmtitle[100] ; sprintf(weightmtitle,"weightm%s%dn",axis.Data(),ran) ; m2dweight->SetName(weightmtitle) ;
    
  // Fill histo
  for(Int_t indexX = xminBin ; indexX <= xmaxBin ; indexX++)
    {
      for(Int_t indexY = yminBin ; indexY <= ymaxBin ; indexY++)
	{
	  for(Int_t indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
	    {
	      // Get the content of the 3d cell
	      Double_t cellcontent = mRatio->GetBinContent(indexX,indexY,indexZ) ;

	      // Get error of this cell
	      Double_t cellContentError = mRatio->GetBinError(indexX,indexY,indexZ) ;
			  
	      // Get the content of the 3d cell using the fit function
	      TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
				mRatio->GetYaxis()->GetBinCenter(indexY),
				mRatio->GetZaxis()->GetBinCenter(indexZ));

	      Double_t cellContentFromFit = mFitter->TheoreticalValue(&position,mFitParameters) ;
	      
	      // Fill the 2d histos
	      Int_t proBin = -100 ;
			  
	      if (axis == "x")        
		{ 
		  proBin = m2dpro->FindBin(mRatio->GetYaxis()->GetBinCenter(indexY),mRatio->GetZaxis()->GetBinCenter(indexZ)) ; 
		}  
	      else if (axis == "y")   
		{ 
		  proBin = m2dpro->FindBin(mRatio->GetZaxis()->GetBinCenter(indexZ),mRatio->GetXaxis()->GetBinCenter(indexX)) ; 
		}
	      else if (axis == "z")   
		{ 
		  proBin = m2dpro->FindBin(mRatio->GetXaxis()->GetBinCenter(indexX),mRatio->GetYaxis()->GetBinCenter(indexY)) ; 
		}
	      
	      if( cellcontent>0.5 && cellcontent <1.5 && cellContentError > 0.0 )
		{
		  // Fill projection 2d histo 
		  m2dpro->AddBinContent(proBin,cellcontent/cellContentError) ;
				    
		  // Fill fit 2d histo
		  m2dfit->AddBinContent(proBin,cellContentFromFit/cellContentError) ;

		  // Fill 1/error histogram to normalize
		  m2dweight->AddBinContent(proBin,1.0/cellContentError) ;
		}
	    }
	}
    }

  // normalize fit and projection histos
  m2dpro->Divide(m2dweight) ;
  m2dfit->Divide(m2dweight) ;
    
  // return projection
  return m2dpro ;    
}
//___________________
 void T49HbtFitBase::doFit()
{
  // everything is set up: start minimization !

  // Errorflag
  Int_t ierflg = 0 ;
  
  // check whether the TMinuit-object was initialized or not
  if (mMinuit->GetNumPars() == 0)
    {
      cout << "TMinuit not intitialized => calling standard initialization." << endl;
      InitMinuit();
    }

  // if normalization is derived elsewhere we simply fix it here
  if(mfixNormalization)
    {
      mMinuit->mnparm(0, "Normalization",  1.0, 0.01,  0,0,ierflg) ;
      mMinuit->FixParameter(0) ;
    }
  
  // set function to call
  mMinuit->SetObjectFit(this) ;
  mMinuit->SetFCN(gfcn) ;
  mMinuit->Migrad() ;
}


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.