Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:13:10 

0001 double factorial(int n)
0002 {
0003   double result = 1;
0004   for(int i = n; i>0; i--)
0005     result*=i;
0006   
0007   return result;
0008 }
0009 
0010 double  integral(int nfg,int nbg,float lowerLimit, float upperLimit)
0011 {
0012   // See AN 195 by Mike Tannenbaum, equation 9, for where this comes from
0013   // x is the signal, ie. FG-BG
0014   // lowerLimit and upperLimit are just the min and max signal to be considered
0015 
0016   //double stepsize=0.5;
0017 
0018   int steps = 200;
0019   double intgrl = 0;
0020   
0021   double dx=(upperLimit-lowerLimit)/(double)steps;
0022 
0023   for(int k = 0; k<=nfg; k++) {
0024     double factorial_up = factorial(nfg+nbg-k);
0025     double factorial_down = factorial(nbg)*factorial(nfg-k)*factorial(k);
0026     double tmp = factorial_up/factorial_down*0.5*pow(0.5, nfg+nbg-k);
0027     for(int itmp=0;itmp<=steps;itmp++)
0028       {
0029     double x=lowerLimit+dx*(float)itmp;
0030     intgrl += tmp*pow(x,k)*exp(-x)*dx;
0031       }
0032   }
0033   return intgrl;
0034 }
0035 
0036 void error_fg_bg_pair(int nfg, int nbg, double& err_up,double& err_down)
0037 {
0038   //nfg: # of fg in the bin
0039   //nbg: # of bg in the bin
0040   //err_up: upside error
0041   //err_down: downside error
0042 
0043   if(nfg+nbg > 170)
0044     {
0045       cout << " error_fg_bg_pair: factorials too large, return symmetric error " << endl;
0046       err_up = sqrt(nfg+nbg);
0047       err_down = err_up;
0048       return;
0049     }
0050 
0051   int nmax = 50;
0052   if(nfg-nbg>nmax) 
0053     {
0054       //cout<<" error_fg_bg_pair: Signal > " << nmax << " counts - use symmetric errors "<<endl;
0055       err_up = sqrt(nfg+nbg);
0056       err_down = err_up;
0057       return;
0058     }
0059 
0060   
0061   if(nfg<=nbg) 
0062     {
0063       //cout<<" error_fg_bg_pair: This code does not work in the case of nfg<=nbg"<<endl;
0064       cout << "nfg = " << nfg << " nbg = " << nbg << endl;
0065       err_up=sqrt(nfg+nbg);
0066       err_down = err_up;
0067       return;
0068     }
0069   
0070   // Establish the upper signal boundary for integration - Poisson distribution
0071   // should be negligible by the time the signal reaches this number
0072   // Used only for getting the total integral
0073 
0074   float maxsignal=0;
0075   if(nfg+nbg < 10)
0076     maxsignal = (nfg+nbg)* 4;
0077   else if(nfg+nbg < 30)
0078     maxsignal = (nfg+nbg)* 3;
0079   else
0080     maxsignal = (nfg+nbg)* 2;
0081 
0082   float total=integral(nfg,nbg,0,maxsignal);
0083 
0084   float binSize=0.1;
0085 
0086   float cl_up, cl_down;
0087   int i = 0;
0088   bool is_cl_down = false;
0089   //.. calculate the lower error bars using 16% in lower tail
0090   while(1) {
0091     if(!is_cl_down) err_down = i*binSize;
0092     if(err_down<0) err_down = 0;
0093     cl_down = integral(nfg,nbg,0, err_down)/total;
0094     if(cl_down>=0.159) {
0095       is_cl_down = true;
0096       break;
0097     }
0098     i++;
0099   }
0100   
0101   i = 0;
0102   while(1) {
0103     err_up = err_down+i*binSize;
0104     cl_up = cl_down+integral(nfg,nbg,err_down, err_up)/total;
0105     if(cl_up>=0.841) break;
0106     i++;
0107   }
0108 
0109   //cout<< " nfg = " << nfg << " n bg = " << nbg << " signal = " << nfg-nbg << " signal+err_up ="<<err_up
0110   //    <<" signal-err_down = "<<err_down<<endl;
0111   //cout<<" result =" << nfg-nbg << "+" << err_up-(nfg-nbg) << "-" << (nfg-nbg)-err_down << endl;
0112 
0113   err_down = (nfg-nbg)-err_down;
0114   err_up = err_up - (nfg-nbg); 
0115 }
0116 
0117 
0118 double CBFunction(double *x, double *p)
0119 {
0120   double norm = p[0];
0121   double alpha = p[1];  // tail start
0122   double n = p[2];      // tail shape 
0123   double mu = p[3];     // upsilon mass
0124   double sigma = p[4];  // upsilon width
0125 
0126   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0127   double B = n/fabs(alpha) - fabs(alpha);
0128   double k = (x[0]-mu)/sigma;
0129 
0130   double val;
0131   if( k > -alpha )
0132     val = norm*TMath::Exp(-0.5*pow(k,2));
0133   else
0134     val = norm*A*pow(B-k,-n);
0135 
0136   if( TMath::IsNaN(val) ) val = 0.0;
0137 
0138   return val;
0139 }
0140 
0141 double TripleCBFunction(double *x, double *p)
0142 {
0143   double norm1 = p[0];  // amplitude of Y(1S) peak
0144   double alpha = p[1];  // tail start
0145   double n = p[2];      // tail shape 
0146   double mu1 = p[3];    // upsilon (1S) mass
0147   double sigma = p[4];  // upsilon width
0148   double norm2 = p[5];
0149   double norm3 = p[6];
0150   double mu2 = mu1*1.0595;
0151   double mu3 = mu1*1.0946;
0152 
0153   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0154   double B = n/fabs(alpha) - fabs(alpha);
0155   double k1 = (x[0]-mu1)/sigma;
0156   double k2 = (x[0]-mu2)/sigma;
0157   double k3 = (x[0]-mu3)/sigma;
0158 
0159   double val,val1,val2,val3;
0160 
0161   if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
0162   else { val1 = norm1*A*pow(B-k1,-n); }
0163   if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
0164   else { val2 = norm2*A*pow(B-k2,-n); }
0165   if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
0166   else { val3 = norm3*A*pow(B-k3,-n); }
0167 
0168   val = val1 + val2 + val3;
0169 
0170   if( TMath::IsNaN(val) ) val = 0.0;
0171 
0172   return val;
0173 }
0174 
0175 double SandB_CBFunction(double *x, double *p)
0176 {
0177   double norm1 = p[0];  // amplitude of Y(1S) peak
0178   double alpha = p[1];  // tail start
0179   double n = p[2];      // tail shape 
0180   double mu1 = p[3];     // upsilon (1S) mass
0181   double sigma = p[4];  // upsilon width
0182   double norm2 = p[5];
0183   double norm3 = p[6];
0184   double mu2 = mu1*1.0595;
0185   double mu3 = mu1*1.0946;
0186 
0187   double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
0188   double B = n/fabs(alpha) - fabs(alpha);
0189   double k1 = (x[0]-mu1)/sigma;
0190   double k2 = (x[0]-mu2)/sigma;
0191   double k3 = (x[0]-mu3)/sigma;
0192 
0193   double val,val1,val2,val3;
0194 
0195   if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
0196   else { val1 = norm1*A*pow(B-k1,-n); }
0197   if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
0198   else { val2 = norm2*A*pow(B-k2,-n); }
0199   if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
0200   else { val3 = norm3*A*pow(B-k3,-n); }
0201 
0202   double bgnorm1 = p[7];
0203   double bgslope1 = p[8];
0204 
0205   double bg = exp(bgnorm1+x[0]*bgslope1);
0206 
0207   val = val1 + val2 + val3 + bg;
0208   if( TMath::IsNaN(val) ) val = 0.0;
0209 
0210   return val;
0211 }
0212 
0213 //----------------------------------------------------------------
0214 //----------------------------------------------------------------
0215 //----------------------------------------------------------------
0216 
0217 void centrality_newplotbg_newsupp() {
0218 
0219 //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
0220 //SetsPhenixStyle();
0221 
0222 gStyle->SetOptStat(0);
0223 gStyle->SetOptFit(0);
0224 
0225 TRandom* myrandom = new TRandom3();
0226 const int nbins = 15;
0227 double eideff = 0.9;
0228 string times = "*";
0229 TLatex* tl[15];
0230 char tlchar[999];
0231 
0232 // Updates for 2020 BUP from Jamie (by ADF, August 14, 2020)
0233 // 3 year run plan for 24 (28) cryoweeks / year:
0234 //    AuAu year 1 25B (39B) MB events
0235 //    AuAu year 3 85B (103B) MB evemts
0236 //    Total AuAu 110B (142) MB events (i.e. 11B (14.2B) in 10% central) 
0237 //    old: pp year 2 73 /Pb (101 /Pb) or 3T (4T) sampled events
0238 //    updated: pp year (28 wks) 2 62 /Pb  or 2.6T  sampled events
0239 // 5 year run plan
0240 //    AuAu year 5 adds 206B MB events for a total of 316B (348B) (i.e. 31.6B (34.8B) in 10% central)
0241 //    old: pp year 4 adds 130 /Pb sampled or 5.5T sampled events for a total of 8.5T (9.6T)
0242 //   updated:  pp year 4 adds 80 /Pb sampled or 3.4T sampled events for a total of 5.9T 
0243 
0244 // statscale factor is relative to 100B events for AuAu, and relative to 7.35T events in pp
0245  double set_statscale[4] = {3.48, 3.16, 1.42, 1.1};  
0246  double set_statscalepp[4] = {0.8, 1.16, 0.35, 0.41};
0247  std::string run_plan[4] = {"5 years/28 cryo-weeks", "5 years/24 cryo-weeks","3 years/28 cryo-weeks","3 years/24 cryo-weeks"};
0248  int scenario = 0;  // scenario = {0,1,2,3} = {5yrs_28wks, 5yrs_24wks, 3yrs_28wks, 3yrs_24wks}  (BUP 2020 uses 2 and 0)
0249  double statscale = set_statscale[scenario]; 
0250  double statscalepp = set_statscalepp[scenario]; 
0251 
0252  int auauevts = (int) (statscale * 10.0);  // 10% most central events = statscale * 10% of 100B
0253  
0254  double statscale_lowlim = 7.0;
0255  double statscale_uplim = 14.0;
0256 
0257   TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
0258   TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
0259   TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
0260   TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
0261   TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
0262   TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
0263   TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
0264   TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
0265   TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
0266   TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
0267   TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
0268   TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
0269 
0270 //---------------------------------------------------------
0271 // Upsilons
0272 //---------------------------------------------------------
0273 
0274 // Need different raa at each centrality
0275 
0276   // this is for the Y(2S) and Y(3S) centrality dependence binning
0277   static const int NCENT = 7;
0278   double centlow[NCENT] = {0,5,10,20,30,40,60};  
0279   double centhigh[NCENT] = {5, 10,20,30,40,60,92};
0280   double Ncoll[NCENT] = {1067, 858, 610, 378, 224, 94.2, 14.5};
0281 
0282   double raacent_ups1s[NCENT] = {0.527483, 0.544815, 0.581766, 0.644214, 0.721109, 0.853888, 0.998349}; 
0283   double raacent_ups2s[NCENT] = {0.165314, 0.179333, 0.210351, 0.26748, 0.34664, 0.510807, 0.979886}; 
0284   double raacent_ups3s[NCENT] = {0.0302562, 0.0362929, 0.047993, 0.0706653, 0.108092, 0.231598, 0.961714};
0285 
0286 
0287 /*
0288   // not used
0289   static const int NCENT = 4;
0290   double centlow[NCENT] = {0,20,40,60};
0291   double centhigh[NCENT] = {20,40,60,92};
0292   double Ncoll[NCENT] = {783, 301, 94.2, 14.5};
0293   double raacent_ups3s[4] = {0.038949, 0.0838922, 0.220475, 0.924073};
0294   double raacent_ups1s[NCENT] = {0.527483, 0.544815, 0.581766, 0.644214}; // dummy 
0295   double raacent_ups2s[NCENT] = {0.165314, 0.179333, 0.210351, 0.26748};  // dummy 
0296 */
0297 
0298 /*
0299   // this is to get the 0-10% pT dependence 
0300   static const int NCENT = 6;
0301   double centlow[NCENT] = {0,10,20,30,40,60};
0302   double centhigh[NCENT] = {10,20,30,40,60,92};
0303   double Ncoll[NCENT] = {962, 610, 378, 224, 94.2, 14.5};
0304 */
0305 
0306 /*
0307   // this is for the Y(3S) centrality dependence binning 
0308   static const int NCENT = 3;
0309   double centlow[NCENT] = {0,30,60};
0310   double centhigh[NCENT] = {30,60,92};
0311   double Ncoll[NCENT] = {648, 137, 14.5};
0312   double raacent_ups1s[NCENT] = {0.527483, 0.544815, 0.644214}; // dummy 
0313   double raacent_ups2s[NCENT] = {0.165314, 0.179333, 0.26748};  // dummy 
0314   double raacent_ups3s[3] = {0.0458138, 0.167527, 0.924073};
0315 */
0316 
0317   double centwidth[NCENT];
0318   for(int i=0; i<NCENT; ++i)
0319     {
0320       centwidth[i] = (centhigh[i] - centlow[i]) / 100.0;
0321       cout << " cent bin " << i << " has width " << centwidth[i] << endl;
0322     }
0323 
0324   int icent = 6;
0325 
0326   double raapt[9],raa1s[9],raa2s[9],raa3s[9];
0327 raapt[0] = 1.5;
0328 raapt[1] = 4.5;
0329 raapt[2] = 7.5;
0330 raapt[3] = 10.5;
0331 raapt[4] = 13.5;
0332 raa1s[0] = raacent_ups1s[icent];
0333 raa1s[1] = raacent_ups1s[icent];
0334 raa1s[2] = raacent_ups1s[icent];
0335 raa1s[3] = raacent_ups1s[icent];
0336 raa1s[4] = raacent_ups1s[icent];
0337 raa2s[0] = raacent_ups2s[icent];
0338 raa2s[1] = raacent_ups2s[icent];
0339 raa2s[2] = raacent_ups2s[icent];
0340 raa2s[3] = raacent_ups2s[icent];
0341 raa2s[4] = raacent_ups2s[icent];
0342 raa3s[0] = raacent_ups3s[icent];
0343 raa3s[1] = raacent_ups3s[icent];
0344 raa3s[2] = raacent_ups3s[icent];
0345 raa3s[3] = raacent_ups3s[icent];
0346 raa3s[4] = raacent_ups3s[icent];
0347 
0348 /*
0349 raa1s[0] = 0.535;
0350 raa1s[1] = 0.535;
0351 raa1s[2] = 0.535;
0352 raa1s[3] = 0.535;
0353 raa1s[4] = 0.535;
0354 raa2s[0] = 0.170;
0355 raa2s[1] = 0.170;
0356 raa2s[2] = 0.170;
0357 raa2s[3] = 0.170;
0358 raa2s[4] = 0.170;
0359 raa1s[0] = 0.4960;
0360 raa1s[1] = 0.4960;
0361 raa1s[2] = 0.4955;
0362 raa1s[3] = 0.4968;
0363 raa1s[4] = 0.4743;
0364 raa2s[0] = 0.1710;
0365 raa2s[1] = 0.1629;
0366 raa2s[2] = 0.1326;
0367 raa2s[3] = 0.1232;
0368 raa2s[4] = 0.0928;
0369 raa3s[0] = 0.035;
0370 raa3s[1] = 0.035;
0371 raa3s[2] = 0.035;
0372 raa3s[3] = 0.035;
0373 raa3s[4] = 0.035;
0374 */
0375 
0376 TGraph* grRAA1S = new TGraph(5,raapt,raa1s);
0377 TGraph* grRAA2S = new TGraph(5,raapt,raa2s);
0378 TGraph* grRAA3S = new TGraph(5,raapt,raa3s);
0379 
0380 double NNN = statscale * 7780./0.49;      // from sPHENIX proposal for 100B minbias Au+Au events
0381 double NNNpp = statscalepp * 12130./0.90; // Upsilons in p+p from sPHENIX proposal for 7350B p+p events
0382 double ppauaurat = NNNpp/NNN;
0383 
0384 // NNN is the unsuppressed number of Upsilons in 0-10% most central events for 100B MB AuAu collisions 
0385 // That is for Ncoll = 962, for any other centrality, scale it by Ncoll
0386 // Also need to scale by number of events, i.e. centrality bin width
0387  double Ncoll_ref = (Ncoll[0] + Ncoll[1]) / 2.0;  // reference is 0-10% centrality, so average 0-5 and 5-10
0388  double ncollfact = (Ncoll[icent] / Ncoll_ref);
0389  double centwidth_ref = (centwidth[0]+centwidth[1]);  // reference is 0-10% centrality, so sum 0-5 and 5-10
0390  double centwidthfact = (centwidth[icent] / centwidth_ref);
0391  NNN *= ncollfact*centwidthfact;   // unsuppressed Upsilons for this centrality bin
0392 
0393 double frac[3];  // upsilons fractions
0394   frac[0] = 0.7117;
0395   frac[1] = 0.1851;
0396   frac[2] = 0.1032;
0397 double scale[3]; // mass scaling
0398   scale[0] = 1.0;
0399   scale[1] = 1.0595;
0400   scale[2] = 1.0946;
0401 //double supcor[3]; // suppression
0402 //  supcor[0]=0.535; 
0403 //  supcor[1]=0.17; 
0404 //  supcor[2]=0.035;
0405 //int Nups1 = int(NNN*eideff*eideff*frac[0]*supcor[0]);
0406 //int Nups2 = int(NNN*eideff*eideff*frac[1]*supcor[1]);
0407 //int Nups3 = int(NNN*eideff*eideff*frac[2]*supcor[2]);
0408 int Nups1nosup = int(NNN*eideff*eideff*frac[0]);
0409 int Nups2nosup = int(NNN*eideff*eideff*frac[1]);
0410 int Nups3nosup = int(NNN*eideff*eideff*frac[2]);
0411 int Nups1pp = int(NNNpp*0.90*frac[0]); // always 95% eideff in p+p
0412 int Nups2pp = int(NNNpp*0.90*frac[1]);
0413 int Nups3pp = int(NNNpp*0.90*frac[2]);
0414 
0415 int nchan=400;
0416 double start=0.0;
0417 double stop=20.0;
0418 
0419 string str_UpsilonPt = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))";
0420 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
0421 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
0422 double upsnorm = fUpsilonPt->Integral(0.,20.);
0423 
0424 // count all upsilons with suppression
0425 double tmpsum1=0.;
0426 double tmpsum2=0.;
0427 double tmpsum3=0.;
0428 for(int j=0; j<nbins; j++) {
0429   double s1 = j*1.0;
0430   double s2 = s1 + 1.0;
0431   tmpsum1 += Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
0432   tmpsum2 += Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
0433   tmpsum3 += Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
0434 }
0435 int Nups1 = int(tmpsum1);
0436 int Nups2 = int(tmpsum2);
0437 int Nups3 = int(tmpsum3);
0438 
0439   cout << "Total number of ALL Upsilons in AuAu = " << NNN << endl;
0440   cout << "Total number of ALL Upsilons in pp = " << NNNpp << endl;
0441   cout << "Number of upsilons in pp in acceptance = " << NNNpp*frac[0] << " " << NNNpp*frac[1] << " " << NNNpp*frac[2] << endl;
0442   cout << "Number of upsilons in AuAu in acceptance = " << NNN*frac[0] << " " << NNN*frac[1] << " " << NNN*frac[2] << endl;
0443   cout << "Number of upsilons in AuAu after eID efficiency = " << Nups1nosup << " " << Nups2nosup << " " << Nups3nosup << endl;
0444   cout << "Number of upsilons in AuAu plot = " << Nups1 << " " << Nups2 << " " << Nups3 << endl;
0445   cout << "Number of upsilons in pp plot = " << Nups1pp << " " << Nups2pp << " " << Nups3pp << endl;
0446 
0447 //====================================================================================================
0448 /*
0449 f=new TFile("UpsilonsInHijing.root");
0450   cout << endl << endl << "Number of entries in Upsilon NTuple = " << ntpups->GetEntries() << endl;
0451 
0452   TH1D* hforfit = new TH1D("hforfit","",nchan,start,stop);
0453   TH1D* hUps1 = new TH1D("hUps1","",nchan,start,stop);
0454   TH1D* hUps2 = new TH1D("hUps2","",nchan,start,stop);
0455   TH1D* hUps3 = new TH1D("hUps3","",nchan,start,stop);
0456   TH1D* hUps1pp = new TH1D("hUps1pp","",nchan,start,stop);
0457   TH1D* hUps2pp = new TH1D("hUps2pp","",nchan,start,stop);
0458   TH1D* hUps3pp = new TH1D("hUps3pp","",nchan,start,stop);
0459     hforfit->Sumw2();
0460     hUps1->Sumw2();
0461     hUps2->Sumw2();
0462     hUps3->Sumw2();
0463     hUps1pp->Sumw2();
0464     hUps2pp->Sumw2();
0465     hUps3pp->Sumw2();
0466 
0467 //------ Additional smearing to reproduce mass resolution -----
0468 //double smear = 0.098; // 127
0469 //double smear = 0.090;   // 119
0470 //double smear = 0.067;   // 100 MeV benchmark
0471 double smear = 0.046;  // 86 MeV resolution, Aalysis Note result 
0472 float mass;
0473 TChain* tree = new TChain("ntpups");
0474 tree->Add("UpsilonsInHijing.root");
0475 tree->SetBranchAddress("mass",&mass);
0476 int countups=0;
0477   for(int j=0; j<tree->GetEntries(); j++){
0478     tree->GetEvent(j);
0479     double newmass = myrandom->Gaus(mass,smear);
0480     hforfit->Fill(newmass);
0481     if(countups<=Nups1) hUps1->Fill(newmass);
0482     if(countups<=Nups2) hUps2->Fill(newmass*1.0595);
0483     if(countups<=Nups3) hUps3->Fill(newmass*1.0946);
0484     if(countups<=Nups1pp) hUps1pp->Fill(newmass);
0485     if(countups<=Nups2pp) hUps2pp->Fill(newmass*1.0595);
0486     if(countups<=Nups3pp) hUps3pp->Fill(newmass*1.0946);
0487     countups++;
0488   }
0489 delete tree;
0490 //---------------------------------------------------------------
0491 
0492   hforfit->SetDirectory(gROOT);
0493   hUps1->SetDirectory(gROOT);
0494   hUps2->SetDirectory(gROOT);
0495   hUps3->SetDirectory(gROOT);
0496   hUps1pp->SetDirectory(gROOT);
0497   hUps2pp->SetDirectory(gROOT);
0498   hUps3pp->SetDirectory(gROOT);
0499 
0500 f->Close();
0501 
0502 TH1F* hUps = (TH1F*)hUps1->Clone("hUps");
0503 hUps->Add(hUps2);
0504 hUps->Add(hUps3);
0505 hUps->SetLineWidth(2);
0506 TH1F* hUpspp = (TH1F*)hUps1pp->Clone("hUpspp");
0507 hUpspp->Add(hUps2pp);
0508 hUpspp->Add(hUps3pp);
0509 hUpspp->SetLineWidth(2);
0510 
0511 TCanvas* cups = new TCanvas("cups","Upsilons (all pT)",100,100,1000,500);
0512 cups->Divide(2,1);
0513 
0514   cups_1->cd();
0515   cups_1->SetLogy();
0516   fCB->SetParameter(0,5000.);
0517   fCB->SetParameter(1,1.5); 
0518   fCB->SetParameter(2,1.2); 
0519   fCB->SetParameter(3,9.448);
0520   fCB->SetParameter(4,0.070);
0521   hforfit->Fit("fCB","rl","",5.,9.7);
0522   hforfit->SetAxisRange(5.,11.);
0523   hforfit->SetLineWidth(2);
0524   hforfit->Draw();
0525   double myupsmass = fCB->GetParameter(3);
0526   double upswidth = fCB->GetParameter(4);
0527 
0528   cups_2->cd();
0529   cups_2->SetLogy();
0530   fTCB->SetParameter(0,2000.);
0531   fTCB->SetParameter(1,1.5);
0532   fTCB->SetParameter(2,1.25);
0533   fTCB->SetParameter(3,9.448);
0534   fTCB->FixParameter(4,upswidth); // fix width from the single peak fit
0535   fTCB->SetParameter(5,500.);
0536   fTCB->SetParameter(6,100.);
0537   hUps->Fit(fTCB,"rl","",6.,11.);
0538   double upswidth3 = fTCB->GetParameter(4);
0539   double upswidth3err = fTCB->GetParError(4);
0540   hUps->SetAxisRange(7.0,11.0);
0541   hUps->Draw();
0542 
0543 cout << "MASS RESOLUTION = " << upswidth << " " << upswidth3 << " +- " << upswidth3err << endl;
0544 */
0545 //====================================================================================
0546 //====================================================================================
0547 //====================================================================================
0548 
0549 //
0550 // these histograms (hhups*) are randomly generated using the fit above to single 
0551 // peak (fCB function) and used for the RAA uncertainty calculation
0552 //
0553 // FORCE specific RESOLUTION and tail shape
0554   double tonypar1 = 0.98;  // Tony's 100 pion simulation April 2019
0555   double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
0556   //double tonypar3 = 9.437;  // Tony's 100 pion simulation April 2019
0557   double tonypar3 = 9.448;  // benchmark
0558   double tonypar4 = 0.100;  // benchmark
0559   double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
0560   fCBpp->SetParameter(0,1000.); 
0561   fCBpp->SetParameter(1,tonypar1); 
0562   fCBpp->SetParameter(2,tonypar2); 
0563   fCBpp->SetParameter(3,tonypar3); 
0564   fCBpp->SetParameter(4,tonypar4pp);  
0565   fCBauau->SetParameter(0,1000.); 
0566   fCBauau->SetParameter(1,tonypar1); 
0567   fCBauau->SetParameter(2,tonypar2); 
0568   fCBauau->SetParameter(3,tonypar3); 
0569   fCBauau->SetParameter(4,tonypar4);  
0570 
0571 
0572 char hhname[999];
0573 TH1D* hhups[nbins+1];
0574 TH1D* hhups1[nbins+1];
0575 TH1D* hhups2[nbins+1];
0576 TH1D* hhups3[nbins+1];
0577 TH1D* hhupspp[nbins+1];
0578 TH1D* hhups1pp[nbins+1];
0579 TH1D* hhups2pp[nbins+1];
0580 TH1D* hhups3pp[nbins+1];
0581 for(int i=0; i<nbins+1; i++) {
0582   sprintf(hhname,"hhups_%d",i);
0583   hhups[i] = new TH1D(hhname,"",nchan,start,stop);
0584   hhups[i]->Sumw2();
0585   sprintf(hhname,"hhups1_%d",i);
0586   hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
0587   hhups1[i]->Sumw2();
0588   sprintf(hhname,"hhups2_%d",i);
0589   hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
0590   hhups2[i]->Sumw2();
0591   sprintf(hhname,"hhups3_%d",i);
0592   hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
0593   hhups3[i]->Sumw2();
0594   hhups[i]->SetLineWidth(2);
0595   hhups1[i]->SetLineWidth(2);
0596   hhups2[i]->SetLineWidth(2);
0597   hhups3[i]->SetLineWidth(2);
0598   sprintf(hhname,"hhupspp_%d",i);
0599   hhupspp[i] = new TH1D(hhname,"",nchan,start,stop);
0600   hhupspp[i]->Sumw2();
0601   sprintf(hhname,"hhups1pp_%d",i);
0602   hhups1pp[i] = new TH1D(hhname,"",nchan,start,stop);
0603   hhups1pp[i]->Sumw2();
0604   sprintf(hhname,"hhups2pp_%d",i);
0605   hhups2pp[i] = new TH1D(hhname,"",nchan,start,stop);
0606   hhups2pp[i]->Sumw2();
0607   sprintf(hhname,"hhups3pp_%d",i);
0608   hhups3pp[i] = new TH1D(hhname,"",nchan,start,stop);
0609   hhups3pp[i]->Sumw2();
0610   hhupspp[i]->SetLineWidth(2);
0611   hhups1pp[i]->SetLineWidth(2);
0612   hhups2pp[i]->SetLineWidth(2);
0613   hhups3pp[i]->SetLineWidth(2);
0614 }
0615 
0616 for(int j=0; j<nbins; j++) {
0617 
0618   double s1 = j*1.0;
0619   double s2 = s1 + 1.0;
0620   double tmpnups1 = Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
0621   double tmpnups2 = Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
0622   double tmpnups3 = Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
0623   //cout << "Nups1["<<j<<"]="<<tmpnups1<<";"<<endl;
0624   //cout << "Nups2["<<j<<"]="<<tmpnups2<<";"<<endl;
0625   //cout << "Nups3["<<j<<"]="<<tmpnups3<<";"<<endl;
0626   fCBauau->SetParameter(3,tonypar3);
0627     for(int i=0; i<int(tmpnups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0628   fCBauau->SetParameter(3,tonypar3*scale[1]);
0629     for(int i=0; i<int(tmpnups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0630   fCBauau->SetParameter(3,tonypar3*scale[2]);
0631     for(int i=0; i<int(tmpnups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0632 
0633   double tmpnups1pp = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0634   double tmpnups2pp = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0635   double tmpnups3pp = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
0636   cout << "Nups1pp["<<j<<"]="<<tmpnups1pp<<";"<<endl;
0637   cout << "Nups2pp["<<j<<"]="<<tmpnups2pp<<";"<<endl;
0638   cout << "Nups3pp["<<j<<"]="<<tmpnups3pp<<";"<<endl;
0639   fCBpp->SetParameter(3,tonypar3);
0640     for(int i=0; i<int(tmpnups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0641   fCBpp->SetParameter(3,tonypar3*scale[1]);
0642     for(int i=0; i<int(tmpnups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0643   fCBpp->SetParameter(3,tonypar3*scale[2]);
0644     for(int i=0; i<int(tmpnups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
0645 
0646 } // end loop over pT bins
0647 
0648 // all pT
0649 
0650   fCBpp->SetParameter(3,tonypar3);
0651   fCBauau->SetParameter(3,tonypar3);
0652     for(int i=0; i<int(Nups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0653     for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0654   fCBpp->SetParameter(3,tonypar3*scale[1]);
0655   fCBauau->SetParameter(3,tonypar3*scale[1]);
0656     for(int i=0; i<int(Nups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0657     for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0658   fCBpp->SetParameter(3,tonypar3*scale[2]);
0659   fCBauau->SetParameter(3,tonypar3*scale[2]);
0660     for(int i=0; i<int(Nups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
0661     for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
0662 
0663 // fit and draw pp no bg
0664 
0665 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
0666   fTCBpp->SetParameter(0,2000.);
0667   fTCBpp->FixParameter(1,tonypar1);
0668   fTCBpp->FixParameter(2,tonypar2);
0669   fTCBpp->SetParameter(3,tonypar3);
0670   fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
0671   fTCBpp->SetParameter(5,500.);
0672   fTCBpp->SetParameter(6,100.);
0673   hhupspp[nbins]->Fit(fTCBpp,"rl","",7.,11.);
0674   hhupspp[nbins]->SetAxisRange(7.,11.);
0675   hhupspp[nbins]->SetMarkerSize(1.0);
0676   hhupspp[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0677   hhupspp[nbins]->GetXaxis()->SetTitleOffset(1.0);
0678   double tmpamp1 = hhupspp[nbins]->GetFunction("fTCBpp")->GetParameter(0);
0679   double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0680   double tmpamp6 = tmpamp1*frac[2]/frac[0];
0681   hhupspp[nbins]->Draw();
0682 
0683   fCB1s->SetLineColor(kBlue);
0684   fCB1s->SetLineWidth(1);
0685     fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
0686     fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
0687     fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
0688     fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
0689     fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
0690   fCB2s->SetLineColor(kRed);
0691   fCB2s->SetLineWidth(1);
0692     fCB2s->SetParameter(0,tmpamp5);
0693     fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
0694     fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
0695     fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
0696     fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
0697   fCB3s->SetLineColor(kGreen+2);
0698   fCB3s->SetLineWidth(1);
0699     fCB3s->SetParameter(0,tmpamp6);
0700     fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
0701     fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
0702     fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
0703     fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
0704   fCB1s->Draw("same");
0705   fCB2s->Draw("same");
0706   fCB3s->Draw("same");
0707 
0708 // fit and draw AuAu no bg
0709 
0710 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Au+Au",100,100,600,600);
0711   fTCBauau->SetParameter(0,2000.);
0712   fTCBauau->FixParameter(1,tonypar1);
0713   fTCBauau->FixParameter(2,tonypar2);
0714   fTCBauau->SetParameter(3,tonypar3);
0715   fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
0716   fTCBauau->SetParameter(5,500.);
0717   fTCBauau->SetParameter(6,100.);
0718   hhups[nbins]->Fit(fTCBauau,"rl","",7.,11.);
0719   hhups[nbins]->SetAxisRange(7.,11.);
0720   hhups[nbins]->SetMarkerSize(1.0);
0721   hhups[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0722   hhups[nbins]->GetXaxis()->SetTitleOffset(1.0);
0723   tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
0724   tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0725   tmpamp6 = tmpamp1*frac[2]/frac[0];
0726   hhups[nbins]->Draw();
0727 
0728 TCanvas* cupsvspt = new TCanvas("cupsvspt","Upsilons vs. p_{T}",100,100,1200,900);
0729 cupsvspt->Divide(4,3);
0730 for(int i=0; i<12; i++) {
0731 if(i==0)  {cupsvspt->cd(1);}
0732 if(i==1)  {cupsvspt->cd(2);}
0733 if(i==2)  {cupsvspt->cd(3);}
0734 if(i==3)  {cupsvspt->cd(4);}
0735 if(i==4)  {cupsvspt->cd(5);}
0736 if(i==5)  {cupsvspt->cd(6);}
0737 if(i==6)  {cupsvspt->cd(7);}
0738 if(i==7)  {cupsvspt->cd(8);}
0739 if(i==8)  {cupsvspt->cd(9);}
0740 if(i==9)  {cupsvspt->cd(10);}
0741 if(i==10) {cupsvspt->cd(11);}
0742 if(i==11) {cupsvspt->cd(12);}
0743 hhups[i]->SetAxisRange(7.0,11.0); hhups[i]->Draw();
0744   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhups[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0745 }
0746 
0747 //----------------------------------------------------
0748 //  Backgrounds
0749 //----------------------------------------------------
0750 
0751 TH1D* hhall[nbins+1];
0752 TH1D* hhall_scaled[nbins+1];
0753 
0754 TH1D* hhtotbg[nbins+1]; 
0755 TH1D* hhtotbg_scaled[nbins+1]; 
0756 TH1D* hhcombbg[nbins+1];
0757 TH1D* hhcombbg_scaled[nbins+1];
0758 TH1D* hhfakefake[nbins+1];
0759 TH1D* hhfakehf[nbins+1];
0760 TH1D* hhbottom[nbins+1];
0761 TH1D* hhcharm[nbins+1];
0762 TH1D* hhdy[nbins+1];
0763 TH1D* hhcorrbg[nbins+1];
0764 TH1D* hhcorrbg_scaled[nbins+1];
0765 TH1D* hhfit[nbins+1];
0766 char tmpname[999];
0767 
0768 //----------------------------------------------------------------------------------------
0769 // Correlated BG calculated for 10B central AuAu events
0770 //----------------------------------------------------------------------------------------
0771 
0772 // The correlated background should scale with Ncoll, right?
0773 
0774 double corrbgfitpar0;
0775 double corrbgfitpar1;
0776 
0777 TFile* f=new TFile("ccbb_eideff09.root");
0778 for(int i=0; i<nbins+1; i++) {
0779   sprintf(tmpname,"hhbottom_%d",i);
0780   hhbottom[i] = (TH1D*)f->Get(tmpname);
0781   hhbottom[i]->SetDirectory(gROOT);
0782   sprintf(tmpname,"hhcharm_%d",i);
0783   hhcharm[i] = (TH1D*)f->Get(tmpname);
0784   hhcharm[i]->SetDirectory(gROOT);
0785   sprintf(tmpname,"hhdy_%d",i);
0786   hhdy[i] = (TH1D*)f->Get(tmpname);
0787   hhdy[i]->SetDirectory(gROOT);
0788   sprintf(tmpname,"hhcorrbg_%d",i);
0789   hhcorrbg[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0790   hhcorrbg[i]->Add(hhcharm[i]);
0791   hhcorrbg[i]->Add(hhdy[i]);
0792   sprintf(tmpname,"hhcorrbg_scaled_%d",i);
0793   hhcorrbg_scaled[i] = (TH1D*)hhcorrbg[i]->Clone(tmpname);
0794   hhcorrbg[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0795   hhbottom[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0796   hhdy[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0797     if(i==nbins) {
0798        corrbgfitpar0 = hhcorrbg[i]->GetFunction("expo")->GetParameter(0);
0799        corrbgfitpar1 = hhcorrbg[i]->GetFunction("expo")->GetParameter(1);
0800     }
0801     cout << "bgpar0["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(0)+TMath::Log(statscale)<<";"<< endl; 
0802     cout << "bgpar1["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(1)<<";"<< endl; 
0803     for(int k=1; k<=hhcorrbg[i]->GetNbinsX(); k++) {
0804       if(hhcorrbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[i]->GetBinLowEdge(k)+hhcorrbg[i]->GetBinWidth(k))>statscale_uplim) {
0805         hhcorrbg_scaled[i]->SetBinContent(k,0.); 
0806         hhcorrbg_scaled[i]->SetBinError(k,0.);
0807       }
0808       else {
0809     // assumes corr bg scales with Ncoll
0810         double tmp = ncollfact * centwidthfact * statscale * hhcorrbg[i]->GetFunction("expo")->Eval(hhcorrbg[i]->GetBinCenter(k));
0811         double tmprnd = myrandom->Poisson(tmp); 
0812         if(tmprnd<0.) { tmprnd=0.; }
0813         hhcorrbg_scaled[i]->SetBinContent(k,tmprnd); 
0814         hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0815       }
0816     }
0817   hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0818   hhcorrbg[i]->SetDirectory(gROOT);
0819   hhcorrbg_scaled[i]->SetDirectory(gROOT);
0820 }
0821 f->Close();
0822 
0823 TCanvas* c2 = new TCanvas("c2","Correlated BG (all p_{T})",100,100,600,600);
0824   
0825 hhbottom[nbins]->SetAxisRange(7.0,14.0);
0826 hhbottom[nbins]->SetLineColor(kBlue);
0827 hhbottom[nbins]->SetLineWidth(2);
0828 
0829 hhbottom[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0830 hhbottom[nbins]->GetXaxis()->SetTitleOffset(1.0);
0831 hhbottom[nbins]->GetXaxis()->SetTitleColor(1);
0832 hhbottom[nbins]->GetXaxis()->SetTitleSize(0.040);
0833 hhbottom[nbins]->GetXaxis()->SetLabelSize(0.040);
0834 //hhbottom[nbins]->GetYaxis()->SetTitle("Correlated background");
0835 hhbottom[nbins]->GetYaxis()->SetTitleOffset(1.3);
0836 hhbottom[nbins]->GetYaxis()->SetTitleSize(0.040);
0837 hhbottom[nbins]->GetYaxis()->SetLabelSize(0.040);
0838 
0839 hhcorrbg[nbins]->SetAxisRange(7.0,14.0);
0840 hhcorrbg[nbins]->SetMinimum(0.1);
0841 hhcorrbg[nbins]->SetLineColor(kBlack);
0842 hhcorrbg[nbins]->SetLineWidth(2);
0843 hhcorrbg[nbins]->Draw("hist");
0844 
0845 hhbottom[nbins]->SetMarkerColor(kBlue);
0846 hhbottom[nbins]->SetLineColor(kBlue);
0847 hhbottom[nbins]->Draw("same");
0848 
0849 hhdy[nbins]->SetMarkerColor(kGreen+2);
0850 hhdy[nbins]->SetLineColor(kGreen+2);
0851 hhdy[nbins]->SetLineWidth(2);
0852 hhdy[nbins]->Draw("same");
0853 
0854 hhcharm[nbins]->SetMarkerColor(kRed);
0855 hhcharm[nbins]->SetLineColor(kRed);
0856 hhcharm[nbins]->SetLineWidth(2);
0857 hhcharm[nbins]->Draw("same");
0858 
0859 TCanvas* c0 = new TCanvas("c0","Correlated BG vs. p_{T} 10B events",100,100,1200,900);
0860 c0->Divide(4,3);
0861 for(int i=0; i<nbins; i++) {
0862 if(i==0)  {c0->cd(1);}
0863 if(i==1)  {c0->cd(2);}
0864 if(i==2)  {c0->cd(3);}
0865 if(i==3)  {c0->cd(4);}
0866 if(i==4)  {c0->cd(5);}
0867 if(i==5)  {c0->cd(6);}
0868 if(i==6)  {c0->cd(7);}
0869 if(i==7)  {c0->cd(8);}
0870 if(i==8)  {c0->cd(9);}
0871 if(i==9)  {c0->cd(10);}
0872 if(i==10) {c0->cd(11);}
0873 if(i==11) {c0->cd(12);}
0874   hhcorrbg[i]->SetAxisRange(7.,14.);
0875   hhcorrbg[i]->GetFunction("expo")->SetLineColor(kBlack); 
0876   hhbottom[i]->GetFunction("expo")->SetLineColor(kBlue); 
0877   hhdy[i]->GetFunction("expo")->SetLineColor(kGreen+2); 
0878   hhcorrbg[i]->SetLineColor(kBlack); 
0879   hhbottom[i]->SetLineColor(kBlue); 
0880   hhdy[i]->SetLineColor(kGreen+2); 
0881   hhcorrbg[i]->Draw();
0882   hhbottom[i]->Draw("same");
0883   hhdy[i]->Draw("same");
0884   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0885 }
0886 
0887 TCanvas* c0scaled = new TCanvas("c0scaled","SCALED Correlated BG vs. p_{T}",100,100,1200,900);
0888 c0scaled->Divide(4,3);
0889 for(int i=0; i<nbins; i++) {
0890 if(i==0)  {c0scaled->cd(1);}
0891 if(i==1)  {c0scaled->cd(2);}
0892 if(i==2)  {c0scaled->cd(3);}
0893 if(i==3)  {c0scaled->cd(4);}
0894 if(i==4)  {c0scaled->cd(5);}
0895 if(i==5)  {c0scaled->cd(6);}
0896 if(i==6)  {c0scaled->cd(7);}
0897 if(i==7)  {c0scaled->cd(8);}
0898 if(i==8)  {c0scaled->cd(9);}
0899 if(i==9)  {c0scaled->cd(10);}
0900 if(i==10) {c0scaled->cd(11);}
0901 if(i==11) {c0scaled->cd(12);}
0902 hhcorrbg_scaled[i]->SetAxisRange(7.0,14.0); hhcorrbg_scaled[i]->Draw();
0903   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
0904 }
0905 
0906 //-----------------------------------------------------------------------------------------
0907 //  Correlated bg in p+p
0908 //-----------------------------------------------------------------------------------------
0909 
0910 TH1D* hhbottom_pp[nbins+1];
0911 TH1D* hhdy_pp[nbins+1];
0912 TH1D* hhcorrbg_pp[nbins+1];
0913 TH1D* hhall_pp[nbins+1];
0914 
0915 double ppcorr = 8300./10./962.;
0916 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
0917 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
0918 
0919 for(int i=0; i<nbins+1; i++) {
0920 
0921   sprintf(tmpname,"hhbottom_pp_%d",i);
0922   hhbottom_pp[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
0923   for(int k=1; k<=hhbottom_pp[i]->GetNbinsX(); k++) {
0924     if(hhbottom_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp[i]->GetBinLowEdge(k)+hhbottom_pp[i]->GetBinWidth(k))>statscale_uplim) {
0925       hhbottom_pp[i]->SetBinContent(k,0.);
0926       hhbottom_pp[i]->SetBinError(k,0.);
0927     }
0928     else {
0929       double tmp = ppcorr * fbottom_nosup_corr->Eval(hhbottom[i]->GetBinCenter(k)) * hhbottom[i]->GetFunction("expo")->Eval(hhbottom[i]->GetBinCenter(k));
0930       double tmprnd = myrandom->Poisson(tmp);
0931       if(tmprnd<0.) { tmprnd=0.; }
0932       hhbottom_pp[i]->SetBinContent(k,tmprnd);
0933       hhbottom_pp[i]->SetBinError(k,sqrt(tmprnd));
0934     }
0935   }
0936 
0937   sprintf(tmpname,"hhdy_pp_%d",i);
0938   hhdy_pp[i] = (TH1D*)hhdy[i]->Clone(tmpname);
0939   for(int k=1; k<=hhdy_pp[i]->GetNbinsX(); k++) {
0940     if(hhdy_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp[i]->GetBinLowEdge(k)+hhdy_pp[i]->GetBinWidth(k))>statscale_uplim) {
0941       hhdy_pp[i]->SetBinContent(k,0.);
0942       hhdy_pp[i]->SetBinError(k,0.);
0943     }
0944     else {
0945       double tmp = ppcorr * hhdy[i]->GetFunction("expo")->Eval(hhdy[i]->GetBinCenter(k));
0946       double tmprnd = myrandom->Poisson(tmp);
0947       if(tmprnd<0.) { tmprnd=0.; }
0948       hhdy_pp[i]->SetBinContent(k,tmprnd);
0949       hhdy_pp[i]->SetBinError(k,sqrt(tmprnd));
0950     }
0951   }
0952 
0953   sprintf(tmpname,"hhcorrbg_pp_%d",i);
0954   hhcorrbg_pp[i] = (TH1D*)hhbottom_pp[i]->Clone(tmpname);
0955   hhcorrbg_pp[i]->Add(hhdy_pp[i]);
0956     hhcorrbg_pp[i]->SetMarkerColor(kBlack);
0957     hhcorrbg_pp[i]->SetLineColor(kBlack);
0958     hhbottom_pp[i]->SetLineColor(kBlue);
0959     hhdy_pp[i]->SetLineColor(kGreen+2);
0960   sprintf(tmpname,"hhall_pp_%d",i);
0961   hhall_pp[i] = (TH1D*)hhcorrbg_pp[i]->Clone(tmpname);
0962   hhall_pp[i]->Add(hhupspp[i]);
0963     hhall_pp[i]->SetLineColor(kMagenta);
0964     hhall_pp[i]->SetMarkerColor(kMagenta);
0965 
0966   hhcorrbg_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0967   hhcorrbg_pp[i]->GetFunction("expo")->SetLineColor(kBlack);
0968   hhbottom_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0969   hhbottom_pp[i]->GetFunction("expo")->SetLineColor(kBlue);
0970   hhdy_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0971   hhdy_pp[i]->GetFunction("expo")->SetLineColor(kGreen+2);
0972 
0973 } // end loop over pT bins
0974 
0975 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
0976   hhall_pp[nbins]->SetAxisRange(7.,12.);
0977   hhcorrbg_pp[nbins]->Draw("pehist");
0978   hhbottom_pp[nbins]->Draw("histsame");
0979   hhdy_pp[nbins]->Draw("histsame");
0980 
0981 
0982 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
0983   hhall_pp[nbins]->SetAxisRange(7.,12.);
0984   hhall_pp[nbins]->Draw("pehist");
0985   hhcorrbg_pp[nbins]->Draw("pesame");
0986   hhbottom_pp[nbins]->Draw("same");
0987   hhdy_pp[nbins]->Draw("same");
0988 
0989 TCanvas* cpp_vspt = new TCanvas("cpp_vspt","corr bg + sig vs pt in pp",50,50,1200,900);
0990 cpp_vspt->Divide(4,3);
0991 for(int i=0; i<nbins; i++) {
0992 if(i==0)  {cpp_vspt->cd(1);}
0993 if(i==1)  {cpp_vspt->cd(2);}
0994 if(i==2)  {cpp_vspt->cd(3);}
0995 if(i==3)  {cpp_vspt->cd(4);}
0996 if(i==4)  {cpp_vspt->cd(5);}
0997 if(i==5)  {cpp_vspt->cd(6);}
0998 if(i==6)  {cpp_vspt->cd(7);}
0999 if(i==7)  {cpp_vspt->cd(8);}
1000 if(i==8)  {cpp_vspt->cd(9);}
1001 if(i==9)  {cpp_vspt->cd(10);}
1002 if(i==10) {cpp_vspt->cd(11);}
1003 if(i==11) {cpp_vspt->cd(12);}
1004 hhall_pp[i]->SetAxisRange(7.0,12.0); hhall_pp[i]->Draw("hist"); hhcorrbg_pp[i]->Draw("histsame");
1005   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcorrbg_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
1006 }
1007 
1008 //-----------------------------------------------------------------------------------------
1009 // Combinatorial BG calculated for 10B central AuAu events
1010 //-----------------------------------------------------------------------------------------
1011 TCanvas* dummy = new TCanvas("dummy","dummy",0,0,500,500);
1012 
1013 f = new TFile("fakee_eideff09.root");
1014 for(int i=0; i<nbins+1; i++) {
1015   sprintf(tmpname,"hhfakefake_%d",i);
1016   hhfakefake[i] = (TH1D*)f->Get(tmpname);
1017   hhfakefake[i]->SetDirectory(gROOT);
1018 }
1019 f->Close();
1020 
1021 f = new TFile("crossterms_eideff09.root");
1022 for(int i=0; i<nbins+1; i++) {
1023   sprintf(tmpname,"hhfakehf_%d",i);
1024   hhfakehf[i] = (TH1D*)f->Get(tmpname);
1025   hhfakehf[i]->SetDirectory(gROOT);
1026 }
1027 f->Close();
1028 
1029 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
1030 fbg->SetParameters(10., -1.0, 4., -0.1);
1031 fbg->SetParLimits(1.,-999.,0.);
1032 fbg->SetParLimits(3.,-999.,0.);
1033 
1034 for(int i=0; i<nbins+1; i++) {
1035   sprintf(tmpname,"hhcombbg_%d",i);
1036   hhcombbg[i] = (TH1D*)hhfakefake[i]->Clone(tmpname);
1037   hhcombbg[i]->Add(hhfakehf[i]);
1038   sprintf(tmpname,"hhcombbg_scaled_%d",i);
1039   hhcombbg_scaled[i] = (TH1D*)hhcombbg[i]->Clone(tmpname);
1040   if(i==nbins) { fbg->SetParameters(10., -1.0, 4., -0.1); }
1041   hhcombbg[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
1042 
1043     for(int k=1; k<=hhcombbg[i]->GetNbinsX(); k++) {
1044       if(hhcombbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[i]->GetBinLowEdge(k)+hhcombbg[i]->GetBinWidth(k))>statscale_uplim) {
1045         hhcombbg_scaled[i]->SetBinContent(k,0.);
1046         hhcombbg_scaled[i]->SetBinError(k,0.);
1047       }
1048       else {
1049     // assumes comb bg scales with Ncoll^2
1050         double tmp = ncollfact * ncollfact * centwidthfact * statscale * hhcombbg[i]->GetFunction("fbg")->Eval(hhcombbg[i]->GetBinCenter(k));
1051         double tmprnd = myrandom->Poisson(tmp);
1052         if(tmprnd<0.) { tmprnd=0.; }
1053         hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
1054         hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
1055       }
1056     }
1057   hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
1058 }
1059 
1060 delete dummy;
1061 
1062 TCanvas* C1 = new TCanvas("C1","Combinatorial BG (ALL p_{T})",100,100,600,600);
1063 C1->SetLogy();
1064   hhfakefake[nbins]->SetAxisRange(7.0,14.0);
1065   hhfakefake[nbins]->SetMinimum(0.1);
1066   hhfakefake[nbins]->SetMaximum(5000.);
1067   hhfakefake[nbins]->SetLineColor(kGreen+2);
1068   hhfakefake[nbins]->SetLineWidth(2);
1069   hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1070   hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
1071   hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
1072   hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
1073   hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
1074   hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
1075   hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
1076   hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
1077   hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
1078   hhfakefake[nbins]->Draw("e");
1079 
1080   hhfakehf[nbins]->SetLineColor(kOrange+4);
1081   hhfakehf[nbins]->SetLineWidth(2);
1082   hhfakehf[nbins]->Draw("esame");
1083 
1084   hhcombbg[nbins]->SetLineColor(kBlack);
1085   hhcombbg[nbins]->SetLineWidth(2);
1086   hhcombbg[nbins]->Draw("esame");
1087 
1088 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG (ALL p_{T})",100,100,600,600);
1089 C1sc->SetLogy(); 
1090   hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
1091   hhcombbg_scaled[nbins]->Draw("esame");
1092 
1093 TCanvas* c00 = new TCanvas("c00","Combinatorial BG vs. p_{T}",150,150,1200,900);
1094 c00->Divide(4,3);
1095 
1096 for(int i=0; i<nbins; i++) {
1097 if(i==0)  {c00->cd(1);}
1098 if(i==1)  {c00->cd(2);}
1099 if(i==2)  {c00->cd(3);}
1100 if(i==3)  {c00->cd(4);}
1101 if(i==4)  {c00->cd(5);}
1102 if(i==5)  {c00->cd(6);}
1103 if(i==6)  {c00->cd(7);}
1104 if(i==7)  {c00->cd(8);}
1105 if(i==8)  {c00->cd(9);}
1106 if(i==9)  {c00->cd(10);}
1107 if(i==10) {c00->cd(11);}
1108 if(i>=11) {c00->cd(12);}
1109 hhcombbg[i]->SetAxisRange(7.0,14.0); hhcombbg[i]->Draw();
1110   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
1111 }
1112 
1113 TCanvas* c00scaled = new TCanvas("c00scaled","SCALED Combinatorial BG vs. p_{T}",150,150,1200,900);
1114 c00scaled->Divide(4,3);
1115 
1116 for(int i=0; i<nbins; i++) {
1117 if(i==0)  {c00scaled->cd(1);}
1118 if(i==1)  {c00scaled->cd(2);}
1119 if(i==2)  {c00scaled->cd(3);}
1120 if(i==3)  {c00scaled->cd(4);}
1121 if(i==4)  {c00scaled->cd(5);}
1122 if(i==5)  {c00scaled->cd(6);}
1123 if(i==6)  {c00scaled->cd(7);}
1124 if(i==7)  {c00scaled->cd(8);}
1125 if(i==8)  {c00scaled->cd(9);}
1126 if(i==9)  {c00scaled->cd(10);}
1127 if(i==10) {c00scaled->cd(11);}
1128 if(i>=11) {c00scaled->cd(12);}
1129 hhcombbg_scaled[i]->SetAxisRange(statscale_lowlim,statscale_uplim); hhcombbg_scaled[i]->Draw();
1130   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhcombbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
1131 }
1132 
1133 //---------------------------------------------
1134 //  Sibnal + BG
1135 //---------------------------------------------
1136 
1137 for(int i=0; i<nbins+1; i++) {
1138   sprintf(tmpname,"hhtotbg_scaled_%d",i);
1139   hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
1140   hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
1141 }
1142 
1143 for(int i=0; i<nbins+1; i++) {
1144 sprintf(tmpname,"hhall_scaled_%d",i);
1145 hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
1146 hhall_scaled[i]->Add(hhups[i]);
1147 }
1148 
1149 TCanvas* c000 = new TCanvas("c000","Signal + Background vs. p_{T}",200,200,1200,900);
1150 c000->Divide(4,3);
1151 for(int i=0; i<nbins; i++) {
1152 if(i==0)  {c000->cd(1);}
1153 if(i==1)  {c000->cd(2);}
1154 if(i==2)  {c000->cd(3);}
1155 if(i==3)  {c000->cd(4);}
1156 if(i==4)  {c000->cd(5);}
1157 if(i==5)  {c000->cd(6);}
1158 if(i==6)  {c000->cd(7);}
1159 if(i==7)  {c000->cd(8);}
1160 if(i==8)  {c000->cd(9);}
1161 if(i==9)  {c000->cd(10);}
1162 if(i==10) {c000->cd(11);}
1163 if(i==11) {c000->cd(12);}
1164 hhall_scaled[i]->SetAxisRange(7.0,14.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
1165   sprintf(tlchar,"%d-%d GeV",i,i+1);   tl[i] = new TLatex(8.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
1166 }
1167 
1168 // Fit Signal + Correlated BG
1169 
1170 // hhfit[] is signal + correlated bg
1171 for(int i=0; i<nbins+1; i++) {
1172   sprintf(tmpname,"hhfit_%d",i);
1173   hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
1174   for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
1175     hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
1176     hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
1177   }
1178 }
1179 
1180 //---------------------------------------------------------------------
1181 // plot signal + correlated bg for all pT
1182 //---------------------------------------------------------------------
1183 
1184   double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
1185   double tmppar1 = corrbgfitpar1;
1186 cout << "###### " << tmppar0 << " " << tmppar1 << endl;
1187 
1188 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
1189 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
1190 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
1191 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
1192 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
1193 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
1194 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
1195 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
1196 fSandBpp->SetParameter(7,tmppar0);
1197 fSandBpp->SetParameter(8,tmppar1);
1198 fSandBpp->SetLineColor(kBlue);
1199 fSandBpp->SetLineStyle(2);
1200 
1201 cout << "######### " << fTCBauau->GetParameter(0) << " " << fTCBauau->GetParameter(5) << " " << fTCBauau->GetParameter(6) << endl;
1202 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
1203 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
1204 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
1205 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
1206 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
1207 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
1208 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
1209 fSandBauau->SetParameter(7,tmppar0);
1210 fSandBauau->SetParameter(8,tmppar1);
1211 fSandBauau->SetLineColor(kRed);
1212 
1213 
1214 TCanvas* cfitall = new TCanvas("cfitall","FIT all pT",270,270,600,600);
1215   hhfit[nbins]->SetAxisRange(7.0,14.);
1216   hhfit[nbins]->GetXaxis()->CenterTitle();
1217   hhfit[nbins]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
1218   hhfit[nbins]->GetXaxis()->SetTitleOffset(1.1);
1219   hhfit[nbins]->GetXaxis()->SetLabelSize(0.045);
1220   hhfit[nbins]->GetYaxis()->CenterTitle();
1221   hhfit[nbins]->GetYaxis()->SetLabelSize(0.045);
1222   hhfit[nbins]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
1223   hhfit[nbins]->GetYaxis()->SetTitleOffset(1.5);
1224   hhfit[nbins]->Draw("pehist");
1225   fSandBpp->Draw("same");
1226   fSandBauau->Draw("same");
1227     //hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1228     //hhcorrbg_scaled[nbins]->Scale(0.82);
1229     //hhcorrbg_scaled[nbins]->Scale(0.70);
1230     //hhcorrbg_scaled[nbins]->Draw("histsame");
1231 
1232   TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
1233   fmycorrbg->SetParameters(tmppar0,tmppar1);
1234   fmycorrbg->SetLineStyle(2);
1235   fmycorrbg->SetLineColor(kRed);
1236   fmycorrbg->Draw("same");
1237 
1238 
1239 double myheight = fTCBauau->GetParameter(0);
1240 TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
1241 ld1->SetTextSize(0.035);
1242 ld1->Draw();
1243 TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
1244 ld2->SetTextSize(0.035);
1245 ld2->Draw();
1246 TLatex* ld3;
1247 
1248  char evstr[500];
1249  sprintf(evstr,"%i billion events",auauevts);
1250  ld3 = new TLatex(10.1,myheight-200.,evstr);
1251 ld3->SetTextSize(0.035);
1252 ld3->Draw();
1253 
1254 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
1255 TH1D* hhfit_tmp = (TH1D*)hhfit[nbins]->Clone("hhfit_tmp");
1256 hhfit_tmp->SetAxisRange(8.0,11.);
1257 hhfit_tmp->Draw("pehist");
1258 fSandBauau->Draw("same");
1259 fmycorrbg->Draw("same");
1260 
1261 
1262 //----------------------------------------------------------------------
1263 
1264 // plot signal + both backgrounds for all pT
1265 
1266 TCanvas* callpt = new TCanvas("callpt","Signal+BG (all p_{T})",300,300,600,600);
1267 
1268   hhall_scaled[nbins]->GetXaxis()->SetTitle("Invariant mass GeV/c");
1269   hhall_scaled[nbins]->SetLineColor(kBlack);
1270   hhall_scaled[nbins]->SetMarkerColor(kBlack);
1271   hhall_scaled[nbins]->SetMarkerStyle(20);
1272   hhall_scaled[nbins]->SetAxisRange(8.0,10.8);
1273   hhall_scaled[nbins]->Draw("pehist");
1274     hhcombbg_scaled[nbins]->SetLineColor(kBlue);
1275     hhcombbg_scaled[nbins]->Draw("histsame");
1276       hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1277       hhcorrbg_scaled[nbins]->Draw("histsame");
1278 
1279 //----------------------------------------------------------------
1280 //  Calculate RAA
1281 //----------------------------------------------------------------
1282 
1283 double u1start = 9.25;
1284 double u1stop  = 9.65;
1285 double u2start = 9.80;
1286 double u2stop  = 10.20;
1287 double u3start = 10.20;
1288 double u3stop  = 10.55;
1289 
1290  double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1], erraa3_up[nbins+1], erraa3_dn[nbins+1];
1291   for(int i=0; i<nbins; i++) { 
1292     //raa1[i] = supcor[0]; raa2[i] = supcor[1]; raa3[i] = supcor[2]; 
1293     raa1[i] = grRAA1S->Eval(0.5+i*1.0);
1294     raa2[i] = grRAA2S->Eval(0.5+i*1.0);
1295     raa3[i] = grRAA3S->Eval(0.5+i*1.0);
1296   }
1297   raa1[nbins] = raa1[0];
1298   raa2[nbins] = raa2[0];
1299   raa3[nbins] = raa3[0];
1300 
1301   int fbin1 = hhall_scaled[nbins]->FindBin(u1start + 0.001);
1302   int lbin1 = hhall_scaled[nbins]->FindBin(u1stop  - 0.001);
1303   int fbin2 = hhall_scaled[nbins]->FindBin(u2start + 0.001);
1304   int lbin2 = hhall_scaled[nbins]->FindBin(u2stop  - 0.001);
1305   int fbin3 = hhall_scaled[nbins]->FindBin(u3start + 0.001);
1306   int lbin3 = hhall_scaled[nbins]->FindBin(u3stop  - 0.001);
1307   cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
1308   cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
1309   cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
1310   cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
1311   cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
1312   cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
1313 
1314   double sum1[99]   = {0.};
1315   double truesum1[99]   = {0.};
1316   double ersum1[99] = {0.};
1317   double sumpp1[99]   = {0.};
1318   double ersumpp1[99] = {0.};
1319   double sum2[99]   = {0.};
1320   double truesum2[99]   = {0.};
1321   double ersum2[99] = {0.};
1322   double sumpp2[99]   = {0.};
1323   double ersumpp2[99] = {0.};
1324   double sum3[99]   = {0.};
1325   double truesum3[99]   = {0.};
1326   double allsum3[99]   = {0.};
1327   double ersum3[99] = {0.};
1328   double ersum3_up[99] = {0.};
1329   double ersum3_dn[99] = {0.};
1330   double sumpp3[99]   = {0.};
1331   double ersumpp3[99] = {0.};
1332 
1333   double sumsum1[99]   = {0.};
1334   double sumsum2[99]   = {0.};
1335   double sumsum3[99]   = {0.};
1336   double sumsum1pp[99]   = {0.};
1337   double sumsum2pp[99]   = {0.};
1338   double sumsum3pp[99]   = {0.};
1339 
1340   for(int i=0; i<nbins+1; i++) {
1341 
1342     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1343 
1344     for(int j=fbin1; j<=lbin1; j++) {
1345       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)));
1346       truesum1[i] += hhups1[i]->GetBinContent(j);
1347       ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1348       sumpp1[i]   += hhups1pp[i]->GetBinContent(j);
1349       ersumpp1[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1350     }
1351       sumsum1[i]     = truesum1[i];                   // direct count in mass range
1352       sumsum1pp[i]   = sumpp1[i];       
1353 
1354       //sumsum1[i]     = hhups1[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1355       //sumsum1pp[i]   = hhups1pp[i]->GetEntries();
1356       //sumsum1[i]     = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1357       //sumsum1pp[i]   = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1358 
1359         if(sumsum1[i]>0. && sumsum1pp[i]>0.) {
1360           erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1[i]/sumsum1pp[i]/sumsum1pp[i]);
1361       cout << "i " << i << " raa1 " << raa1[i] << " erraa1  " << erraa1[i] << endl;
1362         } else {raa1[i]=-1.0; erraa1[i] = 999.; }
1363 
1364     for(int j=fbin2; j<=lbin2; j++) {
1365       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)));
1366       truesum2[i] += hhups2[i]->GetBinContent(j);
1367       ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1368       sumpp2[i]   += hhups2pp[i]->GetBinContent(j);
1369       ersumpp2[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1370     }
1371       sumsum2[i]     = truesum2[i];                   // direct count in mass range
1372       sumsum2pp[i]   = sumpp2[i];       
1373 
1374       //sumsum2[i]     = hhups2[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1375       //sumsum2pp[i]   = hhups2pp[i]->GetEntries();
1376       //sumsum2[i]     = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1377       //sumsum2pp[i]   = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1378       
1379         if(sumsum2[i]>0. && sumsum2pp[i]>0.) {
1380           erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2[i]/sumsum2pp[i]/sumsum2pp[i]);
1381       cout << "i " << i << " raa2 " << raa2[i] << " erraa2  " << erraa2[i] << endl;
1382         } else {raa2[i]=-1.0; erraa2[i] = 999.; }
1383 
1384     for(int j=fbin3; j<=lbin3; j++) {
1385       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)));
1386       truesum3[i] += hhups3[i]->GetBinContent(j);
1387       ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1388       allsum3[i] +=  hhall_scaled[i]->GetBinContent(j);
1389       sumpp3[i]   += hhups3pp[i]->GetBinContent(j);
1390       ersumpp3[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1391     }
1392       sumsum3[i]     = truesum3[i];                   // direct count in mass range
1393       sumsum3pp[i]   = sumpp3[i];       
1394       // get Poisson errors on truesum3
1395       //double err3_up, err3_dn;
1396       //error_fg_bg_pair((int) allsum3[i], (int) (allsum3[i]-truesum3[i]), err3_up, err3_dn);
1397       //ersum3_up[i] = err3_up * err3_up;
1398       //ersum3_dn[i] = err3_dn * err3_dn;
1399       //sumsum3[i]     = hhups3[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1400       //sumsum3pp[i]   = hhups3pp[i]->GetEntries();
1401       //sumsum3[i]     = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;   // total number of upsilons in pT bin
1402       //sumsum3pp[i]   = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1403       
1404         if(truesum3[i]>0. && sumpp3[i]>0.) {
1405       erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
1406           //erraa3_up[i] = raa3[i]*sqrt(ersum3_up[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
1407           //erraa3_dn[i] = raa3[i]*sqrt(ersum3_dn[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
1408       cout << "i " << i << " raa3 " << raa3[i] << " erraa3_up  " << erraa3_up[i] << " erraa3_dn " << erraa3_dn[i] << endl;
1409         } else {raa3[i]=-1.0; erraa3[i] = 999;}
1410   
1411   }
1412 
1413 cout << "====== Y(1S):" << endl;
1414   for(int i=0; i<nbins+1; i++) {
1415     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1416     cout << "   " << i << " " << truesum1[i] << "(" << Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1417          << sqrt(ersum1[i]) << " \t\t pp: " << sumpp1[i] << " +- " << sqrt(ersumpp1[i]) << " raa " << raa1[i] << " erraa " << erraa1[i] << endl;
1418   }
1419 cout << "====== Y(2S):" << endl;
1420   for(int i=0; i<nbins+1; i++) {
1421     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1422     cout << "   " << i << " " << truesum2[i] << "(" << Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- " 
1423          << sqrt(ersum2[i]) << " \t\t pp: " << sumpp2[i] << " +- " << sqrt(ersumpp2[i]) << " raa " << raa2[i] << " erraa " << erraa2[i] << endl;
1424   }
1425 cout << "====== Y(3S):" << endl;
1426   for(int i=0; i<nbins+1; i++) {
1427 
1428     double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1429     cout << "   " << i << " " << truesum3[i] << "(" << Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " (all " << allsum3[i] << ") "<< " +- " 
1430       //     << sqrt(ersum3_up[i]) << " - " << sqrt(ersum3_dn[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << " raa " << raa3[i] << " erraa " << erraa3[i] << endl;
1431       << sqrt(ersum3[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << " raa " << raa3[i] << " erraa " << erraa3[i] << endl;
1432   }
1433 
1434 //-------------------------------------------------
1435 //  Plot RAA vs pT
1436 //-------------------------------------------------
1437 
1438 int npts1 = 9;
1439  int npts2 = 7;
1440 int npts3 = 7;
1441 
1442 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,600,600);
1443 TH2F* hh2 = new TH2F("hh2"," ",10,0.,float(npts1+1),10,0.,1.5);
1444 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1445 hh2->GetXaxis()->SetTitleOffset(1.0);
1446 hh2->GetXaxis()->SetTitleColor(1);
1447 hh2->GetXaxis()->SetTitleSize(0.040);
1448 hh2->GetXaxis()->SetLabelSize(0.040);
1449 hh2->GetYaxis()->SetTitle("R_{AA}");
1450 hh2->GetYaxis()->SetTitleOffset(1.3);
1451 hh2->GetYaxis()->SetTitleSize(0.040);
1452 hh2->GetYaxis()->SetLabelSize(0.040);
1453 hh2->Draw();
1454 
1455 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.5 + double(i);}
1456 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.43 + double(i);}
1457 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.57 + double(i);}
1458 
1459 TGraphErrors* gr1 = new TGraphErrors(npts1,xx1,raa1,0,erraa1);
1460 gr1->SetMarkerStyle(20);
1461 gr1->SetMarkerColor(kBlack);
1462 gr1->SetLineColor(kBlack);
1463 gr1->SetLineWidth(2);
1464 gr1->SetMarkerSize(1.5);
1465 gr1->SetName("gr1");
1466 gr1->Draw("p");
1467 
1468 erraa2[7] = erraa2[7]*0.85;
1469 
1470 TGraphErrors* gr2 = new TGraphErrors(npts2,xx2,raa2,0,erraa2);
1471 gr2->SetMarkerStyle(20);
1472 gr2->SetMarkerColor(kRed);
1473 gr2->SetLineColor(kRed);
1474 gr2->SetLineWidth(2);
1475 gr2->SetMarkerSize(1.5);
1476 gr2->SetName("gr2");
1477 gr2->Draw("p");
1478 
1479 /*
1480 //--- 3S state -------------------
1481  double dummy[9];
1482  for(int i=0; i<9; i++) {dummy[i]=-99.;}
1483  TGraphErrors* gr3 = new TGraphErrors(6,xx3,dummy,0,erraa3);
1484  gr3->SetMarkerStyle(20);
1485  gr3->SetMarkerColor(kBlue);
1486  gr3->SetLineColor(kBlue);
1487  gr3->SetLineWidth(2);
1488  gr3->SetMarkerSize(1.5);
1489  gr3->SetName("gr3");
1490 // gr3->Draw("p");
1491 
1492 TArrow* aa[9];
1493 TLine* ll[9];
1494 //erraa3[3] = erraa3[3]*1.5;
1495 erraa3[4] = erraa3[4]*0.6;
1496 //erraa3[4] = erraa3[4]*1.5;
1497 erraa3[6] = erraa3[6]*1.2;
1498 
1499 for(int i=0; i<npts3; i++) {
1500   aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1501   aa[i]->SetLineColor(kBlue);
1502   aa[i]->SetLineWidth(2);
1503   aa[i]->Draw();
1504   ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1505   ll[i]->SetLineColor(kBlue);
1506   ll[i]->SetLineWidth(2);
1507   ll[i]->Draw();
1508 }
1509 //--- end 3S state --------------
1510 */
1511 
1512 TLegend *leg = new TLegend(0.14,0.75,0.34,0.86);
1513 leg->SetFillColor(10);
1514 leg->SetFillStyle(1001);
1515 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1516 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1517 //TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1518 leg->Draw();
1519 
1520 TLatex* l1 = new TLatex(3.35,1.33,"sPHENIX Simulation");
1521  l1->SetTextSize(0.042);
1522  l1->Draw();
1523 
1524 TLatex* lum = new TLatex(3.35,1.23,"0-10% Au+Au #sqrt{s} = 200 GeV");
1525  lum->SetTextSize(0.042);
1526  lum->Draw();
1527 
1528  TLatex* run = new TLatex(3.35,1.13,run_plan[scenario].c_str());
1529  run->SetTextSize(0.042);
1530  run->Draw();
1531 
1532 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1533 lll->SetLineColor(kBlue);
1534 lll->SetLineWidth(2);
1535 //lll->Draw();
1536 
1537 //==================================================================================
1538 
1539 /*
1540 TFile *fout = new TFile("test.root","recreate");
1541 for(int i=0; i<nbins+1; i++) {
1542   sprintf(tmpname,"hhfit_%d",i);
1543   hhfit[i]->Write();
1544   hhups1pp[i]->Write();
1545   hhups2pp[i]->Write();
1546   hhups3pp[i]->Write();
1547   hhupspp[i]->Write();
1548   hhups1[i]->Write();
1549   hhups2[i]->Write();
1550   hhups3[i]->Write();
1551   hhups[i]->Write();
1552   hhall_pp[i]->Write();
1553   hhcorrbg_pp[i]->Write();
1554   hhcorrbg_scaled[i]->Write();
1555   hhcombbg_scaled[i]->Write();
1556 }
1557 fout->Write();
1558 fout->Close();
1559 */
1560 
1561 }
1562
[ 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!