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