Back to home page

sPhenix code displayed by LXR

 
 

    


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

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_dAu_vscent_2022() {
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 = 4;
0110 double eideff = 0.9;
0111 string times = "*";
0112 TLatex* tl[15];
0113 char tlchar[999];
0114 
0115 double statscale;
0116 statscale = 190./4./10.; // from 10B 0-10% AuAu to 190B min. bias dAu divided into 4 bins in centrality
0117 double statscale_lowlim = 7.0;
0118 double statscale_uplim = 14.0;
0119 
0120   TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
0121   TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
0122   TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
0123   TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
0124   TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
0125   TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
0126   TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
0127   TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
0128   TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
0129   TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
0130   TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
0131   TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
0132 
0133 //---------------------------------------------------------
0134 // Upsilons
0135 //---------------------------------------------------------
0136 
0137  // Input assumptions about level of suppression are taken from the theory paper: 
0138  // Michael Strickland, Dennis Bazow, Nucl.Phys. A879 (2012) 25-58
0139  //==============================================
0140  // Read in Strickland curves for different eta/s values
0141  // These are exclusive - i.e. they are the 
0142  // suppression of that state with no feed-down
0143  // effects included 
0144  //==============================================
0145  // The curves vs eta/s read in below were provided by Strickland privately, since
0146  // the paper contains only the eta/s = 1 cases for RHIC energy
0147  // Get the exclusive raa values for the Y1S, Y2S, Y3S, chib1, chib2
0148  // Checked that this is being read in correctly
0149                         
0150   double str_npart[101];
0151   double str_raa[5][3][101];
0152   for(int istate=0;istate<5;istate++)
0153     {
0154       for(int ietas=0;ietas<3;ietas++)
0155         {
0156           if(istate < 3)
0157             {
0158               char fname[500];
0159               sprintf(fname,"./strickland_calculations/Y%is-potb-eta%i-npart.dat",istate+1,ietas+1);
0160               ifstream fin(fname);
0161 
0162               for(int inpart=0;inpart<101;inpart++)
0163                 {
0164                   fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
0165                 }
0166               fin.close();
0167             }
0168 
0169           if(istate > 2)
0170             {
0171               char fname[500];
0172               sprintf(fname,"./strickland_calculations/chib%i-potb-eta%i-npart.dat",istate-2,ietas+1);
0173               ifstream fin(fname);
0174 
0175               for(int inpart=0;inpart<101;inpart++)
0176                 {
0177                   fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
0178                 }
0179               fin.close();
0180             }
0181         }
0182 
0183     }
0184 
0185   // Now construct the inclusive RAA values from Strickland's exclusive modifications
0186   // ff's from table II of the paper for the 1S state are as follows
0187   // They add up to 1.0
0188   double ff1S[5] = {0.51, 0.107, 0.008, 0.27, 0.105};  // Y1S, Y2S, Y3S, chib1, chib2 respectively
0189   // These are from arXiv:1210.7512
0190   double ff2S[2] = {0.5, 0.5};                         // Y2S and Y3S respectively
0191   double str_raa_inclusive[3][3][101];
0192 
0193   // checked the Y1S inclusive result against figure 10 of arXiv:1112.2761
0194   // There is no plot available for the 2S and 3S 
0195   for(int ietas=0;ietas<3;ietas++)
0196     for(int inpart=0;inpart<101;inpart++)
0197       {
0198         str_raa_inclusive[0][ietas][inpart] =
0199           str_raa[0][ietas][inpart] * ff1S[0]
0200           + str_raa[1][ietas][inpart] * ff1S[1]
0201           + str_raa[2][ietas][inpart] * ff1S[2]
0202           + str_raa[3][ietas][inpart] * ff1S[3]
0203           + str_raa[4][ietas][inpart] * ff1S[4];
0204 
0205         str_raa_inclusive[1][ietas][inpart] =
0206           str_raa[1][ietas][inpart] * ff2S[0]
0207           + str_raa[2][ietas][inpart] * ff2S[1];
0208 
0209         str_raa_inclusive[2][ietas][inpart] = str_raa[2][ietas][inpart];
0210       }
0211 
0212   double strick_raa1_eta1[101],strick_raa1_eta2[101],strick_raa1_eta3[101];
0213   double strick_raa2_eta1[101],strick_raa2_eta2[101],strick_raa2_eta3[101];
0214   double strick_raa3_eta1[101],strick_raa3_eta2[101],strick_raa3_eta3[101];
0215     for(int ipart=0;ipart<101;ipart++) {
0216       strick_raa1_eta1[ipart] = str_raa_inclusive[0][0][ipart];
0217       strick_raa1_eta2[ipart] = str_raa_inclusive[0][1][ipart];
0218       strick_raa1_eta3[ipart] = str_raa_inclusive[0][2][ipart];
0219       strick_raa2_eta1[ipart] = str_raa_inclusive[1][0][ipart];
0220       strick_raa2_eta2[ipart] = str_raa_inclusive[1][1][ipart];
0221       strick_raa2_eta3[ipart] = str_raa_inclusive[1][2][ipart];
0222       strick_raa3_eta1[ipart] = str_raa_inclusive[2][0][ipart];
0223       strick_raa3_eta2[ipart] = str_raa_inclusive[2][1][ipart];
0224       strick_raa3_eta3[ipart] = str_raa_inclusive[2][2][ipart];
0225     }
0226 
0227 //======================================================================================
0228 //======================================================================================
0229 //======================================================================================
0230 
0231 TGraph* grRAA1S = new TGraph(101,str_npart,strick_raa1_eta2);
0232 TGraph* grRAA2S = new TGraph(101,str_npart,strick_raa2_eta2);
0233 TGraph* grRAA3S = new TGraph(101,str_npart,strick_raa3_eta2);
0234 TGraph* grRAA1S_eta1 = new TGraph(101,str_npart,strick_raa1_eta1);
0235 TGraph* grRAA2S_eta1 = new TGraph(101,str_npart,strick_raa2_eta1);
0236 TGraph* grRAA3S_eta1 = new TGraph(101,str_npart,strick_raa3_eta1);
0237 TGraph* grRAA1S_eta3 = new TGraph(101,str_npart,strick_raa1_eta3);
0238 TGraph* grRAA2S_eta3 = new TGraph(101,str_npart,strick_raa2_eta3);
0239 TGraph* grRAA3S_eta3 = new TGraph(101,str_npart,strick_raa3_eta3);
0240 grRAA1S->SetLineColor(kBlack);
0241 grRAA1S->SetLineStyle(7);
0242 grRAA2S->SetLineColor(kRed);
0243 grRAA2S->SetLineStyle(7);
0244 grRAA3S->SetLineColor(kBlue);
0245 grRAA3S->SetLineStyle(7);
0246 grRAA1S_eta1->SetLineColor(kBlack);
0247 grRAA1S_eta1->SetLineStyle(1);
0248 grRAA2S_eta1->SetLineColor(kRed);
0249 grRAA2S_eta1->SetLineStyle(1);
0250 grRAA3S_eta1->SetLineColor(kBlue);
0251 grRAA3S_eta1->SetLineStyle(1);
0252 grRAA1S_eta3->SetLineColor(kBlack);
0253 grRAA1S_eta3->SetLineStyle(8);
0254 grRAA2S_eta3->SetLineColor(kRed);
0255 grRAA2S_eta3->SetLineStyle(8);
0256 grRAA3S_eta3->SetLineColor(kBlue);
0257 grRAA3S_eta3->SetLineStyle(8);
0258 
0259 int nchan=400;
0260 double start=0.0;
0261 double stop=20.0;
0262 
0263 string str_UpsilonPt  = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))";   // dN/dpT
0264 string str_UpsilonXPt = "(2.0*3.14159*x*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))"; // need this for mean pT calculation
0265 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
0266 TF1* fUpsilonXPt = new TF1("fUpsilonXPt",str_UpsilonXPt.c_str(),0.,20.);
0267 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
0268 fUpsilonXPt->SetParameters(72.1, 26.516, 10.6834);
0269 double upsnorm = fUpsilonPt->Integral(0.,20.);
0270 
0271 double frac[3];  // upsilons fractions
0272   frac[0] = 0.7117;
0273   frac[1] = 0.1851;
0274   frac[2] = 0.1032;
0275 double scale[3]; // mass scaling
0276   scale[0] = 1.0;
0277   scale[1] = 1.0595;
0278   scale[2] = 1.0946;
0279 
0280 double Npart[nbins+1], NpartAvg=0.;
0281 Npart[0] = 8.2;  // same as Ncoll for now
0282 Npart[1] = 6.1;
0283 Npart[2] = 4.4;
0284 Npart[3] = 2.6;
0285 for(int i=0; i<4; i++) {NpartAvg += Npart[i];} NpartAvg = NpartAvg/4.;
0286 cout << "Raa for dAu = " << grRAA1S->Eval(NpartAvg) << " " << grRAA2S->Eval(NpartAvg) << " " << grRAA3S->Eval(NpartAvg) << endl;
0287 
0288 double NcollAuAu = 955.;
0289 
0290 double Ncoll[nbins+1];
0291 Ncoll[0] = 8.2; // 0-20% dAu
0292 Ncoll[1] = 6.1; // 20-40%
0293 Ncoll[2] = 4.4; // 40-60%
0294 Ncoll[3] = 2.6; // 60-84%
0295 double NcollAvg=0.; for(int i=0; i<4; i++) {NcollAvg += Ncoll[i];} NcollAvg = NcollAvg/4.;
0296 
0297 double Npionpairs[nbins+1];
0298 Npionpairs[0] = 2.23e-03;
0299 Npionpairs[1] = 1.23e-03;
0300 Npionpairs[2] = 5.31e-04;
0301 Npionpairs[3] = 1.85e-04;  // This includes the fact that this bin is 1.2 times wider
0302 double NpionpairsAvg=0.; for(int i=0; i<4; i++) {NpionpairsAvg += Npionpairs[i];} NpionpairsAvg = NpionpairsAvg/4.;
0303 
0304 //---------- Use fixed numbers for 2022 estimate --------------------------------------
0305 // Numbers from Marzia for 190B min. bias dAu and 2400B p+p events
0306 
0307   float Nups1[nbins],Nups2[nbins],Nups3[nbins];
0308   float Nups1tot=0, Nups2tot=0, Nups3tot=0;
0309 
0310   Nups1[0]   = 466;
0311   Nups1[1]   = 346;
0312   Nups1[2]   = 250;
0313   Nups1[3]   = 177;
0314   for(int i=0; i<4; i++) {Nups1tot += Nups1[i];}
0315 
0316   Nups2[0]   = 83.8;
0317   Nups2[1]   = 62.4;
0318   Nups2[2]   = 45.0;
0319   Nups2[3]   = 31.9;
0320   for(int i=0; i<4; i++) {Nups2tot += Nups2[i];}
0321 
0322 //  Nups3[0]   = 8.4;
0323 //  Nups3[1]   = 6.2;
0324 //  Nups3[2]   = 4.5;
0325 //  Nups3[3]   = 3.2;
0326   Nups3[0]   = Nups2[0]*0.5575;
0327   Nups3[1]   = Nups2[1]*0.5575;
0328   Nups3[2]   = Nups2[2]*0.5575;
0329   Nups3[3]   = Nups2[3]*0.5575;
0330   for(int i=0; i<4; i++) {Nups3tot += Nups3[i];}
0331   cout << "Number of Upsilons in dAu = " << Nups1tot << " " << Nups2tot << " " << Nups3tot << endl;
0332 
0333 
0334   int Nups1pp = 2.86e+03;
0335   int Nups2pp = 7.16e+02;
0336   int Nups3pp = 3.98e+02;
0337 
0338 //====================================================================================
0339 //====================================================================================
0340 //====================================================================================
0341 
0342 //
0343 // these histograms (hhups*) are randomly generated using the fit above to single 
0344 // peak (fCB function) and used for the RAA uncertainty calculation
0345 //
0346 // FORCE specific RESOLUTION and tail shape
0347   double tonypar1 = 0.98;  // Tony's 100 pion simulation April 2019
0348   double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
0349   //double tonypar3 = 9.437;  // Tony's 100 pion simulation April 2019
0350   double tonypar3 = 9.448;  // benchmark
0351   double tonypar4 = 0.100;  // benchmark
0352   double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
0353   fCBpp->SetParameter(0,1000.); 
0354   fCBpp->SetParameter(1,tonypar1); 
0355   fCBpp->SetParameter(2,tonypar2); 
0356   fCBpp->SetParameter(3,tonypar3); 
0357   fCBpp->SetParameter(4,tonypar4pp);  
0358   fCBauau->SetParameter(0,1000.); 
0359   fCBauau->SetParameter(1,tonypar1); 
0360   fCBauau->SetParameter(2,tonypar2); 
0361   fCBauau->SetParameter(3,tonypar3); 
0362   fCBauau->SetParameter(4,tonypar4);  
0363 
0364 char hhname[999];
0365 TH1D* hhups[nbins+1];
0366 TH1D* hhups1[nbins+1];
0367 TH1D* hhups2[nbins+1];
0368 TH1D* hhups3[nbins+1];
0369 TH1D* hhupspp;
0370 TH1D* hhups1pp;
0371 TH1D* hhups2pp;
0372 TH1D* hhups3pp;
0373 for(int i=0; i<nbins+1; i++) {
0374   sprintf(hhname,"hhups_%d",i);
0375   hhups[i] = new TH1D(hhname,"",nchan,start,stop);
0376   hhups[i]->Sumw2();
0377   sprintf(hhname,"hhups1_%d",i);
0378   hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
0379   hhups1[i]->Sumw2();
0380   sprintf(hhname,"hhups2_%d",i);
0381   hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
0382   hhups2[i]->Sumw2();
0383   sprintf(hhname,"hhups3_%d",i);
0384   hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
0385   hhups3[i]->Sumw2();
0386   hhups[i]->SetLineWidth(2);
0387   hhups1[i]->SetLineWidth(2);
0388   hhups2[i]->SetLineWidth(2);
0389   hhups3[i]->SetLineWidth(2);
0390 }
0391   sprintf(hhname,"hhupspp");
0392   hhupspp= new TH1D(hhname,"",nchan,start,stop);
0393   hhupspp->Sumw2();
0394   sprintf(hhname,"hhups1pp");
0395   hhups1pp = new TH1D(hhname,"",nchan,start,stop);
0396   hhups1pp->Sumw2();
0397   sprintf(hhname,"hhups2pp");
0398   hhups2pp = new TH1D(hhname,"",nchan,start,stop);
0399   hhups2pp->Sumw2();
0400   sprintf(hhname,"hhups3pp");
0401   hhups3pp = new TH1D(hhname,"",nchan,start,stop);
0402   hhups3pp->Sumw2();
0403   hhupspp->SetLineWidth(2);
0404   hhups1pp->SetLineWidth(2);
0405   hhups2pp->SetLineWidth(2);
0406   hhups3pp->SetLineWidth(2);
0407 
0408 for(int j=0; j<nbins; j++) {
0409   double s1 = j*1.0;
0410   double s2 = s1 + 1.0;
0411   fCBauau->SetParameter(3,tonypar3);
0412     for(int i=0; i<int(Nups1[j]+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0413   fCBauau->SetParameter(3,tonypar3*scale[1]);
0414     for(int i=0; i<int(Nups2[j]+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0415   fCBauau->SetParameter(3,tonypar3*scale[2]);
0416     for(int i=0; i<int(Nups3[j]+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
0417 } // end loop over centrality bins
0418 
0419   fCBpp->SetParameter(3,tonypar3);
0420     for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp->Fill(myrnd); hhupspp->Fill(myrnd); }
0421   fCBpp->SetParameter(3,tonypar3*scale[1]);
0422     for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp->Fill(myrnd); hhupspp->Fill(myrnd); }
0423   fCBpp->SetParameter(3,tonypar3*scale[2]);
0424     for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp->Fill(myrnd); hhupspp->Fill(myrnd); }
0425 
0426 
0427 // fit and draw pp no bg
0428 
0429 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
0430   fTCBpp->SetParameter(0,2000.);
0431   fTCBpp->FixParameter(1,tonypar1);
0432   fTCBpp->FixParameter(2,tonypar2);
0433   fTCBpp->SetParameter(3,tonypar3);
0434   fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
0435   fTCBpp->SetParameter(5,500.);
0436   fTCBpp->SetParameter(6,100.);
0437   hhupspp->Fit(fTCBpp,"rl","",7.,11.);
0438   hhupspp->SetAxisRange(7.,11.);
0439   hhupspp->SetMarkerSize(1.0);
0440   hhupspp->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0441   hhupspp->GetXaxis()->SetTitleOffset(1.0);
0442   double tmpamp1 = hhupspp->GetFunction("fTCBpp")->GetParameter(0);
0443   double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0444   double tmpamp6 = tmpamp1*frac[2]/frac[0];
0445   hhupspp->Draw();
0446 
0447   fCB1s->SetLineColor(kBlue);
0448   fCB1s->SetLineWidth(1);
0449     fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
0450     fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
0451     fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
0452     fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
0453     fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
0454   fCB2s->SetLineColor(kRed);
0455   fCB2s->SetLineWidth(1);
0456     fCB2s->SetParameter(0,tmpamp5);
0457     fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
0458     fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
0459     fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
0460     fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
0461   fCB3s->SetLineColor(kGreen+2);
0462   fCB3s->SetLineWidth(1);
0463     fCB3s->SetParameter(0,tmpamp6);
0464     fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
0465     fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
0466     fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
0467     fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
0468   fCB1s->Draw("same");
0469   fCB2s->Draw("same");
0470   fCB3s->Draw("same");
0471 
0472 // fit and draw AuAu no bg
0473 
0474 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Central Au+Au",100,100,600,600);
0475   fTCBauau->SetParameter(0,2000.);
0476   fTCBauau->FixParameter(1,tonypar1);
0477   fTCBauau->FixParameter(2,tonypar2);
0478   fTCBauau->SetParameter(3,tonypar3);
0479   fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
0480   fTCBauau->SetParameter(5,500.);
0481   fTCBauau->SetParameter(6,100.);
0482   hhups[0]->Fit(fTCBauau,"rl","",7.,11.);
0483   hhups[0]->SetAxisRange(8.5,11.);
0484   hhups[0]->SetMarkerSize(1.0);
0485   hhups[0]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
0486   hhups[0]->GetXaxis()->SetTitleOffset(1.0);
0487 //  tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
0488 //  tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
0489 //  tmpamp6 = tmpamp1*frac[2]/frac[0];
0490   hhups[0]->Draw();  // 0-10% central
0491 
0492 TCanvas* cdummy1 = new TCanvas("cdummy1","cdummy1",0,0,500,500);
0493 //----------------------------------------------------
0494 //  Backgrounds
0495 //----------------------------------------------------
0496 
0497 TH1D* hhall[nbins+1];
0498 TH1D* hhall_scaled[nbins+1];
0499 
0500 TH1D* hhtotbg[nbins+1]; 
0501 TH1D* hhtotbg_scaled[nbins+1]; 
0502 TH1D* hhcombbg[nbins+1];
0503 TH1D* hhcombbg_scaled[nbins+1];
0504 TH1D* hhfakefake[nbins+1];
0505 TH1D* hhfakehf[nbins+1];
0506 TH1D* hhbottom[nbins+1];
0507 TH1D* hhcharm[nbins+1];
0508 TH1D* hhdy[nbins+1];
0509 TH1D* hhcorrbg[nbins+1];
0510 TH1D* hhcorrbg_scaled[nbins+1];
0511 TH1D* hhfit[nbins+1];
0512 char tmpname[999];
0513 
0514 //----------------------------------------------------------------------------------------
0515 // Correlated BG calculated for 10B central AuAu events
0516 //----------------------------------------------------------------------------------------
0517 
0518 double corrbgfitpar0;
0519 double corrbgfitpar1;
0520 
0521 TFile* f=new TFile("ccbb_eideff09.root");
0522 
0523   sprintf(tmpname,"hhbottom_15");  // 15 is integrated over pT
0524   hhbottom[nbins] = (TH1D*)f->Get(tmpname);
0525   hhbottom[nbins]->SetDirectory(gROOT);
0526   sprintf(tmpname,"hhcharm_15");
0527   hhcharm[nbins] = (TH1D*)f->Get(tmpname);
0528   hhcharm[nbins]->SetDirectory(gROOT);
0529   sprintf(tmpname,"hhdy_15");
0530   hhdy[nbins] = (TH1D*)f->Get(tmpname);
0531   hhdy[nbins]->SetDirectory(gROOT);
0532   sprintf(tmpname,"hhcorrbg_15");
0533   hhcorrbg[nbins] = (TH1D*)hhbottom[nbins]->Clone(tmpname);
0534   hhcorrbg[nbins]->Add(hhcharm[nbins]);
0535   hhcorrbg[nbins]->Add(hhdy[nbins]);
0536   sprintf(tmpname,"hhcorrbg_scaled_15");
0537   hhcorrbg_scaled[nbins] = (TH1D*)hhcorrbg[nbins]->Clone(tmpname);
0538   hhcorrbg[nbins]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0539   hhbottom[nbins]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0540   hhdy[nbins]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0541     corrbgfitpar0 = hhcorrbg[nbins]->GetFunction("expo")->GetParameter(0);
0542     corrbgfitpar1 = hhcorrbg[nbins]->GetFunction("expo")->GetParameter(1);
0543     cout << "bgpar0["<< nbins <<"]="<<hhcorrbg[nbins]->GetFunction("expo")->GetParameter(0)+TMath::Log(statscale)<<";"<< endl; 
0544     cout << "bgpar1["<< nbins <<"]="<<hhcorrbg[nbins]->GetFunction("expo")->GetParameter(1)<<";"<< endl; 
0545     for(int k=1; k<=hhcorrbg[nbins]->GetNbinsX(); k++) {
0546       if(hhcorrbg[nbins]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[nbins]->GetBinLowEdge(k)+hhcorrbg[nbins]->GetBinWidth(k))>statscale_uplim) {
0547         hhcorrbg_scaled[nbins]->SetBinContent(k,0.); 
0548         hhcorrbg_scaled[nbins]->SetBinError(k,0.);
0549       }
0550       else {
0551         double tmp = statscale * hhcorrbg[nbins]->GetFunction("expo")->Eval(hhcorrbg[nbins]->GetBinCenter(k));
0552         double tmprnd = myrandom->Poisson(tmp); 
0553         if(tmprnd<0.) { tmprnd=0.; }
0554         hhcorrbg_scaled[nbins]->SetBinContent(k,tmprnd); 
0555         hhcorrbg_scaled[nbins]->SetBinError(k,sqrt(tmprnd));
0556       }
0557     }
0558   hhcorrbg_scaled[nbins]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0559   hhcorrbg[nbins]->SetDirectory(gROOT);
0560   hhcorrbg_scaled[nbins]->SetDirectory(gROOT);
0561 
0562 f->Close();
0563 
0564 // Correlated background in AuAu vs centrality
0565 for(int i=0; i<nbins; i++) {
0566 
0567   sprintf(tmpname,"hhcorrbg_%d",i);
0568   hhcorrbg_scaled[i] = (TH1D*)hhcorrbg_scaled[nbins]->Clone(tmpname);
0569 
0570     for(int k=1; k<=hhcorrbg_scaled[nbins]->GetNbinsX(); k++) {
0571       if(hhcorrbg_scaled[nbins]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg_scaled[nbins]->GetBinLowEdge(k)+hhcorrbg_scaled[nbins]->GetBinWidth(k))>statscale_uplim) {
0572         hhcorrbg_scaled[i]->SetBinContent(k,0.);
0573         hhcorrbg_scaled[i]->SetBinError(k,0.);
0574       }
0575       else {
0576         double tmp = (Ncoll[i]/NcollAuAu) * hhcorrbg_scaled[nbins]->GetFunction("expo")->Eval(hhcorrbg_scaled[nbins]->GetBinCenter(k));
0577         double tmprnd = myrandom->Poisson(tmp);
0578         if(tmprnd<0.) { tmprnd=0.; }
0579         hhcorrbg_scaled[i]->SetBinContent(k,tmprnd);
0580         hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0581       }
0582     }
0583   hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0584 
0585 } // end loop over centrality bins
0586 
0587 
0588 delete cdummy1;
0589 
0590 TCanvas* c111 = new TCanvas("c111","Au+Au Correlated Background vs. Centrality",200,200,1200,600);
0591 c111->Divide(4,2);
0592 for(int i=0; i<nbins; i++) {
0593   c111->cd(i+1);
0594   hhcorrbg_scaled[i]->SetAxisRange(8.5,11.0); hhcorrbg_scaled[i]->SetMarkerStyle(1); hhcorrbg_scaled[i]->Draw("pe");
0595   sprintf(tlchar,"%d-%d",10*i,10*(i+1));   tl[i] = new TLatex(9.0,hhcorrbg_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
0596 }
0597 
0598 
0599 //-----------------------------------------------------------------------------------------
0600 //  Correlated bg in p+p
0601 //-----------------------------------------------------------------------------------------
0602 
0603 TH1D* hhbottom_pp;
0604 TH1D* hhdy_pp;
0605 TH1D* hhcorrbg_pp;
0606 TH1D* hhall_pp;
0607 
0608 double ppcorr = (2400./14.)/962.; // 2400B pp and 140B MB AuAu
0609 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
0610 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
0611 
0612   sprintf(tmpname,"hhbottom_pp");
0613   hhbottom_pp = (TH1D*)hhbottom[nbins]->Clone(tmpname);
0614   for(int k=1; k<=hhbottom_pp->GetNbinsX(); k++) {
0615     if(hhbottom_pp->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp->GetBinLowEdge(k)+hhbottom_pp->GetBinWidth(k))>statscale_uplim) {
0616       hhbottom_pp->SetBinContent(k,0.);
0617       hhbottom_pp->SetBinError(k,0.);
0618     }
0619     else {
0620       double tmp = ppcorr * fbottom_nosup_corr->Eval(hhbottom[nbins]->GetBinCenter(k)) * hhbottom[nbins]->GetFunction("expo")->Eval(hhbottom[nbins]->GetBinCenter(k));
0621       double tmprnd = myrandom->Poisson(tmp);
0622       if(tmprnd<0.) { tmprnd=0.; }
0623       hhbottom_pp->SetBinContent(k,tmprnd);
0624       hhbottom_pp->SetBinError(k,sqrt(tmprnd));
0625     }
0626   }
0627 
0628   sprintf(tmpname,"hhdy_pp");
0629   hhdy_pp = (TH1D*)hhdy[nbins]->Clone(tmpname);
0630   for(int k=1; k<=hhdy_pp->GetNbinsX(); k++) {
0631     if(hhdy_pp->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp->GetBinLowEdge(k)+hhdy_pp->GetBinWidth(k))>statscale_uplim) {
0632       hhdy_pp->SetBinContent(k,0.);
0633       hhdy_pp->SetBinError(k,0.);
0634     }
0635     else {
0636       double tmp = ppcorr * hhdy[nbins]->GetFunction("expo")->Eval(hhdy[nbins]->GetBinCenter(k));
0637       double tmprnd = myrandom->Poisson(tmp);
0638       if(tmprnd<0.) { tmprnd=0.; }
0639       hhdy_pp->SetBinContent(k,tmprnd);
0640       hhdy_pp->SetBinError(k,sqrt(tmprnd));
0641     }
0642   }
0643 
0644   sprintf(tmpname,"hhcorrbg_pp");
0645   hhcorrbg_pp = (TH1D*)hhbottom_pp->Clone(tmpname);
0646   hhcorrbg_pp->Add(hhdy_pp);
0647     hhcorrbg_pp->SetMarkerColor(kBlack);
0648     hhcorrbg_pp->SetLineColor(kBlack);
0649     hhbottom_pp->SetLineColor(kBlue);
0650     hhdy_pp->SetLineColor(kGreen+2);
0651   sprintf(tmpname,"hhall_pp");
0652   hhall_pp = (TH1D*)hhcorrbg_pp->Clone(tmpname);
0653   hhall_pp->Add(hhupspp);
0654     hhall_pp->SetLineColor(kMagenta);
0655     hhall_pp->SetMarkerColor(kMagenta);
0656 
0657 
0658 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
0659 
0660   hhcorrbg_pp->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0661   hhcorrbg_pp->GetFunction("expo")->SetLineColor(kBlack);
0662   hhbottom_pp->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0663   hhbottom_pp->GetFunction("expo")->SetLineColor(kBlue);
0664   hhdy_pp->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
0665   hhdy_pp->GetFunction("expo")->SetLineColor(kGreen+2);
0666 
0667   hhall_pp->SetAxisRange(7.,12.);
0668   hhcorrbg_pp->Draw("pehist");
0669   hhbottom_pp->Draw("histsame");
0670   hhdy_pp->Draw("histsame");
0671 
0672 
0673 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
0674   hhall_pp->SetAxisRange(7.,12.);
0675   hhall_pp->Draw("pehist");
0676   hhcorrbg_pp->Draw("pesame");
0677   hhbottom_pp->Draw("same");
0678   hhdy_pp->Draw("same");
0679 
0680 //-----------------------------------------------------------------------------------------
0681 // Combinatorial BG calculated for 10B central AuAu events
0682 //-----------------------------------------------------------------------------------------
0683 
0684 TCanvas* cdummy = new TCanvas("cdummy","cdummy",0,0,500,500);
0685 
0686 f = new TFile("fakee_eideff09.root");
0687   sprintf(tmpname,"hhfakefake_15"); // 15 is integrated over pT
0688   hhfakefake[nbins] = (TH1D*)f->Get(tmpname);
0689   hhfakefake[nbins]->SetDirectory(gROOT);
0690 f->Close();
0691 
0692 f = new TFile("crossterms_eideff09.root");
0693   sprintf(tmpname,"hhfakehf_15");
0694   hhfakehf[nbins] = (TH1D*)f->Get(tmpname);
0695   hhfakehf[nbins]->SetDirectory(gROOT);
0696 f->Close();
0697 
0698 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
0699 fbg->SetParameters(10., -1.0, 4., -0.1);
0700 fbg->SetParLimits(1.,-999.,0.);
0701 fbg->SetParLimits(3.,-999.,0.);
0702 
0703   sprintf(tmpname,"hhcombbg_15");
0704   hhcombbg[nbins] = (TH1D*)hhfakefake[nbins]->Clone(tmpname);
0705   hhcombbg[nbins]->Add(hhfakehf[nbins]);
0706   sprintf(tmpname,"hhcombbg_scaled_15");
0707   hhcombbg_scaled[nbins] = (TH1D*)hhcombbg[nbins]->Clone(tmpname);
0708   fbg->SetParameters(10., -1.0, 4., -0.1); 
0709   hhcombbg[nbins]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0710 
0711     for(int k=1; k<=hhcombbg[nbins]->GetNbinsX(); k++) {
0712       if(hhcombbg[nbins]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[nbins]->GetBinLowEdge(k)+hhcombbg[nbins]->GetBinWidth(k))>statscale_uplim) {
0713         hhcombbg_scaled[nbins]->SetBinContent(k,0.);
0714         hhcombbg_scaled[nbins]->SetBinError(k,0.);
0715       }
0716       else {
0717         double tmp = statscale * hhcombbg[nbins]->GetFunction("fbg")->Eval(hhcombbg[nbins]->GetBinCenter(k));
0718         double tmprnd = myrandom->Poisson(tmp);
0719         if(tmprnd<0.) { tmprnd=0.; }
0720         hhcombbg_scaled[nbins]->SetBinContent(k,tmprnd);
0721         hhcombbg_scaled[nbins]->SetBinError(k,sqrt(tmprnd));
0722       }
0723     }
0724   hhcombbg_scaled[nbins]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0725 
0726 delete cdummy;
0727 
0728 TCanvas* C1 = new TCanvas("C1","Combinatorial BG Central Au+Au",100,100,600,600);
0729 C1->SetLogy();
0730   hhfakefake[nbins]->SetAxisRange(7.0,14.0);
0731   hhfakefake[nbins]->SetMinimum(0.1);
0732   hhfakefake[nbins]->SetMaximum(5000.);
0733   hhfakefake[nbins]->SetLineColor(kGreen+2);
0734   hhfakefake[nbins]->SetLineWidth(2);
0735   hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
0736   hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
0737   hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
0738   hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
0739   hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
0740   hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
0741   hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
0742   hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
0743   hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
0744   hhfakefake[nbins]->Draw("e");
0745 
0746   hhfakehf[nbins]->SetLineColor(kOrange+4);
0747   hhfakehf[nbins]->SetLineWidth(2);
0748   hhfakehf[nbins]->Draw("esame");
0749 
0750   hhcombbg[nbins]->SetLineColor(kBlack);
0751   hhcombbg[nbins]->SetLineWidth(2);
0752   hhcombbg[nbins]->Draw("esame");
0753 
0754 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG Central Au+Au",100,100,600,600);
0755 C1sc->SetLogy(); 
0756   hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
0757   hhcombbg_scaled[nbins]->Draw("esame");
0758 
0759 // Scaled Combinatorial background vs centrality
0760 
0761 for(int i=0; i<nbins; i++) {
0762 
0763   sprintf(tmpname,"hhcombbg_%d",i);
0764   hhcombbg_scaled[i] = (TH1D*)hhcombbg_scaled[nbins]->Clone(tmpname);
0765 
0766     for(int k=1; k<=hhcombbg_scaled[nbins]->GetNbinsX(); k++) {
0767       if(hhcombbg_scaled[nbins]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg_scaled[nbins]->GetBinLowEdge(k)+hhcombbg_scaled[nbins]->GetBinWidth(k))>statscale_uplim) {
0768         hhcombbg_scaled[i]->SetBinContent(k,0.);
0769         hhcombbg_scaled[i]->SetBinError(k,0.);
0770       }
0771       else {
0772         double tmp = Npionpairs[i] * hhcombbg_scaled[nbins]->GetFunction("fbg")->Eval(hhcombbg_scaled[nbins]->GetBinCenter(k));
0773         double tmprnd = myrandom->Poisson(tmp);
0774         if(tmprnd<0.) { tmprnd=0.; }
0775         hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
0776         hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
0777       }
0778     }
0779   hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
0780 
0781 } // end over centrality bins
0782 
0783 TCanvas* c_comb_scaled = new TCanvas("c_comb_scaled","Combinatorial Background vs. Centrality",200,100,1200,600);
0784 c_comb_scaled->Divide(4,2);
0785 for(int i=0; i<nbins; i++) {
0786   c_comb_scaled->cd(i+1);
0787   hhcombbg_scaled[i]->SetAxisRange(8.5,11.0); hhcombbg_scaled[i]->SetMarkerStyle(1); hhcombbg_scaled[i]->Draw("pe");
0788   sprintf(tlchar,"%d-%d",10*i,10*(i+1));   tl[i] = new TLatex(9.0,hhcombbg_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
0789 }
0790 
0791 
0792 //---------------------------------------------
0793 //  Au+Au Signal + all BG
0794 //---------------------------------------------
0795 
0796 for(int i=0; i<nbins; i++) {
0797   sprintf(tmpname,"hhtotbg_scaled_%d",i);
0798   hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
0799   hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
0800 }
0801 
0802 for(int i=0; i<nbins; i++) {
0803   sprintf(tmpname,"hhall_scaled_%d",i);
0804   hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
0805   hhall_scaled[i]->Add(hhups[i]);
0806 }
0807 
0808 TCanvas* c000 = new TCanvas("c000","Au+Au Signal + All Background vs. Centrality",200,200,1200,600);
0809 c000->Divide(4,2);
0810 for(int i=0; i<nbins; i++) {
0811   c000->cd(i+1);
0812   hhall_scaled[i]->SetAxisRange(8.5,11.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
0813   sprintf(tlchar,"%d-%d",10*i,10*(i+1));   tl[i] = new TLatex(9.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
0814 }
0815 
0816 /*
0817 // Fit Signal + Correlated BG
0818 
0819 // hhfit[] is signal + correlated bg
0820 for(int i=0; i<nbins; i++) {
0821   sprintf(tmpname,"hhfit_%d",i);
0822   hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
0823   for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
0824     hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
0825     hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
0826   }
0827 }
0828 
0829 //---------------------------------------------------------------------
0830 // plot signal + correlated bg for all pT
0831 //---------------------------------------------------------------------
0832 
0833   double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
0834   double tmppar1 = corrbgfitpar1;
0835 cout << "###### " << tmppar0 << " " << tmppar1 << endl;
0836 
0837 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
0838 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
0839 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
0840 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
0841 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
0842 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
0843 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
0844 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
0845 fSandBpp->SetParameter(7,tmppar0);
0846 fSandBpp->SetParameter(8,tmppar1);
0847 fSandBpp->SetLineColor(kBlue);
0848 fSandBpp->SetLineStyle(2);
0849 
0850 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
0851 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
0852 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
0853 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
0854 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
0855 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
0856 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
0857 fSandBauau->SetParameter(7,tmppar0);
0858 fSandBauau->SetParameter(8,tmppar1);
0859 fSandBauau->SetLineColor(kRed);
0860 
0861 TCanvas* cfitall = new TCanvas("cfitall","Fit Central",270,270,600,600);
0862   hhfit[0]->SetAxisRange(7.0,14.);
0863   hhfit[0]->GetXaxis()->CenterTitle();
0864   hhfit[0]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
0865   hhfit[0]->GetXaxis()->SetTitleOffset(1.1);
0866   hhfit[0]->GetXaxis()->SetLabelSize(0.045);
0867   hhfit[0]->GetYaxis()->CenterTitle();
0868   hhfit[0]->GetYaxis()->SetLabelSize(0.045);
0869   hhfit[0]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
0870   hhfit[0]->GetYaxis()->SetTitleOffset(1.5);
0871   hhfit[0]->Draw("pehist");
0872 //  fSandBpp->Draw("same");
0873 //  fSandBauau->Draw("same");
0874     //hhcorrbg_scaled[0]->SetLineColor(kRed);
0875     //hhcorrbg_scaled[0]->Scale(0.82);
0876     //hhcorrbg_scaled[0]->Scale(0.70);
0877     //hhcorrbg_scaled[0]->Draw("histsame");
0878 
0879   TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
0880   fmycorrbg->SetParameters(tmppar0,tmppar1);
0881   fmycorrbg->SetLineStyle(2);
0882   fmycorrbg->SetLineColor(kRed);
0883   fmycorrbg->Draw("same");
0884 
0885 //double myheight = fTCBauau->GetParameter(0);
0886 //TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
0887 //ld1->SetTextSize(0.035);
0888 //ld1->Draw();
0889 //TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
0890 //ld2->SetTextSize(0.035);
0891 //ld2->Draw();
0892 
0893 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
0894 TH1D* hhfit_tmp = (TH1D*)hhfit[0]->Clone("hhfit_tmp");
0895 hhfit_tmp->SetAxisRange(8.0,11.);
0896 hhfit_tmp->Draw("pehist");
0897 //fSandBauau->Draw("same");
0898 fmycorrbg->Draw("same");
0899 */
0900 
0901 //----------------------------------------------------------------------
0902 
0903 // plot signal + both backgrounds for central Au+Au
0904 
0905 TCanvas* callpt = new TCanvas("callpt","Signal + All BG Central Au+Au",300,300,600,600);
0906 
0907   hhall_scaled[0]->GetXaxis()->SetTitle("Invariant mass GeV/c");
0908   hhall_scaled[0]->SetLineColor(kBlack);
0909   hhall_scaled[0]->SetMarkerColor(kBlack);
0910   hhall_scaled[0]->SetMarkerStyle(20);
0911   hhall_scaled[0]->SetAxisRange(8.0,10.8);
0912   hhall_scaled[0]->Draw("pehist");
0913     hhcombbg_scaled[0]->SetLineColor(kBlue);
0914     hhcombbg_scaled[0]->Draw("histsame");
0915       hhcorrbg_scaled[0]->SetLineColor(kRed);
0916       hhcorrbg_scaled[0]->Draw("histsame");
0917 
0918 
0919 
0920 //----------------------------------------------------------------
0921 //  Calculate RAA
0922 //----------------------------------------------------------------
0923 
0924 double u1start = 9.25;
0925 double u1stop  = 9.65;
0926 double u2start = 9.80;
0927 double u2stop  = 10.20;
0928 double u3start = 10.20;
0929 double u3stop  = 10.55;
0930 cout << "kuku2" << endl;
0931   double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
0932   for(int i=0; i<nbins; i++) { 
0933     //raa1[i] = raa1s[i];
0934     //raa2[i] = raa2s[i];
0935     //raa3[i] = raa3s[i];
0936     raa1[i] = grRAA1S->Eval(Npart[i]);
0937     raa2[i] = grRAA2S->Eval(Npart[i]);
0938     raa3[i] = grRAA3S->Eval(Npart[i]); 
0939     //raa3[i] = -1.; 
0940     cout << "Npart, Raa = " << Npart[i] << " " << raa1[i] << " " << raa2[i] << " " << raa3[i] << endl;
0941   }
0942   int fbin1 = hhall_scaled[0]->FindBin(u1start + 0.001);
0943   int lbin1 = hhall_scaled[0]->FindBin(u1stop  - 0.001);
0944   int fbin2 = hhall_scaled[0]->FindBin(u2start + 0.001);
0945   int lbin2 = hhall_scaled[0]->FindBin(u2stop  - 0.001);
0946   int fbin3 = hhall_scaled[0]->FindBin(u3start + 0.001);
0947   int lbin3 = hhall_scaled[0]->FindBin(u3stop  - 0.001);
0948   cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
0949   cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
0950   cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
0951   cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
0952   cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
0953   cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
0954 
0955   double sum1[99]   = {0.};
0956   double truesum1[99]   = {0.};
0957   double ersum1[99] = {0.};
0958   double sumpp1   = 0.;
0959   double ersumpp1 = 0.;
0960   double sum2[99]   = {0.};
0961   double truesum2[99]   = {0.};
0962   double ersum2[99] = {0.};
0963   double sumpp2   = 0.;
0964   double ersumpp2 = 0.;
0965   double sum3[99]   = {0.};
0966   double truesum3[99]   = {0.};
0967   double ersum3[99] = {0.};
0968   double sumpp3   = 0.;
0969   double ersumpp3 = 0.;
0970 
0971   double sumsum1[99]   = {0.};
0972   double sumsum2[99]   = {0.};
0973   double sumsum3[99]   = {0.};
0974   double sumsum1pp   = 0.;
0975   double sumsum2pp   = 0.;
0976   double sumsum3pp   = 0.;
0977 
0978   for(int j=fbin1; j<=lbin1; j++) {
0979     sumpp1   += hhups1pp->GetBinContent(j);
0980     ersumpp1 += hhupspp->GetBinError(j)*hhupspp->GetBinError(j);
0981   }
0982   for(int j=fbin2; j<=lbin2; j++) {
0983     sumpp2   += hhups2pp->GetBinContent(j);
0984     ersumpp2 += hhupspp->GetBinError(j)*hhupspp->GetBinError(j);
0985   }
0986   for(int j=fbin3; j<=lbin3; j++) {
0987     sumpp3   += hhups3pp->GetBinContent(j);
0988     ersumpp3 += hhupspp->GetBinError(j)*hhupspp->GetBinError(j);
0989   }
0990 
0991   for(int i=0; i<nbins; i++) {
0992 
0993     //double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
0994 
0995     for(int j=fbin1; j<=lbin1; j++) {
0996       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)));
0997       truesum1[i] += hhups1[i]->GetBinContent(j);
0998       ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
0999     }
1000       sumsum1[i]     = truesum1[i];                   // direct count in mass range
1001       sumsum1pp   = sumpp1;       
1002       //sumsum1[i]     = hhups1[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1003       //sumsum1pp   = hhups1pp->GetEntries();
1004       //sumsum1[i]     = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1005       //sumsum1pp   = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1006 
1007         if(sumsum1[i]>0. && sumsum1pp>0.) {
1008           erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1/sumsum1pp/sumsum1pp);
1009         } else {raa1[i]=-1.0; erraa1[i] = 999.; }
1010 
1011     for(int j=fbin2; j<=lbin2; j++) {
1012       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)));
1013       truesum2[i] += hhups2[i]->GetBinContent(j);
1014       ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1015     }
1016       sumsum2[i]     = truesum2[i];                   // direct count in mass range
1017       sumsum2pp   = sumpp2;       
1018       //sumsum2[i]     = hhups2[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1019       //sumsum2pp   = hhups2pp->GetEntries();
1020       //sumsum2[i]     = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;  // total number of upsilons in pT bin
1021       //sumsum2pp   = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1022       
1023         if(sumsum2[i]>0. && sumsum2pp>0.) {
1024           erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2/sumsum2pp/sumsum2pp);
1025         } else {raa2[i]=-1.0; erraa2[i] = 999.; }
1026 
1027     for(int j=fbin3; j<=lbin3; j++) {
1028       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)));
1029       truesum3[i] += hhups3[i]->GetBinContent(j);
1030       ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1031     }
1032       sumsum3[i]     = truesum3[i];                   // direct count in mass range
1033       sumsum3pp   = sumpp3;       
1034       //sumsum3[i]     = hhups3[i]->GetEntries();       // total number of upsilons in pT bin (rounded up)
1035       //sumsum3pp   = hhups3pp->GetEntries();
1036       //sumsum3[i]     = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;   // total number of upsilons in pT bin
1037       //sumsum3pp   = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1038       
1039         if(truesum3[i]>0. && sumpp3>0.) {
1040           erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3/sumsum3pp/sumsum3pp);
1041         } else {raa3[i]=-1.0; erraa3[i] = 999.; }
1042   
1043   } // end loop over centrality
1044 
1045   erraa2[2] = erraa2[2]*0.70;
1046 
1047   for(int i=0; i<nbins; i++) {
1048     cout << "Npart, Raa = " << Npart[i] << " " << raa1[i] << " " << raa2[i] << " " << raa3[i] << endl;
1049   }
1050 
1051 cout << "====== Y(1S):" << endl;
1052   for(int i=0; i<nbins; i++) {
1053     cout << "   " << i << " " << sumsum1[i] << "(" << Nups1[i] << ")" << " +- " << sqrt(ersum1[i]) 
1054          << " \t\t pp: " << sumsum1pp << " +- " << sqrt(ersumpp1) << endl;
1055   }
1056 cout << "====== Y(2S):" << endl;
1057   for(int i=0; i<nbins; i++) {
1058     cout << "   " << i << " " << sumsum2[i] << "(" << Nups2[i] << ")" << " +- " << sqrt(ersum2[i]) 
1059          << " \t\t pp: " << sumsum2pp << " +- " << sqrt(ersumpp2) << endl;
1060   }
1061 cout << "====== Y(3S):" << endl;
1062   for(int i=0; i<nbins; i++) {
1063     cout << "   " << i << " " << sumsum3[i] << "(" << Nups3[i] << ")" << " +- " << sqrt(ersum3[i]) 
1064          << " \t\t pp: " << sumsum3pp << " +- " << sqrt(ersumpp3) << endl;
1065   }
1066 
1067 
1068 double raa3_rebin[9],raapt_rebin3[9],erraa3_rebin[9];
1069 double sum3_rebin[9],ersum3_rebin[9],sum3pp_rebin[9],ersumpp3_rebin[9];
1070 
1071 // Rebin Y(3S) by 4
1072 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin3[0]);
1073 raa3_rebin[1] = grRAA3S->Eval(raapt_rebin3[1]);
1074 sum3_rebin[0] = truesum3[0]+truesum3[1];
1075 sum3_rebin[1] = truesum3[2]+truesum3[3];
1076 ersum3_rebin[0] = ersum3[0]+ersum3[1];
1077 ersum3_rebin[1] = ersum3[2]+ersum3[3];
1078 sum3pp_rebin[0] = sumsum3pp;
1079 sum3pp_rebin[1] = sumsum3pp;
1080 ersumpp3_rebin[0] = ersumpp3;
1081 ersumpp3_rebin[1] = ersumpp3;
1082   erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1083   erraa3_rebin[1] = raa3[1]*sqrt(ersum3_rebin[1]/sum3_rebin[1]/sum3_rebin[1] + ersumpp3_rebin[1]/sum3pp_rebin[1]/sum3pp_rebin[1]);
1084 
1085 //-------------------------------------------------
1086 //  Plot RAA
1087 //-------------------------------------------------
1088 
1089 int npts1 = nbins;
1090 int npts2 = nbins;
1091 int npts3 = nbins;
1092 int npts2_rebin = 4;
1093 int npts3_rebin = 2;
1094 
1095 TCanvas* craa = new TCanvas("craa","R_{dAu}",120,120,800,600);
1096 TH2F* hh2 = new TH2F("hh2"," ",10,0.,10.,10,0.,1.5);
1097 hh2->GetXaxis()->SetTitle("N_{coll}");
1098 hh2->GetXaxis()->SetTitleOffset(0.9);
1099 hh2->GetXaxis()->SetTitleColor(1);
1100 hh2->GetXaxis()->SetTitleSize(0.050);
1101 hh2->GetXaxis()->SetLabelSize(0.040);
1102 hh2->GetYaxis()->SetTitle("R_{dAu}");
1103 hh2->GetYaxis()->SetTitleOffset(0.9);
1104 hh2->GetYaxis()->SetTitleSize(0.050);
1105 hh2->GetYaxis()->SetLabelSize(0.040);
1106 hh2->Draw();
1107 
1108 double xx1[nbins+1]; for(int i=0; i<nbins; i++) {xx1[i] = Npart[i];}
1109 double xx2[nbins+1]; for(int i=0; i<nbins; i++) {xx2[i] = Npart[i] - 0.1;}
1110 double xx3[nbins+1]; for(int i=0; i<nbins; i++) {xx3[i] = Npart[i] + 0.1;}
1111 double xx3_rebin[nbins+1]; 
1112 xx3_rebin[0] = (Ncoll[0]+Ncoll[1])/2.;
1113 xx3_rebin[1] = (Ncoll[2]+1.2*Ncoll[3])/2.; // last bin wider
1114 
1115 TGraphErrors* gr1 = new TGraphErrors(npts1,xx1,raa1,0,erraa1);
1116 gr1->SetMarkerStyle(kFullCircle);
1117 gr1->SetMarkerColor(kBlack);
1118 gr1->SetLineColor(kBlack);
1119 gr1->SetLineWidth(2);
1120 gr1->SetMarkerSize(1.5);
1121 gr1->SetName("gr1");
1122 gr1->Draw("p");
1123 
1124 TGraphErrors* gr2 = new TGraphErrors(npts2,xx2,raa2,0,erraa2);
1125 gr2->SetMarkerStyle(kFullSquare);
1126 gr2->SetMarkerColor(kRed);
1127 gr2->SetLineColor(kRed);
1128 gr2->SetLineWidth(2);
1129 gr2->SetMarkerSize(1.5);
1130 gr2->SetName("gr2");
1131 gr2->Draw("p");
1132 
1133 //TGraphErrors* gr2_rebin = new TGraphErrors(npts2_rebin,xx2_rebin,raa2_rebin,0,erraa2_rebin);
1134 //gr2_rebin->SetMarkerStyle(20);
1135 //gr2_rebin->SetMarkerColor(kRed);
1136 //gr2_rebin->SetLineColor(kRed);
1137 //gr2_rebin->SetLineWidth(2);
1138 //gr2_rebin->SetMarkerSize(1.5);
1139 //gr2_rebin->SetName("gr2");
1140 //gr2_rebin->Draw("p");
1141 
1142 //--- 3S state -------------------
1143 // double dummy[9]; for(int i=0; i<9; i++) {dummy[i]=-99.;}
1144  TGraphErrors* gr3 = new TGraphErrors(2,xx3_rebin,raa3,0,erraa3_rebin);
1145  gr3->SetMarkerStyle(kFullDiamond);
1146  gr3->SetMarkerColor(kBlue);
1147  gr3->SetLineColor(kBlue);
1148  gr3->SetLineWidth(2);
1149  gr3->SetMarkerSize(2.5);
1150  gr3->SetName("gr3");
1151  gr3->Draw("p");
1152 /*
1153  TGraphErrors* gr3_rebin = new TGraphErrors(npts3_rebin,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1154  gr3_rebin->SetMarkerStyle(20);
1155  gr3_rebin->SetMarkerColor(kBlue);
1156  gr3_rebin->SetLineColor(kBlue);
1157  gr3_rebin->SetLineWidth(2);
1158  gr3_rebin->SetMarkerSize(1.5);
1159  gr3_rebin->SetName("gr3");
1160  gr3_rebin->Draw("p");
1161 */
1162 /*
1163  TGraphErrors* gr3_rebin = new TGraphErrors(1,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1164  gr3_rebin->SetMarkerStyle(20);
1165  gr3_rebin->SetMarkerColor(kBlue);
1166  gr3_rebin->SetLineColor(kBlue);
1167  gr3_rebin->SetLineWidth(2);
1168  gr3_rebin->SetMarkerSize(1.5);
1169  gr3_rebin->SetName("gr3");
1170  gr3_rebin->Draw("p");
1171 */
1172 //TArrow* aa[9];
1173 //TLine* ll[9];
1174 //erraa3[3] = erraa3[3]*1.5;
1175 //erraa3[4] = erraa3[4]*0.6;
1176 //erraa3[4] = erraa3[4]*1.5;
1177 //erraa3[6] = erraa3[6]*1.2;
1178 /*
1179 for(int i=0; i<npts3; i++) {
1180   aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1181   aa[i]->SetLineColor(kBlue);
1182   aa[i]->SetLineWidth(2);
1183   aa[i]->Draw();
1184   ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1185   ll[i]->SetLineColor(kBlue);
1186   ll[i]->SetLineWidth(2);
1187   ll[i]->Draw();
1188 }
1189 */
1190 //--- end 3S state --------------
1191 
1192 TLegend *leg = new TLegend(0.70,0.20,0.88,0.38);
1193 leg->SetBorderSize(0);
1194 leg->SetFillColor(10);
1195 leg->SetFillStyle(1001);
1196 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1197 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1198 //TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1199 TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","p");
1200 leg->Draw();
1201 
1202 TLatex* l1 = new TLatex(1.,0.40,"#font[72]{sPHENIX} Projection"); l1->SetTextFont(42); l1->Draw();
1203 //TLatex* l11 = new TLatex(20.,0.90,"0-10% cent. Au+Au, Years 1-3"); l11->SetTextFont(42); l11->Draw();
1204 TLatex* l2 = new TLatex(1.,0.30,"190B rec. d+Au events"); l2->SetTextFont(42); l2->Draw();
1205 TLatex* l3 = new TLatex(1.,0.20,"62 pb^{-1} samp. #it{p+p}"); l3->SetTextFont(42); l3->Draw();
1206 
1207 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1208 lll->SetLineColor(kBlue);
1209 lll->SetLineWidth(2);
1210 //lll->Draw();
1211 
1212 /*
1213 grRAA1S->Draw("l");
1214 grRAA1S_eta1->Draw("l");
1215 grRAA1S_eta3->Draw("l");
1216 grRAA2S->Draw("l");
1217 grRAA2S_eta1->Draw("l");
1218 grRAA2S_eta3->Draw("l");
1219 grRAA3S->Draw("l");
1220 grRAA3S_eta1->Draw("l");
1221 grRAA3S_eta3->Draw("l");
1222 */
1223 
1224 //==================================================================================
1225 
1226 /*
1227 fout = new TFile("test.root","recreate");
1228 for(int i=0; i<nbins+1; i++) {
1229   sprintf(tmpname,"hhfit_%d",i);
1230   hhfit[i]->Write();
1231   hhups1pp[i]->Write();
1232   hhups2pp[i]->Write();
1233   hhups3pp[i]->Write();
1234   hhupspp[i]->Write();
1235   hhups1[i]->Write();
1236   hhups2[i]->Write();
1237   hhups3[i]->Write();
1238   hhups[i]->Write();
1239   hhall_pp[i]->Write();
1240   hhcorrbg_pp[i]->Write();
1241   hhcorrbg_scaled[i]->Write();
1242 }
1243 fout->Write();
1244 fout->Close();
1245 */
1246 
1247 }
1248