PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitHandler.cxx

   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