analysis/FrawleyTrackEvalMacros/CBfitter_jpsi_1state.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 <TGraphErrors.h>
0009 #include <TLegend.h>
0010 #include <TStyle.h>
0011 #include <TROOT.h>
0012 #include <TLatex.h>
0013
0014 TH1D *recomass;
0015 TH1D *recomassFG;
0016 TH1D *recomassBG;
0017
0018 Double_t CBcalc(Double_t *xx, Double_t *par)
0019 {
0020   // Crystal Ball lineshape plus exponential background
0021   // Used for chisq fit to mixed event background data
0022
0023   double f;
0024   double x = xx[0];
0025
0026   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0027   double alpha = par[0];
0028   double n = par[1];
0029   double x_mean = par[2];
0030   double sigma = par[3];
0031   double N = par[4];
0032
0033   // we need:
0034   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0035   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0036
0037   // The Crystal Ball function is:
0038   if( (x-x_mean)/sigma > -alpha)
0039     {
0040       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0041     }
0042   else
0043     {
0044       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0045     }
0046
0047   return f;
0048 }
0049
0050 Double_t CBcalc_exp(Double_t *xx, Double_t *par)
0051 {
0052   // Crystal Ball lineshape plus exponential background
0053   // Used for chisq fit to mixed event background data
0054
0055   double f;
0056   double x = xx[0];
0057
0058   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0059   double alpha = par[0];
0060   double n = par[1];
0061   double x_mean = par[2];
0062   double sigma = par[3];
0063   double N = par[4];
0064   // and we add an exponential background
0065   double Nexp = par[5];
0066   double slope = par[6];
0067
0068   // we need:
0069   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0070   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0071
0072   // The Crystal Ball function is:
0073   if( (x-x_mean)/sigma > -alpha)
0074     {
0075       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0076     }
0077   else
0078     {
0079       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0080     }
0081
0082   f += Nexp * exp(slope * x);
0083
0084   return f;
0085 }
0086
0087 Double_t CBcalc_LL(Double_t *xx, Double_t *par)
0088 {
0089   // Crystal Ball lineshape plus exponential background plus mixed event combinatorial background
0090   // Used for LL fit to foreground data
0091   double f;
0092   double x = xx[0];
0093
0094   // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
0095   double alpha = par[0];
0096   double n = par[1];
0097   double x_mean = par[2];
0098   double sigma = par[3];
0099   double N = par[4];
0100   // and we add an exponential background
0101   double Nexp = par[5];
0102   double slope = par[6];
0103
0104   /*
0105   TF1 *fcbexp = new TF1("fcbexp",CBcalc_exp,2.0,4.0,7);
0106   fcbexp->SetParameter(0,alpha);
0107   fcbexp->SetParameter(1,n);
0108   fcbexp->SetParameter(2,x_mean);
0109   fcbexp->SetParameter(3,sigma);
0110   fcbexp->SetParameter(4, N);
0111   fcbexp->SetParameter(5,Nexp);
0112   fcbexp->SetParameter(6,slope);
0113   f = fcbexp->Eval(x);
0114   */
0115
0116   // we need:
0117   double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
0118   double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
0119
0120   // The Crystal Ball function is:
0121   if( (x-x_mean)/sigma > -alpha)
0122     {
0123       f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
0124     }
0125   else
0126     {
0127       f = N * A * pow(B - (x-x_mean)/sigma, -n);
0128     }
0129
0130   f += Nexp * exp(slope * x);
0131
0132   // Add the mixed event background to the fit function
0133   int ibin = recomassBG->FindBin(x);
0134
0135   double bg12 = recomassBG->GetBinContent(ibin);
0136   f += bg12;
0137
0138
0139   return f;
0140 }
0141
0142 Double_t Upscalc(Double_t *xx, Double_t *par)
0143 {
0144   double x = xx[0];
0145
0146   // The input parameters are:
0147   double alpha1s = par[0];
0148   double n1s = par[1];
0149   double sigma1s = par[2];
0150   double m1s = par[3];
0151   double m2s = par[4];
0152   double m3s = par[5];
0153   double N1s = par[6];
0154   double N2s = par[7];
0155   double N3s = par[8];
0156   double Nexp = par[9];
0157   double slope = par[10];
0158
0159   // Construct the Y(1S) CB shape
0160
0161   TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
0162   f1->SetParameter(0,alpha1s);
0163   f1->SetParameter(1,n1s);
0164   f1->SetParameter(2,m1s);
0165   f1->SetParameter(3,sigma1s);
0166   f1->SetParameter(4,N1s);
0167
0168   TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
0169   f2->SetParameter(0,alpha1s);
0170   f2->SetParameter(1,n1s);
0171   f2->SetParameter(2,m2s);
0172   f2->SetParameter(3,sigma1s);
0173   f2->SetParameter(4,N2s);
0174
0175   TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
0176   f3->SetParameter(0,alpha1s);
0177   f3->SetParameter(1,n1s);
0178   f3->SetParameter(2,m3s);
0179   f3->SetParameter(3,sigma1s);
0180   f3->SetParameter(4,N3s);
0181
0182   TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
0183   fexp->SetParameter(0,Nexp);
0184   fexp->SetParameter(1,slope);
0185
0186
0187   double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
0188
0189   return mass;
0190 }
0191
0192 void CBfitter_jpsi_1state()
0193 {
0194   gROOT->SetStyle("Plain");
0195   gStyle->SetOptStat(0);
0196   gStyle->SetOptFit(0);
0197   gStyle->SetOptTitle(0);
0198
0199   bool LL_fit = true;
0200
0201   double maxmass = 3.7;
0202   double minmass = 2.4;
0203
0204   TFile *file1S;
0205   file1S = new TFile("Runcuts_DIMU_MUTRONLY_Run15pp200_JPSISIM_NOEMBED_C_PTBIN0_0_NONE.root");
0206   if(!file1S)
0207     {
0208       cout << " Failed to open input root file" << endl;
0209       return;
0210     }
0211
0212   TH2D *recomass2D[2];
0213   TH1D *recomass[2][4];
0214   for(int iarm = 0; iarm < 2; ++iarm)
0215     {
0216       // Get the 2D histos
0217       char hname[500];
0218       sprintf(hname, "mass_y_same_arm%i_chg0", iarm);
0219       recomass2D[iarm] = (TH2D *)file1S->Get(hname);
0220
0221       if(!recomass2D[iarm])
0222     cout << "Failed to get 2D histogram " << endl;
0223
0224       // and project the four rapidity bins
0225       for(int irap =0; irap < 4; ++irap)
0226     {
0227       char h1dname[500];
0228       sprintf(h1dname, "recomass%i%i", iarm, irap);
0229       recomass[iarm][irap] = recomass2D[iarm]->ProjectionX(h1dname, irap+1, irap+1);
0230     }
0231     }
0232
0233   TF1 *f1S[2][4];
0234   TCanvas *cups = new TCanvas("cups","cups",5,5,1600,800);
0235   cups->Divide(4,2);
0236   for(int iarm = 0; iarm<2; ++iarm)
0237     {
0238       for(int irap = 0; irap < 4; ++irap)
0239     {
0240       cout << "Fit for iarm = " << iarm << " irap = " << irap << endl;
0241
0242       cups->cd(irap + iarm*4 + 1);
0243       recomass[iarm][irap]->Draw();
0244
0245       char fname[500];
0246       sprintf(fname, "f1S%i%i", iarm, irap);
0247       f1S[iarm][irap] = new TF1(fname,CBcalc, minmass, maxmass, 5);
0248
0249       f1S[iarm][irap]->SetParameter(0, 1.5);     // alpha
0250       f1S[iarm][irap]->SetParameter(1, 100.0);      // n
0251       f1S[iarm][irap]->SetParameter(2, 3.12);     // mass
0252       f1S[iarm][irap]->SetParameter(3, 0.13);     // sigma
0253       f1S[iarm][irap]->SetParameter(4, 300);    // N
0254       f1S[iarm][irap]->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
0255       f1S[iarm][irap]->SetLineColor(kBlue);
0256       f1S[iarm][irap]->SetLineWidth(3);
0257       f1S[iarm][irap]->SetLineStyle(kDashed);
0258
0259       recomass[iarm][irap]->GetXaxis()->SetRangeUser(2.2,3.8);
0260       recomass[iarm][irap]->Fit(f1S[iarm][irap]);
0261       //recomass[iarm][irap]->Fit(f1S[iarm][irap],"L");
0262       recomass[iarm][irap]->DrawCopy();
0263       f1S[iarm][irap]->Draw("same");
0264     }
0265     }
0266
0267   // Get the mass and width so we can plot them
0268
0269   TCanvas *cgr = new TCanvas("cgr", "cgr", 5, 5, 1600, 800);
0270   cgr->Divide(2,2);
0271
0272   TLatex *lab[2];
0273   double rap[2][4] = {-2.075, -1.825, -1.575, -1.325, 1.325, 1.575, 1.825, 2.075};
0274   double rap_err[4] = {0};
0275   //TGraph *grmean[2];
0276   //TGraph *grsigma[2];
0277   TGraphErrors *grmean[2];
0278   TGraphErrors *grsigma[2];
0279   double mean[2][4];
0280   double sigma[2][4];
0281   double mean_err[2][4];
0282   double sigma_err[2][4];
0283   for(int iarm = 0; iarm < 2; ++iarm)
0284     {
0285       for(int irap = 0; irap < 4; ++irap)
0286     {
0287       mean[iarm][irap] =  f1S[iarm][irap]->GetParameter(2);
0288       mean_err[iarm][irap] =  f1S[iarm][irap]->GetParError(2);
0289       sigma[iarm][irap] =  f1S[iarm][irap]->GetParameter(3);
0290       sigma_err[iarm][irap] =  f1S[iarm][irap]->GetParError(3);
0291     }
0292       grmean[iarm] = new TGraphErrors(4, rap[iarm], mean[iarm], rap_err, mean_err[iarm]);
0293       grsigma[iarm] = new TGraphErrors(4, rap[iarm], sigma[iarm], rap_err, sigma_err[iarm]);
0294
0295       cgr->cd(1+iarm);
0296       grmean[iarm]->SetMarkerStyle(20);
0297       grmean[iarm]->SetMarkerSize(2);
0298       grmean[iarm]->GetYaxis()->SetTitle("mass (GeV/c^{2})");
0299       grmean[iarm]->GetYaxis()->SetTitleSize(0.05);
0300       grmean[iarm]->GetXaxis()->SetTitle("y");
0301       grmean[iarm]->GetXaxis()->SetTitleSize(0.05);
0302       grmean[iarm]->Draw();
0303
0304       char arm[500];
0305       sprintf(arm,"arm %i",iarm);
0306       lab[iarm] = new TLatex(0.2, 0.350, arm);
0307       lab[iarm]->SetNDC();
0308       lab[iarm]->SetTextSize(0.1);
0309       lab[iarm]->Draw();
0310
0311       cgr->cd(3+iarm);
0312       grsigma[iarm]->SetMarkerStyle(20);
0313       grsigma[iarm]->SetMarkerSize(2);
0314       grsigma[iarm]->GetYaxis()->SetTitle("width (GeV/c^{2})");
0315       grsigma[iarm]->GetYaxis()->SetTitleSize(0.05);
0316       grsigma[iarm]->GetXaxis()->SetTitle("y");
0317       grsigma[iarm]->GetYaxis()->SetTitleSize(0.05);
0318       grsigma[iarm]->Draw();
0319       lab[iarm]->Draw();
0320     }
0321
0322
0323
0324 }
0325
0326
0327
0328
0329