Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
crossterms.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file crossterms.cc
1 #include "TFile.h"
2 #include "TF1.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 //#include "TNtuple.h"
6 #include "TRandom3.h"
7 #include "TLorentzVector.h"
8 
9 #include <cstdlib>
10 #include <cstdio>
11 #include <iostream>
12 #include <iomanip>
13 #include <fstream>
14 #include <vector>
15 
16 using namespace std;
17 
18 int crossterms(double, string, int);
19 
20 double BFunction(double *x, double *p) {
21  double b0 = p[0];
22  double b1 = p[1];
23  double b2 = p[2];
24  double b3 = p[3];
25  double bb0 = p[4];
26  double bb1 = p[5];
27  double bb2 = p[6];
28  double bb3 = p[7];
29  double eideff = p[8];
30  double xx = x[0];
31  double val=0.;
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); }
34  return val;
35 }
36 double CFunction(double *x, double *p) {
37  double b0 = p[0];
38  double b1 = p[1];
39  double b2 = p[2];
40  double b3 = p[3];
41  double bb0 = p[4];
42  double bb1 = p[5];
43  double bb2 = p[6];
44  double bb3 = p[7];
45  double eideff = p[8];
46  double xx = x[0];
47  double val=0.;
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); }
50  return val;
51 }
52 
53 
54 // This program will generate invariant mass distribution from fake
55 // electrons (mis-identified charged hadrons) in 10B 0-10% central AuAu collisions
56 
57 int main(int argc, char* argv[]) {
58  double eideff = 0.7;
59  string ofname="crossterms.root";
60  int scalefactor = 1;
61  if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl;
62  //if(argc==2) {centr = atoi(argv[1]);}
63  //if(argc==3) {centr = atoi(argv[1]); ofname=argv[2];}
64  return crossterms(eideff, ofname, scalefactor);
65 }
66 
67 
68 int crossterms(double eideff, string ofname, int scalefactor) {
69 
70 bool useoldhf=false;
71 TF1* fcharm;
72 TF1* fbottom;
73 double ptcut = 2.0;
74 double start = ptcut;
75 double stop = 20.0;
76 
77 if(useoldhf) {
78 // OLD STUFF
79 // Functions below are dN/dPt per event for Au+Au events
80 // Charm/bottom functions are from ccbar.C and bbbar.C macros
81  char ccharm[999];
82  char cbottom[999];
83  if(eideff==1.0) {
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);
86  } else {
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);
89  }
90  fcharm = new TF1("fcharm", ccharm, start,stop);
91  fbottom = new TF1("fbottom",cbottom,start,stop);
92 }
93 // END OLD STUFF
94 else {
95 
96 // Functions below are dN/dPt per event for Au+Au events
97 // Charm/bottom functions are from fitsingles.C macro
98  double c0,c1,c2,c3,cc0,cc1,cc2,cc3;
99  c0 = 4.27243e+09;
100  c1 = 2.17776;
101  c2 = -6.70255;
102  c3 = -0.700129;
103  cc0 = 553246;
104  cc1 = 0.788563;
105  cc2 = 14.975;
106  cc3 = -36.3373;
107  fcharm = new TF1("fcharm",CFunction,ptcut,stop,9);
108  fcharm->SetParameters(c0,c1,c2,c3,cc0,cc1,cc2,cc3,eideff);
109 
110  double b0,b1,b2,b3,bb0,bb1,bb2,bb3;
111  b0 = 5.98862e+07;
112  b1 = 7.02853;
113  b2 = -14.4721;
114  b3 = -0.215377;
115  bb0 = 5.47498e+09;
116  bb1 = 1.50573;
117  bb2 = 0.168816;
118  bb3 = -6.8832;
119  fbottom = new TF1("fbottom",BFunction,ptcut,stop,9);
120  fbottom->SetParameters(b0,b1,b2,b3,bb0,bb1,bb2,bb3,eideff);
121 }
122 
123 // fake electrons (hadron spectra corrected for rejection factor)
124 // from compare.C macro
125 
126 string fake_str;
127 cout << "eID efficiency = " << eideff << endl;
128 
129 if(eideff==0.9) {
130 // eid efficiency 90%
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))";
132 }
133 else if(eideff==0.7) {
134 // eid efficiency 70%
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))";
136 }
137 else if(eideff==1.0) { // constant pion rejection factor = 90
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.";
140 }
141 else {cerr << "ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl;}
142 
143 TF1* ffake = new TF1("ffake",fake_str.c_str(),ptcut,stop);
144 
145 
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.0e+09*nc*nfake);
153 long int Nbottom = (long int) (10.0e+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;
156 
157 //TRandom2* rnd = new TRandom2(52835623);
158 TRandom* rnd = new TRandom3();
159 rnd->SetSeed(0);
160 
161 TFile* fout = new TFile(ofname.c_str(),"RECREATE");
162 
163 const int nbins = 15;
164 double binlim[nbins+1];
165 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);}
166 
167 // last bin = integrated over pT
168 TH1D* hhmassc[nbins+1];
169 TH1D* hhmassb[nbins+1];
170 char hhname[999];
171 int nchan = 400;
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();
179 }
180 
181 double fscalefactor = double(scalefactor);
182 
183 
184 // Generate events for charm/fake pairs
185 cout << "Generating " << Ncharm*scalefactor << " charm/fake events..." << endl;
186 
187 
188 for(int i=0; i<Ncharm*scalefactor; i++) {
189 
190  if(i%500000==0) cout << "processing event # " << i << endl;
191 
192 // fake electron
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);
201  double ee1 = pp1;
202  TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
203 
204 // charm positron
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);
213  double ee2 = pp2;
214  TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2);
215 
216  if(pt1<ptcut || pt2<ptcut) continue;
217 
218  TLorentzVector pair = vec1 + vec2;
219  double invmass = pair.M();
220  double ptpair = pair.Pt();
221 
222  (hhmassc[nbins])->Fill(invmass);
223  int mybin = int(ptpair);
224  if(mybin>=0 && mybin<nbins) { (hhmassc[mybin])->Fill(invmass); }
225 
226 } // end loop over charm/fake events
227 
228 
229 // Generate events for bottom/fake pairs
230 cout << "Generating " << Nbottom*scalefactor << " bottom/fake events..." << endl;
231 
232 for(int i=0; i<Nbottom*scalefactor; i++) {
233 
234  if(i%500000==0) cout << "processing event # " << i << endl;
235 
236 // fake electron
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);
245  double ee1 = pp1;
246  TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
247 
248 // bottom positron
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);
257  double ee2 = pp2;
258  TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2);
259 
260  if(pt1<ptcut || pt2<ptcut) continue;
261 
262  TLorentzVector pair = vec1 + vec2;
263  double invmass = pair.M();
264  double ptpair = pair.Pt();
265 
266  (hhmassb[nbins])->Fill(invmass);
267  int mybin = int(ptpair);
268  if(mybin>=0 && mybin<nbins) { (hhmassb[mybin])->Fill(invmass); }
269 
270 } // end loop over charm/fake events
271 
272 
273  for(int i=0; i<nbins+1; i++) { (hhmassc[i])->Scale(1./fscalefactor); (hhmassb[i])->Scale(1./fscalefactor); }
274 
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]);
280  }
281 
282  fout->Write();
283  fout->Close();
284 
285  return 0;
286 }
287 
288 
289