[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(49),
	fPrintStatus(kFALSE),
        fParamStartValues(),
	fUp(0),
	fX(0x0),
	fM(0x0),
	fPt(0x0),
	fType(0x0),
        fLikely(0x0),
	fNcand(0),
	fContPlot1(0x0),
	fContPlot2(0x0),
	fContPlot3(0x0),
	fitter(0)
{
	//
	// default constructor
	//
}
//_________________________________________________________________________________________________
AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, 
		Double_t* invariantmass, Double_t *pt,Int_t *type, Int_t ncand) :
	fIsParamFixed(49),
	fPrintStatus(kFALSE),
        fParamStartValues(),
	fUp(0),
	fX(decaytime),
	fM(invariantmass),
	fPt(pt),
        fType(type),
	fLikely(0x0),
	fNcand(ncand),
	fContPlot1(0x0),
	fContPlot2(0x0),
	fContPlot3(0x0),
	fitter(0)
{
	//
	// constructor
	//
	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("Parameter 45 ----> fOneOvLamSym1 (additional parameter)");
        AliInfo("Parameter 46 ----> fSym1 (additional parameter)");
        AliInfo("Parameter 47 ----> fPolyn4 (additional parameter)");
        AliInfo("Parameter 48 ----> fPolyn5 (additional parameter)");

	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;
                for (Int_t i=0; i<49; ++i) fParamStartValues[i]=c.fParamStartValues[i];
		fUp           = c.fUp;
		fX            = c.fX;
		fM            = c.fM;
		fPt           = c.fPt;
		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),
        fParamStartValues(),
	fUp(c.fUp),
	fX(c.fX),
	fM(c.fM),
	fPt(c.fPt),
	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<49; ++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){ 
		fitter = (TFitter*)TVirtualFitter::Fitter(this,49);
		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, -1.e+04, 1.e+04);
		fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, -1.e+04, 1.e+04);
		fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, -1.e+04, 1.e+04);
		fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, -1.e+04, 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);
                fitter->SetParameter(45,"fOneOvLamSym1",fParamStartValues[45], 1.e-10, 0.0000001, 5.e+01);
                fitter->SetParameter(46,"fSym1",fParamStartValues[46], 1.e-08, 0., 1.e+06);
                fitter->SetParameter(47,"fBkgInvMassPol4",fParamStartValues[47], 1.e-08, 0.0000001, 5.e+01);
                fitter->SetParameter(48,"fBkgInvMassPol5",fParamStartValues[48], 1.e-08, 0., 5.e+06);
                //(TMinuit*)fitter->GetMinuit()->SetErrorDef(fUp); 
                }

	for(UInt_t indexparam = 0; indexparam < 49; indexparam++){
		if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam); 
		else fitter->ReleaseParameter((Int_t)indexparam);
	}
	Double_t arglist[2]={10000,0.1}; Int_t iret = 0;
	if(step == 2) {Int_t  iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
	if(step == 0) { fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.); 
                         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();
                c2->SaveAs("contourPlot.root");	
	}

	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,fPt, 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[49])
{
	for(Int_t index=0; index < 49; 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::SetExponentialFunction(Bool_t okExp)
{
        //
        // Sets the CB as the parametrization for the signal invariant mass spectrum 
        // (otherwise Landau is chosen)
        //
        fLikely->SetExponentialFunction(okExp);
}


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