Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
compare.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file compare.C
1 void compare() {
2 
3 // From ppg30
4 
5 //Differential Invariant Cross Section
6 //Rapididty coverage measured: |Dy| = 0.7
7 //The data have been corrected to |Dy|=1 making a
8 //flat extrapolation from 0.7 to 1.
9 //-----------------------------------------------------------------
10 //E d3sigma/ dp3 = 1/2pi 1/pt d2sigma/dpt dy UNITS [mb/GeV^2]
11 //-----------------------------------------------------------------
12 // If you need to convert this to yield-per-event you must divide by sigma_inel(sqrt_s = 200GeV) = 42 mb
13 
14 float x[18],y[18],ex[18],ey[18];
15 float x1[18],y1[18],ey1[18];
16 float x2[18],y2[18],ey2[18];
17 
18 //piplus
19 x1[0]=0.35, y1[0]=25.6287, ex[0]=0, ey1[0]=0.246673;
20 x1[1]=0.45, y1[1]=14.6531, ex[1]=0, ey1[1]=0.14773;
21 x1[2]=0.55, y1[2]=8.62473, ex[2]=0, ey1[2]=0.0943008;
22 x1[3]=0.65, y1[3]=4.80858, ex[3]=0, ey1[3]=0.0579757;
23 x1[4]=0.75, y1[4]=2.85599, ex[4]=0, ey1[4]=0.03862;
24 x1[5]=0.85, y1[5]=1.71022, ex[5]=0, ey1[5]=0.0258761;
25 x1[6]=0.95, y1[6]=1.08056, ex[6]=0, ey1[6]=0.0183583;
26 x1[7]=1.05, y1[7]=0.68761, ex[7]=0, ey1[7]=0.013222;
27 x1[8]=1.15, y1[8]=0.431937, ex[8]=0, ey1[8]=0.00938769;
28 x1[9]=1.25, y1[9]=0.308209, ex[9]=0, ey1[9]=0.00732418;
29 x1[10]=1.35, y1[10]=0.221919, ex[10]=0, ey1[10]=0.00582217;
30 x1[11]=1.45, y1[11]=0.147961, ex[11]=0, ey1[11]=0.00437927;
31 x1[12]=1.6, y1[12]=0.0825923, ex[12]=0, ey1[12]=0.00205954;
32 x1[13]=1.8, y1[13]=0.0410802, ex[13]=0, ey1[13]=0.00131989;
33 x1[14]=2, y1[14]=0.0211058, ex[14]=0, ey1[14]=0.000876262;
34 x1[15]=2.2, y1[15]=0.0116705, ex[15]=0, ey1[15]=0.000638873;
35 x1[16]=2.4, y1[16]=0.00787318, ex[16]=0, ey1[16]=0.000557839;
36 x1[17]=2.6, y1[17]=0.00447392, ex[17]=0, ey1[17]=0.000458999;
37 
38 // piminus
39 x2[0]=0.35, y2[0]=25.8453, ex[0]=0, ey2[0]=0.21199;
40 x2[1]=0.45, y2[1]=14.5864, ex[1]=0, ey2[1]=0.127527;
41 x2[2]=0.55, y2[2]=8.27867, ex[2]=0, ey2[2]=0.0789266;
42 x2[3]=0.65, y2[3]=4.72545, ex[3]=0, ey2[3]=0.0498388;
43 x2[4]=0.75, y2[4]=2.756, ex[4]=0, ey2[4]=0.032014;
44 x2[5]=0.85, y2[5]=1.66385, ex[5]=0, ey2[5]=0.0215601;
45 x2[6]=0.95, y2[6]=1.05331, ex[6]=0, ey2[6]=0.0152551;
46 x2[7]=1.05, y2[7]=0.648492, ex[7]=0, ey2[7]=0.0105207;
47 x2[8]=1.15, y2[8]=0.431804, ex[8]=0, ey2[8]=0.0079342;
48 x2[9]=1.25, y2[9]=0.286134, ex[9]=0, ey2[9]=0.00585618;
49 x2[10]=1.35, y2[10]=0.199094, ex[10]=0, ey2[10]=0.004624;
50 x2[11]=1.45, y2[11]=0.133006, ex[11]=0, ey2[11]=0.003509;
51 x2[12]=1.6, y2[12]=0.0795068, ex[12]=0, ey2[12]=0.001787;
52 x2[13]=1.8, y2[13]=0.0422071, ex[13]=0, ey2[13]=0.001195;
53 x2[14]=2, y2[14]=0.0208344, ex[14]=0, ey2[14]=0.000781061;
54 x2[15]=2.2, y2[15]=0.0131348, ex[15]=0, ey2[15]=0.000625054;
55 x2[16]=2.4, y2[16]=0.0070146, ex[16]=0, ey2[16]=0.000479789;
56 x2[17]=2.6, y2[17]=0.00425537, ex[17]=0, ey2[17]=0.000408;
57 
58 // Kplus
59 float x3[13],y3[13],ey3[13];
60 x3[0]=0.45, y3[0]=1.95293, ex[0]=0, ey3[0]=0.0724524;
61 x3[1]=0.55, y3[1]=1.40404, ex[1]=0, ey3[1]=0.0466304;
62 x3[2]=0.65, y3[2]=0.925078, ex[2]=0, ey3[2]=0.0296139;
63 x3[3]=0.75, y3[3]=0.624339, ex[3]=0, ey3[3]=0.020472;
64 x3[4]=0.85, y3[4]=0.403909, ex[4]=0, ey3[4]=0.0142644;
65 x3[5]=0.95, y3[5]=0.281567, ex[5]=0, ey3[5]=0.0107128;
66 x3[6]=1.05, y3[6]=0.199669, ex[6]=0, ey3[6]=0.00826777;
67 x3[7]=1.15, y3[7]=0.137853, ex[7]=0, ey3[7]=0.00620796;
68 x3[8]=1.25, y3[8]=0.101317, ex[8]=0, ey3[8]=0.00483626;
69 x3[9]=1.35, y3[9]=0.0694105, ex[9]=0, ey3[9]=0.00362463;
70 x3[10]=1.45, y3[10]=0.0526523, ex[10]=0, ey3[10]=0.00301769;
71 x3[11]=1.6, y3[11]=0.0324873, ex[11]=0, ey3[11]=0.00151623;
72 x3[12]=1.8, y3[12]=0.0179392, ex[12]=0, ey3[12]=0.00100964;
73 
74 // Kminus
75 float x4[13],y4[13],ey4[13];
76 x4[0]=0.45, y4[0]=1.66952, ex[0]=0, ey4[0]=0.0568774;
77 x4[1]=0.55, y4[1]=1.25238, ex[1]=0, ey4[1]=0.0374863;
78 x4[2]=0.65, y4[2]=0.828785, ex[2]=0, ey4[2]=0.0243554;
79 x4[3]=0.75, y4[3]=0.604526, ex[3]=0, ey4[3]=0.0177572;
80 x4[4]=0.85, y4[4]=0.406696, ex[4]=0, ey4[4]=0.0125949;
81 x4[5]=0.95, y4[5]=0.261138, ex[5]=0, ey4[5]=0.00878815;
82 x4[6]=1.05, y4[6]=0.189842, ex[6]=0, ey4[6]=0.00685834;
83 x4[7]=1.15, y4[7]=0.118611, ex[7]=0, ey4[7]=0.00483604;
84 x4[8]=1.25, y4[8]=0.0977945, ex[8]=0, ey4[8]=0.00407271;
85 x4[9]=1.35, y4[9]=0.0673309, ex[9]=0, ey4[9]=0.00315524;
86 x4[10]=1.45, y4[10]=0.0439384, ex[10]=0, ey4[10]=0.00239946;
87 x4[11]=1.6, y4[11]=0.0283855, ex[11]=0, ey4[11]=0.00125486;
88 x4[12]=1.8, y4[12]=0.0181311, ex[12]=0, ey4[12]=0.000911162;
89 
90 // protons
91 float x5[15],y5[15],ey5[15];
92 x5[0]=0.65, y5[0]=0.467596, ex[0]=0, ey5[0]=0.0124957;
93 x5[1]=0.75, y5[1]=0.361103, ex[1]=0, ey5[1]=0.0100168;
94 x5[2]=0.85, y5[2]=0.250201, ex[2]=0, ey5[2]=0.00749516;
95 x5[3]=0.95, y5[3]=0.181417, ex[3]=0, ey5[3]=0.00593336;
96 x5[4]=1.05, y5[4]=0.126166, ex[4]=0, ey5[4]=0.00445335;
97 x5[5]=1.15, y5[5]=0.0902262, ex[5]=0, ey5[5]=0.00351448;
98 x5[6]=1.25, y5[6]=0.0651584, ex[6]=0, ey5[6]=0.00275751;
99 x5[7]=1.35, y5[7]=0.0448966, ex[7]=0, ey5[7]=0.00214243;
100 x5[8]=1.45, y5[8]=0.0328418, ex[8]=0, ey5[8]=0.00172412;
101 x5[9]=1.6, y5[9]=0.0203516, ex[9]=0, ey5[9]=0.000881656;
102 x5[10]=1.8, y5[10]=0.0111527, ex[10]=0, ey5[10]=0.000593141;
103 x5[11]=2, y5[11]=0.00477154, ex[11]=0, ey5[11]=0.000353859;
104 x5[12]=2.2, y5[12]=0.00309521, ex[12]=0, ey5[12]=0.000265977;
105 x5[13]=2.4, y5[13]=0.0018804, ex[13]=0, ey5[13]=0.0001962;
106 x5[14]=2.6, y5[14]=0.00122164, ex[14]=0, ey5[14]=0.000150502;
107 
108 // anti-protons
109 float x6[15],y6[15],ey6[15];
110 x6[0]=0.65, y6[0]=0.306682, ex[0]=0, ey6[0]=0.00840118;
111 x6[1]=0.75, y6[1]=0.241129, ex[1]=0, ey6[1]=0.00669015;
112 x6[2]=0.85, y6[2]=0.175801, ex[2]=0, ey6[2]=0.00512706;
113 x6[3]=0.95, y6[3]=0.124961, ex[3]=0, ey6[3]=0.00389379;
114 x6[4]=1.05, y6[4]=0.0902856, ex[4]=0, ey6[4]=0.00301377;
115 x6[5]=1.15, y6[5]=0.0675941, ex[5]=0, ey6[5]=0.00250251;
116 x6[6]=1.25, y6[6]=0.0451574, ex[6]=0, ey6[6]=0.00186547;
117 x6[7]=1.35, y6[7]=0.0329953, ex[7]=0, ey6[7]=0.00152812;
118 x6[8]=1.45, y6[8]=0.0219945, ex[8]=0, ey6[8]=0.00118843;
119 x6[9]=1.6, y6[9]=0.0151914, ex[9]=0, ey6[9]=0.000644838;
120 x6[10]=1.8, y6[10]=0.00740469, ex[10]=0, ey6[10]=0.000417221;
121 x6[11]=2, y6[11]=0.0037345, ex[11]=0, ey6[11]=0.000274817;
122 x6[12]=2.2, y6[12]=0.00251093, ex[12]=0, ey6[12]=0.000213118;
123 x6[13]=2.4, y6[13]=0.00116624, ex[13]=0, ey6[13]=0.000138554;
124 x6[14]=2.6, y6[14]=0.000707305, ex[14]=0, ey6[14]=0.000103754;
125 
126 
127 
128 
129 
130 
131 // calculate average and divide by cross-section, by 2pi and by 2 units of rapidity, and by pT
132 
133 float yk[99],eyk[99];
134 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]; }
135 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]; }
136 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]; }
137 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]; }
138 
139 
140 gStyle->SetOptStat(0);
141 TCanvas* c1 = new TCanvas("c1","",20,20,600,600);
142 c1->SetLogy();
143 
144 TH1F* hh = new TH1F("hh"," ",20,0.,5.);
145 hh->GetXaxis()->SetTitle("p_{T} (GeV/c)");
146 hh->SetYTitle("dN/dP_{T}");
147 hh->GetYaxis()->SetTitleOffset(1.0);
148 hh->GetXaxis()->SetTitleOffset(1.0);
149 hh->SetMinimum(0.00005);
150 hh->SetMaximum(10.);
151 hh->Draw();
152 
153 TGraphErrors* gr1 = new TGraphErrors(18,x1,y,0,ey);
154 gr1->SetMarkerStyle(20);
155 gr1->SetMarkerColor(kBlue);
156 gr1->SetLineColor(kBlue);
157 gr1->SetMarkerSize(1);
158 
159 TGraphErrors* gr2 = new TGraphErrors(13,x3,yk,0,eyk);
160 gr2->SetMarkerStyle(20);
161 gr2->SetMarkerColor(kMagenta);
162 gr2->SetLineColor(kMagenta);
163 gr2->SetMarkerSize(1);
164 
165 TGraphErrors* gr3 = new TGraphErrors(15,x5,y5,0,ey5);
166 gr3->SetMarkerStyle(20);
167 gr3->SetMarkerColor(kGreen+2);
168 gr3->SetLineColor(kGreen+2);
169 gr3->SetMarkerSize(1);
170 
171 TGraphErrors* gr4 = new TGraphErrors(15,x6,y6,0,ey6);
172 gr4->SetMarkerStyle(20);
173 gr4->SetMarkerColor(kBlack);
174 gr4->SetLineColor(kBlack);
175 gr4->SetMarkerSize(1);
176 
177 
178 
179 float start = 1.;
180 float stop = 20.0;
181 float ncoll = 955.; // number of binary collisions for 0-10% Au+Au
182 
183 // Suppression in AuAu
184 // this is taken from suppressionAuAu.C macro
185 string str_supppi = "(0.157654+1.947663*pow(x,-2.581256))";
186 TF1 supppi = TF1("supppi",str_supppi.c_str(),start,stop);
187 string str_suppp = "(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)";
188 TF1 suppp = TF1("suppp",str_suppp.c_str(),start,stop);
189 string str_suppk = str_supppi; // kaon suppression in AuAu same as pion
190 TF1 suppk = TF1("suppk",str_suppk.c_str(),start,stop);
191 
192 char srejpi[999],srejk[999],srejp[999],srejap[999];
193 
194 TF1 frejpi,frejk,frejp,frejantiprot;
195 double eideff=0.7;
196 string str_srejpi;
197 string str_srejk;
198 string str_srejp;
199 string str_srejap;
200 
201 if(eideff==0.9) {
202 // Hadron inverse rejection factors for eid efficiency 90%
203 
204 float pirejpar0 = -3.17086e-03;
205 float pirejpar1 = 2.02062e-01;
206 float pirejpar2 = -2.41365e+00;
207 float pirejpar3 = 5.81442e-04;
208 sprintf(srejpi,"(%f+%f*pow(x,%f)+%f*x)",pirejpar0,pirejpar1,pirejpar2,pirejpar3);
209 str_srejpi = srejpi;
210 frejpi = TF1("frejpi",str_srejpi.c_str(),start,stop);
211 
212 float krejpar0 = -1.78362e-03;
213 float krejpar1 = 2.77532e-01;
214 float krejpar2 = -2.72674e+00;
215 float krejpar3 = 2.19957e-04;
216 sprintf(srejk,"(%f+%f*pow(x,%f)+%f*x)",krejpar0,krejpar1,krejpar2,krejpar3);
217 str_srejk = srejk;
218 frejk = TF1("frejk",str_srejk.c_str(),start,stop);
219 
220 float prejpar0 = -1.24742e-04;
221 float prejpar1 = 1.82736e-01;
222 float prejpar2 = -3.66287e+00;
223 float prejpar3 = 2.30891e-05;
224 sprintf(srejp,"(%f+%f*pow(x,%f)+%f*x)",prejpar0,prejpar1,prejpar2,prejpar3);
225 string str_srejp = srejp;
226 frejp = TF1("frejp",str_srejp.c_str(),start,stop);
227 
228 float aprejpar0 = 3.04653e-03;
229 float aprejpar1 = 7.17479e-01;
230 float aprejpar2 = -1.26446e+00;
231 sprintf(srejap,"(%f+%f*exp(x/%f))",aprejpar0,aprejpar1,aprejpar2);
232 str_srejap = srejap;
233 frejantiprot = TF1("frejantiprot",str_srejap.c_str(),start,stop);
234 
235 }
236 else if (eideff==0.7) {
237 // Hadron inverse rejection factors for eid efficiency 90%
238 
239 float pirejpar0_eid70 = -4.71224e-04;
240 float pirejpar1_eid70 = 1.07319e-01;
241 float pirejpar2_eid70 = -2.72607e+00;
242 float pirejpar3_eid70 = 1.51839e-04;
243 sprintf(srejpi,"(%f+%f*pow(x,%f)+%f*x)",pirejpar0_eid70,pirejpar1_eid70,pirejpar2_eid70,pirejpar3_eid70);
244 str_srejpi = srejpi;
245 frejpi = TF1("frejpi",str_srejpi.c_str(),start,stop);
246 
247 // just divide by 2 for lack of statistics
248 float krejpar0 = -1.78362e-03;
249 float krejpar1 = 2.77532e-01;
250 float krejpar2 = -2.72674e+00;
251 float krejpar3 = 2.19957e-04;
252 sprintf(srejk,"(%f+%f*pow(x,%f)+%f*x)/2.",krejpar0,krejpar1,krejpar2,krejpar3);
253 str_srejk = srejk;
254 frejk = TF1("frejk",str_srejk.c_str(),start,stop);
255 
256 // just divide by 2 for lack of statistics
257 float prejpar0 = -1.24742e-04;
258 float prejpar1 = 1.82736e-01;
259 float prejpar2 = -3.66287e+00;
260 float prejpar3 = 2.30891e-05;
261 sprintf(srejp,"(%f+%f*pow(x,%f)+%f*x)/2.",prejpar0,prejpar1,prejpar2,prejpar3);
262 str_srejp = srejp;
263 frejp = TF1("frejp",str_srejp.c_str(),start,stop);
264 
265 float aprejpar0_eid70 = 1.08465e-03;
266 float aprejpar1_eid70 = 5.22870e-01;
267 float aprejpar2_eid70 = -1.11407e+00;
268 sprintf(srejap,"(%f+%f*exp(x/%f))",aprejpar0_eid70,aprejpar1_eid70,aprejpar2_eid70);
269 str_srejap = srejap;
270 frejantiprot = TF1("frejantiprot",str_srejap.c_str(),start,stop);
271 
272 } else {
273  cerr << "ERROR: eid efficiency must be 0.9 or 0.7 !!!" << endl;
274 }
275 
276 
277 string multiply = "*";
278 string divide = "/";
279 string myplus = "+";
280 
281 //================ pi-zero ============================================================
282 
283 
284 // Invariant cross-section of pi-zero in p+p from ppg063
285 // The units are mb * GeV^-2 * c^3
286 char invdistr[999];
287 float a0 = 3.13729e+00;
288 float a1 = 1.45571e+00;
289 float a2 = -3.02276e+00;
290 float a3 = -2.35733e+00;
291 sprintf(invdistr,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f))",a0,a3,a1,a1,a2);
292 TF1 invptdistr = TF1("invptdistr",invdistr,start,stop);
293 
294 // Multiply by pT, divide by p+p total crossection (40mb), and multiply by 2pi
295 // This should be dN/dpT per event per unit of rapidity
296 // Then multiply by 2 since we are looking at +-1 unit of rapidity
297 
298 a0 = a0 * 2.*3.141592654 / 40. * 2.;
299 
300 char piondistribution[999];
301 sprintf(piondistribution,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)",a0,a3,a1,a1,a2);
302 TF1* pionptdistr = new TF1("pionptdistr",piondistribution,start,stop); // this is pion dN/dpT in p+p
303 string str_pionptdistr = piondistribution;
304 
305 // convert to AuAu
306 
307 char piondistributionAuAu[999];
308 float a0auau = a0 * ncoll; // scale to central AuAu
309 
310 sprintf(piondistributionAuAu,"(%f*pow(x,%f)*pow((1.+(x/%f)*(x/%f)),%f)*x)",a0auau,a3,a1,a1,a2);
311 string str_piondistributionAuAu = piondistributionAuAu;
312 str_piondistributionAuAu = str_piondistributionAuAu + multiply + str_supppi; // pion suppression in AuAu
313 TF1* pionptdistrAuAu_norej = new TF1("pionptdistrAuAu_norej",str_piondistributionAuAu.c_str(),start,stop);
314 
315 string str_piondistributionAuAu_withrej = str_piondistributionAuAu + multiply + str_srejpi; // multiply by inverse pion rejection factor
316 
317 TF1* pionptdistrAuAu = new TF1("pionptdistrAuAu",str_piondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from pions in AuAu
318 
319 char piondistributionAuAu_rej90[999];
320 string str_piondistributionAuAu_rej90 = str_piondistributionAuAu + divide + "90.";
321 TF1* pionptdistrAuAu_rej90 = new TF1("pionptdistrAuAu_rej90",str_piondistributionAuAu_rej90.c_str(),start,stop);
322 
323 cout << "CONSTANT REJECTION = 90:" << endl;
324 cout << str_piondistributionAuAu_rej90 << endl << endl;
325 
326 
327 
328 //============ kaons ============================================
329 
330 char kaondistribution[999];
331 float scalekaon = 0.29;
332 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);
333 string str_kaondistribution = kaondistribution;
334 TF1* kaonptdistr = new TF1("kaonptdistr",str_kaondistribution.c_str(),start,stop);
335 
336 char kaondistributionAuAu[999];
337 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);
338 string str_kaondistributionAuAu = kaondistributionAuAu;
339 str_kaondistributionAuAu = str_kaondistributionAuAu + multiply + str_suppk; // kaon suppression in AuAu
340 TF1* kaonptdistrAuAu_norej = new TF1("kaonptdistrAuAu_norej",str_kaondistributionAuAu.c_str(),start,stop);
341 
342 string str_kaondistributionAuAu_withrej = str_kaondistributionAuAu + multiply + str_srejk; // multiply by inverse kaon rejection factor
343 
344 TF1* kaonptdistrAuAu = new TF1("kaonptdistrAuAu",str_kaondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from Kaons in AuAu
345 
346 
347 //=========== protons ============================================
348 
349 char protondistribution[999];
350 float scalep = 0.14;
351 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);
352 string str_protondistribution = protondistribution;
353 TF1* protonptdistr = new TF1("protonptdistr",str_protondistribution.c_str(),start,stop);
354 
355 char protondistributionAuAu[999];
356 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);
357 string str_protondistributionAuAu = protondistributionAuAu;
358 str_protondistributionAuAu = str_protondistributionAuAu + multiply + str_suppp; // proton suppression in AuAu
359 TF1* protonptdistrAuAu_norej = new TF1("protonptdistrAuAu_norej",str_protondistributionAuAu.c_str(),start,stop);
360 
361 string str_protondistributionAuAu_withrej = str_protondistributionAuAu + multiply + str_srejp; // multiply by inverse proton rejection factor
362 
363 TF1* protonptdistrAuAu = new TF1("protonptdistrAuAu",str_protondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from protons in AuAu
364 
365 
366 //=========== anti-protons =======================================
367 
368 char aprotondistribution[999];
369 float scaleap = 0.09;
370 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);
371 string str_aprotondistribution = aprotondistribution;
372 TF1* aprotonptdistr = new TF1("aprotonptdistr",str_aprotondistribution.c_str(),start,stop);
373 
374 char aprotondistributionAuAu[999];
375 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);
376 string str_aprotondistributionAuAu = aprotondistributionAuAu;
377 str_aprotondistributionAuAu = str_aprotondistributionAuAu + multiply + str_suppp; // proton suppression in AuAu
378 TF1* aprotonptdistrAuAu_norej = new TF1("aprotonptdistrAuAu_norej",str_aprotondistributionAuAu.c_str(),start,stop);
379 
380 string str_aprotondistributionAuAu_withrej = str_aprotondistributionAuAu + multiply + str_srejap; // multiply by inverse anti-proton rejection factor
381 
382 TF1* aprotonptdistrAuAu = new TF1("aprotonptdistrAuAu",str_aprotondistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from anti-protons in AuAu
383 
384 
385 
386 gr1->Draw("P");
387 gr2->Draw("P");
388 gr3->Draw("P");
389 gr4->Draw("P");
390 
391 pionptdistr->SetLineColor(kBlue);
392 pionptdistr->Draw("same");
393 kaonptdistr->SetLineColor(kMagenta);
394 kaonptdistr->Draw("same");
395 protonptdistr->SetLineColor(kGreen+2);
396 protonptdistr->Draw("same");
397 aprotonptdistr->SetLineColor(kBlack);
398 aprotonptdistr->Draw("same");
399 
400 
401 
402 
403 TCanvas* c2 = new TCanvas("c2","",120,120,600,600);
404 c2->SetLogy();
405 
406 TH1F* hh2 = new TH1F("hh2"," ",20,0.,5.);
407 hh2->GetXaxis()->SetTitle("p_{T} (GeV/c)");
408 hh2->SetYTitle("dN/dP_{T}");
409 hh2->GetYaxis()->SetTitleOffset(1.0);
410 hh2->GetXaxis()->SetTitleOffset(1.0);
411 hh2->SetMinimum(0.00005);
412 hh2->SetMaximum(500.);
413 hh2->Draw();
414 
415 pionptdistrAuAu->SetLineColor(kBlue);
416 pionptdistrAuAu->Draw();
417 kaonptdistrAuAu->SetLineColor(kMagenta);
418 kaonptdistrAuAu->Draw("same");
419 protonptdistrAuAu->SetLineColor(kGreen+2);
420 protonptdistrAuAu->Draw("same");
421 aprotonptdistrAuAu->SetLineColor(kBlack);
422 aprotonptdistrAuAu->Draw("same");
423 
424 pionptdistrAuAu_rej90->SetLineColor(kBlue);
425 pionptdistrAuAu_rej90->SetLineStyle(2);
426 //pionptdistrAuAu_rej90->Draw("same");
427 
428 //
429 // Now add up everything
430 //
431 
432 string str_plusdistributionAuAu_withrej = str_piondistributionAuAu_withrej + myplus + str_kaondistributionAuAu_withrej + myplus + str_protondistributionAuAu_withrej;
433 TF1* plusptdistrAuAu = new TF1("plusptdistrAuAu",str_plusdistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake positrons from all sources in AuAu
434 plusptdistrAuAu->SetLineColor(kRed);
435 plusptdistrAuAu->Draw("same");
436 
437 cout << str_plusdistributionAuAu_withrej << endl;
438 
439 
440 string str_minusdistributionAuAu_withrej = str_piondistributionAuAu_withrej + myplus + str_kaondistributionAuAu_withrej + myplus + str_aprotondistributionAuAu_withrej;
441 TF1* minusptdistrAuAu = new TF1("minusptdistrAuAu",str_minusdistributionAuAu_withrej.c_str(),start,stop); // this is dN/dpT for fake electrons from all sources in AuAu
442 minusptdistrAuAu->SetLineColor(kRed+2);
443 minusptdistrAuAu->Draw("same");
444 
445 cout << str_minusdistributionAuAu_withrej << endl;
446 
447 
448 /*
449 TCanvas* c3 = new TCanvas("c3","",220,220,600,600);
450 c3->SetLogy();
451 hh3 = new TH1F("hh2"," ",20,0.,5.);
452 hh3->SetMinimum(0.00005);
453 hh3->SetMaximum(500.);
454 hh3->Draw();
455 
456 pionptdistrAuAu_norej->SetLineColor(kBlue);
457 pionptdistrAuAu_norej->Draw();
458 kaonptdistrAuAu_norej->SetLineColor(kMagenta);
459 kaonptdistrAuAu_norej->Draw("same");
460 protonptdistrAuAu_norej->SetLineColor(kGreen+2);
461 protonptdistrAuAu_norej->Draw("same");
462 aprotonptdistrAuAu_norej->SetLineColor(kBlack);
463 aprotonptdistrAuAu_norej->Draw("same");
464 */
465 
466 
467 }