QuarkoniaBG/compiled/crossterms.cc

0001 #include "TFile.h"
0002 #include "TF1.h"
0003 #include "TH1.h"
0004 #include "TH2.h"
0005 //#include "TNtuple.h"
0006 #include "TRandom3.h"
0007 #include "TLorentzVector.h"
0008
0009 #include <cstdlib>
0010 #include <cstdio>
0011 #include <iostream>
0012 #include <iomanip>
0013 #include <fstream>
0014 #include <vector>
0015
0016 using namespace std;
0017
0018 int crossterms(double, string, int);
0019
0020 double BFunction(double *x, double *p) {
0021   double b0 = p[0];
0022   double b1 = p[1];
0023   double b2 = p[2];
0024   double b3 = p[3];
0025   double bb0 = p[4];
0026   double bb1 = p[5];
0027   double bb2 = p[6];
0028   double bb3 = p[7];
0029   double eideff = p[8];
0030   double xx = x[0];
0031   double val=0.;
0032   if(xx<4.8) { val = 1.0e-10*eideff*b0*pow(xx,b3)*pow((1.+(xx/b1)*(xx/b1)),b2); }
0033     else { val     = 1.0e-10*eideff*bb0*pow(xx,bb3)*pow((1.+(xx/bb1)*(xx/bb1)),bb2); }
0034   return val;
0035 }
0036 double CFunction(double *x, double *p) {
0037   double b0 = p[0];
0038   double b1 = p[1];
0039   double b2 = p[2];
0040   double b3 = p[3];
0041   double bb0 = p[4];
0042   double bb1 = p[5];
0043   double bb2 = p[6];
0044   double bb3 = p[7];
0045   double eideff = p[8];
0046   double xx = x[0];
0047   double val=0.;
0048   if(xx<2.5) { val = 1.0e-10*eideff*b0*pow(xx,b3)*pow((1.+(xx/b1)*(xx/b1)),b2); }
0049     else { val     = 1.0e-10*eideff*bb0*pow(xx,bb3)*pow((1.+(xx/bb1)*(xx/bb1)),bb2); }
0050   return val;
0051 }
0052
0053
0054 // This program will generate invariant mass distribution from fake
0055 // electrons (mis-identified charged hadrons) in 10B 0-10% central AuAu collisions
0056
0057 int main(int argc, char* argv[]) {
0058   double eideff = 0.7;
0059   string ofname="crossterms.root";
0060   int scalefactor  = 1;
0061   if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl;
0062   //if(argc==2) {centr = atoi(argv[1]);}
0063   //if(argc==3) {centr = atoi(argv[1]); ofname=argv[2];}
0064   return crossterms(eideff, ofname, scalefactor);
0065 }
0066
0067
0068 int crossterms(double eideff, string ofname, int scalefactor) {
0069
0070 bool useoldhf=false;
0071 TF1* fcharm;
0072 TF1* fbottom;
0073 double ptcut = 2.0;
0074 double start = ptcut;
0075 double stop = 20.0;
0076
0077 if(useoldhf) {
0078 // OLD STUFF
0079 // Functions below are dN/dPt per event for Au+Au events
0080 // Charm/bottom functions are from ccbar.C and bbbar.C macros
0081   char ccharm[999];
0082   char cbottom[999];
0083   if(eideff==1.0) {
0084     sprintf(ccharm,"%f*0.030311*pow(x,2.071907)*pow((1.+(x/2.047743)*(x/2.047743)),-5.699277)",0.7);
0085     sprintf(cbottom,"%f*0.002977*pow(x,0.244158)*pow((1.+(x/3.453862)*(x/3.453862)),-5.024055)",0.7);
0086   } else {
0087     sprintf(ccharm,"%f*0.030311*pow(x,2.071907)*pow((1.+(x/2.047743)*(x/2.047743)),-5.699277)",eideff);
0088     sprintf(cbottom,"%f*0.002977*pow(x,0.244158)*pow((1.+(x/3.453862)*(x/3.453862)),-5.024055)",eideff);
0089   }
0090   fcharm =  new TF1("fcharm", ccharm, start,stop);
0091   fbottom = new TF1("fbottom",cbottom,start,stop);
0092 }
0093 // END OLD STUFF
0094 else {
0095
0096 // Functions below are dN/dPt per event for Au+Au events
0097 // Charm/bottom functions are from fitsingles.C macro
0098   double c0,c1,c2,c3,cc0,cc1,cc2,cc3;
0099   c0 = 4.27243e+09;
0100   c1 = 2.17776;
0101   c2 = -6.70255;
0102   c3 = -0.700129;
0103   cc0 = 553246;
0104   cc1 = 0.788563;
0105   cc2 = 14.975;
0106   cc3 = -36.3373;
0107   fcharm = new TF1("fcharm",CFunction,ptcut,stop,9);
0108   fcharm->SetParameters(c0,c1,c2,c3,cc0,cc1,cc2,cc3,eideff);
0109
0110   double b0,b1,b2,b3,bb0,bb1,bb2,bb3;
0111   b0 = 5.98862e+07;
0112   b1 = 7.02853;
0113   b2 = -14.4721;
0114   b3 = -0.215377;
0115   bb0 = 5.47498e+09;
0116   bb1 = 1.50573;
0117   bb2 = 0.168816;
0118   bb3 = -6.8832;
0119   fbottom = new TF1("fbottom",BFunction,ptcut,stop,9);
0120   fbottom->SetParameters(b0,b1,b2,b3,bb0,bb1,bb2,bb3,eideff);
0121 }
0122
0123 // fake electrons (hadron spectra corrected for rejection factor)
0124 // from compare.C macro
0125
0126 string fake_str;
0127 cout << "eID efficiency = " << eideff << endl;
0128
0129 if(eideff==0.9) {
0130 // eid efficiency 90%
0131 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))";
0132 }
0133 else if(eideff==0.7) {
0134 // eid efficiency 70%
0135 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))";
0136 }
0137 else if(eideff==1.0) { // constant pion rejection factor = 90
0138 cout << "Using constant pion rejection factor = 90." << endl;
0139 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.";
0140 }
0141 else {cerr << "ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl;}
0142
0143 TF1* ffake =  new TF1("ffake",fake_str.c_str(),ptcut,stop);
0144
0145
0146 double nc    = fcharm ->Integral(ptcut,stop);
0147 double nb    = fbottom->Integral(ptcut,stop);
0148 double nfake = ffake  ->Integral(ptcut,stop);
0149 cout << "nfake = " << nfake << endl;
0150 cout << "nc = " << nc << "     one event = " << 1./nc/nfake << " Au+Au events." << endl;
0151 cout << "nb = " << nb << "     one event = " << 1./nb/nfake << " Au+Au events." << endl;
0152 long int Ncharm = (long int) (10.0e+09*nc*nfake);
0153 long int Nbottom = (long int) (10.0e+09*nb*nfake);
0154 cout << " Number of events to generate for charm: " << Ncharm << endl;
0155 cout << " Number of events to generate for bottom: " << Nbottom << endl;
0156
0157 //TRandom2* rnd = new TRandom2(52835623);
0158 TRandom* rnd = new TRandom3();
0159 rnd->SetSeed(0);
0160
0161 TFile* fout = new TFile(ofname.c_str(),"RECREATE");
0162
0163 const int nbins = 15;
0164 double binlim[nbins+1];
0165 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);}
0166
0167 // last bin = integrated over pT
0168 TH1D* hhmassc[nbins+1];
0169 TH1D* hhmassb[nbins+1];
0170 char hhname[999];
0171 int nchan = 400;
0172 for(int i=0; i<nbins+1; i++) {
0173   sprintf(hhname,"hhmassc_%d",i);
0174   hhmassc[i] = new TH1D(hhname,"",nchan,0.,20.);
0175   sprintf(hhname,"hhmassb_%d",i);
0176   hhmassb[i] = new TH1D(hhname,"",nchan,0.,20.);
0177   (hhmassc[i])->Sumw2();
0178   (hhmassb[i])->Sumw2();
0179 }
0180
0181 double fscalefactor = double(scalefactor);
0182
0183
0184 // Generate events for charm/fake pairs
0185 cout << "Generating " << Ncharm*scalefactor << " charm/fake events..." << endl;
0186
0187
0188 for(int i=0; i<Ncharm*scalefactor; i++) {
0189
0190   if(i%500000==0) cout << "processing event # " << i << endl;
0191
0192 // fake electron
0193   double pt1 = ffake->GetRandom();
0194   double eta1 = rnd->Uniform(-1.0,1.0);
0195   double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.);
0196   double phi1 = rnd->Uniform(-3.141592654,3.141592654);
0197   double px1 = pt1*cos(phi1);
0198   double py1 = pt1*sin(phi1);
0199   double pp1 = pt1/sin(theta1);
0200   double pz1 = pp1*cos(theta1);
0201   double ee1 = pp1;
0202   TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
0203
0204 // charm positron
0205   double pt2 = fcharm->GetRandom();
0206   double eta2 = rnd->Uniform(-1.0,1.0);
0207   double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.);
0208   double phi2 = rnd->Uniform(-3.141592654,3.141592654);
0209   double px2 = pt2*cos(phi2);
0210   double py2 = pt2*sin(phi2);
0211   double pp2 = pt2/sin(theta2);
0212   double pz2 = pp2*cos(theta2);
0213   double ee2 = pp2;
0214   TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2);
0215
0216   if(pt1<ptcut || pt2<ptcut) continue;
0217
0218     TLorentzVector pair = vec1 + vec2;
0219     double invmass = pair.M();
0220     double ptpair = pair.Pt();
0221
0222     (hhmassc[nbins])->Fill(invmass);
0223       int mybin = int(ptpair);
0224       if(mybin>=0 && mybin<nbins) { (hhmassc[mybin])->Fill(invmass); }
0225
0226 } // end loop over charm/fake events
0227
0228
0229 // Generate events for bottom/fake pairs
0230 cout << "Generating " << Nbottom*scalefactor << " bottom/fake events..." << endl;
0231
0232 for(int i=0; i<Nbottom*scalefactor; i++) {
0233
0234   if(i%500000==0) cout << "processing event # " << i << endl;
0235
0236 // fake electron
0237   double pt1 = ffake->GetRandom();
0238   double eta1 = rnd->Uniform(-1.0,1.0);
0239   double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.);
0240   double phi1 = rnd->Uniform(-3.141592654,3.141592654);
0241   double px1 = pt1*cos(phi1);
0242   double py1 = pt1*sin(phi1);
0243   double pp1 = pt1/sin(theta1);
0244   double pz1 = pp1*cos(theta1);
0245   double ee1 = pp1;
0246   TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
0247
0248 // bottom positron
0249   double pt2 = fbottom->GetRandom();
0250   double eta2 = rnd->Uniform(-1.0,1.0);
0251   double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.);
0252   double phi2 = rnd->Uniform(-3.141592654,3.141592654);
0253   double px2 = pt2*cos(phi2);
0254   double py2 = pt2*sin(phi2);
0255   double pp2 = pt2/sin(theta2);
0256   double pz2 = pp2*cos(theta2);
0257   double ee2 = pp2;
0258   TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2);
0259
0260   if(pt1<ptcut || pt2<ptcut) continue;
0261
0262     TLorentzVector pair = vec1 + vec2;
0263     double invmass = pair.M();
0264     double ptpair = pair.Pt();
0265
0266     (hhmassb[nbins])->Fill(invmass);
0267       int mybin = int(ptpair);
0268       if(mybin>=0 && mybin<nbins) { (hhmassb[mybin])->Fill(invmass); }
0269
0270 } // end loop over charm/fake events
0271
0272
0273   for(int i=0; i<nbins+1; i++) { (hhmassc[i])->Scale(1./fscalefactor); (hhmassb[i])->Scale(1./fscalefactor); }
0274
0275   TH1D* hhfakehf[nbins+1];
0276   for(int i=0; i<nbins+1; i++) {
0277     sprintf(hhname,"hhfakehf_%d",i);
0278     hhfakehf[i]  = (TH1D*)(hhmassc[i])->Clone(hhname);
0279     (hhfakehf[i])->Add(hhmassb[i]);
0280   }
0281
0282   fout->Write();
0283   fout->Close();
0284
0285   return 0;
0286 }
0287
0288
0289