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(45),
  50         fPrintStatus(kFALSE),
  51         fUp(0),
  52         fX(0x0),
  53         fM(0x0),
  54         fType(0x0),
  55         fLikely(0x0),
  56         fNcand(0),
  57         fContPlot1(0x0),
  58         fContPlot2(0x0),
  59         fContPlot3(0x0),
  60         fitter(0)
  61 {
  62         //
  63         // default constructor
  64         //
  65         for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
  66 }
  67 //_________________________________________________________________________________________________
  68 AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime,
  69                 Double_t* invariantmass, Int_t *type, Int_t ncand) :
  70         fIsParamFixed(45),
  71         fPrintStatus(kFALSE),
  72         fUp(0),
  73         fX(decaytime),
  74         fM(invariantmass),
  75         fType(type),
  76         fLikely(0x0),
  77         fNcand(ncand),
  78         fContPlot1(0x0),
  79         fContPlot2(0x0),
  80         fContPlot3(0x0),
  81         fitter(0)
  82 {
  83         //
  84         // constructor
  85         //
  86         for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
  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");
 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)");
 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;
 150                 fType         = c.fType;
 151                 fLikely       = c.fLikely;
 152                 fNcand        = c.fNcand;
 153                 fContPlot1    = c.fContPlot1;
 154                 fContPlot2    = c.fContPlot2;
 155                 fContPlot3    = c.fContPlot3;
 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),
 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)
 175 {
 176         //
 177         // Copy Constructor
 178         //
 179         for (Int_t i=0; i<45; ++i) fParamStartValues[i]=c.fParamStartValues[i];
 180 }
 181 //_______________________________________________________________________________________
 182 AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
 183 {
 184         //
 185         //destructor
 186         //
 187         delete fLikely;
 188 }
 189 //_______________________________________________________________________________________
 190 Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
 191 {
 192         //
 193         // performs the minimization
 194         //
 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                 }
 245
 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;}
 252         Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
 253         fitter->PrintResults(4,0);
 254
 255         if(step == 3) {
 256
 257                 TMinuit* minuitPoint = fitter->GetMinuit();
 258
 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");
 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
 311         f = fLikely->EvaluateLikelihood(fX,fM,fType,fNcand);
 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 //_______________________________________________________________________________________
 321 void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[45])
 322 {
 323         for(Int_t index=0; index < 45; index++) fParamStartValues[index] = inputparamvalues[index];
 324 }
 325 //_______________________________________________________________________________________
 326 void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
 327 {
 328         //
 329         // Sets constants for the resolution function
 330         //
 331         fLikely->SetResolutionConstants(resolutionConst,type);
 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