Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-03 08:15:29

0001 double CBFunction(double *x, double *p)
0002 {
0003   double norm = p[0];
0004   double alpha = p[1];  // tail start
0005   double n = p[2];      // tail shape 
0006   double mu = p[3];     // upsilon mass
0007   double sigma = p[4];  // upsilon width
0008 
0009   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0010   double B = n/fabs(alpha) - fabs(alpha);
0011   double k = (x[0]-mu)/sigma;
0012 
0013   double val;
0014   if( k > -alpha )
0015     val = norm*TMath::Exp(-0.5*pow(k,2));
0016   else
0017     val = norm*A*pow(B-k,-n);
0018 
0019   if( TMath::IsNaN(val) ) val = 0.0;
0020 
0021   return val;
0022 }
0023 
0024 double TripleCBFunction(double *x, double *p)
0025 {
0026   double norm1 = p[0];  // amplitude of Y(1S) peak
0027   double alpha = p[1];  // tail start
0028   double n = p[2];      // tail shape 
0029   double mu1 = p[3];    // upsilon (1S) mass
0030   double sigma = p[4];  // upsilon width
0031   double norm2 = p[5];
0032   double norm3 = p[6];
0033   double mu2 = mu1*1.0595;
0034   double mu3 = mu1*1.0946;
0035 
0036   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0037   double B = n/fabs(alpha) - fabs(alpha);
0038   double k1 = (x[0]-mu1)/sigma;
0039   double k2 = (x[0]-mu2)/sigma;
0040   double k3 = (x[0]-mu3)/sigma;
0041 
0042   double val,val1,val2,val3;
0043 
0044   if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
0045   else { val1 = norm1*A*pow(B-k1,-n); }
0046   if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
0047   else { val2 = norm2*A*pow(B-k2,-n); }
0048   if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
0049   else { val3 = norm3*A*pow(B-k3,-n); }
0050 
0051   val = val1 + val2 + val3;
0052 
0053   if( TMath::IsNaN(val) ) val = 0.0;
0054 
0055   return val;
0056 }
0057 
0058 double SandB_CBFunction(double *x, double *p)
0059 {
0060   double norm1 = p[0];  // amplitude of Y(1S) peak
0061   double alpha = p[1];  // tail start
0062   double n = p[2];      // tail shape 
0063   double mu1 = p[3];     // upsilon (1S) mass
0064   double sigma = p[4];  // upsilon width
0065   double norm2 = p[5];
0066   double norm3 = p[6];
0067   double mu2 = mu1*1.0595;
0068   double mu3 = mu1*1.0946;
0069 
0070   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0071   double B = n/fabs(alpha) - fabs(alpha);
0072   double k1 = (x[0]-mu1)/sigma;
0073   double k2 = (x[0]-mu2)/sigma;
0074   double k3 = (x[0]-mu3)/sigma;
0075 
0076   double val,val1,val2,val3;
0077 
0078   if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
0079   else { val1 = norm1*A*pow(B-k1,-n); }
0080   if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
0081   else { val2 = norm2*A*pow(B-k2,-n); }
0082   if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
0083   else { val3 = norm3*A*pow(B-k3,-n); }
0084 
0085   double bgnorm1 = p[7];
0086   double bgslope1 = p[8];
0087 
0088   double bg = exp(bgnorm1+x[0]*bgslope1);
0089 
0090   val = val1 + val2 + val3 + bg;
0091   if( TMath::IsNaN(val) ) val = 0.0;
0092 
0093   return val;
0094 }
0095 
0096 //----------------------------------------------------------------
0097 //----------------------------------------------------------------
0098 //----------------------------------------------------------------
0099 
0100 void newplotbg_newsupp_2022_60pct() {
0101 
0102 //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
0103 //SetsPhenixStyle();
0104 
0105 gStyle->SetOptStat(0);
0106 gStyle->SetOptFit(0);
0107 
0108 TRandom* myrandom = new TRandom3();
0109 const int nbins = 15;
0110 double eideff = 0.9;
0111 string times = "*";
0112 TLatex* tl[15];
0113 char tlchar[999];
0114 
0115 double statscale = 1.4; // from 10B to 14B
0116 double statscale_lowlim = 7.0;
0117 double statscale_uplim = 14.0;
0118 
0119 double ppcorr = (2400./14.)/962.; // 2400B pp and 140B MB AuAu
0120 
0121   TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
0122   TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
0123   TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
0124   TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
0125   TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
0126   TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
0127   TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
0128   TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
0129   TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
0130   TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
0131   TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
0132   TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
0133 
0134 //---------------------------------------------------------
0135 // Upsilons
0136 //---------------------------------------------------------
0137 
0138 int whichsuppression = 1; // 0=CMS suppression, 1=Strickland suppression
0139 
0140 int npts1 = 10;
0141 int npts2 = 8;
0142 int npts3 = 7;
0143 int npts3_rebin = 1;
0144 if(whichsuppression==0) {npts3_rebin=4;}
0145 
0146 string str_UpsilonPt  = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))";   // dN/dpT
0147 string str_UpsilonXPt = "(2.0*3.14159*x*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))"; // need this for mean pT calculation
0148 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
0149 TF1* fUpsilonXPt = new TF1("fUpsilonXPt",str_UpsilonXPt.c_str(),0.,20.);
0150 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
0151 fUpsilonXPt->SetParameters(72.1, 26.516, 10.6834);
0152 double upsnorm = fUpsilonPt->Integral(0.,20.);
0153 
0154 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.5 + double(i);}
0155 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.5 + double(i);}
0156 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.5 + double(i);}
0157 double xx3_rebin[nbins+1];
0158 if(npts3_rebin==2) {
0159   xx3_rebin[0] = fUpsilonXPt->Integral(0.,4.)/fUpsilonPt->Integral(0.,4.);
0160   xx3_rebin[1] = fUpsilonXPt->Integral(4.,10.)/fUpsilonPt->Integral(4.,10.);
0161 }  else if(npts3_rebin==4) {
0162   xx3_rebin[0] = fUpsilonXPt->Integral(0.,2.)/fUpsilonPt->Integral(0.,2.);
0163   xx3_rebin[1] = fUpsilonXPt->Integral(2.,4.)/fUpsilonPt->Integral(2.,4.);
0164   xx3_rebin[2] = fUpsilonXPt->Integral(4.,6.)/fUpsilonPt->Integral(4.,6.);
0165   xx3_rebin[3] = fUpsilonXPt->Integral(6.,10.)/fUpsilonPt->Integral(6.,10.);
0166 } else if(npts3_rebin==1) {
0167   xx3_rebin[0] = fUpsilonXPt->Integral(0.,10.)/fUpsilonPt->Integral(0.,10.);
0168 } else {cerr << "ERROR: Wrong Y(3S) binning." << endl; return;}
0169 
0170 
0171 // Raa for 0-60% centrality bin
0172 double raa1s[nbins+1];
0173 double raa2s[nbins+1];
0174 double raa3s[nbins+1];
0175 for(int i=0; i<nbins+1; i++) {raa1s[i]=0.654; raa2s[i]=0.277; raa3s[i]=0.075;}
0176 if(whichsuppression==0) {for(int i=0; i<nbins+1; i++) {raa3s[i]=0.138;}}
0177 
0178 TGraph* grRAA1S = new TGraph(nbins,xx1,raa1s);
0179 TGraph* grRAA2S = new TGraph(nbins,xx2,raa2s);
0180 TGraph* grRAA3S = new TGraph(nbins,xx3,raa3s);
0181 
0182 int nchan=400;
0183 double start=0.0;
0184 double stop=20.0;
0185 
0186 double frac[3];  // upsilons fractions
0187   frac[0] = 0.7117;
0188   frac[1] = 0.1851;
0189   frac[2] = 0.1032;
0190 double scale[3]; // mass scaling
0191   scale[0] = 1.0;
0192   scale[1] = 1.0595;
0193   scale[2] = 1.0946;
0194 
0195 //Npart for 0-60% centrality = 160.2
0196 //Ncoll for 0-60% centrality = 388.8
0197 //Ncoll for 0-10% centrality = 955.
0198 //Npionpairs for 0-60% centrality = 0.406667
0199 //Number of Upsilons in 0-60% centrality = 19405 1920 549
0200 double combbg_corr_60pct = 0.406667 * 6.;
0201 double corrbg_corr_60pct = 388.8/955. * 6.;
0202 
0203 //---------- Use fixed numbers for 2022 estimate --------------------------------------
0204 // Numbers from Marzia for 140 B MB Au+Au and 2400 B p+p events
0205 // RAA for Y(3S) is 8.5% Number of Upsilons for 0-60% centrality bin
0206   int Nups1   = 19405;
0207   int Nups2   = 1920;
0208   int Nups3   = 549; // with CMS suppression
0209   if(whichsuppression==1) {Nups3   = 549*(3.5/8.5);} // with Strickland's suppression
0210   int Nups23 = Nups2+Nups3;
0211   int Nups1pp = 2.86e+03;
0212   int Nups2pp = 7.16e+02;
0213   int Nups3pp = 3.98e+02;
0214   int Nups23pp = Nups2pp+Nups3pp;
0215 
0216 //-------------------------------------------------------------------------------------
0217 
0218 //====================================================================================
0219 //====================================================================================
0220 //====================================================================================
0221 
0222 //
0223 // these histograms (hhups*) are randomly generated using the fit above to single 
0224 // peak (fCB function) and used for the RAA uncertainty calculation
0225 //
0226 // FORCE specific RESOLUTION and tail shape
0227   double tonypar1 = 0.98;  // Tony's 100 pion simulation April 2019
0228   double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
0229   //double tonypar3 = 9.437;  // Tony's 100 pion simulation April 2019
0230   double tonypar3 = 9.448;  // benchmark
0231   double tonypar4 = 0.100;  // benchmark
0232   double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
0233   fCBpp->SetParameter(0,1000.); 
0234   fCBpp->SetParameter(1,tonypar1); 
0235   fCBpp->SetParameter(2,tonypar2); 
0236   fCBpp->SetParameter(3,tonypar3); 
0237   fCBpp->SetParameter(4,tonypar4pp);  
0238   fCBauau->SetParameter(0,1000.); 
0239   fCBauau->SetParameter(1,tonypar1); 
0240   fCBauau->SetParameter(2,tonypar2); 
0241   fCBauau->SetParameter(3,tonypar3); 
0242   fCBauau->SetParameter(4,tonypar4);  
0243 
0244 
0245 char hhname[999];
0246 TH1D* hhups[nbins+1];
0247 TH1D* hhups1[nbins+1];
0248 TH1D* hhups2[nbins+1];
0249 TH1D* hhups3[nbins+1];
0250 TH1D* hhupspp[nbins+1];
0251 TH1D* hhups1pp[nbins+1];
0252 TH1D* hhups2pp[nbins+1];
0253 TH1D* hhups3pp[nbins+1];
0254 for(int i=0; i<nbins+1; i++) {
0255   sprintf(hhname,"hhups_%d",i);
0256   hhups[i] = new TH1D(hhname,"",nchan,start,stop);
0257   hhups[i]->Sumw2();
0258   sprintf(hhname,"hhups1_%d",i);
0259   hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
0260   hhups1[i]->Sumw2();
0261   sprintf(hhname,"hhups2_%d",i);
0262   hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
0263   hhups2[i]->Sumw2();
0264   sprintf(hhname,"hhups3_%d",i);
0265   hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
0266   hhups3[i]->Sumw2();
0267   hhups[i]->SetLineWidth(2);
0268   hhups1[i]->SetLineWidth(2);
0269   hhups2[i]->SetLineWidth(2);
0270   hhups3[i]->SetLineWidth(2);
0271   sprintf(hhname,"hhupspp_%d",i);
0272   hhupspp[i] = new TH1D(hhname,"",nchan,start,stop);
0273   hhupspp[i]->Sumw2();
0274   sprintf(hhname,"hhups1pp_%d",i);
0275   hhups1pp[i] = new TH1D(hhname,"",nchan,start,stop);
0276   hhups1pp[i]->Sumw2();
0277   sprintf(hhname,"hhups2pp_%d",i);
0278   hhups2pp[i] = new TH1D(hhname,"",nchan,start,stop);
0279   hhups2pp[i]->Sumw2();
0280   sprintf(hhname,"hhups3pp_%d",i);
0281   hhups3pp[i] = new TH1D(hhname,"",nchan,start,stop);
0282   hhups3pp[i]->Sumw2();
0283   hhupspp[i]->SetLineWidth(2);
0284   hhups1pp[i]->SetLineWidth(2);
0285   hhups2pp[i]->SetLineWidth(2);
0286   hhups3pp[i]->SetLineWidth(2);
0287 }
0288 
0289 for(int j=0; j<nbins; j++) {
0290 
0291   double s1 = j*1.0;
0292   double s2 = s1 + 1.0;
0293   double tmpnups1 = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;
0294   double tmpnups2 = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;
0295   double tmpnups3 = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;
0296   fCBauau->SetParameter(3,tonypar3);
0297     for(int i=0; i<int(tmpnups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0298   fCBauau->SetParameter(3,tonypar3*scale[1]);
0299     for(int i=0; i<int(tmpnups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0300   fCBauau->SetParameter(3,tonypar3*scale[2]);
0301     for(int i=0; i<int(tmpnups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0302 
0303   double tmpnups1pp = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0304   double tmpnups2pp = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0305   double tmpnups3pp = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0306   fCBpp->SetParameter(3,tonypar3);
0307     for(int i=0; i<int(tmpnups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0308   fCBpp->SetParameter(3,tonypar3*scale[1]);
0309     for(int i=0; i<int(tmpnups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0310   fCBpp->SetParameter(3,tonypar3*scale[2]);
0311     for(int i=0; i<int(tmpnups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0312 
0313 } // end loop over pT bins
0314 
0315 // all pT
0316 
0317   fCBpp->SetParameter(3,tonypar3);
0318   fCBauau->SetParameter(3,tonypar3);
0319     for(int i=0; i<int(Nups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0320     for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0321   fCBpp->SetParameter(3,tonypar3*scale[1]);
0322   fCBauau->SetParameter(3,tonypar3*scale[1]);
0323     for(int i=0; i<int(Nups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0324     for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0325   fCBpp->SetParameter(3,tonypar3*scale[2]);
0326   fCBauau->SetParameter(3,tonypar3*scale[2]);
0327     for(int i=0; i<int(Nups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0328     for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0329 
0330 // fit and draw pp no bg
0331 
0332 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
0333   fTCBpp->SetParameter(0,2000.);
0334   fTCBpp->FixParameter(1,tonypar1);
0335   fTCBpp->FixParameter(2,tonypar2);
0336   fTCBpp->SetParameter(3,tonypar3);
0337   fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
0338   fTCBpp->SetParameter(5,500.);
0339   fTCBpp->SetParameter(6,100.);
0340   hhupspp[nbins]->Fit(fTCBpp,"rl","",7.,11.);
0341   hhupspp[nbins]->SetAxisRange(7.,11.);
0342   hhupspp[nbins]->SetMarkerSize(1.0);
0343   hhupspp[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0344   hhupspp[nbins]->GetXaxis()->SetTitleOffset(1.0);
0345   double tmpamp1 = hhupspp[nbins]->GetFunction("fTCBpp")->GetParameter(0);
0346   double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0347   double tmpamp6 = tmpamp1*frac[2]/frac[0];
0348   hhupspp[nbins]->Draw();
0349 
0350   fCB1s->SetLineColor(kBlue);
0351   fCB1s->SetLineWidth(1);
0352     fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
0353     fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
0354     fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
0355     fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
0356     fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
0357   fCB2s->SetLineColor(kRed);
0358   fCB2s->SetLineWidth(1);
0359     fCB2s->SetParameter(0,tmpamp5);
0360     fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
0361     fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
0362     fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
0363     fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
0364   fCB3s->SetLineColor(kGreen+2);
0365   fCB3s->SetLineWidth(1);
0366     fCB3s->SetParameter(0,tmpamp6);
0367     fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
0368     fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
0369     fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
0370     fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
0371   fCB1s->Draw("same");
0372   fCB2s->Draw("same");
0373   fCB3s->Draw("same");
0374 
0375 // fit and draw AuAu no bg
0376 
0377 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Au+Au",100,100,600,600);
0378   fTCBauau->SetParameter(0,2000.);
0379   fTCBauau->FixParameter(1,tonypar1);
0380   fTCBauau->FixParameter(2,tonypar2);
0381   fTCBauau->SetParameter(3,tonypar3);
0382   fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
0383   fTCBauau->SetParameter(5,500.);
0384   fTCBauau->SetParameter(6,100.);
0385   hhups[nbins]->Fit(fTCBauau,"rl","",7.,11.);
0386   hhups[nbins]->SetAxisRange(7.,11.);
0387   hhups[nbins]->SetMarkerSize(1.0);
0388   hhups[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0389   hhups[nbins]->GetXaxis()->SetTitleOffset(1.0);
0390   tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
0391   tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0392   tmpamp6 = tmpamp1*frac[2]/frac[0];
0393   hhups[nbins]->Draw();
0394 
0395 TCanvas* cupsvspt = new TCanvas("cupsvspt","Upsilons vs. p_{T}",100,100,1200,900);
0396 cupsvspt->Divide(4,3);
0397 for(int i=0; i<12; i++) {
0398   cupsvspt->cd(i+1);
0399   hhups[i]->SetAxisRange(7.0,11.0); hhups[i]->Draw();
0400   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhups[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0401 }
0402 
0403 //----------------------------------------------------
0404 //  Backgrounds
0405 //----------------------------------------------------
0406 
0407 TH1D* hhall[nbins+1];
0408 TH1D* hhall_scaled[nbins+1];
0409 
0410 TH1D* hhtotbg[nbins+1]; 
0411 TH1D* hhtotbg_scaled[nbins+1]; 
0412 TH1D* hhcombbg[nbins+1];
0413 TH1D* hhcombbg_scaled[nbins+1];
0414 TH1D* hhfakefake[nbins+1];
0415 TH1D* hhfakehf[nbins+1];
0416 TH1D* hhbottom[nbins+1];
0417 TH1D* hhcharm[nbins+1];
0418 TH1D* hhdy[nbins+1];
0419 TH1D* hhcorrbg[nbins+1];
0420 TH1D* hhcorrbg_scaled[nbins+1];
0421 TH1D* hhfit[nbins+1];
0422 char tmpname[999];
0423 
0424 //----------------------------------------------------------------------------------------
0425 // Correlated BG calculated for 10B central AuAu events
0426 //----------------------------------------------------------------------------------------
0427 
0428 double corrbgfitpar0;
0429 double corrbgfitpar1;
0430 
0431 TFile* f=new TFile("ccbb_eideff09.root");
0432 for(int i=0; i<nbins+1; i++) {
0433   sprintf(tmpname,"hhbottom_%d",i);
0434   hhbottom[i] = (TH1D*)f->Get(tmpname);
0435   hhbottom[i]->SetDirectory(gROOT);
0436   sprintf(tmpname,"hhcharm_%d",i);
0437   hhcharm[i] = (TH1D*)f->Get(tmpname);
0438   hhcharm[i]->SetDirectory(gROOT);
0439   sprintf(tmpname,"hhdy_%d",i);
0440   hhdy[i] = (TH1D*)f->Get(tmpname);
0441   hhdy[i]->SetDirectory(gROOT);
0442   sprintf(tmpname,"hhcorrbg_%d",i);
0443   hhcorrbg[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0444   hhcorrbg[i]->Add(hhcharm[i]);
0445   hhcorrbg[i]->Add(hhdy[i]);
0446   sprintf(tmpname,"hhcorrbg_scaled_%d",i);
0447   hhcorrbg_scaled[i] = (TH1D*)hhcorrbg[i]->Clone(tmpname);
0448   hhcorrbg[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0449   hhbottom[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0450   hhdy[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0451     if(i==nbins) {
0452        corrbgfitpar0 = hhcorrbg[i]->GetFunction("expo")->GetParameter(0);
0453        corrbgfitpar1 = hhcorrbg[i]->GetFunction("expo")->GetParameter(1);
0454     }
0455     for(int k=1; k<=hhcorrbg[i]->GetNbinsX(); k++) {
0456       if(hhcorrbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[i]->GetBinLowEdge(k)+hhcorrbg[i]->GetBinWidth(k))>statscale_uplim) {
0457         hhcorrbg_scaled[i]->SetBinContent(k,0.); 
0458         hhcorrbg_scaled[i]->SetBinError(k,0.);
0459       }
0460       else {
0461         double tmp = corrbg_corr_60pct * statscale * hhcorrbg[i]->GetFunction("expo")->Eval(hhcorrbg[i]->GetBinCenter(k));
0462         double tmprnd = myrandom->Poisson(tmp); 
0463         if(tmprnd<0.) { tmprnd=0.; }
0464         hhcorrbg_scaled[i]->SetBinContent(k,tmprnd); 
0465         hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0466       }
0467     }
0468   hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0469   hhcorrbg[i]->SetDirectory(gROOT);
0470   hhcorrbg_scaled[i]->SetDirectory(gROOT);
0471 }
0472 f->Close();
0473 
0474 TCanvas* c2 = new TCanvas("c2","Correlated BG (all p_{T})",100,100,600,600);
0475   
0476 hhbottom[nbins]->SetAxisRange(7.0,14.0);
0477 hhbottom[nbins]->SetLineColor(kBlue);
0478 hhbottom[nbins]->SetLineWidth(2);
0479 
0480 hhbottom[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0481 hhbottom[nbins]->GetXaxis()->SetTitleOffset(1.0);
0482 hhbottom[nbins]->GetXaxis()->SetTitleColor(1);
0483 hhbottom[nbins]->GetXaxis()->SetTitleSize(0.040);
0484 hhbottom[nbins]->GetXaxis()->SetLabelSize(0.040);
0485 //hhbottom[nbins]->GetYaxis()->SetTitle("Correlated background");
0486 hhbottom[nbins]->GetYaxis()->SetTitleOffset(1.3);
0487 hhbottom[nbins]->GetYaxis()->SetTitleSize(0.040);
0488 hhbottom[nbins]->GetYaxis()->SetLabelSize(0.040);
0489 
0490 hhcorrbg[nbins]->SetAxisRange(7.0,14.0);
0491 hhcorrbg[nbins]->SetMinimum(0.1);
0492 hhcorrbg[nbins]->SetLineColor(kBlack);
0493 hhcorrbg[nbins]->SetLineWidth(2);
0494 hhcorrbg[nbins]->Draw("hist");
0495 
0496 hhbottom[nbins]->SetMarkerColor(kBlue);
0497 hhbottom[nbins]->SetLineColor(kBlue);
0498 hhbottom[nbins]->Draw("same");
0499 
0500 hhdy[nbins]->SetMarkerColor(kGreen+2);
0501 hhdy[nbins]->SetLineColor(kGreen+2);
0502 hhdy[nbins]->SetLineWidth(2);
0503 hhdy[nbins]->Draw("same");
0504 
0505 hhcharm[nbins]->SetMarkerColor(kRed);
0506 hhcharm[nbins]->SetLineColor(kRed);
0507 hhcharm[nbins]->SetLineWidth(2);
0508 hhcharm[nbins]->Draw("same");
0509 
0510 TCanvas* c0 = new TCanvas("c0","Correlated BG vs. p_{T} 10B events",100,100,1200,900);
0511 c0->Divide(4,3);
0512 for(int i=0; i<nbins; i++) {
0513   c0->cd(i+1);
0514   hhcorrbg[i]->SetAxisRange(7.,14.);
0515   hhcorrbg[i]->GetFunction("expo")->SetLineColor(kBlack); 
0516   hhbottom[i]->GetFunction("expo")->SetLineColor(kBlue); 
0517   hhdy[i]->GetFunction("expo")->SetLineColor(kGreen+2); 
0518   hhcorrbg[i]->SetLineColor(kBlack); 
0519   hhbottom[i]->SetLineColor(kBlue); 
0520   hhdy[i]->SetLineColor(kGreen+2); 
0521   hhcorrbg[i]->Draw();
0522   hhbottom[i]->Draw("same");
0523   hhdy[i]->Draw("same");
0524   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0525 }
0526 
0527 TCanvas* c0scaled = new TCanvas("c0scaled","SCALED Correlated BG vs. p_{T}",100,100,1200,900);
0528 c0scaled->Divide(4,3);
0529 for(int i=0; i<nbins; i++) {
0530   c0scaled->cd(i+1);
0531   hhcorrbg_scaled[i]->SetAxisRange(7.0,14.0); hhcorrbg_scaled[i]->Draw();
0532   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0533 }
0534 
0535 //-----------------------------------------------------------------------------------------
0536 //  Correlated bg in p+p
0537 //-----------------------------------------------------------------------------------------
0538 
0539 TH1D* hhbottom_pp[nbins+1];
0540 TH1D* hhdy_pp[nbins+1];
0541 TH1D* hhcorrbg_pp[nbins+1];
0542 TH1D* hhall_pp[nbins+1];
0543 
0544 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
0545 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
0546 
0547 for(int i=0; i<nbins+1; i++) {
0548 
0549   sprintf(tmpname,"hhbottom_pp_%d",i);
0550   hhbottom_pp[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0551   for(int k=1; k<=hhbottom_pp[i]->GetNbinsX(); k++) {
0552     if(hhbottom_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp[i]->GetBinLowEdge(k)+hhbottom_pp[i]->GetBinWidth(k))>statscale_uplim) {
0553       hhbottom_pp[i]->SetBinContent(k,0.);
0554       hhbottom_pp[i]->SetBinError(k,0.);
0555     }
0556     else {
0557       double tmp = corrbg_corr_60pct * ppcorr * fbottom_nosup_corr->Eval(hhbottom[i]->GetBinCenter(k)) * hhbottom[i]->GetFunction("expo")->Eval(hhbottom[i]->GetBinCenter(k));
0558       double tmprnd = myrandom->Poisson(tmp);
0559       if(tmprnd<0.) { tmprnd=0.; }
0560       hhbottom_pp[i]->SetBinContent(k,tmprnd);
0561       hhbottom_pp[i]->SetBinError(k,sqrt(tmprnd));
0562     }
0563   }
0564 
0565   sprintf(tmpname,"hhdy_pp_%d",i);
0566   hhdy_pp[i] = (TH1D*)hhdy[i]->Clone(tmpname);
0567   for(int k=1; k<=hhdy_pp[i]->GetNbinsX(); k++) {
0568     if(hhdy_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp[i]->GetBinLowEdge(k)+hhdy_pp[i]->GetBinWidth(k))>statscale_uplim) {
0569       hhdy_pp[i]->SetBinContent(k,0.);
0570       hhdy_pp[i]->SetBinError(k,0.);
0571     }
0572     else {
0573       double tmp = ppcorr * hhdy[i]->GetFunction("expo")->Eval(hhdy[i]->GetBinCenter(k));
0574       double tmprnd = myrandom->Poisson(tmp);
0575       if(tmprnd<0.) { tmprnd=0.; }
0576       hhdy_pp[i]->SetBinContent(k,tmprnd);
0577       hhdy_pp[i]->SetBinError(k,sqrt(tmprnd));
0578     }
0579   }
0580 
0581   sprintf(tmpname,"hhcorrbg_pp_%d",i);
0582   hhcorrbg_pp[i] = (TH1D*)hhbottom_pp[i]->Clone(tmpname);
0583   hhcorrbg_pp[i]->Add(hhdy_pp[i]);
0584     hhcorrbg_pp[i]->SetMarkerColor(kBlack);
0585     hhcorrbg_pp[i]->SetLineColor(kBlack);
0586     hhbottom_pp[i]->SetLineColor(kBlue);
0587     hhdy_pp[i]->SetLineColor(kGreen+2);
0588   sprintf(tmpname,"hhall_pp_%d",i);
0589   hhall_pp[i] = (TH1D*)hhcorrbg_pp[i]->Clone(tmpname);
0590   hhall_pp[i]->Add(hhupspp[i]);
0591     hhall_pp[i]->SetLineColor(kMagenta);
0592     hhall_pp[i]->SetMarkerColor(kMagenta);
0593 
0594   hhcorrbg_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0595   hhcorrbg_pp[i]->GetFunction("expo")->SetLineColor(kBlack);
0596   hhbottom_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0597   hhbottom_pp[i]->GetFunction("expo")->SetLineColor(kBlue);
0598   hhdy_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0599   hhdy_pp[i]->GetFunction("expo")->SetLineColor(kGreen+2);
0600 
0601 } // end loop over pT bins
0602 
0603 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
0604   hhall_pp[nbins]->SetAxisRange(7.,12.);
0605   hhcorrbg_pp[nbins]->Draw("pehist");
0606   hhbottom_pp[nbins]->Draw("histsame");
0607   hhdy_pp[nbins]->Draw("histsame");
0608 
0609 
0610 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
0611   hhall_pp[nbins]->SetAxisRange(7.,12.);
0612   hhall_pp[nbins]->Draw("pehist");
0613   hhcorrbg_pp[nbins]->Draw("pesame");
0614   hhbottom_pp[nbins]->Draw("same");
0615   hhdy_pp[nbins]->Draw("same");
0616 
0617 TCanvas* cpp_vspt = new TCanvas("cpp_vspt","corr bg + sig vs pt in pp",50,50,1200,900);
0618 cpp_vspt->Divide(4,3);
0619 for(int i=0; i<nbins; i++) {
0620   cpp_vspt->cd(i+1);
0621   hhall_pp[i]->SetAxisRange(7.0,12.0); hhall_pp[i]->Draw("hist"); hhcorrbg_pp[i]->Draw("histsame");
0622   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0623 }
0624 
0625 //-----------------------------------------------------------------------------------------
0626 // Combinatorial BG calculated for 10B central AuAu events
0627 //-----------------------------------------------------------------------------------------
0628 TCanvas* cdummy = new TCanvas("cdummy","cdummy",0,0,500,500);
0629 
0630 f = new TFile("fakee_eideff09.root");
0631 for(int i=0; i<nbins+1; i++) {
0632   sprintf(tmpname,"hhfakefake_%d",i);
0633   hhfakefake[i] = (TH1D*)f->Get(tmpname);
0634   hhfakefake[i]->SetDirectory(gROOT);
0635 }
0636 f->Close();
0637 
0638 f = new TFile("crossterms_eideff09.root");
0639 for(int i=0; i<nbins+1; i++) {
0640   sprintf(tmpname,"hhfakehf_%d",i);
0641   hhfakehf[i] = (TH1D*)f->Get(tmpname);
0642   hhfakehf[i]->SetDirectory(gROOT);
0643 }
0644 f->Close();
0645 
0646 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
0647 fbg->SetParameters(10., -1.0, 4., -0.1);
0648 fbg->SetParLimits(1.,-999.,0.);
0649 fbg->SetParLimits(3.,-999.,0.);
0650 
0651 for(int i=0; i<nbins+1; i++) {
0652   sprintf(tmpname,"hhcombbg_%d",i);
0653   hhcombbg[i] = (TH1D*)hhfakefake[i]->Clone(tmpname);
0654   hhcombbg[i]->Add(hhfakehf[i]);
0655   sprintf(tmpname,"hhcombbg_scaled_%d",i);
0656   hhcombbg_scaled[i] = (TH1D*)hhcombbg[i]->Clone(tmpname);
0657   if(i==nbins) { fbg->SetParameters(10., -1.0, 4., -0.1); }
0658   hhcombbg[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0659 
0660     for(int k=1; k<=hhcombbg[i]->GetNbinsX(); k++) {
0661       if(hhcombbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[i]->GetBinLowEdge(k)+hhcombbg[i]->GetBinWidth(k))>statscale_uplim) {
0662         hhcombbg_scaled[i]->SetBinContent(k,0.);
0663         hhcombbg_scaled[i]->SetBinError(k,0.);
0664       }
0665       else {
0666         double tmp = combbg_corr_60pct * statscale * hhcombbg[i]->GetFunction("fbg")->Eval(hhcombbg[i]->GetBinCenter(k));
0667         double tmprnd = myrandom->Poisson(tmp);
0668         if(tmprnd<0.) { tmprnd=0.; }
0669         hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
0670         hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0671       }
0672     }
0673   hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0674 }
0675 
0676 delete cdummy;
0677 
0678 TCanvas* C1 = new TCanvas("C1","Combinatorial BG (ALL p_{T})",100,100,600,600);
0679 C1->SetLogy();
0680   hhfakefake[nbins]->SetAxisRange(7.0,14.0);
0681   hhfakefake[nbins]->SetMinimum(0.1);
0682   hhfakefake[nbins]->SetMaximum(5000.);
0683   hhfakefake[nbins]->SetLineColor(kGreen+2);
0684   hhfakefake[nbins]->SetLineWidth(2);
0685   hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
0686   hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
0687   hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
0688   hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
0689   hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
0690   hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
0691   hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
0692   hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
0693   hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
0694   hhfakefake[nbins]->Draw("e");
0695 
0696   hhfakehf[nbins]->SetLineColor(kOrange+4);
0697   hhfakehf[nbins]->SetLineWidth(2);
0698   hhfakehf[nbins]->Draw("esame");
0699 
0700   hhcombbg[nbins]->SetLineColor(kBlack);
0701   hhcombbg[nbins]->SetLineWidth(2);
0702   hhcombbg[nbins]->Draw("esame");
0703 
0704 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG (ALL p_{T})",100,100,600,600);
0705 C1sc->SetLogy(); 
0706   hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
0707   hhcombbg_scaled[nbins]->Draw("esame");
0708 
0709 TCanvas* c00 = new TCanvas("c00","Combinatorial BG vs. p_{T}",150,150,1200,900);
0710 c00->Divide(4,3);
0711 
0712 for(int i=0; i<nbins; i++) {
0713   c00->cd(i+1);
0714   hhcombbg[i]->SetAxisRange(7.0,14.0); hhcombbg[i]->Draw();
0715   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0716 }
0717 
0718 TCanvas* c00scaled = new TCanvas("c00scaled","SCALED Combinatorial BG vs. p_{T}",150,150,1200,900);
0719 c00scaled->Divide(4,3);
0720 
0721 for(int i=0; i<nbins; i++) {
0722   c00scaled->cd(i+1);
0723   hhcombbg_scaled[i]->SetAxisRange(statscale_lowlim,statscale_uplim); hhcombbg_scaled[i]->Draw();
0724   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0725 }
0726 
0727 //---------------------------------------------
0728 //  Signal + BG
0729 //---------------------------------------------
0730 
0731 for(int i=0; i<nbins+1; i++) {
0732   sprintf(tmpname,"hhtotbg_scaled_%d",i);
0733   hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
0734   hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
0735 }
0736 
0737 for(int i=0; i<nbins+1; i++) {
0738 sprintf(tmpname,"hhall_scaled_%d",i);
0739 hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
0740 hhall_scaled[i]->Add(hhups[i]);
0741 }
0742 
0743 TCanvas* c000 = new TCanvas("c000","Signal + Background vs. p_{T}",200,200,1200,900);
0744 c000->Divide(4,3);
0745 for(int i=0; i<nbins; i++) {
0746   c000->cd(i+1);
0747   hhall_scaled[i]->SetAxisRange(7.0,14.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
0748   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
0749 }
0750 
0751 // Fit Signal + Correlated BG
0752 
0753 // hhfit[] is signal + correlated bg
0754 for(int i=0; i<nbins+1; i++) {
0755   sprintf(tmpname,"hhfit_%d",i);
0756   hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
0757   for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
0758     hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
0759     hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
0760   }
0761 }
0762 
0763 //---------------------------------------------------------------------
0764 // plot signal + correlated bg for all pT
0765 //---------------------------------------------------------------------
0766 
0767   double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
0768   double tmppar1 = corrbgfitpar1;
0769 //cout << "###### " << tmppar0 << " " << tmppar1 << endl;
0770 
0771 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
0772 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
0773 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
0774 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
0775 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
0776 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
0777 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
0778 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
0779 fSandBpp->SetParameter(7,tmppar0);
0780 fSandBpp->SetParameter(8,tmppar1);
0781 fSandBpp->SetLineColor(kBlue);
0782 fSandBpp->SetLineStyle(2);
0783 
0784 //cout << "######### " << fTCBauau->GetParameter(0) << " " << fTCBauau->GetParameter(5) << " " << fTCBauau->GetParameter(6) << endl;
0785 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
0786 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
0787 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
0788 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
0789 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
0790 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
0791 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
0792 fSandBauau->SetParameter(7,tmppar0);
0793 fSandBauau->SetParameter(8,tmppar1);
0794 fSandBauau->SetLineColor(kRed);
0795 
0796 
0797 TCanvas* cfitall = new TCanvas("cfitall","FIT all pT",270,270,600,600);
0798   hhfit[nbins]->SetAxisRange(7.0,14.);
0799   hhfit[nbins]->GetXaxis()->CenterTitle();
0800   hhfit[nbins]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
0801   hhfit[nbins]->GetXaxis()->SetTitleOffset(1.1);
0802   hhfit[nbins]->GetXaxis()->SetLabelSize(0.045);
0803   hhfit[nbins]->GetYaxis()->CenterTitle();
0804   hhfit[nbins]->GetYaxis()->SetLabelSize(0.045);
0805   hhfit[nbins]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
0806   hhfit[nbins]->GetYaxis()->SetTitleOffset(1.5);
0807   hhfit[nbins]->Draw("pehist");
0808   fSandBpp->Draw("same");
0809   fSandBauau->Draw("same");
0810     //hhcorrbg_scaled[nbins]->SetLineColor(kRed);
0811     //hhcorrbg_scaled[nbins]->Scale(0.82);
0812     //hhcorrbg_scaled[nbins]->Scale(0.70);
0813     //hhcorrbg_scaled[nbins]->Draw("histsame");
0814 
0815   TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
0816   fmycorrbg->SetParameters(tmppar0,tmppar1);
0817   fmycorrbg->SetLineStyle(2);
0818   fmycorrbg->SetLineColor(kRed);
0819   fmycorrbg->Draw("same");
0820 
0821 
0822 double myheight = fTCBauau->GetParameter(0);
0823 TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
0824 ld1->SetTextSize(0.035);
0825 ld1->Draw();
0826 TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
0827 ld2->SetTextSize(0.035);
0828 ld2->Draw();
0829 
0830 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
0831 TH1D* hhfit_tmp = (TH1D*)hhfit[nbins]->Clone("hhfit_tmp");
0832 hhfit_tmp->SetAxisRange(8.0,11.);
0833 hhfit_tmp->Draw("pehist");
0834 fSandBauau->Draw("same");
0835 fmycorrbg->Draw("same");
0836 
0837 
0838 //----------------------------------------------------------------------
0839 
0840 // plot signal + both backgrounds for all pT
0841 
0842 TCanvas* callpt = new TCanvas("callpt","Signal+BG (all p_{T})",300,300,600,600);
0843 
0844   hhall_scaled[nbins]->GetXaxis()->SetTitle("Invariant mass GeV/c");
0845   hhall_scaled[nbins]->SetLineColor(kBlack);
0846   hhall_scaled[nbins]->SetMarkerColor(kBlack);
0847   hhall_scaled[nbins]->SetMarkerStyle(20);
0848   hhall_scaled[nbins]->SetAxisRange(8.0,10.8);
0849   hhall_scaled[nbins]->Draw("pehist");
0850     hhcombbg_scaled[nbins]->SetLineColor(kBlue);
0851     hhcombbg_scaled[nbins]->Draw("histsame");
0852       hhcorrbg_scaled[nbins]->SetLineColor(kRed);
0853       hhcorrbg_scaled[nbins]->Draw("histsame");
0854 
0855 //----------------------------------------------------------------
0856 //  Calculate RAA
0857 //----------------------------------------------------------------
0858 
0859 double u1start = 9.25;
0860 double u1stop  = 9.65;
0861 double u2start = 9.80;
0862 double u2stop  = 10.20;
0863 double u3start = 10.20;
0864 double u3stop  = 10.55;
0865 
0866   double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
0867   for(int i=0; i<nbins; i++) { 
0868     raa1[i] = grRAA1S->Eval(0.5+i*1.0);
0869     raa2[i] = grRAA2S->Eval(0.5+i*1.0);
0870     raa3[i] = grRAA3S->Eval(0.5+i*1.0);
0871   }
0872   int fbin1 = hhall_scaled[nbins]->FindBin(u1start + 0.001);
0873   int lbin1 = hhall_scaled[nbins]->FindBin(u1stop  - 0.001);
0874   int fbin2 = hhall_scaled[nbins]->FindBin(u2start + 0.001);
0875   int lbin2 = hhall_scaled[nbins]->FindBin(u2stop  - 0.001);
0876   int fbin3 = hhall_scaled[nbins]->FindBin(u3start + 0.001);
0877   int lbin3 = hhall_scaled[nbins]->FindBin(u3stop  - 0.001);
0878   cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
0879   cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
0880   cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
0881   cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
0882   cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
0883   cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
0884 
0885   double sum1[99]   = {0.};
0886   double truesum1[99]   = {0.};
0887   double ersum1[99] = {0.};
0888   double sumpp1[99]   = {0.};
0889   double ersumpp1[99] = {0.};
0890   double sum2[99]   = {0.};
0891   double truesum2[99]   = {0.};
0892   double ersum2[99] = {0.};
0893   double sumpp2[99]   = {0.};
0894   double ersumpp2[99] = {0.};
0895   double sum3[99]   = {0.};
0896   double truesum3[99]   = {0.};
0897   double ersum3[99] = {0.};
0898   double sumpp3[99]   = {0.};
0899   double ersumpp3[99] = {0.};
0900 
0901   double sumsum1[99]   = {0.};
0902   double sumsum2[99]   = {0.};
0903   double sumsum3[99]   = {0.};
0904   double sumsum1pp[99]   = {0.};
0905   double sumsum2pp[99]   = {0.};
0906   double sumsum3pp[99]   = {0.};
0907 
0908   for(int i=0; i<nbins+1; i++) {
0909 
0910     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
0911 
0912     for(int j=fbin1; j<=lbin1; j++) {
0913       sum1[i]   += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
0914       truesum1[i] += hhups1[i]->GetBinContent(j);
0915       ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
0916       sumpp1[i]   += hhups1pp[i]->GetBinContent(j);
0917       ersumpp1[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
0918     }
0919       sumsum1[i]     = truesum1[i];                   // direct count in mass range
0920       sumsum1pp[i]   = sumpp1[i];       
0921       //sumsum1[i]     = hhups1[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
0922       //sumsum1pp[i]   = hhups1pp[i]->GetEntries();
0923       //sumsum1[i]     = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
0924       //sumsum1pp[i]   = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0925 
0926         if(sumsum1[i]>0. && sumsum1pp[i]>0.) {
0927           erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1[i]/sumsum1pp[i]/sumsum1pp[i]);
0928         } else {raa1[i]=-1.0; erraa1[i] = 999.; }
0929 
0930     for(int j=fbin2; j<=lbin2; j++) {
0931       sum2[i]   += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
0932       truesum2[i] += hhups2[i]->GetBinContent(j);
0933       ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
0934       sumpp2[i]   += hhups2pp[i]->GetBinContent(j);
0935       ersumpp2[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
0936     }
0937       sumsum2[i]     = truesum2[i];                   // direct count in mass range
0938       sumsum2pp[i]   = sumpp2[i];       
0939       //sumsum2[i]     = hhups2[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
0940       //sumsum2pp[i]   = hhups2pp[i]->GetEntries();
0941       //sumsum2[i]     = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
0942       //sumsum2pp[i]   = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0943       
0944         if(sumsum2[i]>0. && sumsum2pp[i]>0.) {
0945           erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2[i]/sumsum2pp[i]/sumsum2pp[i]);
0946         } else {raa2[i]=-1.0; erraa2[i] = 999.; }
0947 
0948     for(int j=fbin3; j<=lbin3; j++) {
0949       sum3[i]   += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
0950       truesum3[i] += hhups3[i]->GetBinContent(j);
0951       ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
0952       sumpp3[i]   += hhups3pp[i]->GetBinContent(j);
0953       ersumpp3[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
0954     }
0955       sumsum3[i]     = truesum3[i];                   // direct count in mass range
0956       sumsum3pp[i]   = sumpp3[i];       
0957       //sumsum3[i]     = hhups3[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
0958       //sumsum3pp[i]   = hhups3pp[i]->GetEntries();
0959       //sumsum3[i]     = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;   // total number of upsilons in pT bin
0960       //sumsum3pp[i]   = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0961       
0962         if(truesum3[i]>0. && sumpp3[i]>0.) {
0963           erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
0964         } else {raa3[i]=-1.0; erraa3[i] = 999.; }
0965   
0966   }
0967 
0968 double raa2_rebin[9],raapt_rebin2[9],erraa2_rebin[9];
0969 double raa3_rebin[9],raapt_rebin3[9],erraa3_rebin[9];
0970 double sum2_rebin[9],ersum2_rebin[9],sum2pp_rebin[9],ersumpp2_rebin[9];
0971 double sum3_rebin[9],ersum3_rebin[9],sum3pp_rebin[9],ersumpp3_rebin[9];
0972 
0973 // Rebin Y(3S) by 2
0974 if(npts3_rebin==4) {
0975 raapt_rebin3[0] = 1.;
0976 raapt_rebin3[1] = 3.;
0977 raapt_rebin3[2] = 5.;
0978 raapt_rebin3[3] = 7.;
0979 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin2[0]);
0980 raa3_rebin[1] = grRAA3S->Eval(raapt_rebin2[1]);
0981 raa3_rebin[2] = grRAA3S->Eval(raapt_rebin2[2]);
0982 raa3_rebin[3] = grRAA3S->Eval(raapt_rebin2[3]);
0983 sum3_rebin[0] = truesum3[0]+truesum3[1];
0984 sum3_rebin[1] = truesum3[2]+truesum3[3];
0985 sum3_rebin[2] = truesum3[4]+truesum3[5];
0986 sum3_rebin[3] = truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
0987 ersum3_rebin[0] = ersum3[0]+ersum3[1];
0988 ersum3_rebin[1] = ersum3[2]+ersum3[3];
0989 ersum3_rebin[2] = ersum3[4]+ersum3[5];
0990 ersum3_rebin[3] = ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
0991 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1];
0992 sum3pp_rebin[1] = sumpp3[2]+sumpp3[3];
0993 sum3pp_rebin[2] = sumpp3[4]+sumpp3[5];
0994 sum3pp_rebin[3] = sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
0995 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1];
0996 ersumpp3_rebin[1] = ersumpp3[2]+ersumpp3[3];
0997 ersumpp3_rebin[2] = ersumpp3[4]+ersumpp3[5];
0998 ersumpp3_rebin[3] = ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
0999   erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1000   erraa3_rebin[1] = raa3[1]*sqrt(ersum3_rebin[1]/sum3_rebin[1]/sum3_rebin[1] + ersumpp3_rebin[1]/sum3pp_rebin[1]/sum3pp_rebin[1]);
1001   erraa3_rebin[2] = raa3[2]*sqrt(ersum3_rebin[2]/sum3_rebin[2]/sum3_rebin[2] + ersumpp3_rebin[2]/sum3pp_rebin[2]/sum3pp_rebin[2]);
1002   erraa3_rebin[3] = raa3[3]*sqrt(ersum3_rebin[3]/sum3_rebin[3]/sum3_rebin[3] + ersumpp3_rebin[3]/sum3pp_rebin[3]/sum3pp_rebin[3]);
1003 }
1004 else if(npts3_rebin==2) {
1005 // Rebin Y(3S) by 4
1006 raapt_rebin3[0] = 2.;
1007 raapt_rebin3[1] = 6.;
1008 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin3[0]);
1009 raa3_rebin[1] = grRAA3S->Eval(raapt_rebin3[1]);
1010 sum3_rebin[0] = truesum3[0]+truesum3[1]+truesum3[2]+truesum3[3];
1011 sum3_rebin[1] = truesum3[4]+truesum3[5]+truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
1012 ersum3_rebin[0] = ersum3[0]+ersum3[1]+ersum3[2]+ersum3[3];
1013 ersum3_rebin[1] = ersum3[4]+ersum3[5]+ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
1014 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1]+sumpp3[3]+sumpp3[3];
1015 sum3pp_rebin[1] = sumpp3[4]+sumpp3[5]+sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
1016 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1]+ersumpp3[2]+ersumpp3[3];
1017 ersumpp3_rebin[1] = ersumpp3[4]+ersumpp3[5]+ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
1018   erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1019   erraa3_rebin[1] = raa3[1]*sqrt(ersum3_rebin[1]/sum3_rebin[1]/sum3_rebin[1] + ersumpp3_rebin[1]/sum3pp_rebin[1]/sum3pp_rebin[1]);
1020 } 
1021 else if(npts3_rebin==1) {
1022 // Rebin Y(3S) by 8
1023 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin3[0]);
1024 sum3_rebin[0] = truesum3[0]+truesum3[1]+truesum3[2]+truesum3[3]+
1025                 truesum3[4]+truesum3[5]+truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
1026 ersum3_rebin[0] = ersum3[0]+ersum3[1]+ersum3[2]+ersum3[3]+
1027                   ersum3[4]+ersum3[5]+ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
1028 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1]+sumpp3[3]+sumpp3[3]+
1029                   sumpp3[4]+sumpp3[5]+sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
1030 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1]+ersumpp3[2]+ersumpp3[3]+
1031                     ersumpp3[4]+ersumpp3[5]+ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
1032 erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1033 } else {cerr << "ERROR: Wring Y(3S) binning." << endl; return;}
1034 
1035 cout << "====== Y(1S):" << endl;
1036   for(int i=0; i<nbins+1; i++) {
1037     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1038     cout << "   " << i << " " << truesum1[i] << "(" << Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1039          << sqrt(ersum1[i]) << " \t\t pp: " << sumpp1[i] << " +- " << sqrt(ersumpp1[i]) << endl;
1040   }
1041 cout << "====== Y(2S):" << endl;
1042   for(int i=0; i<nbins+1; i++) {
1043     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1044     cout << "   " << i << " " << truesum2[i] << "(" << Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1045          << sqrt(ersum2[i]) << " \t\t pp: " << sumpp2[i] << " +- " << sqrt(ersumpp2[i]) << endl;
1046   }
1047 cout << "====== Y(3S):" << endl;
1048   for(int i=0; i<nbins+1; i++) {
1049 
1050     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1051     cout << "   " << i << " " << truesum3[i] << "(" << Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1052          << sqrt(ersum3[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << endl;
1053   }
1054 
1055 //-------------------------------------------------
1056 //  Plot RAA
1057 //-------------------------------------------------
1058 
1059 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,850,600);
1060 TH2F* hh2 = new TH2F("hh2"," ",10,0.,float(npts1),10,0.,1.1);
1061 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1062 hh2->GetXaxis()->SetTitleOffset(1.0);
1063 hh2->GetXaxis()->SetTitleColor(1);
1064 hh2->GetXaxis()->SetTitleSize(0.040);
1065 hh2->GetXaxis()->SetLabelSize(0.040);
1066 hh2->GetYaxis()->SetTitle("R_{AA}");
1067 hh2->GetYaxis()->SetTitleOffset(0.7);
1068 hh2->GetYaxis()->SetTitleSize(0.050);
1069 hh2->GetYaxis()->SetLabelSize(0.040);
1070 hh2->Draw();
1071 
1072 TGraphErrors* gr1 = new TGraphErrors(npts1-1,xx1,raa1,0,erraa1);
1073 gr1->SetMarkerStyle(20);
1074 gr1->SetMarkerColor(kBlack);
1075 gr1->SetLineColor(kBlack);
1076 gr1->SetLineWidth(2);
1077 gr1->SetMarkerSize(1.5);
1078 gr1->SetName("gr1");
1079 gr1->Draw("p");
1080 
1081 TGraphErrors* gr2 = new TGraphErrors(npts2-1,xx2,raa2,0,erraa2);
1082 gr2->SetMarkerStyle(20);
1083 gr2->SetMarkerColor(kRed);
1084 gr2->SetLineColor(kRed);
1085 gr2->SetLineWidth(2);
1086 gr2->SetMarkerSize(1.5);
1087 gr2->SetName("gr2");
1088 gr2->Draw("p");
1089 
1090 //--- 3S state -------------------
1091 TGraphErrors* gr3_rebin;
1092 if(npts3_rebin==4) {
1093  gr3_rebin = new TGraphErrors(npts3_rebin-1,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1094 } else {
1095  gr3_rebin = new TGraphErrors(npts3_rebin,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1096 }
1097  gr3_rebin->SetMarkerStyle(20);
1098  gr3_rebin->SetMarkerColor(kBlue);
1099  gr3_rebin->SetLineColor(kBlue);
1100  gr3_rebin->SetLineWidth(2);
1101  gr3_rebin->SetMarkerSize(1.5);
1102  gr3_rebin->SetName("gr3");
1103  gr3_rebin->Draw("p");
1104 
1105 //TArrow* aa[9];
1106 //TLine* ll[9];
1107 //erraa3[3] = erraa3[3]*1.5;
1108 //erraa3[4] = erraa3[4]*0.6;
1109 //erraa3[4] = erraa3[4]*1.5;
1110 //erraa3[6] = erraa3[6]*1.2;
1111 /*
1112 for(int i=0; i<npts3; i++) {
1113   aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1114   aa[i]->SetLineColor(kBlue);
1115   aa[i]->SetLineWidth(2);
1116   aa[i]->Draw();
1117   ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1118   ll[i]->SetLineColor(kBlue);
1119   ll[i]->SetLineWidth(2);
1120   ll[i]->Draw();
1121 }
1122 */
1123 //--- end 3S state --------------
1124 
1125 
1126 TLegend *leg = new TLegend(0.72,0.70,0.88,0.88);
1127 leg->SetBorderSize(0);
1128 leg->SetFillColor(10);
1129 leg->SetFillStyle(1001);
1130 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1131 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1132 //TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1133 TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","p");
1134 leg->Draw();
1135 
1136 TLatex* l1 = new TLatex(0.5,1.02,"#font[72]{sPHENIX} Projection"); l1->SetTextSize(0.045); l1->SetTextFont(42); l1->Draw();
1137 TLatex* l11 = new TLatex(0.5,0.93,"0-60% cent. Au+Au, Years 1-3"); l11->SetTextSize(0.045); l11->SetTextFont(42); l11->Draw();
1138 TLatex* l2 = new TLatex(0.5,0.84,"21 nb^{-1} rec. Au+Au, 62 pb^{-1} samp. #it{p+p}");    l2->SetTextSize(0.045); l2->SetTextFont(42); l2->Draw();
1139 TLatex* l3 = new TLatex(1.8,0.48,"STAR Prelim. 0-60% Y(2S+3S)");  l3->SetTextColor(kGreen+2);   l3->SetTextSize(0.045); l3->SetTextFont(42); l3->Draw();
1140 TMarker* m1 = new TMarker(1.55,0.50,29); m1->SetMarkerSize(2); m1->SetMarkerColor(kGreen+2); m1->Draw();
1141 
1142 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1143 lll->SetLineColor(kBlue);
1144 lll->SetLineWidth(2);
1145 //lll->Draw();
1146 
1147 //==================================================================================
1148 //  STAR
1149 
1150 double star_x[3],star_erx[3],star[3],star_stater[3],star_syser[3];
1151 star_x[0] = 1.;
1152 star_x[1] = 3.;
1153 star_x[2] = 7.;
1154 star_erx[0] = 0.16;
1155 star_erx[1] = 0.16;
1156 star_erx[2] = 0.16;
1157 star[0] = 0.384;
1158 star[1] = 0.249;
1159 star[2] = 0.449;
1160 star_stater[0] = 0.147;
1161 star_stater[1] = 0.113;
1162 star_stater[2] = 0.164;
1163 star_syser[0] = 0.130;
1164 star_syser[1] = 0.086;
1165 star_syser[2] = 0.128;
1166 
1167 TGraphErrors* gr_star = new TGraphErrors(3,star_x,star,0,star_stater);
1168 gr_star->SetMarkerStyle(kFullStar);
1169 gr_star->SetMarkerColor(kGreen+2);
1170 gr_star->SetLineColor(kGreen+2);
1171 gr_star->SetLineWidth(2);
1172 gr_star->SetMarkerSize(2);
1173 gr_star->SetName("gr_star");
1174 
1175 TGraphErrors* grsys_star = new TGraphErrors(3,star_x,star,star_erx,star_syser);
1176 grsys_star->SetMarkerStyle(kFullStar);
1177 grsys_star->SetMarkerColor(kGreen+2);
1178 grsys_star->SetLineColor(kGreen+2);
1179 grsys_star->SetLineWidth(2);
1180 grsys_star->SetMarkerSize(2);
1181 grsys_star->SetName("gr_star");
1182 grsys_star->Draw("5");
1183 gr_star->Draw("p");
1184  gr3_rebin->Draw("p");
1185 
1186 /*
1187 fout = new TFile("test.root","recreate");
1188 for(int i=0; i<nbins+1; i++) {
1189   sprintf(tmpname,"hhfit_%d",i);
1190   hhfit[i]->Write();
1191   hhups1pp[i]->Write();
1192   hhups2pp[i]->Write();
1193   hhups3pp[i]->Write();
1194   hhupspp[i]->Write();
1195   hhups1[i]->Write();
1196   hhups2[i]->Write();
1197   hhups3[i]->Write();
1198   hhups[i]->Write();
1199   hhall_pp[i]->Write();
1200   hhcorrbg_pp[i]->Write();
1201   hhcorrbg_scaled[i]->Write();
1202 }
1203 fout->Write();
1204 fout->Close();
1205 */
1206 
1207 }
1208