[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronBtoJPSItoEleCDFfitHandler.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
#include <TVirtualFitter.h>
#include <TFitter.h>
#include <TMinuit.h>
#include <TStopwatch.h>
#include <TCanvas.h>
#include "AliLog.h"
#include "AliDielectronBtoJPSItoEleCDFfitHandler.h"
#include "AliDielectronBtoJPSItoEleCDFfitFCN.h"

//-------------------------------------------------------------------------
//                      Class AliDielectronBtoJPSItoEleCDFfitHandler
//            Class to perform unbinned log-likelihood fit
//      
//                        Origin: C. Di Giglio
//     Contact: carmelo.digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it
//-------------------------------------------------------------------------

void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);

ClassImp(AliDielectronBtoJPSItoEleCDFfitHandler)

	//______________________________________________________________________________
void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
	// This function is called by minuit
	// The corresponding member method is called
	// using SetObjectFit/GetObjectFit methods of TMinuit
	AliDielectronBtoJPSItoEleCDFfitHandler* dummy = (AliDielectronBtoJPSItoEleCDFfitHandler *)TVirtualFitter::GetFitter()->GetObjectFit();
	dummy->CdfFCN(npar, gin, f, par, iflag);
}


//_________________________________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler():
	fIsParamFixed(45),
	fPrintStatus(kFALSE),
	fUp(0),
	fX(0x0),
	fM(0x0),
	fType(0x0),
        fLikely(0x0),
	fNcand(0),
	fContPlot1(0x0),
	fContPlot2(0x0),
	fContPlot3(0x0),
	fitter(0)
{
	//
	// default constructor
	//
	for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
}
//_________________________________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, 
		Double_t* invariantmass, Int_t *type, Int_t ncand) :
	fIsParamFixed(45),
	fPrintStatus(kFALSE),
	fUp(0),
	fX(decaytime),
	fM(invariantmass),
        fType(type),
	fLikely(0x0),
	fNcand(ncand),
	fContPlot1(0x0),
	fContPlot2(0x0),
	fContPlot3(0x0),
	fitter(0)
{
	//
	// constructor
	//
        for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
	AliInfo("\n+++\n+++ Minimization object AliDielectronBtoJPSItoEleCDFfitHandler created\n+++\n");
	fLikely = new AliDielectronBtoJPSItoEleCDFfitFCN();
	AliInfo("\n+++\n+++ CDF fit function object AliDielectronBtoJPSItoEleCDFfitFCN created\n+++\n");
	AliInfo("Parameter 0  ----> fWeightRes");
	AliInfo("Parameter 1  ----> fPos");
	AliInfo("Parameter 2  ----> fNeg");
	AliInfo("Parameter 3  ----> fSym");
	AliInfo("Parameter 4  ----> fOneOvLamPlus");
	AliInfo("Parameter 5  ----> fOneOvLamMinus");
	AliInfo("Parameter 6  ----> fOneOvLamSym");
	AliInfo("Parameter 7  ----> fB");
	AliInfo("Parameter 8  ----> fFsig");
	AliInfo("Parameter 9  ----> fMmean");
	AliInfo("Parameter 10 ----> fNexp");
	AliInfo("Parameter 11 ----> fSigma");
	AliInfo("Parameter 12 ----> fAlpha");
	AliInfo("Parameter 13 ----> fNorm");
	AliInfo("Parameter 14 ----> fBkgNorm");
	AliInfo("Parameter 15 ----> fBkgMean");
	AliInfo("Parameter 16 ----> fBkgSlope");
	AliInfo("Parameter 17 ----> fBkgConst");
	AliInfo("Parameter 18 ----> fNormGaus1 (FF)");
	AliInfo("Parameter 19 ----> fNormGaus2 (FF)");
        AliInfo("Parameter 20 ----> fMean1Res (FF)");
        AliInfo("Parameter 21 ----> fsigma1Res (FF)"); 
        AliInfo("Parameter 22 ----> fMean2Res (FF)");
        AliInfo("Parameter 23 ----> fsigma2Res (FF)");
        AliInfo("Parameter 24 ----> fAlfaRes (FF)");
        AliInfo("Parameter 25 ----> fLambdaRes (FF)");
        AliInfo("Parameter 26 ----> fNormResExp (FF)");
        AliInfo("Parameter 27 ----> fNormGaus1 (FS)");
        AliInfo("Parameter 28 ----> fNormGaus2 (FS)");
        AliInfo("Parameter 29 ----> fMean1Res (FS)");
        AliInfo("Parameter 30 ----> fsigma1Res (FS)");
        AliInfo("Parameter 31 ----> fMean2Res (FS)");
        AliInfo("Parameter 32 ----> fsigma2Res (FS)");
        AliInfo("Parameter 33 ----> fAlfaRes (FS)");
        AliInfo("Parameter 34 ----> fLambdaRes (FS)");
        AliInfo("Parameter 35 ----> fNormResExp (FS)");
        AliInfo("Parameter 36 ----> fNormGaus1 (SS)");
        AliInfo("Parameter 37 ----> fNormGaus2 (SS)");
        AliInfo("Parameter 38 ----> fMean1Res (SS)");
        AliInfo("Parameter 39 ----> fsigma1Res (SS)");
        AliInfo("Parameter 40 ----> fMean2Res (SS)");
        AliInfo("Parameter 41 ----> fsigma2Res (SS)");
        AliInfo("Parameter 42 ----> fAlfaRes (SS)");
        AliInfo("Parameter 43 ----> fLambdaRes (SS)");
        AliInfo("Parameter 44 ----> fNormResExp (SS)");

	AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
}
//___________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler& AliDielectronBtoJPSItoEleCDFfitHandler::operator=(const AliDielectronBtoJPSItoEleCDFfitHandler& c)
{
	//
	// Assignment operator
	//
	if (this!=&c) {
		fIsParamFixed = c.fIsParamFixed;
		fPrintStatus  = c.fPrintStatus;
		fUp           = c.fUp;
		fX            = c.fX;
		fM            = c.fM;
		fType         = c.fType;
                fLikely       = c.fLikely;
		fNcand        = c.fNcand;
		fContPlot1    = c.fContPlot1;
		fContPlot2    = c.fContPlot2;
		fContPlot3    = c.fContPlot3;
	}
	return *this;
}

