analysis/FrawleyTrackEvalMacros/CBfitter_1state.C

0001 #include <TF1.h>
0002 #include <TMath.h>
0003 #include <TFile.h>
0004 #include <TCanvas.h>
0005 #include <TH1.h>
0006 #include <TGraph.h>
0007 #include <TLegend.h>
0008 #include <TStyle.h>
0009 #include <TROOT.h>
0010 #include <TLatex.h>
0011
0012 Double_t CBcalc(Double_t *xx, Double_t *par)
0013 {
0014   double f;
0015   double x = xx[0];
0016
0017   // The four parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0018   double alpha = par[0];
0019   double n = par[1];
0020   double x_mean = par[2];
0021   double sigma = par[3];
0022   double N = par[4];
0023
0024   // we need:
0025   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0026   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0027
0028   // The Crystal Ball function is:
0029   if( (x-x_mean)/sigma > -alpha)
0030     {
0031       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0032     }
0033   else
0034     {
0035       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0036     }
0037
0038   double Nexp = par[5];
0039   double slope = par[6];
0040
0041   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0042   fexp->SetParameter(0,Nexp);
0043   fexp->SetParameter(1,slope);
0044
0045   f = f + fexp->Eval(x);
0046   return f;
0047 }
0048
0049 Double_t Upscalc(Double_t *xx, Double_t *par)
0050 {
0051   double x = xx[0];
0052
0053   // The input parameters are:
0054   double alpha1s = par[0];
0055   double n1s = par[1];
0056   double sigma1s = par[2];
0057   double m1s = par[3];
0058   double m2s = par[4];
0059   double m3s = par[5];
0060   double N1s = par[6];
0061   double N2s = par[7];
0062   double N3s = par[8];
0063   double Nexp = par[9];
0064   double slope = par[10];
0065
0066   // Construct the Y(1S) CB shape
0067
0068   TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
0069   f1->SetParameter(0,alpha1s);
0070   f1->SetParameter(1,n1s);
0071   f1->SetParameter(2,m1s);
0072   f1->SetParameter(3,sigma1s);
0073   f1->SetParameter(4,N1s);
0074
0075   TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
0076   f2->SetParameter(0,alpha1s);
0077   f2->SetParameter(1,n1s);
0078   f2->SetParameter(2,m2s);
0079   f2->SetParameter(3,sigma1s);
0080   f2->SetParameter(4,N2s);
0081
0082   TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
0083   f3->SetParameter(0,alpha1s);
0084   f3->SetParameter(1,n1s);
0085   f3->SetParameter(2,m3s);
0086   f3->SetParameter(3,sigma1s);
0087   f3->SetParameter(4,N3s);
0088
0089   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0090   fexp->SetParameter(0,Nexp);
0091   fexp->SetParameter(1,slope);
0092
0093
0094
0095   double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
0096
0097   return mass;
0098 }
0099
0100
0101 void CBfitter_1state()
0102 {
0103   gROOT->SetStyle("Plain");
0104   gStyle->SetOptStat(0);
0105   gStyle->SetOptFit(0);
0106   gStyle->SetOptTitle(0);
0107
0108   TFile *file1S, *file2S, *file3S;
0109   TH1 *recomass1S, *recomass2S, *recomass3S;
0110
0111   bool draw_gauss = false;
0112
0113   // if all are false, do AuAu 0-10%
0114   bool pp = false;
0115   bool pAu = false;
0116   bool AuAu_20pc = false;
0117   bool AuAu_10pc = true;
0118
0119   //bool do_subtracted = true;
0120
0121   file1S = new TFile("root_files/ntp_quarkonium_out.root");
0122
0123   if(!file1S)
0124     {
0125       cout << "Failed to open file  " << endl;
0126       exit(1);
0127     }
0128
0129   recomass1S = (TH1 *)file1S->Get("recomass0");
0130
0131   TH1 *recomass = (TH1*)recomass1S->Clone("recomass1");
0132   if(!recomass)
0133     {
0134       cout << "Failed to get recomass histogram  " << endl;
0135       exit(1);
0136     }
0137
0138   recomass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
0139   recomass->Sumw2();
0140
0141   int nrebin = 1;  // set to 2 to match background histo binning, but this worsens resolution slightly. Use 1 for signal fit
0142   //int nrebin = 2;  // set to 2 to match background histo binning, but this worsens resolution slightly. Use 1 for signal fit
0143   recomass->Rebin(nrebin);
0144
0145   TCanvas *cups = new TCanvas("cups","cups",5,5,800,800);
0146   recomass->SetTitle("Y(1S,2S,3S) #rightarrow e^{+}e^{-}");
0147   recomass->SetMarkerStyle(20);
0148   recomass->SetMarkerSize(1);
0149   recomass->SetLineStyle(kSolid);
0150   recomass->SetLineWidth(2);
0151   //  recomass->SetMaximum(700);
0152   recomass->DrawCopy("p");
0153
0154
0155   TF1 *f1S = new TF1("f1S",CBcalc,7,11,7);
0156   f1S->SetParameter(0, 1.0);     // alpha
0157   f1S->SetParameter(1, 1.0);      // n
0158   f1S->SetParameter(2, 9.46);      // xmean
0159   f1S->SetParameter(3, 0.08);     // sigma
0160   f1S->SetParameter(4, 2000.0);    // N
0161   //f1S->SetParameter(4, 200.0);    // N
0162   f1S->SetParameter(5,3.5);
0163   f1S->SetParameter(6,0.05);
0164
0165   f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
0166   f1S->SetLineColor(kBlue);
0167   f1S->SetLineWidth(3);
0168   f1S->SetLineStyle(kDashed);
0169
0170   recomass->Fit(f1S);
0171   f1S->Draw("same");
0172   cout << "f1S pars " <<  f1S->GetParameter(3) << "   " << f1S->GetParError(3) << endl;
0173
0174   char resstr[500];
0175   sprintf(resstr,"#sigma_{1S} = %.1f #pm %.1f MeV", f1S->GetParameter(3)*1000, f1S->GetParError(3)*1000);
0176   TLatex *res = new TLatex(0.13,0.55,resstr);
0177   res->SetNDC();
0178   res->SetTextSize(0.05);
0179   res->Draw();
0180
0181
0182   double binw = recomass->GetBinWidth(1);
0183   double renorm = 1.0/binw;   // (1 / (bin_width of data in GeV) )
0184   cout << "renorm = " << renorm << endl;
0185
0186   cout << "Area of f1S is " << renorm * f1S->Integral(7,11) << endl;
0187
0188   // Extract ratio of yield in central gaussian to total
0189
0190   TF1 *fgauss = new TF1("fgauss","gaus(0)",7,11);
0191   fgauss->SetParameter(0, f1S->GetParameter(4));
0192   fgauss->SetParameter(1, f1S->GetParameter(2));
0193   fgauss->SetParameter(2, f1S->GetParameter(3));
0194   fgauss->SetLineColor(kRed);
0195   if(draw_gauss) fgauss->Draw("same");
0196
0197   // calculate fraction of yield in gaussian
0198   double area_fgauss =  fgauss->Integral(7,11) * renorm;
0199   double area_f1S = f1S->Integral(7,11) * renorm;
0200   double fraction = area_fgauss / area_f1S;
0201
0202
0203   cout << "Parameters of fgauss = " << fgauss->GetParameter(0) << "  " << fgauss->GetParameter(1) << "  " << fgauss->GetParameter(2) << " Area of fgauss is " << renorm * fgauss->Integral(7,11) << " fraction in fgauss " << area_fgauss / area_f1S << endl;
0204
0205   char labfrac[500];
0206   sprintf(labfrac, "Gauss fraction %.2f", fraction);
0207   TLatex *lab = new TLatex(0.13,0.75,labfrac);
0208   lab->SetNDC();
0209   lab->SetTextSize(0.05);
0210   if(draw_gauss)  lab->Draw();
0211
0212
0213   /*
0214   TLatex *lab;
0215
0216   if(pp)
0217     lab = new TLatex(0.15,0.75,"p+p, 10 weeks");
0218   else if(pAu)
0219     lab = new TLatex(0.20,0.75,"#splitline{p+Au 0-20%}{10 weeks}");
0220   else if(AuAu_20pc)
0221     lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{20B events}");
0222   else
0223     lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-10%}{10B events}");
0224
0225   lab->SetNDC();
0226   lab->SetTextSize(0.05);
0227   lab->Draw();
0228
0229   char resstr[500];
0230   sprintf(resstr,"#sigma_{1S} = %.0f #pm %.0f MeV",sigma1s,sigma1s_err);
0231   TLatex *res = new TLatex(0.13,0.55,resstr);
0232   res->SetNDC();
0233   res->SetTextSize(0.05);
0234   res->Draw();
0235
0236   */
0237   /*
0238   if (do_subtracted)
0239     {
0240       //========================================================
0241       // Now make the background_subtracted spectrum and show the
0242       // CB fits on it
0243
0244       //=======================================
0245       // Read in the subtracted mass spectrum
0246
0247       // There is only one file to read, it was written by upsilon_mass_plus_background(_pp).C
0248       TFile *file_all;
0249       if(pp)
0250     file_all = new TFile("SplusBtotal_pp_to80.root");
0251       else
0252     file_all = new TFile("SplusBtotal_10Bevts_0_10pc_rej90_eff70.root");
0253
0254       TH1 *hsplusb, *hsub, *comb_bgd, *corr_bgd;
0255       hsplusb = (TH1 *)file_all->Get("hSplusB");
0256       hsub = (TH1 *)file_all->Get("hsub");
0257       corr_bgd =  (TH1 *)file_all->Get("hdy_new");
0258       if(!pp)
0259     comb_bgd =  (TH1 *)file_all->Get("hfake_new_bgd");
0260
0261       // Fit the correlated background with a function
0262
0263       //TCanvas *c = new TCanvas("c","c",30,30,900,600);
0264       TCanvas *c = new TCanvas("c","c",30,30,800,800);
0265
0266       TF1 *fbgd = new TF1("fbgd","[0]*exp([1]+[2]*x)",7,11);
0267       fbgd->SetLineColor(kRed);
0268       fbgd->SetLineWidth(3);
0269       fbgd->SetLineStyle(kDashed);
0270       corr_bgd->Fit(fbgd);
0271
0272       hsub->SetMarkerStyle(20);
0273       hsub->SetMarkerSize(1);
0274       hsub->GetXaxis()->SetRangeUser(7.8,11.0);
0275       hsub->Draw("p");
0276       f1S->Draw("same");
0277       f2S->Draw("same");
0278       f3S->Draw("same");
0279       //fexp->Draw("same");
0280       fbgd->Draw("same");
0281
0282       // Now want to make the total of all of the fit functions
0283
0284       static const int NSTEP = 200;
0285       double xm[NSTEP];
0286       double ym[NSTEP];
0287
0288       double step = 4.0 / (double) NSTEP;
0289       for(int i=0;i<NSTEP;i++)
0290     {
0291       double x = 7.0 + (double) i * step;
0292       double y = f1S->Eval(x);
0293       y += f2S->Eval(x);
0294       y += f3S->Eval(x);
0295       y += fbgd->Eval(x);
0296       y += fexp->Eval(x);
0297
0298       xm[i] = x;
0299       ym[i] = y;
0300     }
0301       TGraph *gtot = new TGraph(NSTEP,xm,ym);
0302       gtot->SetLineWidth(3);
0303       gtot->SetLineStyle(kSolid);
0304       gtot->SetLineColor(kBlue);
0305       gtot->Draw("l");
0306
0307       TLegend *leg = new TLegend(0.62,0.58,0.9,0.88,"");
0308       leg->SetBorderSize(0);
0309       leg->SetFillColor(0);
0310       leg->SetFillStyle(0);
0311       leg->AddEntry(hsub,"++,-- subtracted","p");
0312       leg->AddEntry(f1S,"Y(1S)","l");
0313       leg->AddEntry(f2S,"Y(2S)","l");
0314       leg->AddEntry(f3S,"Y(3S)","l");
0315       leg->AddEntry(fbgd,"correlated bkg","l");
0316
0317       leg->Draw();
0318
0319       TLatex *lab2;
0320
0321       if(pp)
0322     lab2 = new TLatex(0.15,0.75,"p+p, 10 weeks");
0323       else
0324     lab2 = new TLatex(0.20,0.75,"#splitline{Au+Au 0-10%}{10B events}");
0325
0326       lab2->SetNDC();
0327       lab2->SetTextSize(0.05);
0328       lab2->Draw();
0329
0330     }
0331   */
0332 }