analysis/FrawleyTrackEvalMacros/CBfitter.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   return f;
0039 }
0040
0041 Double_t Upscalc(Double_t *xx, Double_t *par)
0042 {
0043   double x = xx[0];
0044
0045   // The input parameters are:
0046   double alpha1s = par[0];
0047   double n1s = par[1];
0048   double sigma1s = par[2];
0049   double m1s = par[3];
0050   double m2s = par[4];
0051   double m3s = par[5];
0052   double N1s = par[6];
0053   double N2s = par[7];
0054   double N3s = par[8];
0055   double Nexp = par[9];
0056   double slope = par[10];
0057
0058   // Construct the Y(1S) CB shape
0059
0060   TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
0061   f1->SetParameter(0,alpha1s);
0062   f1->SetParameter(1,n1s);
0063   f1->SetParameter(2,m1s);
0064   f1->SetParameter(3,sigma1s);
0065   f1->SetParameter(4,N1s);
0066
0067   TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
0068   f2->SetParameter(0,alpha1s);
0069   f2->SetParameter(1,n1s);
0070   f2->SetParameter(2,m2s);
0071   f2->SetParameter(3,sigma1s);
0072   f2->SetParameter(4,N2s);
0073
0074   TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
0075   f3->SetParameter(0,alpha1s);
0076   f3->SetParameter(1,n1s);
0077   f3->SetParameter(2,m3s);
0078   f3->SetParameter(3,sigma1s);
0079   f3->SetParameter(4,N3s);
0080
0081   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0082   fexp->SetParameter(0,Nexp);
0083   fexp->SetParameter(1,slope);
0084
0085   double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
0086
0087   return mass;
0088 }
0089
0090
0091 void CBfitter()
0092 {
0093   gROOT->SetStyle("Plain");
0094   gStyle->SetOptStat(0);
0095   gStyle->SetOptTitle(1);
0096
0097   bool ups1s = true;
0098   bool ups2s = true;
0099   bool ups3s = true;
0100
0101   TFile *file1S, *file2S, *file3S;
0102   TH1 *recomass1S, *recomass2S, *recomass3S;
0103
0104   if(ups1s)
0105     {
0106       file1S = new TFile("root_files/unsuppressed_auau_24B_ntp_quarkonium_out.root");
0107       recomass1S = (TH1 *)file1S->Get("recomass0");
0108     }
0109
0110   if(ups2s)
0111     {
0112       //file2S = new TFile(" root_files/ups2s_80ns_100pions_pp.root ");
0113       recomass2S = (TH1 *)file1S->Get("recomass1");
0114     }
0115   if(ups3s)
0116     {
0117       //file3S = new TFile(" root_files/ups3s_80ns_100pions_pp.root");
0118       recomass3S = (TH1 *)file1S->Get("recomass2");
0119     }
0120
0121   TFile *massout = new TFile("upsilon_mass_histos_out.root","recreate");
0122
0123   TH1 *recomass = (TH1*)recomass1S->Clone("recomass");
0124   recomass->Sumw2();
0125   if(ups2s)
0126     recomass->Add(recomass2S);
0127   if(ups3s)
0128     recomass->Add(recomass3S);
0129   int nrebin = 1;  // a rebin of 2 is needed to match background histo binning, use 1 if want best resolution
0130   recomass->Rebin(nrebin);
0131
0132   //TCanvas *cups = new TCanvas("cups","cups",5,5,800,600);
0133   TCanvas *cups = new TCanvas("cups","cups",5,5,800,800);
0134   recomass->SetTitle("Y(1S,2S,3S) #rightarrow e^{+}e^{-}");
0135   recomass->GetYaxis()->SetNdivisions(205);
0136   recomass->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
0137   recomass->SetMarkerStyle(20);
0138   recomass->SetMarkerSize(1);
0139   recomass->SetLineStyle(kSolid);
0140   recomass->SetLineWidth(2);
0141   //recomass->SetMaximum(810);
0142   recomass->DrawCopy("p");
0143
0144
0145   TF1 *fitups = new TF1("fitups",Upscalc,7,11,11);
0146   fitups->SetParameter(0, 1.05);
0147   fitups->SetParameter(1, 1.5);
0148   fitups->SetParameter(2, 0.09);
0149   fitups->SetParameter(3, 9.44);
0150   fitups->SetParameter(4, 10.0);
0151   fitups->SetParameter(5, 10.3);
0152   fitups->SetParameter(6, 900.0);
0153   fitups->SetParameter(7, 200.0);
0154   fitups->SetParameter(8, 100.0);
0155   //fitups->SetParameter(9, 10.0);
0156   fitups->FixParameter(9, 0.0);
0157   fitups->FixParameter(10, -0.2);
0158   fitups->SetLineColor(kBlack);
0159   fitups->SetParNames("alpha1S","n1S","sigma1S","m1S","m2S","m3S","N1S","N2S","N3S","Nexp","expSlope");
0160
0161   recomass->Fit(fitups);
0162   //fitups->Draw("same");
0163
0164   double sigma1s = 1000.0 * fitups->GetParameter(2);
0165   double sigma1s_err = 1000.0 * fitups->GetParError(2);
0166   cout << " sigma1s " << sigma1s << " par2 " << fitups->GetParameter(2) << endl;
0167   // Now draw the individual states
0168
0169   TF1 *f1S = new TF1("f1S",CBcalc,7,11,5);
0170   f1S->SetParameter(0,fitups->GetParameter(0));
0171   f1S->SetParameter(1,fitups->GetParameter(1));
0172   f1S->SetParameter(2,fitups->GetParameter(3));
0173   f1S->SetParameter(3,fitups->GetParameter(2));
0174   f1S->SetParameter(4,fitups->GetParameter(6));
0175   f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
0176   f1S->SetLineColor(kBlue);
0177   f1S->SetLineWidth(3);
0178   f1S->SetLineStyle(kDashed);
0179
0180   f1S->Draw("same");
0181
0182   TF1 *f2S = new TF1("f2S",CBcalc,7,11,5);
0183   f2S->SetParameter(0,fitups->GetParameter(0));
0184   f2S->SetParameter(1,fitups->GetParameter(1));
0185   f2S->SetParameter(2,fitups->GetParameter(4));
0186   f2S->SetParameter(3,fitups->GetParameter(2));
0187   f2S->SetParameter(4,fitups->GetParameter(7));
0188   f2S->SetParNames("alpha1S","n1S","m2S","sigma1S","N2S");
0189   f2S->SetLineColor(kRed);
0190   f2S->SetLineWidth(3);
0191   f2S->SetLineStyle(kDashed);
0192
0193   f2S->Draw("same");
0194
0195   TF1 *f3S = new TF1("f3S",CBcalc,7,11,5);
0196   f3S->SetParameter(0,fitups->GetParameter(0));
0197   f3S->SetParameter(1,fitups->GetParameter(1));
0198   f3S->SetParameter(2,fitups->GetParameter(5));
0199   f3S->SetParameter(3,fitups->GetParameter(2));
0200   f3S->SetParameter(4,fitups->GetParameter(8));
0201   f3S->SetParNames("alpha1S","n1S","m3S","sigma1S","N3S");
0202   f3S->SetLineColor(kMagenta);
0203   f3S->SetLineWidth(3);
0204   f3S->SetLineStyle(kDashed);
0205
0206   f3S->Draw("same");
0207
0208   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0209   fexp->SetParameter(0,fitups->GetParameter(9));
0210   fexp->SetParameter(1,fitups->GetParameter(10));
0211
0212   fexp->SetLineColor(kGreen);
0213   fexp->SetLineWidth(3);
0214   fexp->SetLineStyle(kDashed);
0215   fexp->Draw("same");
0216
0217   double binw = recomass->GetBinWidth(1);
0218   double renorm = 1.0/binw;   // (1 / (bin_width of data in GeV) )
0219   //double renorm = 25.0;   // (1 / (bin_width of data in GeV) )
0220   cout << "renorm = " << renorm << endl;
0221
0222   cout << "Area of f1S is " << renorm * f1S->Integral(7,11) << endl;
0223   cout << "Area of f2S is " << renorm * f2S->Integral(7,11) << endl;
0224   cout << "Area of f3S is " << renorm * f3S->Integral(7,11) << endl;
0225   cout << "Area of fexp is " << renorm * fexp->Integral(7,11) << endl;
0226
0227   bool pp = true;
0228   bool AuAu_20pc = false;
0229
0230   bool do_subtracted = false;
0231
0232   TLatex *lab;
0233
0234   if(pp)
0235     lab = new TLatex(0.15,0.75,"#splitline{p+p, 197 pb^{-1}}{signal only}");
0236   else if(AuAu_20pc)
0237     lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{10B events}");
0238   else
0239     lab = new TLatex(0.20,0.75,"#splitline{Au+Au MB}{50B events}");
0240
0241   lab->SetNDC();
0242   lab->SetTextSize(0.05);
0243   //lab->Draw();
0244
0245   char resstr[500];
0246   sprintf(resstr,"#sigma_{1S} = %.0f #pm %.1f MeV",sigma1s,sigma1s_err);
0247   TLatex *res = new TLatex(0.13,0.55,resstr);
0248   res->SetNDC();
0249   res->SetTextSize(0.05);
0250   //res->Draw();
0251
0252   recomass->Write();
0253   f1S->Write();
0254   f2S->Write();
0255   f3S->Write();
0256
0257   if (do_subtracted)
0258     {
0259       //========================================================
0260       // Now make the background_subtracted spectrum and show the
0261       // CB fits on it
0262
0263       //=======================================
0264       // Read in the subtracted mass spectrum
0265
0266       // There is only one file to read, it was written by upsilon_mass_plus_background(_pp).C
0267       TFile *file_all;
0268       if(pp)
0269     file_all = new TFile("SplusBtotal_pp_to80.root");
0270       else
0271     file_all = new TFile("SplusBtotal_rej90_to80.root");
0272
0273       TH1 *hsplusb, *hsub, *comb_bgd, *corr_bgd;
0274       hsplusb = (TH1 *)file_all->Get("hSplusB");
0275       hsub = (TH1 *)file_all->Get("hsub");
0276       corr_bgd =  (TH1 *)file_all->Get("hdy_new");
0277       if(!pp)
0278     comb_bgd =  (TH1 *)file_all->Get("hfake_new_bgd");
0279
0280       // Fit the correlated background with a function
0281
0282       //TCanvas *c = new TCanvas("c","c",30,30,900,600);
0283       TCanvas *c = new TCanvas("c","c",30,30,800,800);
0284
0285       TF1 *fbgd = new TF1("fbgd","[0]*exp([1]+[2]*x)",7,11);
0286       fbgd->SetLineColor(kRed);
0287       fbgd->SetLineWidth(3);
0288       fbgd->SetLineStyle(kDashed);
0289       corr_bgd->Fit(fbgd);
0290
0291       hsub->SetMarkerStyle(20);
0292       hsub->SetMarkerSize(1);
0293       hsub->Draw("p");
0294       f1S->Draw("same");
0295       f2S->Draw("same");
0296       f3S->Draw("same");
0297       //fexp->Draw("same");
0298       fbgd->Draw("same");
0299
0300       // Now want to make the total of all of the fit functions
0301
0302       static const int NSTEP = 200;
0303       double xm[NSTEP];
0304       double ym[NSTEP];
0305
0306       double step = 4.0 / (double) NSTEP;
0307       for(int i=0;i<NSTEP;i++)
0308     {
0309       double x = 7.0 + (double) i * step;
0310       double y = f1S->Eval(x);
0311       y += f2S->Eval(x);
0312       y += f3S->Eval(x);
0313       y += fbgd->Eval(x);
0314       y += fexp->Eval(x);
0315
0316       xm[i] = x;
0317       ym[i] = y;
0318     }
0319       TGraph *gtot = new TGraph(NSTEP,xm,ym);
0320       gtot->SetLineWidth(3);
0321       gtot->SetLineStyle(kSolid);
0322       gtot->SetLineColor(kBlue);
0323       gtot->Draw("l");
0324
0325       TLegend *leg = new TLegend(0.62,0.58,0.9,0.88,"");
0326       leg->SetBorderSize(0);
0327       leg->SetFillColor(0);
0328       leg->SetFillStyle(0);
0329       leg->AddEntry(hsub,"++,-- subtracted","p");
0330       leg->AddEntry(f1S,"Y(1S)","l");
0331       leg->AddEntry(f2S,"Y(2S)","l");
0332       leg->AddEntry(f3S,"Y(3S)","l");
0333       leg->AddEntry(fbgd,"correlated bkg","l");
0334
0335       leg->Draw();
0336
0337       TLatex *lab2;
0338
0339       if(pp)
0340     lab2 = new TLatex(0.15,0.75,"p+p, 10 weeks");
0341       else
0342     lab2 = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{10B events}");
0343
0344       lab2->SetNDC();
0345       lab2->SetTextSize(0.05);
0346       lab2->Draw();
0347
0348     }
0349 }