//_______________________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(const AliDielectronBtoJPSItoEleCDFfitHandler& c) :
	TNamed(c),
	fIsParamFixed(c.fIsParamFixed),
	fPrintStatus(c.fPrintStatus),
	fUp(c.fUp),
	fX(c.fX),
	fM(c.fM),
	fType(c.fType),
        fLikely(c.fLikely),
	fNcand(c.fNcand),
	fContPlot1(c.fContPlot1),
	fContPlot2(c.fContPlot2),
	fContPlot3(c.fContPlot3),
	fitter(c.fitter)
{
	//
	// Copy Constructor
	//
        for (Int_t i=0; i<45; ++i) fParamStartValues[i]=c.fParamStartValues[i];
}
//_______________________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
{
	//
	//destructor
	//
	delete fLikely;
}
//_______________________________________________________________________________________
Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
{
	//
	// performs the minimization
	//
	if(step == 0){ 
		//fitter = TVirtualFitter::Fitter(this,20);
		fitter = (TFitter*)TVirtualFitter::Fitter(this,45);
		fitter->SetFCN(CDFFunction);
		fitter->SetParameter(0,"fWeightRes",fParamStartValues[0], 1.e-08, 0., 1.e+06);
		fitter->SetParameter(1,"fPos",fParamStartValues[1], 1.e-08, 0.,1.e+06);
		fitter->SetParameter(2,"fNeg",fParamStartValues[2], 1.e-08, 0.,1.e+06);
		fitter->SetParameter(3,"fSym",fParamStartValues[3], 1.e-08, 0.,1.e+06);
		fitter->SetParameter(4,"fOneOvLamPlus",fParamStartValues[4], 1.e-10, 0.0000001, 5.e+01);
		fitter->SetParameter(5,"fOneOvLamMinus",fParamStartValues[5], 1.e-10, 0.00000001, 5.e+01);
		fitter->SetParameter(6,"fOneOvLamSym",fParamStartValues[6], 1.e-10, 0.00000001, 5.e+01);
		fitter->SetParameter(7,"fB",fParamStartValues[7], 1.e-10, 0., 1.);
		fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.);
		fitter->SetParameter(9,"fMmean",fParamStartValues[9], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(10,"fNexp",fParamStartValues[10], 1.e-08, 0., 1.e+02);
		fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(14,"fBkgNorm",fParamStartValues[14], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, 0., 1.e+04);
		fitter->SetParameter(18,"fNormGaus1FF",fParamStartValues[18], 1.e-08, 0., 1.e+05);
		fitter->SetParameter(19,"fNormGaus2FF",fParamStartValues[19], 1.e-08, 0., 1.e+05); 
	        fitter->SetParameter(20,"fMean1ResFF",fParamStartValues[20], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(21,"fSigma1ResFF",fParamStartValues[21], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(22,"fMean2ResFF",fParamStartValues[22], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(23,"fSigma2ResFF",fParamStartValues[23], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(24,"fAlfaResFF",fParamStartValues[24], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(25,"fLambdaResFF",fParamStartValues[25], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(26,"fResNormExpFF",fParamStartValues[26], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(27,"fNormGaus1FS",fParamStartValues[27], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(28,"fNormGaus2FS",fParamStartValues[28], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(29,"fMean1ResFS",fParamStartValues[29], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(30,"fSigma1ResFS",fParamStartValues[30], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(31,"fMean2ResFS",fParamStartValues[31], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(32,"fSigma2ResFS",fParamStartValues[32], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(33,"fAlfaResFS",fParamStartValues[33], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(34,"fLambdaResFS",fParamStartValues[34], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(35,"fResNormExpFS",fParamStartValues[35], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(36,"fNormGaus1SS",fParamStartValues[36], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(37,"fNormGaus2SS",fParamStartValues[37], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(38,"fMean1ResSS",fParamStartValues[38], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(39,"fSigma1ResSS",fParamStartValues[39], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(40,"fMean2ResSS",fParamStartValues[40], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(41,"fSigma2ResSS",fParamStartValues[41], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(42,"fAlfaResSS",fParamStartValues[42], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(43,"fLambdaResSS",fParamStartValues[43], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(44,"fResNormExpSS",fParamStartValues[44], 1.e-08, 0., 1.e+05);
                }

	for(UInt_t indexparam = 0; indexparam < 45; indexparam++){
		if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam); 
		else fitter->ReleaseParameter((Int_t)indexparam);
	}
	Double_t arglist[2]={10000,0.1};
	if(step == 2) {Int_t  iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
	Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
	fitter->PrintResults(4,0);

	if(step == 3) {

		TMinuit* minuitPoint = fitter->GetMinuit(); 

		TCanvas *c2 = new TCanvas("c2","contours",800,800);

		//68.27% (1 sigma) confidence level for 2 parameters (fSIG versus fB)   
		minuitPoint->SetErrorDef(1.15); // 2.3/2
		fContPlot1 = (TGraph*)minuitPoint->Contour(100,7,8);
		fContPlot1->GetXaxis()->SetRange(0,1);
                fContPlot1->GetYaxis()->SetRange(0,1);
                fContPlot1->SetLineColor(42);
		fContPlot1->SetLineWidth(3);

		//95% (2 sigma) confidence level for 2 parameters (fSIG versus fB)   
		minuitPoint->SetErrorDef(2.995); // 5.99/2
		fContPlot2 = (TGraph*)minuitPoint->Contour(100,7,8);
                fContPlot2->GetXaxis()->SetRange(0,1);
                fContPlot2->GetYaxis()->SetRange(0,1);
		fContPlot2->SetLineColor(38);
		fContPlot2->SetLineWidth(3);

		//99.73% (3 sigma) confidence level for 2 parameters (fSIG versus fB)   
		minuitPoint->SetErrorDef(5.915); // 11.83/2
		fContPlot3 = (TGraph*)minuitPoint->Contour(100,7,8);
                fContPlot3->GetXaxis()->SetRange(0,1);
                fContPlot3->GetXaxis()->SetTitle("f_{B}");
                fContPlot3->GetYaxis()->SetTitle("f_{Sig}[2.4-4]");
                fContPlot3->GetYaxis()->SetRange(0,1);
  		fContPlot3->SetLineColor(34);
		fContPlot3->SetLineWidth(3);
          
                fContPlot3->Draw("al");
                fContPlot2->Draw("l");
                fContPlot1->Draw("l");
		
                c2->Draw();	
	}

	AliInfo("Minimization procedure finished\n");
	return iret;
}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */, 
		Double_t * /* gin */,Double_t &f,Double_t *par,Int_t /* iflag */)
{
	// 
	// Definition of the FCN to be used by minuit
	//
	fLikely->SetAllParameters(par);
	fLikely->ComputeMassIntegral();
	if(fPrintStatus)fLikely->PrintStatus();

	TStopwatch t;
	t.Start();

	f = fLikely->EvaluateLikelihood(fX,fM,fType,fNcand);

	t.Stop();
	AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime()));
	AliDebug(2,Form("CPU time spent to calculate function == %f \n", t.CpuTime()));
	AliDebug(2,Form("Actual value of the AliDielectronBtoJPSItoEleCDFfitFCN function == %f \n", f));

	return;
}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[45])
{
	for(Int_t index=0; index < 45; index++) fParamStartValues[index] = inputparamvalues[index];
}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
{
	//
	// Sets constants for the resolution function
	//
	fLikely->SetResolutionConstants(resolutionConst,type);

}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB) 
{
	//
	// Sets the CB as the parametrization for the signal invariant mass spectrum 
	// (otherwise Landau is chosen)
	//
	fLikely->SetCrystalBallFunction(okCB);
}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit) 
{ 
	//
	// Sets upper limit for the invariant mass window (under J/PSI mass peak)
	//
	fLikely->SetMassWndHigh(limit);
}
//_______________________________________________________________________________________
void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndLow(Double_t limit)
{
	//
	// Sets lower limit for the invariant mass window (under J/PSI mass peak)
	//
	fLikely->SetMassWndLow(limit);
}
Commit	Line	Data
ba15fdfb	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	#include <TVirtualFitter.h>
5720c765	16	#include <TFitter.h>
5720c765	17	#include <TMinuit.h>
ba15fdfb	18	#include <TStopwatch.h>
5720c765	19	#include <TCanvas.h>
ba15fdfb	20	#include "AliLog.h"
	21	#include "AliDielectronBtoJPSItoEleCDFfitHandler.h"
	22	#include "AliDielectronBtoJPSItoEleCDFfitFCN.h"
	23
	24	//-------------------------------------------------------------------------
	25	// Class AliDielectronBtoJPSItoEleCDFfitHandler
	26	// Class to perform unbinned log-likelihood fit
	27	//
	28	// Origin: C. Di Giglio
	29	// Contact: carmelo.digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it
	30	//-------------------------------------------------------------------------
	31
	32	void CDFFunction(Int_t &npar, Double_t gin, Double_t &f, Double_t x, Int_t iflag);
	33
	34	ClassImp(AliDielectronBtoJPSItoEleCDFfitHandler)
	35
	36	//______________________________________________________________________________
	37	void CDFFunction(Int_t &npar, Double_t gin, Double_t &f, Double_t par, Int_t iflag)
	38	{
	39	// This function is called by minuit
	40	// The corresponding member method is called
	41	// using SetObjectFit/GetObjectFit methods of TMinuit
	42	AliDielectronBtoJPSItoEleCDFfitHandler* dummy = (AliDielectronBtoJPSItoEleCDFfitHandler *)TVirtualFitter::GetFitter()->GetObjectFit();
	43	dummy->CdfFCN(npar, gin, f, par, iflag);
	44	}
	45
	46
	47	//_________________________________________________________________________________________________
	48	AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler():
5720c765	49	fIsParamFixed(45),
ba15fdfb	50	fPrintStatus(kFALSE),
	51	fUp(0),
	52	fX(0x0),
	53	fM(0x0),
5720c765	54	fType(0x0),
	55	fLikely(0x0),
	56	fNcand(0),
	57	fContPlot1(0x0),
	58	fContPlot2(0x0),
	59	fContPlot3(0x0),
	60	fitter(0)
ba15fdfb	61	{
	62	//
	63	// default constructor
	64	//
ac390e40	65	for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
ba15fdfb	66	}
	67	//_________________________________________________________________________________________________
	68	AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime,
5720c765	69	Double_t* invariantmass, Int_t *type, Int_t ncand) :
5720c765	70	fIsParamFixed(45),
ba15fdfb	71	fPrintStatus(kFALSE),
	72	fUp(0),
	73	fX(decaytime),
	74	fM(invariantmass),
5720c765	75	fType(type),
ba15fdfb	76	fLikely(0x0),
5720c765	77	fNcand(ncand),
	78	fContPlot1(0x0),
	79	fContPlot2(0x0),
	80	fContPlot3(0x0),
	81	fitter(0)
ba15fdfb	82	{
	83	//
	84	// constructor
	85	//
ac390e40	86	for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
ba15fdfb	87	AliInfo("\n+++\n+++ Minimization object AliDielectronBtoJPSItoEleCDFfitHandler created\n+++\n");
	88	fLikely = new AliDielectronBtoJPSItoEleCDFfitFCN();
	89	AliInfo("\n+++\n+++ CDF fit function object AliDielectronBtoJPSItoEleCDFfitFCN created\n+++\n");
	90	AliInfo("Parameter 0 ----> fWeightRes");
	91	AliInfo("Parameter 1 ----> fPos");
	92	AliInfo("Parameter 2 ----> fNeg");
	93	AliInfo("Parameter 3 ----> fSym");
	94	AliInfo("Parameter 4 ----> fOneOvLamPlus");
	95	AliInfo("Parameter 5 ----> fOneOvLamMinus");
	96	AliInfo("Parameter 6 ----> fOneOvLamSym");
	97	AliInfo("Parameter 7 ----> fB");
	98	AliInfo("Parameter 8 ----> fFsig");
	99	AliInfo("Parameter 9 ----> fMmean");
	100	AliInfo("Parameter 10 ----> fNexp");
	101	AliInfo("Parameter 11 ----> fSigma");
	102	AliInfo("Parameter 12 ----> fAlpha");
	103	AliInfo("Parameter 13 ----> fNorm");
	104	AliInfo("Parameter 14 ----> fBkgNorm");
	105	AliInfo("Parameter 15 ----> fBkgMean");
	106	AliInfo("Parameter 16 ----> fBkgSlope");
	107	AliInfo("Parameter 17 ----> fBkgConst");
5720c765	108	AliInfo("Parameter 18 ----> fNormGaus1 (FF)");
	109	AliInfo("Parameter 19 ----> fNormGaus2 (FF)");
	110	AliInfo("Parameter 20 ----> fMean1Res (FF)");
	111	AliInfo("Parameter 21 ----> fsigma1Res (FF)");
	112	AliInfo("Parameter 22 ----> fMean2Res (FF)");
	113	AliInfo("Parameter 23 ----> fsigma2Res (FF)");
	114	AliInfo("Parameter 24 ----> fAlfaRes (FF)");
	115	AliInfo("Parameter 25 ----> fLambdaRes (FF)");
	116	AliInfo("Parameter 26 ----> fNormResExp (FF)");
	117	AliInfo("Parameter 27 ----> fNormGaus1 (FS)");
	118	AliInfo("Parameter 28 ----> fNormGaus2 (FS)");
	119	AliInfo("Parameter 29 ----> fMean1Res (FS)");
	120	AliInfo("Parameter 30 ----> fsigma1Res (FS)");
	121	AliInfo("Parameter 31 ----> fMean2Res (FS)");
	122	AliInfo("Parameter 32 ----> fsigma2Res (FS)");
	123	AliInfo("Parameter 33 ----> fAlfaRes (FS)");
	124	AliInfo("Parameter 34 ----> fLambdaRes (FS)");
	125	AliInfo("Parameter 35 ----> fNormResExp (FS)");
	126	AliInfo("Parameter 36 ----> fNormGaus1 (SS)");
	127	AliInfo("Parameter 37 ----> fNormGaus2 (SS)");
	128	AliInfo("Parameter 38 ----> fMean1Res (SS)");
	129	AliInfo("Parameter 39 ----> fsigma1Res (SS)");
	130	AliInfo("Parameter 40 ----> fMean2Res (SS)");
	131	AliInfo("Parameter 41 ----> fsigma2Res (SS)");
	132	AliInfo("Parameter 42 ----> fAlfaRes (SS)");
	133	AliInfo("Parameter 43 ----> fLambdaRes (SS)");
	134	AliInfo("Parameter 44 ----> fNormResExp (SS)");
ba15fdfb	135
	136	AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
	137	}
	138	//___________________________________________________________________________
	139	AliDielectronBtoJPSItoEleCDFfitHandler& AliDielectronBtoJPSItoEleCDFfitHandler::operator=(const AliDielectronBtoJPSItoEleCDFfitHandler& c)
	140	{
	141	//
	142	// Assignment operator
	143	//
	144	if (this!=&c) {
	145	fIsParamFixed = c.fIsParamFixed;
	146	fPrintStatus = c.fPrintStatus;
	147	fUp = c.fUp;
	148	fX = c.fX;
	149	fM = c.fM;
5720c765	150	fType = c.fType;
5720c765	151	fLikely = c.fLikely;
ba15fdfb	152	fNcand = c.fNcand;
5720c765	153	fContPlot1 = c.fContPlot1;
	154	fContPlot2 = c.fContPlot2;
	155	fContPlot3 = c.fContPlot3;
ba15fdfb	156	}
	157	return *this;
	158	}
	159
	160	//_______________________________________________________________________________________
	161	AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(const AliDielectronBtoJPSItoEleCDFfitHandler& c) :
	162	TNamed(c),
	163	fIsParamFixed(c.fIsParamFixed),
	164	fPrintStatus(c.fPrintStatus),
	165	fUp(c.fUp),
	166	fX(c.fX),
	167	fM(c.fM),
5720c765	168	fType(c.fType),
	169	fLikely(c.fLikely),
	170	fNcand(c.fNcand),
	171	fContPlot1(c.fContPlot1),
	172	fContPlot2(c.fContPlot2),
	173	fContPlot3(c.fContPlot3),
	174	fitter(c.fitter)
ba15fdfb	175	{
	176	//
	177	// Copy Constructor
	178	//
ac390e40	179	for (Int_t i=0; i<45; ++i) fParamStartValues[i]=c.fParamStartValues[i];
ba15fdfb	180	}
	181	//_______________________________________________________________________________________
	182	AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
	183	{
	184	//
	185	//destructor
	186	//
	187	delete fLikely;
	188	}
	189	//_______________________________________________________________________________________
5720c765	190	Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
ba15fdfb	191	{
	192	//
	193	// performs the minimization
	194	//
5720c765	195	if(step == 0){
	196	//fitter = TVirtualFitter::Fitter(this,20);
	197	fitter = (TFitter*)TVirtualFitter::Fitter(this,45);
	198	fitter->SetFCN(CDFFunction);
	199	fitter->SetParameter(0,"fWeightRes",fParamStartValues[0], 1.e-08, 0., 1.e+06);
	200	fitter->SetParameter(1,"fPos",fParamStartValues[1], 1.e-08, 0.,1.e+06);
	201	fitter->SetParameter(2,"fNeg",fParamStartValues[2], 1.e-08, 0.,1.e+06);
	202	fitter->SetParameter(3,"fSym",fParamStartValues[3], 1.e-08, 0.,1.e+06);
	203	fitter->SetParameter(4,"fOneOvLamPlus",fParamStartValues[4], 1.e-10, 0.0000001, 5.e+01);
	204	fitter->SetParameter(5,"fOneOvLamMinus",fParamStartValues[5], 1.e-10, 0.00000001, 5.e+01);
	205	fitter->SetParameter(6,"fOneOvLamSym",fParamStartValues[6], 1.e-10, 0.00000001, 5.e+01);
	206	fitter->SetParameter(7,"fB",fParamStartValues[7], 1.e-10, 0., 1.);
	207	fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.);
	208	fitter->SetParameter(9,"fMmean",fParamStartValues[9], 1.e-08, 0., 1.e+04);
	209	fitter->SetParameter(10,"fNexp",fParamStartValues[10], 1.e-08, 0., 1.e+02);
	210	fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04);
	211	fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04);
	212	fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+04);
	213	fitter->SetParameter(14,"fBkgNorm",fParamStartValues[14], 1.e-08, 0., 1.e+04);
	214	fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, 0., 1.e+04);
	215	fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, 0., 1.e+04);
	216	fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, 0., 1.e+04);
	217	fitter->SetParameter(18,"fNormGaus1FF",fParamStartValues[18], 1.e-08, 0., 1.e+05);
	218	fitter->SetParameter(19,"fNormGaus2FF",fParamStartValues[19], 1.e-08, 0., 1.e+05);
	219	fitter->SetParameter(20,"fMean1ResFF",fParamStartValues[20], 1.e-08, 0., 1.e+05);
	220	fitter->SetParameter(21,"fSigma1ResFF",fParamStartValues[21], 1.e-08, 0., 1.e+05);
	221	fitter->SetParameter(22,"fMean2ResFF",fParamStartValues[22], 1.e-08, 0., 1.e+05);
	222	fitter->SetParameter(23,"fSigma2ResFF",fParamStartValues[23], 1.e-08, 0., 1.e+05);
	223	fitter->SetParameter(24,"fAlfaResFF",fParamStartValues[24], 1.e-08, 0., 1.e+05);
	224	fitter->SetParameter(25,"fLambdaResFF",fParamStartValues[25], 1.e-08, 0., 1.e+05);
	225	fitter->SetParameter(26,"fResNormExpFF",fParamStartValues[26], 1.e-08, 0., 1.e+05);
	226	fitter->SetParameter(27,"fNormGaus1FS",fParamStartValues[27], 1.e-08, 0., 1.e+05);
	227	fitter->SetParameter(28,"fNormGaus2FS",fParamStartValues[28], 1.e-08, 0., 1.e+05);
	228	fitter->SetParameter(29,"fMean1ResFS",fParamStartValues[29], 1.e-08, 0., 1.e+05);
	229	fitter->SetParameter(30,"fSigma1ResFS",fParamStartValues[30], 1.e-08, 0., 1.e+05);
	230	fitter->SetParameter(31,"fMean2ResFS",fParamStartValues[31], 1.e-08, 0., 1.e+05);
	231	fitter->SetParameter(32,"fSigma2ResFS",fParamStartValues[32], 1.e-08, 0., 1.e+05);
	232	fitter->SetParameter(33,"fAlfaResFS",fParamStartValues[33], 1.e-08, 0., 1.e+05);
	233	fitter->SetParameter(34,"fLambdaResFS",fParamStartValues[34], 1.e-08, 0., 1.e+05);
	234	fitter->SetParameter(35,"fResNormExpFS",fParamStartValues[35], 1.e-08, 0., 1.e+05);
	235	fitter->SetParameter(36,"fNormGaus1SS",fParamStartValues[36], 1.e-08, 0., 1.e+05);
	236	fitter->SetParameter(37,"fNormGaus2SS",fParamStartValues[37], 1.e-08, 0., 1.e+05);
	237	fitter->SetParameter(38,"fMean1ResSS",fParamStartValues[38], 1.e-08, 0., 1.e+05);
	238	fitter->SetParameter(39,"fSigma1ResSS",fParamStartValues[39], 1.e-08, 0., 1.e+05);
	239	fitter->SetParameter(40,"fMean2ResSS",fParamStartValues[40], 1.e-08, 0., 1.e+05);
	240	fitter->SetParameter(41,"fSigma2ResSS",fParamStartValues[41], 1.e-08, 0., 1.e+05);
	241	fitter->SetParameter(42,"fAlfaResSS",fParamStartValues[42], 1.e-08, 0., 1.e+05);
	242	fitter->SetParameter(43,"fLambdaResSS",fParamStartValues[43], 1.e-08, 0., 1.e+05);
	243	fitter->SetParameter(44,"fResNormExpSS",fParamStartValues[44], 1.e-08, 0., 1.e+05);
	244	}
ba15fdfb	245
5720c765	246	for(UInt_t indexparam = 0; indexparam < 45; indexparam++){
	247	if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam);
	248	else fitter->ReleaseParameter((Int_t)indexparam);
	249	}
	250	Double_t arglist[2]={10000,0.1};
	251	if(step == 2) {Int_t iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
ba15fdfb	252	Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
	253	fitter->PrintResults(4,0);
	254
5720c765	255	if(step == 3) {
	256
	257	TMinuit* minuitPoint = fitter->GetMinuit();
ba15fdfb	258
5720c765	259	TCanvas *c2 = new TCanvas("c2","contours",800,800);
	260
	261	//68.27% (1 sigma) confidence level for 2 parameters (fSIG versus fB)
	262	minuitPoint->SetErrorDef(1.15); // 2.3/2
	263	fContPlot1 = (TGraph*)minuitPoint->Contour(100,7,8);
	264	fContPlot1->GetXaxis()->SetRange(0,1);
	265	fContPlot1->GetYaxis()->SetRange(0,1);
	266	fContPlot1->SetLineColor(42);
	267	fContPlot1->SetLineWidth(3);
	268
	269	//95% (2 sigma) confidence level for 2 parameters (fSIG versus fB)
	270	minuitPoint->SetErrorDef(2.995); // 5.99/2
	271	fContPlot2 = (TGraph*)minuitPoint->Contour(100,7,8);
	272	fContPlot2->GetXaxis()->SetRange(0,1);
	273	fContPlot2->GetYaxis()->SetRange(0,1);
	274	fContPlot2->SetLineColor(38);
	275	fContPlot2->SetLineWidth(3);
	276
	277	//99.73% (3 sigma) confidence level for 2 parameters (fSIG versus fB)
	278	minuitPoint->SetErrorDef(5.915); // 11.83/2
	279	fContPlot3 = (TGraph*)minuitPoint->Contour(100,7,8);
	280	fContPlot3->GetXaxis()->SetRange(0,1);
	281	fContPlot3->GetXaxis()->SetTitle("f_{B}");
	282	fContPlot3->GetYaxis()->SetTitle("f_{Sig}[2.4-4]");
	283	fContPlot3->GetYaxis()->SetRange(0,1);
	284	fContPlot3->SetLineColor(34);
	285	fContPlot3->SetLineWidth(3);
	286
	287	fContPlot3->Draw("al");
	288	fContPlot2->Draw("l");
	289	fContPlot1->Draw("l");
	290
	291	c2->Draw();
	292	}
	293
	294	AliInfo("Minimization procedure finished\n");
ba15fdfb	295	return iret;
	296	}
	297	//_______________________________________________________________________________________
	298	void AliDielectronBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */,
	299	Double_t * /* gin /,Double_t &f,Double_t par,Int_t /* iflag */)
	300	{
	301	//
	302	// Definition of the FCN to be used by minuit
	303	//
	304	fLikely->SetAllParameters(par);
	305	fLikely->ComputeMassIntegral();
	306	if(fPrintStatus)fLikely->PrintStatus();
	307
	308	TStopwatch t;
	309	t.Start();
	310
5720c765	311	f = fLikely->EvaluateLikelihood(fX,fM,fType,fNcand);
ba15fdfb	312
	313	t.Stop();
	314	AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime()));
	315	AliDebug(2,Form("CPU time spent to calculate function == %f \n", t.CpuTime()));
	316	AliDebug(2,Form("Actual value of the AliDielectronBtoJPSItoEleCDFfitFCN function == %f \n", f));
	317
	318	return;
	319	}
	320	//_______________________________________________________________________________________
5720c765	321	void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[45])
ba15fdfb	322	{
5720c765	323	for(Int_t index=0; index < 45; index++) fParamStartValues[index] = inputparamvalues[index];
ba15fdfb	324	}
ba15fdfb	325	//_______________________________________________________________________________________
5720c765	326	void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
ba15fdfb	327	{
	328	//
	329	// Sets constants for the resolution function
	330	//
5720c765	331	fLikely->SetResolutionConstants(resolutionConst,type);
ba15fdfb	332
	333	}
	334	//_______________________________________________________________________________________
	335	void AliDielectronBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB)
	336	{
	337	//
	338	// Sets the CB as the parametrization for the signal invariant mass spectrum
	339	// (otherwise Landau is chosen)
	340	//
	341	fLikely->SetCrystalBallFunction(okCB);
	342	}
	343	//_______________________________________________________________________________________
	344	void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit)
	345	{
	346	//
	347	// Sets upper limit for the invariant mass window (under J/PSI mass peak)
	348	//
	349	fLikely->SetMassWndHigh(limit);
	350	}
	351	//_______________________________________________________________________________________
	352	void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndLow(Double_t limit)
	353	{
	354	//
	355	// Sets lower limit for the invariant mass window (under J/PSI mass peak)
	356	//
	357	fLikely->SetMassWndLow(limit);
	358	}
	359