Back to home page

sPhenix code displayed by LXR
 
 

    

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

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() {
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 = 2.4;  // 5-year plan all correladed bg plots are for 10B
0116 //double statscale = 1.43;  // 3-year plan all correladed bg plots are for 10B
0117 //double statscalepp = 1.78; // old wrong???
0118 double statscalepp = 1.13; // 5-year plan
0119 double statscale_lowlim = 7.0;
0120 double statscale_uplim = 14.0;
0121 
0122   TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
0123   TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
0124   TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
0125   TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
0126   TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
0127   TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
0128   TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
0129   TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
0130   TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
0131   TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
0132   TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
0133   TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
0134 
0135 //---------------------------------------------------------
0136 // Upsilons
0137 //---------------------------------------------------------
0138 
0139 double raapt[9],raa1s[9],raa2s[9],raa3s[9];
0140 raapt[0] = 1.5;
0141 raapt[1] = 4.5;
0142 raapt[2] = 7.5;
0143 raapt[3] = 10.5;
0144 raapt[4] = 13.5;
0145 raa1s[0] = 0.535;
0146 raa1s[1] = 0.535;
0147 raa1s[2] = 0.535;
0148 raa1s[3] = 0.535;
0149 raa1s[4] = 0.535;
0150 raa2s[0] = 0.170;
0151 raa2s[1] = 0.170;
0152 raa2s[2] = 0.170;
0153 raa2s[3] = 0.170;
0154 raa2s[4] = 0.170;
0155 /*
0156 raa1s[0] = 0.4960;
0157 raa1s[1] = 0.4960;
0158 raa1s[2] = 0.4955;
0159 raa1s[3] = 0.4968;
0160 raa1s[4] = 0.4743;
0161 raa2s[0] = 0.1710;
0162 raa2s[1] = 0.1629;
0163 raa2s[2] = 0.1326;
0164 raa2s[3] = 0.1232;
0165 raa2s[4] = 0.0928;
0166 */
0167 raa3s[0] = 0.035;
0168 raa3s[1] = 0.035;
0169 raa3s[2] = 0.035;
0170 raa3s[3] = 0.035;
0171 raa3s[4] = 0.035;
0172 TGraph* grRAA1S = new TGraph(5,raapt,raa1s);
0173 TGraph* grRAA2S = new TGraph(5,raapt,raa2s);
0174 TGraph* grRAA3S = new TGraph(5,raapt,raa3s);
0175 
0176 double NNN = statscale * 7780./0.49;      // from sPHENIX proposal for 100B minbias Au+Au events
0177 double NNNpp = statscalepp * 12130./0.90; // Upsilons in p+p from sPHENIX proposal for 7350B p+p events
0178 double ppauaurat = NNNpp/NNN;
0179 
0180 double frac[3];  // upsilons fractions
0181   frac[0] = 0.7117;
0182   frac[1] = 0.1851;
0183   frac[2] = 0.1032;
0184 double scale[3]; // mass scaling
0185   scale[0] = 1.0;
0186   scale[1] = 1.0595;
0187   scale[2] = 1.0946;
0188 //double supcor[3]; // suppression
0189 //  supcor[0]=0.535; 
0190 //  supcor[1]=0.17; 
0191 //  supcor[2]=0.035;
0192 //int Nups1 = int(NNN*eideff*eideff*frac[0]*supcor[0]);
0193 //int Nups2 = int(NNN*eideff*eideff*frac[1]*supcor[1]);
0194 //int Nups3 = int(NNN*eideff*eideff*frac[2]*supcor[2]);
0195 int Nups1nosup = int(NNN*eideff*eideff*frac[0]);
0196 int Nups2nosup = int(NNN*eideff*eideff*frac[1]);
0197 int Nups3nosup = int(NNN*eideff*eideff*frac[2]);
0198 int Nups1pp = int(NNNpp*0.90*frac[0]); // always 95% eideff in p+p
0199 int Nups2pp = int(NNNpp*0.90*frac[1]);
0200 int Nups3pp = int(NNNpp*0.90*frac[2]);
0201 
0202 int nchan=400;
0203 double start=0.0;
0204 double stop=20.0;
0205 
0206 string str_UpsilonPt = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))";
0207 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
0208 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
0209 double upsnorm = fUpsilonPt->Integral(0.,20.);
0210 
0211 // count all upsilons with suppression
0212 double tmpsum1=0.;
0213 double tmpsum2=0.;
0214 double tmpsum3=0.;
0215 for(int j=0; j<nbins; j++) {
0216   double s1 = j*1.0;
0217   double s2 = s1 + 1.0;
0218   tmpsum1 += Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
0219   tmpsum2 += Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
0220   tmpsum3 += Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
0221 }
0222 int Nups1 = int(tmpsum1);
0223 int Nups2 = int(tmpsum2);
0224 int Nups3 = int(tmpsum3);
0225 
0226   cout << "Total number of ALL Upsilons in AuAu = " << NNN << endl;
0227   cout << "Total number of ALL Upsilons in pp = " << NNNpp << endl;
0228   cout << "Number of upsilons in pp in acceptance = " << NNNpp*frac[0] << " " << NNNpp*frac[1] << " " << NNNpp*frac[2] << endl;
0229   cout << "Number of upsilons in AuAu in acceptance = " << NNN*frac[0] << " " << NNN*frac[1] << " " << NNN*frac[2] << endl;
0230   cout << "Number of upsilons in AuAu after eID efficiency = " << Nups1nosup << " " << Nups2nosup << " " << Nups3nosup << endl;
0231   cout << "Number of upsilons in AuAu plot = " << Nups1 << " " << Nups2 << " " << Nups3 << endl;
0232   cout << "Number of upsilons in pp plot = " << Nups1pp << " " << Nups2pp << " " << Nups3pp << endl;
0233 
0234 //====================================================================================================
0235 /*
0236 f=new TFile("UpsilonsInHijing.root");
0237   cout << endl << endl << "Number of entries in Upsilon NTuple = " << ntpups->GetEntries() << endl;
0238 
0239   TH1D* hforfit = new TH1D("hforfit","",nchan,start,stop);
0240   TH1D* hUps1 = new TH1D("hUps1","",nchan,start,stop);
0241   TH1D* hUps2 = new TH1D("hUps2","",nchan,start,stop);
0242   TH1D* hUps3 = new TH1D("hUps3","",nchan,start,stop);
0243   TH1D* hUps1pp = new TH1D("hUps1pp","",nchan,start,stop);
0244   TH1D* hUps2pp = new TH1D("hUps2pp","",nchan,start,stop);
0245   TH1D* hUps3pp = new TH1D("hUps3pp","",nchan,start,stop);
0246     hforfit->Sumw2();
0247     hUps1->Sumw2();
0248     hUps2->Sumw2();
0249     hUps3->Sumw2();
0250     hUps1pp->Sumw2();
0251     hUps2pp->Sumw2();
0252     hUps3pp->Sumw2();
0253 
0254 //------ Additional smearing to reproduce mass resolution -----
0255 //double smear = 0.098; // 127
0256 //double smear = 0.090;   // 119
0257 //double smear = 0.067;   // 100 MeV benchmark
0258 double smear = 0.046;  // 86 MeV resolution, Aalysis Note result 
0259 float mass;
0260 TChain* tree = new TChain("ntpups");
0261 tree->Add("UpsilonsInHijing.root");
0262 tree->SetBranchAddress("mass",&mass);
0263 int countups=0;
0264   for(int j=0; j<tree->GetEntries(); j++){
0265     tree->GetEvent(j);
0266     double newmass = myrandom->Gaus(mass,smear);
0267     hforfit->Fill(newmass);
0268     if(countups<=Nups1) hUps1->Fill(newmass);
0269     if(countups<=Nups2) hUps2->Fill(newmass*1.0595);
0270     if(countups<=Nups3) hUps3->Fill(newmass*1.0946);
0271     if(countups<=Nups1pp) hUps1pp->Fill(newmass);
0272     if(countups<=Nups2pp) hUps2pp->Fill(newmass*1.0595);
0273     if(countups<=Nups3pp) hUps3pp->Fill(newmass*1.0946);
0274     countups++;
0275   }
0276 delete tree;
0277 //---------------------------------------------------------------
0278 
0279   hforfit->SetDirectory(gROOT);
0280   hUps1->SetDirectory(gROOT);
0281   hUps2->SetDirectory(gROOT);
0282   hUps3->SetDirectory(gROOT);
0283   hUps1pp->SetDirectory(gROOT);
0284   hUps2pp->SetDirectory(gROOT);
0285   hUps3pp->SetDirectory(gROOT);
0286 
0287 f->Close();
0288 
0289 TH1F* hUps = (TH1F*)hUps1->Clone("hUps");
0290 hUps->Add(hUps2);
0291 hUps->Add(hUps3);
0292 hUps->SetLineWidth(2);
0293 TH1F* hUpspp = (TH1F*)hUps1pp->Clone("hUpspp");
0294 hUpspp->Add(hUps2pp);
0295 hUpspp->Add(hUps3pp);
0296 hUpspp->SetLineWidth(2);
0297 
0298 TCanvas* cups = new TCanvas("cups","Upsilons (all pT)",100,100,1000,500);
0299 cups->Divide(2,1);
0300 
0301   cups_1->cd();
0302   cups_1->SetLogy();
0303   fCB->SetParameter(0,5000.);
0304   fCB->SetParameter(1,1.5); 
0305   fCB->SetParameter(2,1.2); 
0306   fCB->SetParameter(3,9.448);
0307   fCB->SetParameter(4,0.070);
0308   hforfit->Fit("fCB","rl","",5.,9.7);
0309   hforfit->SetAxisRange(5.,11.);
0310   hforfit->SetLineWidth(2);
0311   hforfit->Draw();
0312   double myupsmass = fCB->GetParameter(3);
0313   double upswidth = fCB->GetParameter(4);
0314 
0315   cups_2->cd();
0316   cups_2->SetLogy();
0317   fTCB->SetParameter(0,2000.);
0318   fTCB->SetParameter(1,1.5);
0319   fTCB->SetParameter(2,1.25);
0320   fTCB->SetParameter(3,9.448);
0321   fTCB->FixParameter(4,upswidth); // fix width from the single peak fit
0322   fTCB->SetParameter(5,500.);
0323   fTCB->SetParameter(6,100.);
0324   hUps->Fit(fTCB,"rl","",6.,11.);
0325   double upswidth3 = fTCB->GetParameter(4);
0326   double upswidth3err = fTCB->GetParError(4);
0327   hUps->SetAxisRange(7.0,11.0);
0328   hUps->Draw();
0329 
0330 cout << "MASS RESOLUTION = " << upswidth << " " << upswidth3 << " +- " << upswidth3err << endl;
0331 */
0332 //====================================================================================
0333 //====================================================================================
0334 //====================================================================================
0335 
0336 //
0337 // these histograms (hhups*) are randomly generated using the fit above to single 
0338 // peak (fCB function) and used for the RAA uncertainty calculation
0339 //
0340 // FORCE specific RESOLUTION and tail shape
0341   double tonypar1 = 0.98;  // Tony's 100 pion simulation April 2019
0342   double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
0343   //double tonypar3 = 9.437;  // Tony's 100 pion simulation April 2019
0344   double tonypar3 = 9.448;  // benchmark
0345   double tonypar4 = 0.100;  // benchmark
0346   double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
0347   fCBpp->SetParameter(0,1000.); 
0348   fCBpp->SetParameter(1,tonypar1); 
0349   fCBpp->SetParameter(2,tonypar2); 
0350   fCBpp->SetParameter(3,tonypar3); 
0351   fCBpp->SetParameter(4,tonypar4pp);  
0352   fCBauau->SetParameter(0,1000.); 
0353   fCBauau->SetParameter(1,tonypar1); 
0354   fCBauau->SetParameter(2,tonypar2); 
0355   fCBauau->SetParameter(3,tonypar3); 
0356   fCBauau->SetParameter(4,tonypar4);  
0357 
0358 
0359 char hhname[999];
0360 TH1D* hhups[nbins+1];
0361 TH1D* hhups1[nbins+1];
0362 TH1D* hhups2[nbins+1];
0363 TH1D* hhups3[nbins+1];
0364 TH1D* hhupspp[nbins+1];
0365 TH1D* hhups1pp[nbins+1];
0366 TH1D* hhups2pp[nbins+1];
0367 TH1D* hhups3pp[nbins+1];
0368 for(int i=0; i<nbins+1; i++) {
0369   sprintf(hhname,"hhups_%d",i);
0370   hhups[i] = new TH1D(hhname,"",nchan,start,stop);
0371   hhups[i]->Sumw2();
0372   sprintf(hhname,"hhups1_%d",i);
0373   hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
0374   hhups1[i]->Sumw2();
0375   sprintf(hhname,"hhups2_%d",i);
0376   hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
0377   hhups2[i]->Sumw2();
0378   sprintf(hhname,"hhups3_%d",i);
0379   hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
0380   hhups3[i]->Sumw2();
0381   hhups[i]->SetLineWidth(2);
0382   hhups1[i]->SetLineWidth(2);
0383   hhups2[i]->SetLineWidth(2);
0384   hhups3[i]->SetLineWidth(2);
0385   sprintf(hhname,"hhupspp_%d",i);
0386   hhupspp[i] = new TH1D(hhname,"",nchan,start,stop);
0387   hhupspp[i]->Sumw2();
0388   sprintf(hhname,"hhups1pp_%d",i);
0389   hhups1pp[i] = new TH1D(hhname,"",nchan,start,stop);
0390   hhups1pp[i]->Sumw2();
0391   sprintf(hhname,"hhups2pp_%d",i);
0392   hhups2pp[i] = new TH1D(hhname,"",nchan,start,stop);
0393   hhups2pp[i]->Sumw2();
0394   sprintf(hhname,"hhups3pp_%d",i);
0395   hhups3pp[i] = new TH1D(hhname,"",nchan,start,stop);
0396   hhups3pp[i]->Sumw2();
0397   hhupspp[i]->SetLineWidth(2);
0398   hhups1pp[i]->SetLineWidth(2);
0399   hhups2pp[i]->SetLineWidth(2);
0400   hhups3pp[i]->SetLineWidth(2);
0401 }
0402 
0403 for(int j=0; j<nbins; j++) {
0404 
0405   double s1 = j*1.0;
0406   double s2 = s1 + 1.0;
0407   double tmpnups1 = Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
0408   double tmpnups2 = Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
0409   double tmpnups3 = Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
0410   //cout << "Nups1["<<j<<"]="<<tmpnups1<<";"<<endl;
0411   //cout << "Nups2["<<j<<"]="<<tmpnups2<<";"<<endl;
0412   //cout << "Nups3["<<j<<"]="<<tmpnups3<<";"<<endl;
0413   fCBauau->SetParameter(3,tonypar3);
0414     for(int i=0; i<int(tmpnups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0415   fCBauau->SetParameter(3,tonypar3*scale[1]);
0416     for(int i=0; i<int(tmpnups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0417   fCBauau->SetParameter(3,tonypar3*scale[2]);
0418     for(int i=0; i<int(tmpnups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0419 
0420   double tmpnups1pp = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0421   double tmpnups2pp = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0422   double tmpnups3pp = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0423   cout << "Nups1pp["<<j<<"]="<<tmpnups1pp<<";"<<endl;
0424   cout << "Nups2pp["<<j<<"]="<<tmpnups2pp<<";"<<endl;
0425   cout << "Nups3pp["<<j<<"]="<<tmpnups3pp<<";"<<endl;
0426   fCBpp->SetParameter(3,tonypar3);
0427     for(int i=0; i<int(tmpnups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0428   fCBpp->SetParameter(3,tonypar3*scale[1]);
0429     for(int i=0; i<int(tmpnups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0430   fCBpp->SetParameter(3,tonypar3*scale[2]);
0431     for(int i=0; i<int(tmpnups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0432 
0433 } // end loop over pT bins
0434 
0435 // all pT
0436 
0437   fCBpp->SetParameter(3,tonypar3);
0438   fCBauau->SetParameter(3,tonypar3);
0439     for(int i=0; i<int(Nups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0440     for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0441   fCBpp->SetParameter(3,tonypar3*scale[1]);
0442   fCBauau->SetParameter(3,tonypar3*scale[1]);
0443     for(int i=0; i<int(Nups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0444     for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0445   fCBpp->SetParameter(3,tonypar3*scale[2]);
0446   fCBauau->SetParameter(3,tonypar3*scale[2]);
0447     for(int i=0; i<int(Nups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0448     for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0449 
0450 // fit and draw pp no bg
0451 
0452 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
0453   fTCBpp->SetParameter(0,2000.);
0454   fTCBpp->FixParameter(1,tonypar1);
0455   fTCBpp->FixParameter(2,tonypar2);
0456   fTCBpp->SetParameter(3,tonypar3);
0457   fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
0458   fTCBpp->SetParameter(5,500.);
0459   fTCBpp->SetParameter(6,100.);
0460   hhupspp[nbins]->Fit(fTCBpp,"rl","",7.,11.);
0461   hhupspp[nbins]->SetAxisRange(7.,11.);
0462   hhupspp[nbins]->SetMarkerSize(1.0);
0463   hhupspp[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0464   hhupspp[nbins]->GetXaxis()->SetTitleOffset(1.0);
0465   double tmpamp1 = hhupspp[nbins]->GetFunction("fTCBpp")->GetParameter(0);
0466   double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0467   double tmpamp6 = tmpamp1*frac[2]/frac[0];
0468   hhupspp[nbins]->Draw();
0469 
0470   fCB1s->SetLineColor(kBlue);
0471   fCB1s->SetLineWidth(1);
0472     fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
0473     fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
0474     fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
0475     fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
0476     fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
0477   fCB2s->SetLineColor(kRed);
0478   fCB2s->SetLineWidth(1);
0479     fCB2s->SetParameter(0,tmpamp5);
0480     fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
0481     fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
0482     fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
0483     fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
0484   fCB3s->SetLineColor(kGreen+2);
0485   fCB3s->SetLineWidth(1);
0486     fCB3s->SetParameter(0,tmpamp6);
0487     fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
0488     fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
0489     fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
0490     fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
0491   fCB1s->Draw("same");
0492   fCB2s->Draw("same");
0493   fCB3s->Draw("same");
0494 
0495 // fit and draw AuAu no bg
0496 
0497 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Au+Au",100,100,600,600);
0498   fTCBauau->SetParameter(0,2000.);
0499   fTCBauau->FixParameter(1,tonypar1);
0500   fTCBauau->FixParameter(2,tonypar2);
0501   fTCBauau->SetParameter(3,tonypar3);
0502   fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
0503   fTCBauau->SetParameter(5,500.);
0504   fTCBauau->SetParameter(6,100.);
0505   hhups[nbins]->Fit(fTCBauau,"rl","",7.,11.);
0506   hhups[nbins]->SetAxisRange(7.,11.);
0507   hhups[nbins]->SetMarkerSize(1.0);
0508   hhups[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0509   hhups[nbins]->GetXaxis()->SetTitleOffset(1.0);
0510   tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
0511   tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0512   tmpamp6 = tmpamp1*frac[2]/frac[0];
0513   hhups[nbins]->Draw();
0514 
0515 TCanvas* cupsvspt = new TCanvas("cupsvspt","Upsilons vs. p_{T}",100,100,1200,900);
0516 cupsvspt->Divide(4,3);
0517 for(int i=0; i<12; i++) {
0518 if(i==0)  {cupsvspt->cd(1);}
0519 if(i==1)  {cupsvspt->cd(2);}
0520 if(i==2)  {cupsvspt->cd(3);}
0521 if(i==3)  {cupsvspt->cd(4);}
0522 if(i==4)  {cupsvspt->cd(5);}
0523 if(i==5)  {cupsvspt->cd(6);}
0524 if(i==6)  {cupsvspt->cd(7);}
0525 if(i==7)  {cupsvspt->cd(8);}
0526 if(i==8)  {cupsvspt->cd(9);}
0527 if(i==9)  {cupsvspt->cd(10);}
0528 if(i==10) {cupsvspt->cd(11);}
0529 if(i==11) {cupsvspt->cd(12);}
0530 hhups[i]->SetAxisRange(7.0,11.0); hhups[i]->Draw();
0531   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhups[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0532 }
0533 
0534 //----------------------------------------------------
0535 //  Backgrounds
0536 //----------------------------------------------------
0537 
0538 TH1D* hhall[nbins+1];
0539 TH1D* hhall_scaled[nbins+1];
0540 
0541 TH1D* hhtotbg[nbins+1]; 
0542 TH1D* hhtotbg_scaled[nbins+1]; 
0543 TH1D* hhcombbg[nbins+1];
0544 TH1D* hhcombbg_scaled[nbins+1];
0545 TH1D* hhfakefake[nbins+1];
0546 TH1D* hhfakehf[nbins+1];
0547 TH1D* hhbottom[nbins+1];
0548 TH1D* hhcharm[nbins+1];
0549 TH1D* hhdy[nbins+1];
0550 TH1D* hhcorrbg[nbins+1];
0551 TH1D* hhcorrbg_scaled[nbins+1];
0552 TH1D* hhfit[nbins+1];
0553 char tmpname[999];
0554 
0555 //----------------------------------------------------------------------------------------
0556 // Correlated BG calculated for 10B central AuAu events
0557 //----------------------------------------------------------------------------------------
0558 
0559 double corrbgfitpar0;
0560 double corrbgfitpar1;
0561 
0562 TFile* f=new TFile("ccbb_eideff09.root");
0563 for(int i=0; i<nbins+1; i++) {
0564   sprintf(tmpname,"hhbottom_%d",i);
0565   hhbottom[i] = (TH1D*)f->Get(tmpname);
0566   hhbottom[i]->SetDirectory(gROOT);
0567   sprintf(tmpname,"hhcharm_%d",i);
0568   hhcharm[i] = (TH1D*)f->Get(tmpname);
0569   hhcharm[i]->SetDirectory(gROOT);
0570   sprintf(tmpname,"hhdy_%d",i);
0571   hhdy[i] = (TH1D*)f->Get(tmpname);
0572   hhdy[i]->SetDirectory(gROOT);
0573   sprintf(tmpname,"hhcorrbg_%d",i);
0574   hhcorrbg[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0575   hhcorrbg[i]->Add(hhcharm[i]);
0576   hhcorrbg[i]->Add(hhdy[i]);
0577   sprintf(tmpname,"hhcorrbg_scaled_%d",i);
0578   hhcorrbg_scaled[i] = (TH1D*)hhcorrbg[i]->Clone(tmpname);
0579   hhcorrbg[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0580   hhbottom[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0581   hhdy[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0582     if(i==nbins) {
0583        corrbgfitpar0 = hhcorrbg[i]->GetFunction("expo")->GetParameter(0);
0584        corrbgfitpar1 = hhcorrbg[i]->GetFunction("expo")->GetParameter(1);
0585     }
0586     cout << "bgpar0["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(0)+TMath::Log(statscale)<<";"<< endl; 
0587     cout << "bgpar1["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(1)<<";"<< endl; 
0588     for(int k=1; k<=hhcorrbg[i]->GetNbinsX(); k++) {
0589       if(hhcorrbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[i]->GetBinLowEdge(k)+hhcorrbg[i]->GetBinWidth(k))>statscale_uplim) {
0590         hhcorrbg_scaled[i]->SetBinContent(k,0.); 
0591         hhcorrbg_scaled[i]->SetBinError(k,0.);
0592       }
0593       else {
0594         double tmp = statscale * hhcorrbg[i]->GetFunction("expo")->Eval(hhcorrbg[i]->GetBinCenter(k));
0595         double tmprnd = myrandom->Poisson(tmp); 
0596         if(tmprnd<0.) { tmprnd=0.; }
0597         hhcorrbg_scaled[i]->SetBinContent(k,tmprnd); 
0598         hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0599       }
0600     }
0601   hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0602   hhcorrbg[i]->SetDirectory(gROOT);
0603   hhcorrbg_scaled[i]->SetDirectory(gROOT);
0604 }
0605 f->Close();
0606 
0607 TCanvas* c2 = new TCanvas("c2","Correlated BG (all p_{T})",100,100,600,600);
0608   
0609 hhbottom[nbins]->SetAxisRange(7.0,14.0);
0610 hhbottom[nbins]->SetLineColor(kBlue);
0611 hhbottom[nbins]->SetLineWidth(2);
0612 
0613 hhbottom[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0614 hhbottom[nbins]->GetXaxis()->SetTitleOffset(1.0);
0615 hhbottom[nbins]->GetXaxis()->SetTitleColor(1);
0616 hhbottom[nbins]->GetXaxis()->SetTitleSize(0.040);
0617 hhbottom[nbins]->GetXaxis()->SetLabelSize(0.040);
0618 //hhbottom[nbins]->GetYaxis()->SetTitle("Correlated background");
0619 hhbottom[nbins]->GetYaxis()->SetTitleOffset(1.3);
0620 hhbottom[nbins]->GetYaxis()->SetTitleSize(0.040);
0621 hhbottom[nbins]->GetYaxis()->SetLabelSize(0.040);
0622 
0623 hhcorrbg[nbins]->SetAxisRange(7.0,14.0);
0624 hhcorrbg[nbins]->SetMinimum(0.1);
0625 hhcorrbg[nbins]->SetLineColor(kBlack);
0626 hhcorrbg[nbins]->SetLineWidth(2);
0627 hhcorrbg[nbins]->Draw("hist");
0628 
0629 hhbottom[nbins]->SetMarkerColor(kBlue);
0630 hhbottom[nbins]->SetLineColor(kBlue);
0631 hhbottom[nbins]->Draw("same");
0632 
0633 hhdy[nbins]->SetMarkerColor(kGreen+2);
0634 hhdy[nbins]->SetLineColor(kGreen+2);
0635 hhdy[nbins]->SetLineWidth(2);
0636 hhdy[nbins]->Draw("same");
0637 
0638 hhcharm[nbins]->SetMarkerColor(kRed);
0639 hhcharm[nbins]->SetLineColor(kRed);
0640 hhcharm[nbins]->SetLineWidth(2);
0641 hhcharm[nbins]->Draw("same");
0642 
0643 TCanvas* c0 = new TCanvas("c0","Correlated BG vs. p_{T} 10B events",100,100,1200,900);
0644 c0->Divide(4,3);
0645 for(int i=0; i<nbins; i++) {
0646 if(i==0)  {c0->cd(1);}
0647 if(i==1)  {c0->cd(2);}
0648 if(i==2)  {c0->cd(3);}
0649 if(i==3)  {c0->cd(4);}
0650 if(i==4)  {c0->cd(5);}
0651 if(i==5)  {c0->cd(6);}
0652 if(i==6)  {c0->cd(7);}
0653 if(i==7)  {c0->cd(8);}
0654 if(i==8)  {c0->cd(9);}
0655 if(i==9)  {c0->cd(10);}
0656 if(i==10) {c0->cd(11);}
0657 if(i==11) {c0->cd(12);}
0658   hhcorrbg[i]->SetAxisRange(7.,14.);
0659   hhcorrbg[i]->GetFunction("expo")->SetLineColor(kBlack); 
0660   hhbottom[i]->GetFunction("expo")->SetLineColor(kBlue); 
0661   hhdy[i]->GetFunction("expo")->SetLineColor(kGreen+2); 
0662   hhcorrbg[i]->SetLineColor(kBlack); 
0663   hhbottom[i]->SetLineColor(kBlue); 
0664   hhdy[i]->SetLineColor(kGreen+2); 
0665   hhcorrbg[i]->Draw();
0666   hhbottom[i]->Draw("same");
0667   hhdy[i]->Draw("same");
0668   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0669 }
0670 
0671 TCanvas* c0scaled = new TCanvas("c0scaled","SCALED Correlated BG vs. p_{T}",100,100,1200,900);
0672 c0scaled->Divide(4,3);
0673 for(int i=0; i<nbins; i++) {
0674 if(i==0)  {c0scaled->cd(1);}
0675 if(i==1)  {c0scaled->cd(2);}
0676 if(i==2)  {c0scaled->cd(3);}
0677 if(i==3)  {c0scaled->cd(4);}
0678 if(i==4)  {c0scaled->cd(5);}
0679 if(i==5)  {c0scaled->cd(6);}
0680 if(i==6)  {c0scaled->cd(7);}
0681 if(i==7)  {c0scaled->cd(8);}
0682 if(i==8)  {c0scaled->cd(9);}
0683 if(i==9)  {c0scaled->cd(10);}
0684 if(i==10) {c0scaled->cd(11);}
0685 if(i==11) {c0scaled->cd(12);}
0686 hhcorrbg_scaled[i]->SetAxisRange(7.0,14.0); hhcorrbg_scaled[i]->Draw();
0687   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0688 }
0689 
0690 //-----------------------------------------------------------------------------------------
0691 //  Correlated bg in p+p
0692 //-----------------------------------------------------------------------------------------
0693 
0694 TH1D* hhbottom_pp[nbins+1];
0695 TH1D* hhdy_pp[nbins+1];
0696 TH1D* hhcorrbg_pp[nbins+1];
0697 TH1D* hhall_pp[nbins+1];
0698 
0699 double ppcorr = 8300./10./962.;
0700 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
0701 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
0702 
0703 for(int i=0; i<nbins+1; i++) {
0704 
0705   sprintf(tmpname,"hhbottom_pp_%d",i);
0706   hhbottom_pp[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0707   for(int k=1; k<=hhbottom_pp[i]->GetNbinsX(); k++) {
0708     if(hhbottom_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp[i]->GetBinLowEdge(k)+hhbottom_pp[i]->GetBinWidth(k))>statscale_uplim) {
0709       hhbottom_pp[i]->SetBinContent(k,0.);
0710       hhbottom_pp[i]->SetBinError(k,0.);
0711     }
0712     else {
0713       double tmp = ppcorr * fbottom_nosup_corr->Eval(hhbottom[i]->GetBinCenter(k)) * hhbottom[i]->GetFunction("expo")->Eval(hhbottom[i]->GetBinCenter(k));
0714       double tmprnd = myrandom->Poisson(tmp);
0715       if(tmprnd<0.) { tmprnd=0.; }
0716       hhbottom_pp[i]->SetBinContent(k,tmprnd);
0717       hhbottom_pp[i]->SetBinError(k,sqrt(tmprnd));
0718     }
0719   }
0720 
0721   sprintf(tmpname,"hhdy_pp_%d",i);
0722   hhdy_pp[i] = (TH1D*)hhdy[i]->Clone(tmpname);
0723   for(int k=1; k<=hhdy_pp[i]->GetNbinsX(); k++) {
0724     if(hhdy_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp[i]->GetBinLowEdge(k)+hhdy_pp[i]->GetBinWidth(k))>statscale_uplim) {
0725       hhdy_pp[i]->SetBinContent(k,0.);
0726       hhdy_pp[i]->SetBinError(k,0.);
0727     }
0728     else {
0729       double tmp = ppcorr * hhdy[i]->GetFunction("expo")->Eval(hhdy[i]->GetBinCenter(k));
0730       double tmprnd = myrandom->Poisson(tmp);
0731       if(tmprnd<0.) { tmprnd=0.; }
0732       hhdy_pp[i]->SetBinContent(k,tmprnd);
0733       hhdy_pp[i]->SetBinError(k,sqrt(tmprnd));
0734     }
0735   }
0736 
0737   sprintf(tmpname,"hhcorrbg_pp_%d",i);
0738   hhcorrbg_pp[i] = (TH1D*)hhbottom_pp[i]->Clone(tmpname);
0739   hhcorrbg_pp[i]->Add(hhdy_pp[i]);
0740     hhcorrbg_pp[i]->SetMarkerColor(kBlack);
0741     hhcorrbg_pp[i]->SetLineColor(kBlack);
0742     hhbottom_pp[i]->SetLineColor(kBlue);
0743     hhdy_pp[i]->SetLineColor(kGreen+2);
0744   sprintf(tmpname,"hhall_pp_%d",i);
0745   hhall_pp[i] = (TH1D*)hhcorrbg_pp[i]->Clone(tmpname);
0746   hhall_pp[i]->Add(hhupspp[i]);
0747     hhall_pp[i]->SetLineColor(kMagenta);
0748     hhall_pp[i]->SetMarkerColor(kMagenta);
0749 
0750   hhcorrbg_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0751   hhcorrbg_pp[i]->GetFunction("expo")->SetLineColor(kBlack);
0752   hhbottom_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0753   hhbottom_pp[i]->GetFunction("expo")->SetLineColor(kBlue);
0754   hhdy_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0755   hhdy_pp[i]->GetFunction("expo")->SetLineColor(kGreen+2);
0756 
0757 } // end loop over pT bins
0758 
0759 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
0760   hhall_pp[nbins]->SetAxisRange(7.,12.);
0761   hhcorrbg_pp[nbins]->Draw("pehist");
0762   hhbottom_pp[nbins]->Draw("histsame");
0763   hhdy_pp[nbins]->Draw("histsame");
0764 
0765 
0766 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
0767   hhall_pp[nbins]->SetAxisRange(7.,12.);
0768   hhall_pp[nbins]->Draw("pehist");
0769   hhcorrbg_pp[nbins]->Draw("pesame");
0770   hhbottom_pp[nbins]->Draw("same");
0771   hhdy_pp[nbins]->Draw("same");
0772 
0773 TCanvas* cpp_vspt = new TCanvas("cpp_vspt","corr bg + sig vs pt in pp",50,50,1200,900);
0774 cpp_vspt->Divide(4,3);
0775 for(int i=0; i<nbins; i++) {
0776 if(i==0)  {cpp_vspt->cd(1);}
0777 if(i==1)  {cpp_vspt->cd(2);}
0778 if(i==2)  {cpp_vspt->cd(3);}
0779 if(i==3)  {cpp_vspt->cd(4);}
0780 if(i==4)  {cpp_vspt->cd(5);}
0781 if(i==5)  {cpp_vspt->cd(6);}
0782 if(i==6)  {cpp_vspt->cd(7);}
0783 if(i==7)  {cpp_vspt->cd(8);}
0784 if(i==8)  {cpp_vspt->cd(9);}
0785 if(i==9)  {cpp_vspt->cd(10);}
0786 if(i==10) {cpp_vspt->cd(11);}
0787 if(i==11) {cpp_vspt->cd(12);}
0788 hhall_pp[i]->SetAxisRange(7.0,12.0); hhall_pp[i]->Draw("hist"); hhcorrbg_pp[i]->Draw("histsame");
0789   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0790 }
0791 
0792 //-----------------------------------------------------------------------------------------
0793 // Combinatorial BG calculated for 10B central AuAu events
0794 //-----------------------------------------------------------------------------------------
0795 TCanvas* dummy = new TCanvas("dummy","dummy",0,0,500,500);
0796 
0797 f = new TFile("fakee_eideff09.root");
0798 for(int i=0; i<nbins+1; i++) {
0799   sprintf(tmpname,"hhfakefake_%d",i);
0800   hhfakefake[i] = (TH1D*)f->Get(tmpname);
0801   hhfakefake[i]->SetDirectory(gROOT);
0802 }
0803 f->Close();
0804 
0805 f = new TFile("crossterms_eideff09.root");
0806 for(int i=0; i<nbins+1; i++) {
0807   sprintf(tmpname,"hhfakehf_%d",i);
0808   hhfakehf[i] = (TH1D*)f->Get(tmpname);
0809   hhfakehf[i]->SetDirectory(gROOT);
0810 }
0811 f->Close();
0812 
0813 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
0814 fbg->SetParameters(10., -1.0, 4., -0.1);
0815 fbg->SetParLimits(1.,-999.,0.);
0816 fbg->SetParLimits(3.,-999.,0.);
0817 
0818 for(int i=0; i<nbins+1; i++) {
0819   sprintf(tmpname,"hhcombbg_%d",i);
0820   hhcombbg[i] = (TH1D*)hhfakefake[i]->Clone(tmpname);
0821   hhcombbg[i]->Add(hhfakehf[i]);
0822   sprintf(tmpname,"hhcombbg_scaled_%d",i);
0823   hhcombbg_scaled[i] = (TH1D*)hhcombbg[i]->Clone(tmpname);
0824   if(i==nbins) { fbg->SetParameters(10., -1.0, 4., -0.1); }
0825   hhcombbg[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0826 
0827     for(int k=1; k<=hhcombbg[i]->GetNbinsX(); k++) {
0828       if(hhcombbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[i]->GetBinLowEdge(k)+hhcombbg[i]->GetBinWidth(k))>statscale_uplim) {
0829         hhcombbg_scaled[i]->SetBinContent(k,0.);
0830         hhcombbg_scaled[i]->SetBinError(k,0.);
0831       }
0832       else {
0833         double tmp = statscale * hhcombbg[i]->GetFunction("fbg")->Eval(hhcombbg[i]->GetBinCenter(k));
0834         double tmprnd = myrandom->Poisson(tmp);
0835         if(tmprnd<0.) { tmprnd=0.; }
0836         hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
0837         hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0838       }
0839     }
0840   hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0841 }
0842 
0843 delete dummy;
0844 
0845 TCanvas* C1 = new TCanvas("C1","Combinatorial BG (ALL p_{T})",100,100,600,600);
0846 C1->SetLogy();
0847   hhfakefake[nbins]->SetAxisRange(7.0,14.0);
0848   hhfakefake[nbins]->SetMinimum(0.1);
0849   hhfakefake[nbins]->SetMaximum(5000.);
0850   hhfakefake[nbins]->SetLineColor(kGreen+2);
0851   hhfakefake[nbins]->SetLineWidth(2);
0852   hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
0853   hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
0854   hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
0855   hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
0856   hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
0857   hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
0858   hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
0859   hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
0860   hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
0861   hhfakefake[nbins]->Draw("e");
0862 
0863   hhfakehf[nbins]->SetLineColor(kOrange+4);
0864   hhfakehf[nbins]->SetLineWidth(2);
0865   hhfakehf[nbins]->Draw("esame");
0866 
0867   hhcombbg[nbins]->SetLineColor(kBlack);
0868   hhcombbg[nbins]->SetLineWidth(2);
0869   hhcombbg[nbins]->Draw("esame");
0870 
0871 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG (ALL p_{T})",100,100,600,600);
0872 C1sc->SetLogy(); 
0873   hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
0874   hhcombbg_scaled[nbins]->Draw("esame");
0875 
0876 TCanvas* c00 = new TCanvas("c00","Combinatorial BG vs. p_{T}",150,150,1200,900);
0877 c00->Divide(4,3);
0878 
0879 for(int i=0; i<nbins; i++) {
0880 if(i==0)  {c00->cd(1);}
0881 if(i==1)  {c00->cd(2);}
0882 if(i==2)  {c00->cd(3);}
0883 if(i==3)  {c00->cd(4);}
0884 if(i==4)  {c00->cd(5);}
0885 if(i==5)  {c00->cd(6);}
0886 if(i==6)  {c00->cd(7);}
0887 if(i==7)  {c00->cd(8);}
0888 if(i==8)  {c00->cd(9);}
0889 if(i==9)  {c00->cd(10);}
0890 if(i==10) {c00->cd(11);}
0891 if(i>=11) {c00->cd(12);}
0892 hhcombbg[i]->SetAxisRange(7.0,14.0); hhcombbg[i]->Draw();
0893   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0894 }
0895 
0896 TCanvas* c00scaled = new TCanvas("c00scaled","SCALED Combinatorial BG vs. p_{T}",150,150,1200,900);
0897 c00scaled->Divide(4,3);
0898 
0899 for(int i=0; i<nbins; i++) {
0900 if(i==0)  {c00scaled->cd(1);}
0901 if(i==1)  {c00scaled->cd(2);}
0902 if(i==2)  {c00scaled->cd(3);}
0903 if(i==3)  {c00scaled->cd(4);}
0904 if(i==4)  {c00scaled->cd(5);}
0905 if(i==5)  {c00scaled->cd(6);}
0906 if(i==6)  {c00scaled->cd(7);}
0907 if(i==7)  {c00scaled->cd(8);}
0908 if(i==8)  {c00scaled->cd(9);}
0909 if(i==9)  {c00scaled->cd(10);}
0910 if(i==10) {c00scaled->cd(11);}
0911 if(i>=11) {c00scaled->cd(12);}
0912 hhcombbg_scaled[i]->SetAxisRange(statscale_lowlim,statscale_uplim); hhcombbg_scaled[i]->Draw();
0913   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0914 }
0915 
0916 //---------------------------------------------
0917 //  Sibnal + BG
0918 //---------------------------------------------
0919 
0920 for(int i=0; i<nbins+1; i++) {
0921   sprintf(tmpname,"hhtotbg_scaled_%d",i);
0922   hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
0923   hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
0924 }
0925 
0926 for(int i=0; i<nbins+1; i++) {
0927 sprintf(tmpname,"hhall_scaled_%d",i);
0928 hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
0929 hhall_scaled[i]->Add(hhups[i]);
0930 }
0931 
0932 TCanvas* c000 = new TCanvas("c000","Signal + Background vs. p_{T}",200,200,1200,900);
0933 c000->Divide(4,3);
0934 for(int i=0; i<nbins; i++) {
0935 if(i==0)  {c000->cd(1);}
0936 if(i==1)  {c000->cd(2);}
0937 if(i==2)  {c000->cd(3);}
0938 if(i==3)  {c000->cd(4);}
0939 if(i==4)  {c000->cd(5);}
0940 if(i==5)  {c000->cd(6);}
0941 if(i==6)  {c000->cd(7);}
0942 if(i==7)  {c000->cd(8);}
0943 if(i==8)  {c000->cd(9);}
0944 if(i==9)  {c000->cd(10);}
0945 if(i==10) {c000->cd(11);}
0946 if(i==11) {c000->cd(12);}
0947 hhall_scaled[i]->SetAxisRange(7.0,14.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
0948   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
0949 }
0950 
0951 // Fit Signal + Correlated BG
0952 
0953 // hhfit[] is signal + correlated bg
0954 for(int i=0; i<nbins+1; i++) {
0955   sprintf(tmpname,"hhfit_%d",i);
0956   hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
0957   for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
0958     hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
0959     hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
0960   }
0961 }
0962 
0963 //---------------------------------------------------------------------
0964 // plot signal + correlated bg for all pT
0965 //---------------------------------------------------------------------
0966 
0967   double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
0968   double tmppar1 = corrbgfitpar1;
0969 cout << "###### " << tmppar0 << " " << tmppar1 << endl;
0970 
0971 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
0972 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
0973 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
0974 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
0975 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
0976 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
0977 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
0978 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
0979 fSandBpp->SetParameter(7,tmppar0);
0980 fSandBpp->SetParameter(8,tmppar1);
0981 fSandBpp->SetLineColor(kBlue);
0982 fSandBpp->SetLineStyle(2);
0983 
0984 cout << "######### " << fTCBauau->GetParameter(0) << " " << fTCBauau->GetParameter(5) << " " << fTCBauau->GetParameter(6) << endl;
0985 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
0986 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
0987 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
0988 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
0989 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
0990 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
0991 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
0992 fSandBauau->SetParameter(7,tmppar0);
0993 fSandBauau->SetParameter(8,tmppar1);
0994 fSandBauau->SetLineColor(kRed);
0995 
0996 
0997 TCanvas* cfitall = new TCanvas("cfitall","FIT all pT",270,270,600,600);
0998   hhfit[nbins]->SetAxisRange(7.0,14.);
0999   hhfit[nbins]->GetXaxis()->CenterTitle();
1000   hhfit[nbins]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
1001   hhfit[nbins]->GetXaxis()->SetTitleOffset(1.1);
1002   hhfit[nbins]->GetXaxis()->SetLabelSize(0.045);
1003   hhfit[nbins]->GetYaxis()->CenterTitle();
1004   hhfit[nbins]->GetYaxis()->SetLabelSize(0.045);
1005   hhfit[nbins]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
1006   hhfit[nbins]->GetYaxis()->SetTitleOffset(1.5);
1007   hhfit[nbins]->Draw("pehist");
1008   fSandBpp->Draw("same");
1009   fSandBauau->Draw("same");
1010     //hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1011     //hhcorrbg_scaled[nbins]->Scale(0.82);
1012     //hhcorrbg_scaled[nbins]->Scale(0.70);
1013     //hhcorrbg_scaled[nbins]->Draw("histsame");
1014 
1015   TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
1016   fmycorrbg->SetParameters(tmppar0,tmppar1);
1017   fmycorrbg->SetLineStyle(2);
1018   fmycorrbg->SetLineColor(kRed);
1019   fmycorrbg->Draw("same");
1020 
1021 
1022 double myheight = fTCBauau->GetParameter(0);
1023 TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
1024 ld1->SetTextSize(0.035);
1025 ld1->Draw();
1026 TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
1027 ld2->SetTextSize(0.035);
1028 ld2->Draw();
1029 TLatex* ld3;
1030 if(statscale==2.4) {
1031 ld3 = new TLatex(10.1,myheight-200.,"24 billion events");
1032 }
1033 else {
1034 ld3 = new TLatex(10.1,myheight-200.,"14.3 billion events");
1035 }
1036 ld3->SetTextSize(0.035);
1037 ld3->Draw();
1038 
1039 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
1040 TH1D* hhfit_tmp = (TH1D*)hhfit[nbins]->Clone("hhfit_tmp");
1041 hhfit_tmp->SetAxisRange(8.0,11.);
1042 hhfit_tmp->Draw("pehist");
1043 fSandBauau->Draw("same");
1044 fmycorrbg->Draw("same");
1045 
1046 
1047 //----------------------------------------------------------------------
1048 
1049 // plot signal + both backgrounds for all pT
1050 
1051 TCanvas* callpt = new TCanvas("callpt","Signal+BG (all p_{T})",300,300,600,600);
1052 
1053   hhall_scaled[nbins]->GetXaxis()->SetTitle("Invariant mass GeV/c");
1054   hhall_scaled[nbins]->SetLineColor(kBlack);
1055   hhall_scaled[nbins]->SetMarkerColor(kBlack);
1056   hhall_scaled[nbins]->SetMarkerStyle(20);
1057   hhall_scaled[nbins]->SetAxisRange(8.0,10.8);
1058   hhall_scaled[nbins]->Draw("pehist");
1059     hhcombbg_scaled[nbins]->SetLineColor(kBlue);
1060     hhcombbg_scaled[nbins]->Draw("histsame");
1061       hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1062       hhcorrbg_scaled[nbins]->Draw("histsame");
1063 
1064 //----------------------------------------------------------------
1065 //  Calculate RAA
1066 //----------------------------------------------------------------
1067 
1068 double u1start = 9.25;
1069 double u1stop  = 9.65;
1070 double u2start = 9.80;
1071 double u2stop  = 10.20;
1072 double u3start = 10.20;
1073 double u3stop  = 10.55;
1074 
1075   double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
1076   for(int i=0; i<nbins; i++) { 
1077     //raa1[i] = supcor[0]; raa2[i] = supcor[1]; raa3[i] = supcor[2]; 
1078     raa1[i] = grRAA1S->Eval(0.5+i*1.0);
1079     raa2[i] = grRAA2S->Eval(0.5+i*1.0);
1080     raa3[i] = grRAA3S->Eval(0.5+i*1.0);
1081   }
1082   int fbin1 = hhall_scaled[nbins]->FindBin(u1start + 0.001);
1083   int lbin1 = hhall_scaled[nbins]->FindBin(u1stop  - 0.001);
1084   int fbin2 = hhall_scaled[nbins]->FindBin(u2start + 0.001);
1085   int lbin2 = hhall_scaled[nbins]->FindBin(u2stop  - 0.001);
1086   int fbin3 = hhall_scaled[nbins]->FindBin(u3start + 0.001);
1087   int lbin3 = hhall_scaled[nbins]->FindBin(u3stop  - 0.001);
1088   cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
1089   cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
1090   cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
1091   cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
1092   cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
1093   cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
1094 
1095   double sum1[99]   = {0.};
1096   double truesum1[99]   = {0.};
1097   double ersum1[99] = {0.};
1098   double sumpp1[99]   = {0.};
1099   double ersumpp1[99] = {0.};
1100   double sum2[99]   = {0.};
1101   double truesum2[99]   = {0.};
1102   double ersum2[99] = {0.};
1103   double sumpp2[99]   = {0.};
1104   double ersumpp2[99] = {0.};
1105   double sum3[99]   = {0.};
1106   double truesum3[99]   = {0.};
1107   double ersum3[99] = {0.};
1108   double sumpp3[99]   = {0.};
1109   double ersumpp3[99] = {0.};
1110 
1111   double sumsum1[99]   = {0.};
1112   double sumsum2[99]   = {0.};
1113   double sumsum3[99]   = {0.};
1114   double sumsum1pp[99]   = {0.};
1115   double sumsum2pp[99]   = {0.};
1116   double sumsum3pp[99]   = {0.};
1117 
1118   for(int i=0; i<nbins+1; i++) {
1119 
1120     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1121 
1122     for(int j=fbin1; j<=lbin1; j++) {
1123       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)));
1124       truesum1[i] += hhups1[i]->GetBinContent(j);
1125       ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1126       sumpp1[i]   += hhups1pp[i]->GetBinContent(j);
1127       ersumpp1[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1128     }
1129       sumsum1[i]     = truesum1[i];                   // direct count in mass range
1130       sumsum1pp[i]   = sumpp1[i];       
1131       //sumsum1[i]     = hhups1[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1132       //sumsum1pp[i]   = hhups1pp[i]->GetEntries();
1133       //sumsum1[i]     = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1134       //sumsum1pp[i]   = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1135 
1136         if(sumsum1[i]>0. && sumsum1pp[i]>0.) {
1137           erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1[i]/sumsum1pp[i]/sumsum1pp[i]);
1138         } else {raa1[i]=-1.0; erraa1[i] = 999.; }
1139 
1140     for(int j=fbin2; j<=lbin2; j++) {
1141       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)));
1142       truesum2[i] += hhups2[i]->GetBinContent(j);
1143       ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1144       sumpp2[i]   += hhups2pp[i]->GetBinContent(j);
1145       ersumpp2[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1146     }
1147       sumsum2[i]     = truesum2[i];                   // direct count in mass range
1148       sumsum2pp[i]   = sumpp2[i];       
1149       //sumsum2[i]     = hhups2[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1150       //sumsum2pp[i]   = hhups2pp[i]->GetEntries();
1151       //sumsum2[i]     = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1152       //sumsum2pp[i]   = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1153       
1154         if(sumsum2[i]>0. && sumsum2pp[i]>0.) {
1155           erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2[i]/sumsum2pp[i]/sumsum2pp[i]);
1156         } else {raa2[i]=-1.0; erraa2[i] = 999.; }
1157 
1158     for(int j=fbin3; j<=lbin3; j++) {
1159       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)));
1160       truesum3[i] += hhups3[i]->GetBinContent(j);
1161       ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1162       sumpp3[i]   += hhups3pp[i]->GetBinContent(j);
1163       ersumpp3[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1164     }
1165       sumsum3[i]     = truesum3[i];                   // direct count in mass range
1166       sumsum3pp[i]   = sumpp3[i];       
1167       //sumsum3[i]     = hhups3[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1168       //sumsum3pp[i]   = hhups3pp[i]->GetEntries();
1169       //sumsum3[i]     = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;   // total number of upsilons in pT bin
1170       //sumsum3pp[i]   = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1171       
1172         if(truesum3[i]>0. && sumpp3[i]>0.) {
1173           erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
1174         } else {raa3[i]=-1.0; erraa3[i] = 999.; }
1175   
1176   }
1177 
1178 cout << "====== Y(1S):" << endl;
1179   for(int i=0; i<nbins+1; i++) {
1180     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1181     cout << "   " << i << " " << truesum1[i] << "(" << Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1182          << sqrt(ersum1[i]) << " \t\t pp: " << sumpp1[i] << " +- " << sqrt(ersumpp1[i]) << endl;
1183   }
1184 cout << "====== Y(2S):" << endl;
1185   for(int i=0; i<nbins+1; i++) {
1186     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1187     cout << "   " << i << " " << truesum2[i] << "(" << Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1188          << sqrt(ersum2[i]) << " \t\t pp: " << sumpp2[i] << " +- " << sqrt(ersumpp2[i]) << endl;
1189   }
1190 cout << "====== Y(3S):" << endl;
1191   for(int i=0; i<nbins+1; i++) {
1192 
1193     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1194     cout << "   " << i << " " << truesum3[i] << "(" << Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1195          << sqrt(ersum3[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << endl;
1196   }
1197 
1198 //-------------------------------------------------
1199 //  Plot RAA
1200 //-------------------------------------------------
1201 
1202 int npts1 = 10;
1203 int npts2 = 8;
1204 int npts3 = 7;
1205 
1206 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,600,600);
1207 TH2F* hh2 = new TH2F("hh2"," ",10,0.,float(npts1),10,0.,0.8);
1208 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1209 hh2->GetXaxis()->SetTitleOffset(1.0);
1210 hh2->GetXaxis()->SetTitleColor(1);
1211 hh2->GetXaxis()->SetTitleSize(0.040);
1212 hh2->GetXaxis()->SetLabelSize(0.040);
1213 hh2->GetYaxis()->SetTitle("R_{AA}");
1214 hh2->GetYaxis()->SetTitleOffset(1.3);
1215 hh2->GetYaxis()->SetTitleSize(0.040);
1216 hh2->GetYaxis()->SetLabelSize(0.040);
1217 hh2->Draw();
1218 
1219 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.5 + double(i);}
1220 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.43 + double(i);}
1221 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.57 + double(i);}
1222 
1223 TGraphErrors* gr1 = new TGraphErrors(npts1,xx1,raa1,0,erraa1);
1224 gr1->SetMarkerStyle(20);
1225 gr1->SetMarkerColor(kBlack);
1226 gr1->SetLineColor(kBlack);
1227 gr1->SetLineWidth(2);
1228 gr1->SetMarkerSize(1.5);
1229 gr1->SetName("gr1");
1230 gr1->Draw("p");
1231 
1232 erraa2[7] = erraa2[7]*0.85;
1233 
1234 TGraphErrors* gr2 = new TGraphErrors(npts2,xx2,raa2,0,erraa2);
1235 gr2->SetMarkerStyle(20);
1236 gr2->SetMarkerColor(kRed);
1237 gr2->SetLineColor(kRed);
1238 gr2->SetLineWidth(2);
1239 gr2->SetMarkerSize(1.5);
1240 gr2->SetName("gr2");
1241 gr2->Draw("p");
1242 
1243 /*
1244 //--- 3S state -------------------
1245  double dummy[9];
1246  for(int i=0; i<9; i++) {dummy[i]=-99.;}
1247  TGraphErrors* gr3 = new TGraphErrors(6,xx3,dummy,0,erraa3);
1248  gr3->SetMarkerStyle(20);
1249  gr3->SetMarkerColor(kBlue);
1250  gr3->SetLineColor(kBlue);
1251  gr3->SetLineWidth(2);
1252  gr3->SetMarkerSize(1.5);
1253  gr3->SetName("gr3");
1254 // gr3->Draw("p");
1255 
1256 TArrow* aa[9];
1257 TLine* ll[9];
1258 //erraa3[3] = erraa3[3]*1.5;
1259 erraa3[4] = erraa3[4]*0.6;
1260 //erraa3[4] = erraa3[4]*1.5;
1261 erraa3[6] = erraa3[6]*1.2;
1262 
1263 for(int i=0; i<npts3; i++) {
1264   aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1265   aa[i]->SetLineColor(kBlue);
1266   aa[i]->SetLineWidth(2);
1267   aa[i]->Draw();
1268   ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1269   ll[i]->SetLineColor(kBlue);
1270   ll[i]->SetLineWidth(2);
1271   ll[i]->Draw();
1272 }
1273 //--- end 3S state --------------
1274 */
1275 
1276 TLegend *leg = new TLegend(0.2,0.77,0.5,0.92);
1277 leg->SetFillColor(10);
1278 leg->SetFillStyle(1001);
1279 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1280 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1281 TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1282 leg->Draw();
1283 
1284 TLatex* l1 = new TLatex(3.5,0.73,"sPHENIX Simulation");
1285 l1->Draw();
1286 
1287 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1288 lll->SetLineColor(kBlue);
1289 lll->SetLineWidth(2);
1290 lll->Draw();
1291 
1292 //==================================================================================
1293 
1294 /*
1295 fout = new TFile("test.root","recreate");
1296 for(int i=0; i<nbins+1; i++) {
1297   sprintf(tmpname,"hhfit_%d",i);
1298   hhfit[i]->Write();
1299   hhups1pp[i]->Write();
1300   hhups2pp[i]->Write();
1301   hhups3pp[i]->Write();
1302   hhupspp[i]->Write();
1303   hhups1[i]->Write();
1304   hhups2[i]->Write();
1305   hhups3[i]->Write();
1306   hhups[i]->Write();
1307   hhall_pp[i]->Write();
1308   hhcorrbg_pp[i]->Write();
1309   hhcorrbg_scaled[i]->Write();
1310 }
1311 fout->Write();
1312 fout->Close();
1313 */
1314 
1315 }
1316
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!