[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
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():
bc7a3220	49	fIsParamFixed(49),
ba15fdfb	50	fPrintStatus(kFALSE),
f1ba54ec	51	fParamStartValues(),
ba15fdfb	52	fUp(0),
	53	fX(0x0),
	54	fM(0x0),
bc7a3220	55	fPt(0x0),
5720c765	56	fType(0x0),
	57	fLikely(0x0),
	58	fNcand(0),
	59	fContPlot1(0x0),
	60	fContPlot2(0x0),
	61	fContPlot3(0x0),
	62	fitter(0)
ba15fdfb	63	{
	64	//
	65	// default constructor
	66	//
	67	}
	68	//_________________________________________________________________________________________________
	69	AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime,
bc7a3220	70	Double_t* invariantmass, Double_t pt,Int_t type, Int_t ncand) :
bc7a3220	71	fIsParamFixed(49),
ba15fdfb	72	fPrintStatus(kFALSE),
f1ba54ec	73	fParamStartValues(),
ba15fdfb	74	fUp(0),
	75	fX(decaytime),
	76	fM(invariantmass),
bc7a3220	77	fPt(pt),
5720c765	78	fType(type),
ba15fdfb	79	fLikely(0x0),
5720c765	80	fNcand(ncand),
	81	fContPlot1(0x0),
	82	fContPlot2(0x0),
	83	fContPlot3(0x0),
	84	fitter(0)
ba15fdfb	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");
5720c765	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)");
bc7a3220	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)");
ba15fdfb	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;
f1ba54ec	153	for (Int_t i=0; i<49; ++i) fParamStartValues[i]=c.fParamStartValues[i];
ba15fdfb	154	fUp = c.fUp;
	155	fX = c.fX;
	156	fM = c.fM;
bc7a3220	157	fPt = c.fPt;
5720c765	158	fType = c.fType;
5720c765	159	fLikely = c.fLikely;
ba15fdfb	160	fNcand = c.fNcand;
5720c765	161	fContPlot1 = c.fContPlot1;
	162	fContPlot2 = c.fContPlot2;
	163	fContPlot3 = c.fContPlot3;
ba15fdfb	164	}
	165	return *this;
	166	}
	167
	168	//_______________________________________________________________________________________
	169	AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(const AliDielectronBtoJPSItoEleCDFfitHandler& c) :
	170	TNamed(c),
	171	fIsParamFixed(c.fIsParamFixed),
	172	fPrintStatus(c.fPrintStatus),
f1ba54ec	173	fParamStartValues(),
ba15fdfb	174	fUp(c.fUp),
	175	fX(c.fX),
	176	fM(c.fM),
bc7a3220	177	fPt(c.fPt),
5720c765	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)
ba15fdfb	185	{
	186	//
	187	// Copy Constructor
	188	//
f1ba54ec	189	for (Int_t i=0; i<49; ++i) fParamStartValues[i]=c.fParamStartValues[i];
ba15fdfb	190	}
	191	//_______________________________________________________________________________________
	192	AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
	193	{
	194	//
	195	//destructor
	196	//
	197	delete fLikely;
	198	}
	199	//_______________________________________________________________________________________
5720c765	200	Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
ba15fdfb	201	{
	202	//
	203	// performs the minimization
	204	//
bc7a3220	205	if(step == 0 \|\| !fitter){
bc7a3220	206	fitter = (TFitter*)TVirtualFitter::Fitter(this,49);
5720c765	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);
bc7a3220	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);
5720c765	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);
bc7a3220	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);
5720c765	258	}
ba15fdfb	259
bc7a3220	260	for(UInt_t indexparam = 0; indexparam < 49; indexparam++){
5720c765	261	if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam);
	262	else fitter->ReleaseParameter((Int_t)indexparam);
	263	}
bc7a3220	264	Double_t arglist[2]={10000,0.1}; Int_t iret = 0;
5720c765	265	if(step == 2) {Int_t iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
bc7a3220	266	if(step == 0) { fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.);
bc7a3220	267	iret=fitter->ExecuteCommand("MIGRAD", arglist ,2); }
ba15fdfb	268	fitter->PrintResults(4,0);
ba15fdfb	269
5720c765	270	if(step == 3) {
	271
	272	TMinuit* minuitPoint = fitter->GetMinuit();
ba15fdfb	273
5720c765	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
bc7a3220	306	c2->Draw();
bc7a3220	307	c2->SaveAs("contourPlot.root");
5720c765	308	}
	309
	310	AliInfo("Minimization procedure finished\n");
ba15fdfb	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
bc7a3220	327	f = fLikely->EvaluateLikelihood(fX,fM,fPt, fType,fNcand);
ba15fdfb	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	//_______________________________________________________________________________________
bc7a3220	337	void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[49])
ba15fdfb	338	{
bc7a3220	339	for(Int_t index=0; index < 49; index++) fParamStartValues[index] = inputparamvalues[index];
ba15fdfb	340	}
ba15fdfb	341	//_______________________________________________________________________________________
5720c765	342	void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
ba15fdfb	343	{
	344	//
	345	// Sets constants for the resolution function
	346	//
5720c765	347	fLikely->SetResolutionConstants(resolutionConst,type);
ba15fdfb	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	}
bc7a3220	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
ba15fdfb	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