8 #include "TMCParticle.h"
9 #include "TLorentzVector.h"
29 int main(
int argc,
char* argv[]) {
31 string ofname=
"fakee.root";
32 long int stat=500000000;
33 if(argc==1) cout << argv[0] <<
" is running with standard parameters..." << endl;
49 string plusdistributionAuAu;
50 string minusdistributionAuAu;
53 cout <<
"eID efficiency = " << eideff << endl;
55 plusdistributionAuAu =
"(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.140000)*(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)*(-0.001784+0.277532*pow(x,-2.726740)+0.000220*x)";
56 minusdistributionAuAu =
"(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))";
59 else if(eideff==0.7) {
60 cout <<
"eID efficiency = " << eideff << endl;
62 plusdistributionAuAu =
"(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.140000)*(1.013817*exp(-(x-2.413264)*(x-2.413264)/2./1.423838/1.423838)+0.157654)*(-0.000125+0.182736*pow(x,-3.662870)+0.000023*x)/2.";
63 minusdistributionAuAu =
"(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))";
65 else if(eideff==1.0) {
66 cout <<
"Using constant pion rejection factor = 90." << endl;
67 plusdistributionAuAu =
"(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.";
68 minusdistributionAuAu =
"(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.";
70 else { cerr <<
"ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl; }
72 TF1* plusptdistrAuAu =
new TF1(
"plusptdistrAuAu",plusdistributionAuAu.c_str(),
start,stop);
73 double norm2plus = plusptdistrAuAu->Integral(start,stop);
74 cout <<
"Fake positron multiplicity above " << start <<
"GeV in Au+Au = " << norm2plus << endl;
76 TF1* minusptdistrAuAu =
new TF1(
"minusptdistrAuAu",minusdistributionAuAu.c_str(),
start,stop);
77 double norm2minus = minusptdistrAuAu->Integral(start,stop);
78 cout <<
"Fake electron multiplicity above " << start <<
"GeV in Au+Au = " << norm2minus << endl;
82 TRandom*
rnd =
new TRandom3();
85 TFile*
fout =
new TFile(ofname.c_str(),
"RECREATE");
99 TH1D* hmult1 =
new TH1D(
"hmult1",
"",20,-0.5,19.5);
100 TH1D* hmult2 =
new TH1D(
"hmult2",
"",20,-0.5,19.5);
101 TH1D* hmultpair =
new TH1D(
"hmultpair",
"",20,-0.5,19.5);
113 const int nbins = 15;
114 double binlim[nbins+1];
119 TH1D* hhmass[nbins+1];
122 for(
int i=0;
i<nbins+1;
i++) {
123 sprintf(hhname,
"hhmass_%d",
i);
124 hhmass[
i] =
new TH1D(hhname,
"",nchan,s1,s2);
127 (hhmass[
i])->Sumw2();
133 for(
int ii=0; ii<nstat; ii++) {
134 for(
int i=0;
i<stat;
i++) {
136 if(
i%10000000==0) cout <<
"processing event # " <<
i << endl;
139 float n1smeared = rnd->Gaus(norm2minus,norm2minus*0.25);
140 float n2smeared = rnd->Gaus(norm2plus,norm2plus*0.25);
141 if(n1smeared<=0. || n2smeared<=0.)
continue;
145 int n1 = rnd->Poisson(n1smeared);
146 int n2 = rnd->Poisson(n2smeared);
147 hmult1->Fill(
float(n1));
148 hmult2->Fill(
float(n2));
149 hmultpair->Fill(
float(n1*n2));
150 if(n1<1 || n2<1)
continue;
153 for(
int i1=0; i1<n1; i1++) {
167 double pt11 = minusptdistrAuAu->GetRandom();
168 double eta11 = rnd->Uniform(-1.0,1.0);
169 double theta11 = 3.141592654/2. + (exp(2.*eta11)-1.)/(exp(2.*eta11)+1.);
170 double phi11 = rnd->Uniform(-3.141592654,3.141592654);
171 double px11 = pt11*cos(phi11);
172 double py11 = pt11*sin(phi11);
173 double pp11 = pt11/sin(theta11);
174 double pz11 = pp11*cos(theta11);
176 TLorentzVector vec11 = TLorentzVector(px11, py11, pz11, ee11);
178 for(
int i2=0; i2<n2; i2++) {
191 double pt22 = plusptdistrAuAu->GetRandom();
192 double eta22 = rnd->Uniform(-1.0,1.0);
193 double theta22 = 3.141592654/2. + (exp(2.*eta22)-1.)/(exp(2.*eta22)+1.);
194 double phi22 = rnd->Uniform(-3.141592654,3.141592654);
195 double px22 = pt22*cos(phi22);
196 double py22 = pt22*sin(phi22);
197 double pp22 = pt22/sin(theta22);
198 double pz22 = pp22*cos(theta22);
200 TLorentzVector vec22 = TLorentzVector(px22, py22, pz22, ee22);
220 if(pt11>ptcut && pt22>ptcut) {
221 TLorentzVector upsilon2 = vec11 + vec22;
227 int mybin = int(upsilon2.Pt());
230 if(mybin>=0 && mybin<nbins) { (hhmass[mybin])->
Fill(upsilon2.M()); }
325 TH1D* hhfakefake[
nbins];
326 for(
int i=0;
i<nbins+1;
i++) {
327 sprintf(hhname,
"hhfakefake_%d",
i);
328 hhfakefake[
i] = (TH1D*)((hhmass[
i])->Clone(hhname));