7 #include "TLorentzVector.h"
32 if(xx<4.8) { val = 1.0e-10*eideff*b0*pow(xx,b3)*pow((1.+(xx/b1)*(xx/b1)),b2); }
33 else { val = 1.0e-10*eideff*bb0*pow(xx,bb3)*pow((1.+(xx/bb1)*(xx/bb1)),bb2); }
48 if(xx<2.5) { val = 1.0e-10*eideff*b0*pow(xx,b3)*pow((1.+(xx/b1)*(xx/b1)),b2); }
49 else { val = 1.0e-10*eideff*bb0*pow(xx,bb3)*pow((1.+(xx/bb1)*(xx/bb1)),bb2); }
57 int main(
int argc,
char* argv[]) {
59 string ofname=
"crossterms.root";
61 if(argc==1) cout << argv[0] <<
" is running with standard parameters..." << endl;
64 return crossterms(eideff, ofname, scalefactor);
68 int crossterms(
double eideff,
string ofname,
int scalefactor) {
84 sprintf(ccharm,
"%f*0.030311*pow(x,2.071907)*pow((1.+(x/2.047743)*(x/2.047743)),-5.699277)",0.7);
85 sprintf(cbottom,
"%f*0.002977*pow(x,0.244158)*pow((1.+(x/3.453862)*(x/3.453862)),-5.024055)",0.7);
87 sprintf(ccharm,
"%f*0.030311*pow(x,2.071907)*pow((1.+(x/2.047743)*(x/2.047743)),-5.699277)",eideff);
88 sprintf(cbottom,
"%f*0.002977*pow(x,0.244158)*pow((1.+(x/3.453862)*(x/3.453862)),-5.024055)",eideff);
90 fcharm =
new TF1(
"fcharm", ccharm, start,stop);
91 fbottom =
new TF1(
"fbottom",cbottom,start,stop);
98 double c0,c1,
c2,c3,cc0,cc1,cc2,cc3;
107 fcharm =
new TF1(
"fcharm",
CFunction,ptcut,stop,9);
108 fcharm->SetParameters(c0,c1,c2,c3,cc0,cc1,cc2,cc3,eideff);
110 double b0,
b1,b2,b3,bb0,bb1,bb2,bb3;
119 fbottom =
new TF1(
"fbottom",
BFunction,ptcut,stop,9);
120 fbottom->SetParameters(b0,b1,b2,b3,bb0,bb1,bb2,bb3,eideff);
127 cout <<
"eID efficiency = " << eideff << endl;
131 fake_str =
"(941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*(0.157654+1.947663*pow(x,-2.581256))*(-0.003171+0.202062*pow(x,-2.413650)+0.000581*x)+((941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.4937*0.4937+x*x)*0.290000)*(0.157654+1.947663*pow(x,-2.581256))*(-0.001784+0.277532*pow(x,-2.726740)+0.000220*x)+((941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*0.090000)*(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)*(0.003047+0.717479*exp(x/-1.264460))";
133 else if(eideff==0.7) {
135 fake_str =
"(941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*(0.157654+1.947663*pow(x,-2.581256))*(-0.000471+0.107319*pow(x,-2.726070)+0.000152*x)+((941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.4937*0.4937+x*x)*0.290000)*(0.157654+1.947663*pow(x,-2.581256))*(-0.001784+0.277532*pow(x,-2.726740)+0.000220*x)/2.+((941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*sqrt(0.1396*0.1396+x*x)/sqrt(0.9383*0.9383+x*x)*0.090000)*(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)*(0.001085+0.522870*exp(x/-1.114070))";
137 else if(eideff==1.0) {
138 cout <<
"Using constant pion rejection factor = 90." << endl;
139 fake_str =
"(941.256287*pow(x,-2.357330)*pow((1.+(x/1.455710)*(x/1.455710)),-3.022760)*x)*(0.157654+1.947663*pow(x,-2.581256))/90.";
141 else {cerr <<
"ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl;}
143 TF1* ffake =
new TF1(
"ffake",fake_str.c_str(),ptcut,stop);
146 double nc = fcharm ->Integral(ptcut,stop);
147 double nb = fbottom->Integral(ptcut,stop);
148 double nfake = ffake ->Integral(ptcut,stop);
149 cout <<
"nfake = " << nfake << endl;
150 cout <<
"nc = " << nc <<
" one event = " << 1./nc/nfake <<
" Au+Au events." << endl;
151 cout <<
"nb = " << nb <<
" one event = " << 1./nb/nfake <<
" Au+Au events." << endl;
152 long int Ncharm = (
long int) (10.0
e+09*nc*nfake);
153 long int Nbottom = (
long int) (10.0
e+09*nb*nfake);
154 cout <<
" Number of events to generate for charm: " << Ncharm << endl;
155 cout <<
" Number of events to generate for bottom: " << Nbottom << endl;
158 TRandom*
rnd =
new TRandom3();
161 TFile*
fout =
new TFile(ofname.c_str(),
"RECREATE");
163 const int nbins = 15;
164 double binlim[nbins+1];
168 TH1D* hhmassc[nbins+1];
169 TH1D* hhmassb[nbins+1];
172 for(
int i=0;
i<nbins+1;
i++) {
173 sprintf(hhname,
"hhmassc_%d",
i);
174 hhmassc[
i] =
new TH1D(hhname,
"",nchan,0.,20.);
175 sprintf(hhname,
"hhmassb_%d",
i);
176 hhmassb[
i] =
new TH1D(hhname,
"",nchan,0.,20.);
177 (hhmassc[
i])->Sumw2();
178 (hhmassb[
i])->Sumw2();
181 double fscalefactor =
double(scalefactor);
185 cout <<
"Generating " << Ncharm*scalefactor <<
" charm/fake events..." << endl;
188 for(
int i=0;
i<Ncharm*scalefactor;
i++) {
190 if(
i%500000==0) cout <<
"processing event # " <<
i << endl;
193 double pt1 = ffake->GetRandom();
194 double eta1 = rnd->Uniform(-1.0,1.0);
195 double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.);
196 double phi1 = rnd->Uniform(-3.141592654,3.141592654);
197 double px1 = pt1*cos(phi1);
198 double py1 = pt1*sin(phi1);
199 double pp1 = pt1/sin(theta1);
200 double pz1 = pp1*cos(theta1);
202 TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
205 double pt2 = fcharm->GetRandom();
206 double eta2 = rnd->Uniform(-1.0,1.0);
207 double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.);
208 double phi2 = rnd->Uniform(-3.141592654,3.141592654);
209 double px2 = pt2*cos(phi2);
210 double py2 = pt2*sin(phi2);
211 double pp2 = pt2/sin(theta2);
212 double pz2 = pp2*cos(theta2);
214 TLorentzVector
vec2 = TLorentzVector(px2, py2, pz2, ee2);
216 if(pt1<ptcut || pt2<ptcut)
continue;
218 TLorentzVector pair = vec1 +
vec2;
219 double invmass = pair.M();
220 double ptpair = pair.Pt();
223 int mybin = int(ptpair);
224 if(mybin>=0 && mybin<nbins) { (hhmassc[mybin])->
Fill(invmass); }
230 cout <<
"Generating " << Nbottom*scalefactor <<
" bottom/fake events..." << endl;
232 for(
int i=0;
i<Nbottom*scalefactor;
i++) {
234 if(
i%500000==0) cout <<
"processing event # " <<
i << endl;
237 double pt1 = ffake->GetRandom();
238 double eta1 = rnd->Uniform(-1.0,1.0);
239 double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.);
240 double phi1 = rnd->Uniform(-3.141592654,3.141592654);
241 double px1 = pt1*cos(phi1);
242 double py1 = pt1*sin(phi1);
243 double pp1 = pt1/sin(theta1);
244 double pz1 = pp1*cos(theta1);
246 TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
249 double pt2 = fbottom->GetRandom();
250 double eta2 = rnd->Uniform(-1.0,1.0);
251 double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.);
252 double phi2 = rnd->Uniform(-3.141592654,3.141592654);
253 double px2 = pt2*cos(phi2);
254 double py2 = pt2*sin(phi2);
255 double pp2 = pt2/sin(theta2);
256 double pz2 = pp2*cos(theta2);
258 TLorentzVector
vec2 = TLorentzVector(px2, py2, pz2, ee2);
260 if(pt1<ptcut || pt2<ptcut)
continue;
262 TLorentzVector pair = vec1 +
vec2;
263 double invmass = pair.M();
264 double ptpair = pair.Pt();
267 int mybin = int(ptpair);
268 if(mybin>=0 && mybin<nbins) { (hhmassb[mybin])->
Fill(invmass); }
273 for(
int i=0;
i<nbins+1;
i++) { (hhmassc[
i])->Scale(1./fscalefactor); (hhmassb[
i])->Scale(1./fscalefactor); }
275 TH1D* hhfakehf[nbins+1];
276 for(
int i=0;
i<nbins+1;
i++) {
277 sprintf(hhname,
"hhfakehf_%d",
i);
278 hhfakehf[
i] = (TH1D*)(hhmassc[
i])->Clone(hhname);
279 (hhfakehf[
i])->Add(hhmassb[i]);