analysis/QuarkoniaBG/compare.C

0001 void compare() {
0002
0003 // From ppg30
0004
0005 //Differential Invariant Cross Section
0006 //Rapididty coverage measured: |Dy| = 0.7
0007 //The data have been corrected to |Dy|=1 making a
0008 //flat extrapolation from 0.7 to 1.
0009 //-----------------------------------------------------------------
0010 //E d3sigma/ dp3 = 1/2pi 1/pt d2sigma/dpt dy      UNITS [mb/GeV^2]
0011 //-----------------------------------------------------------------
0012 // If you need to convert this to yield-per-event you must divide by sigma_inel(sqrt_s = 200GeV) = 42 mb
0013
0014 float x[18],y[18],ex[18],ey[18];
0015 float x1[18],y1[18],ey1[18];
0016 float x2[18],y2[18],ey2[18];
0017
0018 //piplus
0019 x1[0]=0.35, y1[0]=25.6287, ex[0]=0, ey1[0]=0.246673;
0020 x1[1]=0.45, y1[1]=14.6531, ex[1]=0, ey1[1]=0.14773;
0021 x1[2]=0.55, y1[2]=8.62473, ex[2]=0, ey1[2]=0.0943008;
0022 x1[3]=0.65, y1[3]=4.80858, ex[3]=0, ey1[3]=0.0579757;
0023 x1[4]=0.75, y1[4]=2.85599, ex[4]=0, ey1[4]=0.03862;
0024 x1[5]=0.85, y1[5]=1.71022, ex[5]=0, ey1[5]=0.0258761;
0025 x1[6]=0.95, y1[6]=1.08056, ex[6]=0, ey1[6]=0.0183583;
0026 x1[7]=1.05, y1[7]=0.68761, ex[7]=0, ey1[7]=0.013222;
0027 x1[8]=1.15, y1[8]=0.431937, ex[8]=0, ey1[8]=0.00938769;
0028 x1[9]=1.25, y1[9]=0.308209, ex[9]=0, ey1[9]=0.00732418;
0029 x1[10]=1.35, y1[10]=0.221919, ex[10]=0, ey1[10]=0.00582217;
0030 x1[11]=1.45, y1[11]=0.147961, ex[11]=0, ey1[11]=0.00437927;
0031 x1[12]=1.6, y1[12]=0.0825923, ex[12]=0, ey1[12]=0.00205954;
0032 x1[13]=1.8, y1[13]=0.0410802, ex[13]=0, ey1[13]=0.00131989;
0033 x1[14]=2, y1[14]=0.0211058, ex[14]=0, ey1[14]=0.000876262;
0034 x1[15]=2.2, y1[15]=0.0116705, ex[15]=0, ey1[15]=0.000638873;
0035 x1[16]=2.4, y1[16]=0.00787318, ex[16]=0, ey1[16]=0.000557839;
0036 x1[17]=2.6, y1[17]=0.00447392, ex[17]=0, ey1[17]=0.000458999;
0037
0038 // piminus
0039 x2[0]=0.35, y2[0]=25.8453, ex[0]=0, ey2[0]=0.21199;
0040 x2[1]=0.45, y2[1]=14.5864, ex[1]=0, ey2[1]=0.127527;
0041 x2[2]=0.55, y2[2]=8.27867, ex[2]=0, ey2[2]=0.0789266;
0042 x2[3]=0.65, y2[3]=4.72545, ex[3]=0, ey2[3]=0.0498388;
0043 x2[4]=0.75, y2[4]=2.756, ex[4]=0, ey2[4]=0.032014;
0044 x2[5]=0.85, y2[5]=1.66385, ex[5]=0, ey2[5]=0.0215601;
0045 x2[6]=0.95, y2[6]=1.05331, ex[6]=0, ey2[6]=0.0152551;
0046 x2[7]=1.05, y2[7]=0.648492, ex[7]=0, ey2[7]=0.0105207;
0047 x2[8]=1.15, y2[8]=0.431804, ex[8]=0, ey2[8]=0.0079342;
0048 x2[9]=1.25, y2[9]=0.286134, ex[9]=0, ey2[9]=0.00585618;
0049 x2[10]=1.35, y2[10]=0.199094, ex[10]=0, ey2[10]=0.004624;
0050 x2[11]=1.45, y2[11]=0.133006, ex[11]=0, ey2[11]=0.003509;
0051 x2[12]=1.6, y2[12]=0.0795068, ex[12]=0, ey2[12]=0.001787;
0052 x2[13]=1.8, y2[13]=0.0422071, ex[13]=0, ey2[13]=0.001195;
0053 x2[14]=2, y2[14]=0.0208344, ex[14]=0, ey2[14]=0.000781061;
0054 x2[15]=2.2, y2[15]=0.0131348, ex[15]=0, ey2[15]=0.000625054;
0055 x2[16]=2.4, y2[16]=0.0070146, ex[16]=0, ey2[16]=0.000479789;
0056 x2[17]=2.6, y2[17]=0.00425537, ex[17]=0, ey2[17]=0.000408;
0057
0058 // Kplus
0059 float x3[13],y3[13],ey3[13];
0060 x3[0]=0.45, y3[0]=1.95293, ex[0]=0, ey3[0]=0.0724524;
0061 x3[1]=0.55, y3[1]=1.40404, ex[1]=0, ey3[1]=0.0466304;
0062 x3[2]=0.65, y3[2]=0.925078, ex[2]=0, ey3[2]=0.0296139;
0063 x3[3]=0.75, y3[3]=0.624339, ex[3]=0, ey3[3]=0.020472;
0064 x3[4]=0.85, y3[4]=0.403909, ex[4]=0, ey3[4]=0.0142644;
0065 x3[5]=0.95, y3[5]=0.281567, ex[5]=0, ey3[5]=0.0107128;
0066 x3[6]=1.05, y3[6]=0.199669, ex[6]=0, ey3[6]=0.00826777;
0067 x3[7]=1.15, y3[7]=0.137853, ex[7]=0, ey3[7]=0.00620796;
0068 x3[8]=1.25, y3[8]=0.101317, ex[8]=0, ey3[8]=0.00483626;
0069 x3[9]=1.35, y3[9]=0.0694105, ex[9]=0, ey3[9]=0.00362463;
0070 x3[10]=1.45, y3[10]=0.0526523, ex[10]=0, ey3[10]=0.00301769;
0071 x3[11]=1.6, y3[11]=0.0324873, ex[11]=0, ey3[11]=0.00151623;
0072 x3[12]=1.8, y3[12]=0.0179392, ex[12]=0, ey3[12]=0.00100964;
0073
0074 // Kminus
0075 float x4[13],y4[13],ey4[13];
0076 x4[0]=0.45, y4[0]=1.66952, ex[0]=0, ey4[0]=0.0568774;
0077 x4[1]=0.55, y4[1]=1.25238, ex[1]=0, ey4[1]=0.0374863;
0078 x4[2]=0.65, y4[2]=0.828785, ex[2]=0, ey4[2]=0.0243554;
0079 x4[3]=0.75, y4[3]=0.604526, ex[3]=0, ey4[3]=0.0177572;
0080 x4[4]=0.85, y4[4]=0.406696, ex[4]=0, ey4[4]=0.0125949;
0081 x4[5]=0.95, y4[5]=0.261138, ex[5]=0, ey4[5]=0.00878815;
0082 x4[6]=1.05, y4[6]=0.189842, ex[6]=0, ey4[6]=0.00685834;
0083 x4[7]=1.15, y4[7]=0.118611, ex[7]=0, ey4[7]=0.00483604;
0084 x4[8]=1.25, y4[8]=0.0977945, ex[8]=0, ey4[8]=0.00407271;
0085 x4[9]=1.35, y4[9]=0.0673309, ex[9]=0, ey4[9]=0.00315524;
0086 x4[10]=1.45, y4[10]=0.0439384, ex[10]=0, ey4[10]=0.00239946;
0087 x4[11]=1.6, y4[11]=0.0283855, ex[11]=0, ey4[11]=0.00125486;
0088 x4[12]=1.8, y4[12]=0.0181311, ex[12]=0, ey4[12]=0.000911162;
0089
0090 // protons
0091 float x5[15],y5[15],ey5[15];
0092 x5[0]=0.65, y5[0]=0.467596, ex[0]=0, ey5[0]=0.0124957;
0093 x5[1]=0.75, y5[1]=0.361103, ex[1]=0, ey5[1]=0.0100168;
0094 x5[2]=0.85, y5[2]=0.250201, ex[2]=0, ey5[2]=0.00749516;
0095 x5[3]=0.95, y5[3]=0.181417, ex[3]=0, ey5[3]=0.00593336;
0096 x5[4]=1.05, y5[4]=0.126166, ex[4]=0, ey5[4]=0.00445335;
0097 x5[5]=1.15, y5[5]=0.0902262, ex[5]=0, ey5[5]=0.00351448;
0098 x5[6]=1.25, y5[6]=0.0651584, ex[6]=0, ey5[6]=0.00275751;
0099 x5[7]=1.35, y5[7]=0.0448966, ex[7]=0, ey5[7]=0.00214243;
0100 x5[8]=1.45, y5[8]=0.0328418, ex[8]=0, ey5[8]=0.00172412;
0101 x5[9]=1.6, y5[9]=0.0203516, ex[9]=0, ey5[9]=0.000881656;
0102 x5[10]=1.8, y5[10]=0.0111527, ex[10]=0, ey5[10]=0.000593141;
0103 x5[11]=2, y5[11]=0.00477154, ex[11]=0, ey5[11]=0.000353859;
0104 x5[12]=2.2, y5[12]=0.00309521, ex[12]=0, ey5[12]=0.000265977;
0105 x5[13]=2.4, y5[13]=0.0018804, ex[13]=0, ey5[13]=0.0001962;
0106 x5[14]=2.6, y5[14]=0.00122164, ex[14]=0, ey5[14]=0.000150502;
0107
0108 // anti-protons
0109 float x6[15],y6[15],ey6[15];
0110 x6[0]=0.65, y6[0]=0.306682, ex[0]=0, ey6[0]=0.00840118;
0111 x6[1]=0.75, y6[1]=0.241129, ex[1]=0, ey6[1]=0.00669015;
0112 x6[2]=0.85, y6[2]=0.175801, ex[2]=0, ey6[2]=0.00512706;
0113 x6[3]=0.95, y6[3]=0.124961, ex[3]=0, ey6[3]=0.00389379;
0114 x6[4]=1.05, y6[4]=0.0902856, ex[4]=0, ey6[4]=0.00301377;
0115 x6[5]=1.15, y6[5]=0.0675941, ex[5]=0, ey6[5]=0.00250251;
0116 x6[6]=1.25, y6[6]=0.0451574, ex[6]=0, ey6[6]=0.00186547;
0117 x6[7]=1.35, y6[7]=0.0329953, ex[7]=0, ey6[7]=0.00152812;
0118 x6[8]=1.45, y6[8]=0.0219945, ex[8]=0, ey6[8]=0.00118843;
0119 x6[9]=1.6, y6[9]=0.0151914, ex[9]=0, ey6[9]=0.000644838;
0120 x6[10]=1.8, y6[10]=0.00740469, ex[10]=0, ey6[10]=0.000417221;
0121 x6[11]=2, y6[11]=0.0037345, ex[11]=0, ey6[11]=0.000274817;
0122 x6[12]=2.2, y6[12]=0.00251093, ex[12]=0, ey6[12]=0.000213118;
0123 x6[13]=2.4, y6[13]=0.00116624, ex[13]=0, ey6[13]=0.000138554;
0124 x6[14]=2.6, y6[14]=0.000707305, ex[14]=0, ey6[14]=0.000103754;
0125
0126
0127
0128
0129
0130
0131 // calculate average and divide by cross-section, by 2pi and by 2 units of rapidity, and by pT
0132
0133 float yk[99],eyk[99];
0134 for(int i=0; i<18; i++) { y[i]=(y1[i]+y2[i])/2./40.*2.*3.141592654*2.*x1[i]; ey[i]=sqrt( (ey1[i]*ey1[i])+(ey2[1]*ey2[i]) )/2./40.*2.*3.141592654*2.*x1[i]; }
0135 for(int i=0; i<13; i++) { yk[i]=(y3[i]+y4[i])/2./40.*2.*3.141592654*2.*x3[i]; eyk[i]=sqrt( (ey3[i]*ey3[i])+(ey4[1]*ey4[i]) )/2./40.*2.*3.141592654*2.*x3[i]; }
0136 for(int i=0; i<15; i++) { y5[i]=y5[i]/2./40.*2.*3.141592654*2.*x5[i]; ey5[i]=ey5[i]/2./40.*2.*3.141592654*2.*x5[i]; }
0137 for(int i=0; i<15; i++) { y6[i]=y6[i]/2./40.*2.*3.141592654*2.*x6[i]; ey6[i]=ey6[i]/2./40.*2.*3.141592654*2.*x6[i]; }
0138
0139
0140 gStyle->SetOptStat(0);
0141 TCanvas* c1 = new TCanvas("c1","",20,20,600,600);
0142 c1->SetLogy();
0143
0144 TH1F* hh = new TH1F("hh"," ",20,0.,5.);
0145 hh->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0146 hh->SetYTitle("dN/dP_{T}");
0147 hh->GetYaxis()->SetTitleOffset(1.0);
0148 hh->GetXaxis()->SetTitleOffset(1.0);
0149 hh->SetMinimum(0.00005);
0150 hh->SetMaximum(10.);
0151 hh->Draw();
0152
0153 TGraphErrors* gr1 = new TGraphErrors(18,x1,y,0,ey);
0154 gr1->SetMarkerStyle(20);
0155 gr1->SetMarkerColor(kBlue);
0156 gr1->SetLineColor(kBlue);
0157 gr1->SetMarkerSize(1);
0158
0159 TGraphErrors* gr2 = new TGraphErrors(13,x3,yk,0,eyk);
0160 gr2->SetMarkerStyle(20);
0161 gr2->SetMarkerColor(kMagenta);
0162 gr2->SetLineColor(kMagenta);
0163 gr2->SetMarkerSize(1);
0164
0165 TGraphErrors* gr3 = new TGraphErrors(15,x5,y5,0,ey5);
0166 gr3->SetMarkerStyle(20);
0167 gr3->SetMarkerColor(kGreen+2);
0168 gr3->SetLineColor(kGreen+2);
0169 gr3->SetMarkerSize(1);
0170
0171 TGraphErrors* gr4 = new TGraphErrors(15,x6,y6,0,ey6);
0172 gr4->SetMarkerStyle(20);
0173 gr4->SetMarkerColor(kBlack);
0174 gr4->SetLineColor(kBlack);
0175 gr4->SetMarkerSize(1);
0176
0177
0178
0179 float start = 1.;
0180 float stop = 20.0;
0181 float  ncoll = 955.; // number of binary collisions for 0-10% Au+Au
0182
0183 // Suppression in AuAu
0184 // this is taken from suppressionAuAu.C macro
0185 string str_supppi = "(0.157654+1.947663*pow(x,-2.581256))";
0186 TF1 supppi = TF1("supppi",str_supppi.c_str(),start,stop);
0187 string str_suppp = "(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)";
0188 TF1 suppp = TF1("suppp",str_suppp.c_str(),start,stop);
0189 string str_suppk = str_supppi; // kaon suppression in AuAu same as pion
0190 TF1 suppk = TF1("suppk",str_suppk.c_str(),start,stop);
0191
0192 char srejpi[999],srejk[999],srejp[999],srejap[999];
0193
0194 TF1 frejpi,frejk,frejp,frejantiprot;
0195 double eideff=0.7;
0196 string str_srejpi;
0197 string str_srejk;
0198 string str_srejp;
0199 string str_srejap;
0200
0201 if(eideff==0.9) {
0202 // Hadron inverse rejection factors for eid efficiency 90%
0203
0204 float pirejpar0 = -3.17086e-03;
0205 float pirejpar1 = 2.02062e-01;
0206 float pirejpar2 = -2.41365e+00;
0207 float pirejpar3 = 5.81442e-04;
0208 sprintf(srejpi,"(%f+%f*pow(x,%f)+%f*x)",pirejpar0,pirejpar1,pirejpar2,pirejpar3);
0209 str_srejpi = srejpi;
0210 frejpi = TF1("frejpi",str_srejpi.c_str(),start,stop);
0211
0212 float krejpar0 = -1.78362e-03;
0213 float krejpar1 = 2.77532e-01;
0214 float krejpar2 = -2.72674e+00;
0215 float krejpar3 = 2.19957e-04;
0216 sprintf(srejk,"(%f+%f*pow(x,%f)+%f*x)",krejpar0,krejpar1,krejpar2,krejpar3);
0217 str_srejk = srejk;
0218 frejk = TF1("frejk",str_srejk.c_str(),start,stop);
0219
0220 float prejpar0 = -1.24742e-04;
0221 float prejpar1 = 1.82736e-01;
0222 float prejpar2 = -3.66287e+00;
0223 float prejpar3 = 2.30891e-05;
0224 sprintf(srejp,"(%f+%f*pow(x,%f)+%f*x)",prejpar0,prejpar1,prejpar2,prejpar3);
0225 string str_srejp = srejp;
0226 frejp = TF1("frejp",str_srejp.c_str(),start,stop);
0227
0228 float aprejpar0 = 3.04653e-03;
0229 float aprejpar1 = 7.17479e-01;
0230 float aprejpar2 = -1.26446e+00;
0231 sprintf(srejap,"(%f+%f*exp(x/%f))",aprejpar0,aprejpar1,aprejpar2);
0232 str_srejap = srejap;
0233 frejantiprot = TF1("frejantiprot",str_srejap.c_str(),start,stop);
0234
0235 }
0236 else if (eideff==0.7) {
0237 // Hadron inverse rejection factors for eid efficiency 90%
0238
0239 float pirejpar0_eid70 = -4.71224e-04;
0240 float pirejpar1_eid70 = 1.07319e-01;
0241 float pirejpar2_eid70 = -2.72607e+00;
0242 float pirejpar3_eid70 = 1.51839e-04;
0243 sprintf(srejpi,"(%f+%f*pow(x,%f)+%f*x)",pirejpar0_eid70,pirejpar1_eid70,pirejpar2_eid70,pirejpar3_eid70);
0244 str_srejpi = srejpi;
0245 frejpi = TF1("frejpi",str_srejpi.c_str(),start,stop);
0246
0247 // just divide by 2 for lack of statistics
0248 float krejpar0 = -1.78362e-03;
0249 float krejpar1 = 2.77532e-01;
0250 float krejpar2 = -2.72674e+00;
0251 float krejpar3 = 2.19957e-04;
0252 sprintf(srejk,"(%f+%f*pow(x,%f)+%f*x)/2.",krejpar0,krejpar1,krejpar2,krejpar3);
0253 str_srejk = srejk;
0254 frejk = TF1("frejk",str_srejk.c_str(),start,stop);
0255
0256 // just divide by 2 for lack of statistics
0257 float prejpar0 = -1.24742e-04;
0258 float prejpar1 = 1.82736e-01;
0259 float prejpar2 = -3.66287e+00;
0260 float prejpar3 = 2.30891e-05;
0261 sprintf(srejp,"(%f+%f*pow(x,%f)+%f*x)/2.",prejpar0,prejpar1,prejpar2,prejpar3);
0262 str_srejp = srejp;
0263 frejp = TF1("frejp",str_srejp.c_str(),start,stop);
0264
0265 float aprejpar0_eid70 = 1.08465e-03;
0266 float aprejpar1_eid70 = 5.22870e-01;
0267 float aprejpar2_eid70 = -1.11407e+00;
0268 sprintf(srejap,"(%f+%f*exp(x/%f))",aprejpar0_eid70,aprejpar1_eid70,aprejpar2_eid70);
0269 str_srejap = srejap;
0270 frejantiprot = TF1("frejantiprot",str_srejap.c_str(),start,stop);
0271
0272 } else {
0273   cerr << "ERROR: eid efficiency must be 0.9 or 0.7 !!!" << endl;
0274 }
0275
0276
0277 string multiply = "*";
0278 string divide = "/";
0279 string myplus = "+";
0280
0281 //================ pi-zero ============================================================
0282
0283
0284 // Invariant cross-section of pi-zero in p+p from ppg063
0285 // The units are mb * GeV^-2 * c^3
0286 char invdistr[999];
0287 float a0 = 3.13729e+00;
0288 float a1 = 1.45571e+00;
0289 float a2 = -3.02276e+00;
0290 float a3 = -2.35733e+00;
0291 sprintf(invdistr,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f))",a0,a3,a1,a1,a2);
0292 TF1 invptdistr = TF1("invptdistr",invdistr,start,stop);
0293
0294 // Multiply by pT, divide by p+p total crossection (40mb), and multiply by 2pi
0295 // This should be dN/dpT per event per unit of rapidity
0296 // Then multiply by 2 since we are looking at +-1 unit of rapidity
0297
0298 a0 = a0 * 2.*3.141592654 / 40. * 2.;
0299
0300 char piondistribution[999];
0301 sprintf(piondistribution,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)",a0,a3,a1,a1,a2);
0302 TF1* pionptdistr = new TF1("pionptdistr",piondistribution,start,stop); // this is pion dN/dpT in p+p
0303 string str_pionptdistr = piondistribution;
0304
0305 // convert to AuAu
0306
0307 char piondistributionAuAu[999];
0308 float a0auau = a0 * ncoll; // scale to central AuAu
0309
0310 sprintf(piondistributionAuAu,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)",a0auau,a3,a1,a1,a2);
0311 string str_piondistributionAuAu = piondistributionAuAu;
0312 str_piondistributionAuAu = str_piondistributionAuAu + multiply + str_supppi; // pion suppression in AuAu
0313 TF1* pionptdistrAuAu_norej = new TF1("pionptdistrAuAu_norej",str_piondistributionAuAu.c_str(),start,stop);
0314
0315 string str_piondistributionAuAu_withrej = str_piondistributionAuAu + multiply + str_srejpi; // multiply by inverse pion rejection factor
0316
0317 TF1* pionptdistrAuAu = new TF1("pionptdistrAuAu",str_piondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from pions in AuAu
0318
0319 char piondistributionAuAu_rej90[999];
0320 string str_piondistributionAuAu_rej90 = str_piondistributionAuAu + divide + "90.";
0321 TF1* pionptdistrAuAu_rej90 = new TF1("pionptdistrAuAu_rej90",str_piondistributionAuAu_rej90.c_str(),start,stop);
0322
0323 cout << "CONSTANT REJECTION = 90:" << endl;
0324 cout << str_piondistributionAuAu_rej90 << endl << endl;
0325
0326
0327
0328 //============ kaons ============================================
0329
0330 char kaondistribution[999];
0331 float scalekaon = 0.29;
0332 sprintf(kaondistribution,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x*sqrt(0.1396*0.1396+x*x)/sqrt(0.4937*0.4937+x*x)*%f)",a0,a3,a1,a1,a2,scalekaon);
0333 string str_kaondistribution = kaondistribution;
0334 TF1* kaonptdistr = new TF1("kaonptdistr",str_kaondistribution.c_str(),start,stop);
0335
0336 char kaondistributionAuAu[999];
0337 sprintf(kaondistributionAuAu,"((%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.4937*0.4937+x*x)*%f)",a0auau,a3,a1,a1,a2,scalekaon);
0338 string str_kaondistributionAuAu = kaondistributionAuAu;
0339 str_kaondistributionAuAu = str_kaondistributionAuAu + multiply + str_suppk; // kaon suppression in AuAu
0340 TF1* kaonptdistrAuAu_norej = new TF1("kaonptdistrAuAu_norej",str_kaondistributionAuAu.c_str(),start,stop);
0341
0342 string str_kaondistributionAuAu_withrej = str_kaondistributionAuAu + multiply + str_srejk; // multiply by inverse kaon rejection factor
0343
0344 TF1* kaonptdistrAuAu = new TF1("kaonptdistrAuAu",str_kaondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from Kaons in AuAu
0345
0346
0347 //=========== protons ============================================
0348
0349 char protondistribution[999];
0350 float scalep = 0.14;
0351 sprintf(protondistribution,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*%f)",a0,a3,a1,a1,a2,scalep);
0352 string str_protondistribution = protondistribution;
0353 TF1* protonptdistr = new TF1("protonptdistr",str_protondistribution.c_str(),start,stop);
0354
0355 char protondistributionAuAu[999];
0356 sprintf(protondistributionAuAu,"((%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*%f)",a0auau,a3,a1,a1,a2,scalep);
0357 string str_protondistributionAuAu = protondistributionAuAu;
0358 str_protondistributionAuAu = str_protondistributionAuAu + multiply + str_suppp; // proton suppression in AuAu
0359 TF1* protonptdistrAuAu_norej = new TF1("protonptdistrAuAu_norej",str_protondistributionAuAu.c_str(),start,stop);
0360
0361 string str_protondistributionAuAu_withrej = str_protondistributionAuAu + multiply + str_srejp; // multiply by inverse proton rejection factor
0362
0363 TF1* protonptdistrAuAu = new TF1("protonptdistrAuAu",str_protondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from protons in AuAu
0364
0365
0366 //=========== anti-protons =======================================
0367
0368 char aprotondistribution[999];
0369 float scaleap = 0.09;
0370 sprintf(aprotondistribution,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*%f)",a0,a3,a1,a1,a2,scaleap);
0371 string str_aprotondistribution = aprotondistribution;
0372 TF1* aprotonptdistr = new TF1("aprotonptdistr",str_aprotondistribution.c_str(),start,stop);
0373
0374 char aprotondistributionAuAu[999];
0375 sprintf(aprotondistributionAuAu,"((%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*%f)",a0auau,a3,a1,a1,a2,scaleap);
0376 string str_aprotondistributionAuAu = aprotondistributionAuAu;
0377 str_aprotondistributionAuAu = str_aprotondistributionAuAu + multiply + str_suppp; // proton suppression in AuAu
0378 TF1* aprotonptdistrAuAu_norej = new TF1("aprotonptdistrAuAu_norej",str_aprotondistributionAuAu.c_str(),start,stop);
0379
0380 string str_aprotondistributionAuAu_withrej = str_aprotondistributionAuAu + multiply + str_srejap; // multiply by inverse anti-proton rejection factor
0381
0382 TF1* aprotonptdistrAuAu = new TF1("aprotonptdistrAuAu",str_aprotondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from anti-protons in AuAu
0383
0384
0385
0386 gr1->Draw("P");
0387 gr2->Draw("P");
0388 gr3->Draw("P");
0389 gr4->Draw("P");
0390
0391 pionptdistr->SetLineColor(kBlue);
0392 pionptdistr->Draw("same");
0393 kaonptdistr->SetLineColor(kMagenta);
0394 kaonptdistr->Draw("same");
0395 protonptdistr->SetLineColor(kGreen+2);
0396 protonptdistr->Draw("same");
0397 aprotonptdistr->SetLineColor(kBlack);
0398 aprotonptdistr->Draw("same");
0399
0400
0401
0402
0403 TCanvas* c2 = new TCanvas("c2","",120,120,600,600);
0404 c2->SetLogy();
0405
0406 TH1F* hh2 = new TH1F("hh2"," ",20,0.,5.);
0407 hh2->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0408 hh2->SetYTitle("dN/dP_{T}");
0409 hh2->GetYaxis()->SetTitleOffset(1.0);
0410 hh2->GetXaxis()->SetTitleOffset(1.0);
0411 hh2->SetMinimum(0.00005);
0412 hh2->SetMaximum(500.);
0413 hh2->Draw();
0414
0415 pionptdistrAuAu->SetLineColor(kBlue);
0416 pionptdistrAuAu->Draw();
0417 kaonptdistrAuAu->SetLineColor(kMagenta);
0418 kaonptdistrAuAu->Draw("same");
0419 protonptdistrAuAu->SetLineColor(kGreen+2);
0420 protonptdistrAuAu->Draw("same");
0421 aprotonptdistrAuAu->SetLineColor(kBlack);
0422 aprotonptdistrAuAu->Draw("same");
0423
0424 pionptdistrAuAu_rej90->SetLineColor(kBlue);
0425 pionptdistrAuAu_rej90->SetLineStyle(2);
0426 //pionptdistrAuAu_rej90->Draw("same");
0427
0428 //
0429 // Now add up everything
0430 //
0431
0432 string str_plusdistributionAuAu_withrej = str_piondistributionAuAu_withrej + myplus + str_kaondistributionAuAu_withrej + myplus + str_protondistributionAuAu_withrej;
0433 TF1* plusptdistrAuAu = new TF1("plusptdistrAuAu",str_plusdistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake positrons from all sources in AuAu
0434 plusptdistrAuAu->SetLineColor(kRed);
0435 plusptdistrAuAu->Draw("same");
0436
0437 cout << str_plusdistributionAuAu_withrej << endl;
0438
0439
0440 string str_minusdistributionAuAu_withrej = str_piondistributionAuAu_withrej + myplus + str_kaondistributionAuAu_withrej + myplus + str_aprotondistributionAuAu_withrej;
0441 TF1* minusptdistrAuAu = new TF1("minusptdistrAuAu",str_minusdistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from all sources in AuAu
0442 minusptdistrAuAu->SetLineColor(kRed+2);
0443 minusptdistrAuAu->Draw("same");
0444
0445 cout << str_minusdistributionAuAu_withrej << endl;
0446
0447
0448 /*
0449 TCanvas* c3 = new TCanvas("c3","",220,220,600,600);
0450 c3->SetLogy();
0451 hh3 = new TH1F("hh2"," ",20,0.,5.);
0452 hh3->SetMinimum(0.00005);
0453 hh3->SetMaximum(500.);
0454 hh3->Draw();
0455
0456 pionptdistrAuAu_norej->SetLineColor(kBlue);
0457 pionptdistrAuAu_norej->Draw();
0458 kaonptdistrAuAu_norej->SetLineColor(kMagenta);
0459 kaonptdistrAuAu_norej->Draw("same");
0460 protonptdistrAuAu_norej->SetLineColor(kGreen+2);
0461 protonptdistrAuAu_norej->Draw("same");
0462 aprotonptdistrAuAu_norej->SetLineColor(kBlack);
0463 aprotonptdistrAuAu_norej->Draw("same");
0464 */
0465
0466
0467 }