analysis/FrawleyTrackEvalMacros/upsilon_raa_sphenix_estimator.C

0001
0002
0003 void upsilon_raa_sphenix_estimator() {
0004
0005   //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
0006   //SetsPhenixStyle();
0007
0008 gStyle->SetOptStat(0);
0009 gStyle->SetOptTitle(0);
0010 gStyle->SetOptFit(0);
0011
0012  bool yrs13 = true;
0013  bool yrs15 = true;
0014  bool ups3s = true;
0015
0016  double offset = 0.0;
0017  if(yrs13 && yrs15)
0018    offset = 3;
0019
0020  // Input assumptions about level of suppression are taken from the theory paper:
0021   // Michael Strickland, Dennis Bazow, Nucl.Phys. A879 (2012) 25-58
0022   //==============================================
0023   // Read in Strickland curves for different eta/s values
0024   // These are exclusive - i.e. they are the
0025   // suppression of that state with no feed-down
0026   // effects included
0027   //==============================================
0028   // The curves vs eta/s read in below were provided by Strickland privately, since
0029   // the paper contains only the eta/s = 1 cases for RHIC energy
0030
0031   // Get the exclusive raa values for the Y1S, Y2S, Y3S, chib1, chib2
0032   // Checked that this is being read in correctly
0033   double str_npart[101];
0034   double str_raa[5][3][101];
0035   for(int istate=0;istate<5;istate++)
0036     {
0037       for(int ietas=0;ietas<3;ietas++)
0038     {
0039       if(istate < 3)
0040         {
0041           char fname[500];
0042           sprintf(fname,"./strickland_calculations/Y%is-potb-eta%i-npart.dat",istate+1,ietas+1);
0043           ifstream fin(fname);
0044
0045           for(int inpart=0;inpart<101;inpart++)
0046         {
0047           fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
0048           //cout << istate << " " << ietas << " " << inpart << " " << str_npart[inpart] << " " << str_raa[istate][ietas][inpart] << endl;
0049         }
0050           fin.close();
0051         }
0052
0053       if(istate > 2)
0054         {
0055           char fname[500];
0056           sprintf(fname,"./strickland_calculations/chib%i-potb-eta%i-npart.dat",istate-2,ietas+1);
0057           ifstream fin(fname);
0058
0059           for(int inpart=0;inpart<101;inpart++)
0060         {
0061           fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
0062           //cout << istate << " " << ietas << " " << inpart << " " << str_npart[inpart] << " " << str_raa[istate][ietas][inpart] << endl;
0063         }
0064           fin.close();
0065         }
0066     }
0067
0068     }
0069
0070   // Now construct the inclusive RAA values from Strickland's exclusive modifications
0071   // ff's from table II of the paper for the 1S state are as follows
0072   // They add up to 1.0
0073   double ff1S[5] = {0.51, 0.107, 0.008, 0.27, 0.105};  // Y1S, Y2S, Y3S, chib1, chib2 respectively
0074   // These are from arXiv:1210.7512
0075   double ff2S[2] = {0.5, 0.5};                         // Y2S and Y3S respectively
0076
0077   double str_raa_inclusive[3][3][101];
0078
0079   // checked the Y1S inclusive result against figure 10 of arXiv:1112.2761
0080   // There is no plot available for the 2S and 3S
0081   for(int ietas=0;ietas<3;ietas++)
0082     for(int inpart=0;inpart<101;inpart++)
0083       {
0084     str_raa_inclusive[0][ietas][inpart] =
0085       str_raa[0][ietas][inpart] * ff1S[0]
0086       + str_raa[1][ietas][inpart] * ff1S[1]
0087       + str_raa[2][ietas][inpart] * ff1S[2]
0088       + str_raa[3][ietas][inpart] * ff1S[3]
0089       + str_raa[4][ietas][inpart] * ff1S[4];
0090
0091     str_raa_inclusive[1][ietas][inpart] =
0092       str_raa[1][ietas][inpart] * ff2S[0]
0093       + str_raa[2][ietas][inpart] * ff2S[1];
0094
0095     str_raa_inclusive[2][ietas][inpart] = str_raa[2][ietas][inpart];
0096       }
0097
0098   TGraph *grth[3][3];
0099   for(int is = 0;is<3;is++)
0100     {
0101       for(int ieta =0; ieta < 3; ++ieta)
0102     {
0103       grth[is][ieta]=new TGraph(101,str_npart,str_raa_inclusive[is][ieta]);
0104     }
0105     }
0106
0107
0108  static const int NCENT = 7;
0109  // static const int NCENT3S = 4;
0110  static const int NCENT3S = 3;
0111
0112  double Ncoll[NCENT] = {1067, 858, 610, 378, 224, 94.2, 14.5};
0113  double Npart[NCENT] = {351, 302, 236, 168, 116, 61.6, 14.4};
0114  double Npart_plot[2][NCENT];
0115  for(int i=0;i<NCENT;++i)
0116    {
0117      Npart_plot[0][i] = Npart[i] + offset;
0118      Npart_plot[1][i] = Npart[i] - offset;
0119    }
0120
0121  // separate binning for Y(3S)
0122  //double Npart3S[NCENT3S] = {281, 142, 62, 14.4};
0123  double Npart3S[NCENT3S] = {243, 80.0, 14.4};
0124  double Npart3S_plot[2][NCENT3S];
0125  for(int i=0;i<NCENT3S;++i)
0126    {
0127      Npart3S_plot[0][i] = Npart3S[i] + offset;
0128      Npart3S_plot[1][i] = Npart3S[i] - offset;
0129    }
0130
0131
0132  double raa1[2][NCENT];
0133  double raa2[2][NCENT];
0134  for(int iset=0;iset<2;++iset)
0135    for(int i=0; i<NCENT; ++i)
0136      {
0137        raa1[iset][i] = grth[0][1]->Eval(Npart_plot[iset][i]);
0138        raa2[iset][i] = grth[1][1]->Eval(Npart_plot[iset][i]);
0139        cout << " icent " << i << " Npart " << Npart_plot[iset][i] << " raa1 " << raa1[iset][i] << " raa2 " << raa2[iset][i] << endl;
0140   }
0141
0142  double raa3[2][NCENT3S];
0143  for(int iset=0;iset<2;++iset)
0144    for(int i=0; i<NCENT3S; ++i)
0145      {
0146        raa3[iset][i] = grth[2][1]->Eval(Npart3S_plot[iset][i]);
0147        cout << " 3S:  icent " << i << " Npart " << Npart3S_plot[iset][i] << " raa1 " << " raa3 " << raa3[iset][i] << endl;
0148   }
0149
0150  // output the RAA values to feed back into the previous macro
0151  cout << " double raa1[" << NCENT << "] = {";
0152  for(int i=0; i< NCENT; ++i)
0153    {
0154      if(i != NCENT-1)
0155        cout << raa1[0][i] << ", ";
0156      else
0157        cout << raa1[0][i] << "}; ";
0158    }
0159  cout << endl;
0160
0161  // output the RAA values to feed back into the previous macro
0162  cout << " double raa2[" << NCENT << "] = {";
0163  for(int i=0; i< NCENT; ++i)
0164    {
0165      if(i != NCENT-1)
0166        cout << raa2[0][i] << ", ";
0167      else
0168        cout << raa2[0][i] << "}; ";
0169    }
0170  cout << endl;
0171
0172  // output the RAA values to feed back into the previous macro
0173  cout << " double raa3[" << NCENT3S << "] = {";
0174  for(int i=0; i< NCENT3S; ++i)
0175    {
0176      if(i != NCENT3S-1)
0177        cout << raa3[0][i] << ", ";
0178      else
0179        cout << raa3[0][i] << "}; ";
0180    }
0181  cout << endl;
0182
0183  // These errors are the only thing that change with different luminosity
0184  // updated 8/26/20
0185  double yrs13_erraa1[NCENT] = {0.021, 0.0227, 0.0211, 0.0251, 0.0315, 0.0365, 0.0665};
0186  double yrs13_erraa2[NCENT] = {0.0327, 0.0424, 0.0335, 0.0401, 0.0537, 0.0662, 0.157};
0187  double yrs13_erraa3[NCENT3S] = {0.042 ,0.094, 0.346};
0188
0189  //double  yrs15_erraa1[NCENT] = {0.012, 0.013, 0.012, 0.015, 0.019, 0.022, 0.042};
0190  double  yrs15_erraa1[NCENT] = {0.0136, 0.0147, 0.0137, 0.0163, 0.0203, 0.0235, 0.0435};
0191  //double yrs15_erraa2[NCENT] = {0.022, 0.024, 0.021, 0.025, 0.034, 0.041, 0.097};
0192  double yrs15_erraa2[NCENT] = {0.0235, 0.0233, 0.0207, 0.0277, 0.0364, 0.0424, 0.103};
0193  double yrs15_erraa3[NCENT3S] = {0.034, 0.0547,  0.208};
0194
0195  double erraa1[2][NCENT];
0196  double erraa2[2][NCENT];
0197  for(int i=0; i<NCENT;++i)
0198    {
0199      erraa1[0][i] = yrs13_erraa1[i];
0200      erraa2[0][i] = yrs13_erraa2[i];
0201      erraa1[1][i] = yrs15_erraa1[i];
0202      erraa2[1][i] = yrs15_erraa2[i];
0203    }
0204  double erraa3[2][NCENT3S];
0205  for(int i=0; i<NCENT3S;++i)
0206    {
0207      erraa3[0][i] = yrs13_erraa3[i];
0208      erraa3[1][i] = yrs15_erraa3[i];
0209    }
0210
0211
0212  TCanvas *craa = new TCanvas("craa","craa",0,0,1000,800);
0213  TH1D *h = new TH1D("h","h",100,-10,400);
0214  h->GetXaxis()->SetTitle("N_{part}");
0215  h->GetYaxis()->SetTitle("R_{AA}");
0216  h->SetMaximum(1.2);
0217  h->SetMinimum(0);
0218
0219  h->Draw();
0220
0221  TGraphErrors *gr1[2];
0222  TGraphErrors *gr2[2];
0223  TGraphErrors *gr3[2];
0224  for(int i=0; i< 2;i++)
0225    {
0226      gr1[i] = new TGraphErrors(NCENT,Npart_plot[i], raa1[i], 0, erraa1[i]);
0227      gr1[i]->SetMarkerStyle(20);
0228      gr1[i]->SetMarkerSize(1.5);
0229
0230      gr2[i] = new TGraphErrors(NCENT,Npart_plot[i], raa2[i], 0, erraa2[i]);
0231      gr2[i]->SetMarkerStyle(21);
0232      gr2[i]->SetMarkerSize(1.5);
0233      gr2[i]->SetMarkerColor(kRed);
0234      gr2[i]->SetLineColor(kRed);
0235
0236      gr3[i] = new TGraphErrors(NCENT3S,Npart3S_plot[i], raa3[i], 0, erraa3[i]);
0237      gr3[i]->SetMarkerStyle(20);
0238      gr3[i]->SetMarkerSize(1.5);
0239      gr3[i]->SetMarkerColor(kBlue);
0240      gr3[i]->SetLineColor(kBlue);
0241    }
0242
0243  if(yrs13)
0244    {
0245      gr1[0]->Draw("p");
0246      gr2[0]->Draw("p");
0247      if(ups3s) gr3[0]->Draw("p");
0248    }
0249
0250  if(yrs15)
0251    {
0252      if(yrs13 && yrs15)
0253      {
0254        gr1[0]->SetMarkerStyle(24);
0255        gr2[0]->SetMarkerStyle(25);
0256        gr3[0]->SetMarkerStyle(25);
0257      }
0258      gr1[1]->Draw("p");
0259      gr2[1]->Draw("p");
0260      if(ups3s) gr3[1]->Draw("p");
0261    }
0262
0263
0264  TLegend *lups = new TLegend(0.38, 0.72, 0.72, 0.87,"#bf{#it{sPHENIX}} Projection");
0265  lups->SetBorderSize(0);
0266  lups->SetFillColor(0);
0267  lups->SetFillStyle(0);
0268  if(yrs13 && yrs15)
0269    {
0270      lups->AddEntry(gr1[1],"Y(1S)","p");
0271      lups->AddEntry(gr2[1],"Y(2S)","p");
0272      if(ups3s) lups->AddEntry(gr3[1],"Y(3S)","p");
0273    }
0274  else
0275    {
0276      lups->AddEntry(gr1[0],"Y(1S)","p");
0277      lups->AddEntry(gr2[0],"Y(2S)","p");
0278      if(ups3s) lups->AddEntry(gr3[0],"Y(3S)","p");
0279    }
0280  lups->Draw();
0281
0282  TLatex *lat1[2];
0283  lat1[0] = new TLatex(0.66, 0.833,"#font[42]{Au+Au, Years 1-3}");
0284  lat1[0]->SetTextSize(0.038);
0285  lat1[0]->SetNDC(1);
0286  lat1[1] = new TLatex(0.66, 0.833,"#font[42]{Au+Au, Years 1-5}");
0287  lat1[1]->SetTextSize(0.038);
0288  lat1[1]->SetNDC(1);
0289
0290  if(yrs13 && !yrs15)
0291    lat1[0]->Draw();
0292  if(yrs15 || (yrs13 && yrs15))
0293    lat1[1]->Draw();
0294
0295  TLatex *lat2[2];
0296  //lat2[0] = new TLatex(0.66, 0.763,"#font[42]{#splitline{62 pb^{-1} trig. p+p}{142B rec. Au+Au}}");
0297  lat2[0] = new TLatex(0.65, 0.763,"#font[42]{#splitline{62 pb^{-1} trig. p+p}{21 nb^{-1} rec. Au+Au}}");
0298  lat2[0]->SetTextSize(0.038);
0299  lat2[0]->SetNDC(1);
0300  //lat2[1] = new TLatex(0.66, 0.763,"#font[42]{#splitline{142 pb^{-1} trig. p+p}{348B rec. Au+Au}}");
0301  lat2[1] = new TLatex(0.65, 0.763,"#font[42]{#splitline{142 pb^{-1} trig. p+p}{51 nb^{-1} rec. Au+Au}}");
0302  lat2[1]->SetTextSize(0.038);
0303  lat2[1]->SetNDC(1);
0304
0305  if(yrs13 && !yrs15)
0306    lat2[0]->Draw();
0307
0308  if(yrs15 || (yrs13 && yrs15))
0309    lat2[1]->Draw();
0310
0311  TLatex *lat4 = new TLatex(0.57,0.625,"#font[42]{#it{open markers: Years 1-3}}");
0312  lat4->SetTextSize(0.038);
0313  lat4->SetNDC(1);
0314  if(yrs15 || (yrs13 && yrs15))
0315    lat4->Draw();
0316
0317   // Plot the theory curves also
0318   double col[3] = {kBlack, kRed, kBlue};
0319   int lstyle[3] = {kDotted, kSolid, kDashDotted};
0320   int nstates = 2;
0321   if(ups3s) nstates = 3;
0322   for(int is = 0;is<nstates;is++)
0323     {
0324       for(int ieta =0; ieta < 3; ++ieta)
0325       {
0326     grth[is][ieta]->SetLineColor(col[is]);
0327     grth[is][ieta]->SetMarkerSize(0.5);
0328     grth[is][ieta]->SetMarkerColor(col[is]);
0329     grth[is][ieta]->SetMarkerStyle(20);
0330     grth[is][ieta]->SetLineStyle(lstyle[ieta]);
0331     grth[is][ieta]->SetLineWidth(2);
0332     grth[is][ieta]->Draw("l");
0333       }
0334     }
0335
0336   TLatex *lat3 = new TLatex(0.4, 0.675,"#font[42]{Strickland,Bazow, N.P. A879 (2012) 25}");
0337   lat3->SetTextSize(0.038);
0338   lat3->SetNDC(1);
0339   lat3->Draw();
0340
0341   TLegend *letas = new TLegend(0.14, 0.14, 0.31, 0.28,"");
0342   letas->SetBorderSize(0);
0343   letas->SetFillColor(0);
0344   letas->SetFillStyle(0);
0345   letas->AddEntry(grth[0][0], "4#pi#eta/S=1","l");
0346   letas->AddEntry(grth[0][1], "4#pi#eta/S=2","l");
0347   letas->AddEntry(grth[0][2], "4#pi#eta/S=3","l");
0348   letas->Draw();
0349 }
0350
0351
0352