analysis/FrawleyTrackEvalMacros/CBfitter_jpsi_LL_single.C

0001 #include <TF1.h>
0002 #include <TMath.h>
0003 #include <TFile.h>
0004 #include <TCanvas.h>
0005 #include <TH1.h>
0006 #include <TH2.h>
0007 #include <TGraph.h>
0008 #include <TLegend.h>
0009 #include <TStyle.h>
0010 #include <TROOT.h>
0011 #include <TLatex.h>
0012
0013 TH1D *recomass;
0014 TH1D *recomassFG;
0015 TH1D *recomassBG;
0016
0017 Double_t CBcalc(Double_t *xx, Double_t *par)
0018 {
0019   // Crystal Ball lineshape plus exponential background
0020   // Used for chisq fit to mixed event background data
0021
0022   double f;
0023   double x = xx[0];
0024
0025   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0026   double alpha = par[0];
0027   double n = par[1];
0028   double x_mean = par[2];
0029   double sigma = par[3];
0030   double N = par[4];
0031
0032   // we need:
0033   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0034   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0035
0036   // The Crystal Ball function is:
0037   if( (x-x_mean)/sigma > -alpha)
0038     {
0039       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0040     }
0041   else
0042     {
0043       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0044     }
0045
0046   return f;
0047 }
0048
0049 Double_t CBcalc_exp(Double_t *xx, Double_t *par)
0050 {
0051   // Crystal Ball lineshape plus exponential background
0052   // Used for chisq fit to mixed event background data
0053
0054   double f;
0055   double x = xx[0];
0056
0057   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0058   double alpha = par[0];
0059   double n = par[1];
0060   double x_mean = par[2];
0061   double sigma = par[3];
0062   double N = par[4];
0063   // and we add an exponential background
0064   double Nexp = par[5];
0065   double slope = par[6];
0066
0067   // we need:
0068   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0069   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0070
0071   // The Crystal Ball function is:
0072   if( (x-x_mean)/sigma > -alpha)
0073     {
0074       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0075     }
0076   else
0077     {
0078       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0079     }
0080
0081   f += Nexp * exp(slope * x);
0082
0083   return f;
0084 }
0085
0086 Double_t CBcalc_LL(Double_t *xx, Double_t *par)
0087 {
0088   // Crystal Ball lineshape plus exponential background plus mixed event combinatorial background
0089   // Used for LL fit to foreground data
0090   double f;
0091   double x = xx[0];
0092
0093   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0094   double alpha = par[0];
0095   double n = par[1];
0096   double x_mean = par[2];
0097   double sigma = par[3];
0098   double N = par[4];
0099   // and we add an exponential background
0100   double Nexp = par[5];
0101   double slope = par[6];
0102
0103   /*
0104   TF1 *fcbexp = new TF1("fcbexp",CBcalc_exp,2.0,4.0,7);
0105   fcbexp->SetParameter(0,alpha);
0106   fcbexp->SetParameter(1,n);
0107   fcbexp->SetParameter(2,x_mean);
0108   fcbexp->SetParameter(3,sigma);
0109   fcbexp->SetParameter(4, N);
0110   fcbexp->SetParameter(5,Nexp);
0111   fcbexp->SetParameter(6,slope);
0112   f = fcbexp->Eval(x);
0113   */
0114
0115   // we need:
0116   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0117   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0118
0119   // The Crystal Ball function is:
0120   if( (x-x_mean)/sigma > -alpha)
0121     {
0122       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0123     }
0124   else
0125     {
0126       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0127     }
0128
0129   f += Nexp * exp(slope * x);
0130
0131   // Add the mixed event background to the fit function
0132   int ibin = recomassBG->FindBin(x);
0133
0134   double bg12 = recomassBG->GetBinContent(ibin);
0135   f += bg12;
0136
0137
0138   return f;
0139 }
0140
0141 Double_t Upscalc(Double_t *xx, Double_t *par)
0142 {
0143   double x = xx[0];
0144
0145   // The input parameters are:
0146   double alpha1s = par[0];
0147   double n1s = par[1];
0148   double sigma1s = par[2];
0149   double m1s = par[3];
0150   double m2s = par[4];
0151   double m3s = par[5];
0152   double N1s = par[6];
0153   double N2s = par[7];
0154   double N3s = par[8];
0155   double Nexp = par[9];
0156   double slope = par[10];
0157
0158   // Construct the Y(1S) CB shape
0159
0160   TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
0161   f1->SetParameter(0,alpha1s);
0162   f1->SetParameter(1,n1s);
0163   f1->SetParameter(2,m1s);
0164   f1->SetParameter(3,sigma1s);
0165   f1->SetParameter(4,N1s);
0166
0167   TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
0168   f2->SetParameter(0,alpha1s);
0169   f2->SetParameter(1,n1s);
0170   f2->SetParameter(2,m2s);
0171   f2->SetParameter(3,sigma1s);
0172   f2->SetParameter(4,N2s);
0173
0174   TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
0175   f3->SetParameter(0,alpha1s);
0176   f3->SetParameter(1,n1s);
0177   f3->SetParameter(2,m3s);
0178   f3->SetParameter(3,sigma1s);
0179   f3->SetParameter(4,N3s);
0180
0181   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0182   fexp->SetParameter(0,Nexp);
0183   fexp->SetParameter(1,slope);
0184
0185
0186   double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
0187
0188   return mass;
0189 }
0190
0191 void CBfitter_jpsi_LL_single()
0192 {
0193   gROOT->SetStyle("Plain");
0194   gStyle->SetOptStat(0);
0195   gStyle->SetOptFit(0);
0196   gStyle->SetOptTitle(0);
0197
0198   bool LL_fit = true;
0199
0200   double maxmass = 3.5;
0201
0202   TFile *file1S;
0203    file1S = new TFile("Final_Hists_mee_meeBG_Type_VetoB0.root");
0204
0205   if(!file1S)
0206     {
0207       cout << " Failed to open input root file" << endl;
0208       return;
0209     }
0210
0211   recomassFG = (TH1D *)file1S->Get("Hunlk");
0212   if(!recomassFG)
0213     cout << "Failed to get Hunlk" << endl;
0214
0215   recomassBG = (TH1D *)file1S->Get("HmxNorm");
0216   if(!recomassBG)
0217     cout << "Failed to get HmxNorm" << endl;
0218
0219   TCanvas *cups = new TCanvas("cups","cups",5,5,1000,800);
0220   recomassFG->Draw();
0221   //recomassBG->Draw("same");
0222
0223   TF1 *f1S;
0224
0225   if(LL_fit)
0226     {
0227       // fit FG with CB+exp+BG12 using LL
0228       f1S = new TF1("f1SLL",CBcalc_LL,1.3,maxmass,7);
0229     }
0230   else
0231     {
0232       // fit FG - BG with CB+exp using chisq
0233       f1S = new TF1("f1S",CBcalc_exp,1.3,maxmass,7);
0234     }
0235
0236   // These are an error weighted average of the four free fits below in the comments
0237   // Use these values to describe the J/psi line shape in fits at all cenralities
0238
0239   f1S->SetParameter(0, -0.5);     // alpha
0240   f1S->SetParameter(1, 2.0);      // n
0241   f1S->SetParameter(2, 3.07);     // mass
0242   f1S->SetParameter(3, 0.11);     // sigma
0243   f1S->SetParameter(4, 850);    // N
0244   f1S->FixParameter(5, 0.0);    // Nexp
0245   f1S->FixParameter(6, -1.00);    // slope
0246   f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S","Nexp","slope");
0247
0248   f1S->SetLineColor(kBlue);
0249   f1S->SetLineWidth(3);
0250   f1S->SetLineStyle(kDashed);
0251
0252   TCanvas *csig = new TCanvas("csig","csig",5,5,1000,800);
0253
0254   TF1 *fitResult= 0;
0255   if(LL_fit)
0256     {
0257       // log likelihood fit with integer data option is "L"
0258       //recomassFG->Fit(f1S,"LQ");
0259       recomassFG->Fit(f1S,"L R");
0260       recomassFG->DrawCopy();
0261       fitResult = recomassFG->GetFunction("f1SLL");
0262     }
0263   else
0264     {
0265       recomassFG->Fit(f1S);
0266       recomassFG->DrawCopy();
0267       fitResult = recomassFG->GetFunction("CBcalc_exp");
0268     }
0269   if(fitResult)
0270     {
0271       cout << " chisquare " <<  fitResult->GetChisquare() << " NDF " <<  fitResult->GetNDF() << endl;
0272     }
0273   else
0274     {
0275       cout << "Failed to get fitResult" << endl;
0276       return;
0277     }
0278
0279   // Plot the fit results
0280   //======================
0281   // Make a function that is only the CB lineshape and plot it
0282   // We want CBcalc for this
0283   TF1 *fcbonly = new TF1("fcbonly",CBcalc,1.3,4.0,5);
0284   fcbonly->SetParameter(0,f1S->GetParameter(0));
0285   fcbonly->SetParameter(1,f1S->GetParameter(1));
0286   fcbonly->SetParameter(2,f1S->GetParameter(2));
0287   fcbonly->SetParameter(3,f1S->GetParameter(3));
0288   fcbonly->SetParameter(4,f1S->GetParameter(4));
0289   fcbonly->SetLineColor(kMagenta);
0290   fcbonly->Draw("same");
0291
0292   // make a function that is only the exponential lineshape and plot it
0293   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",1.3,4.0);
0294   fexp->SetParameter(0,f1S->GetParameter(5));
0295   fexp->SetParameter(1,f1S->GetParameter(6));
0296
0297   fexp->SetLineColor(kGreen);
0298   fexp->SetLineWidth(3);
0299   fexp->SetLineStyle(kDashed);
0300   fexp->Draw("same");
0301
0302   // Finally, plot the mixed event background
0303   if(LL_fit)
0304     {
0305       recomassBG->SetLineColor(kRed);
0306       recomassBG->DrawCopy("same");
0307     }
0308
0309   // calculate yield uncertainty from fit
0310   double frac_yerror = f1S->GetParError(4) / f1S->GetParameter(4);
0311
0312   /*
0313   CB_yield[icent][ipt] = renorm * fcbonly->Integral(2.8,3.4) ;
0314   CB_yield_err[icent][ipt] = CB_yield[icent][ipt] * frac_yerror;
0315
0316   best_CB_yield[ialpha][in1][im1][isig] = CB_yield[icent][ipt];
0317   best_CB_yield_err[ialpha][in1][im1][isig] = CB_yield_err[icent][ipt];
0318   */
0319
0320   // On a separate canvas, plot the FG - BG and the CB+exp
0321   //========================================
0322
0323   recomassFG->DrawCopy("pe");
0324   recomassBG->DrawCopy("same");
0325   fexp->Draw("same");
0326   fcbonly->Draw("same");
0327
0328
0329 }
0330
0331
0332
0333