sPhenix code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-08-05 08:14:51
0001 #include "TFile.h" 0002 #include "TF1.h" 0003 #include "TH1.h" 0004 #include "TH2.h" 0005 #include "TNtuple.h" 0006 #include "TTree.h" 0007 #include "TRandom3.h" 0008 #include "TMCParticle.h" 0009 #include "TLorentzVector.h" 0010 0011 #include <cstdlib> 0012 #include <cstdio> 0013 #include <iostream> 0014 #include <iomanip> 0015 #include <fstream> 0016 #include <vector> 0017 0018 using namespace std; 0019 0020 int fakee_invmass(double, string, long int); 0021 0022 // This program will generate invariant mass distribution from fake 0023 // electrons (mis-identified charged hadrons) in central AuAu collisions 0024 // and scales it up to 5B central events (50B min. bias) 0025 // You can choose 0-10% (centr=0) or 10-20% (centr=1) central events 0026 // You can choose three different pion rejection functions, 0027 // or use constant (Pt-independent) pion rejection factor 0028 0029 int main(int argc, char* argv[]) { 0030 double eideff = 0.7; 0031 string ofname="fakee.root"; 0032 long int stat=500000000; 0033 if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl; 0034 // if(argc==2) {eideff = atoi(argv[1]);} 0035 // if(argc==3) {eideff = atoi(argv[1]); ofname=argv[2]; stat=atoi(argv[3]);} 0036 return fakee_invmass(eideff, ofname, stat); 0037 } 0038 0039 0040 int fakee_invmass(double eideff, string ofname, long int stat) { 0041 0042 int nstat = 20; 0043 0044 // This is dN/dpT for fake electrons from all sources (pions, kaons,protons, anti-protons) 0045 //float start = 1.; 0046 double start = 2.; 0047 double stop = 20.0; 0048 double ptcut=2.0; 0049 string plusdistributionAuAu; 0050 string minusdistributionAuAu; 0051 0052 if(eideff==0.9) { 0053 cout << "eID efficiency = " << eideff << endl; 0054 // eid efficiency 90% 0055 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)"; 0056 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))"; 0057 0058 } 0059 else if(eideff==0.7) { 0060 cout << "eID efficiency = " << eideff << endl; 0061 // eid efficnecy 70% 0062 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."; 0063 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))"; 0064 } 0065 else if(eideff==1.0) { // constant pion rejection factor = 90 0066 cout << "Using constant pion rejection factor = 90." << endl; 0067 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."; 0068 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."; 0069 } 0070 else { cerr << "ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl; } 0071 0072 TF1* plusptdistrAuAu = new TF1("plusptdistrAuAu",plusdistributionAuAu.c_str(),start,stop); 0073 double norm2plus = plusptdistrAuAu->Integral(start,stop); 0074 cout << "Fake positron multiplicity above " << start << "GeV in Au+Au = " << norm2plus << endl; 0075 0076 TF1* minusptdistrAuAu = new TF1("minusptdistrAuAu",minusdistributionAuAu.c_str(),start,stop); 0077 double norm2minus = minusptdistrAuAu->Integral(start,stop); 0078 cout << "Fake electron multiplicity above " << start << "GeV in Au+Au = " << norm2minus << endl; 0079 0080 0081 0082 TRandom* rnd = new TRandom3(); 0083 rnd->SetSeed(0); 0084 0085 TFile* fout = new TFile(ofname.c_str(),"RECREATE"); 0086 int nchan = 400; 0087 double s1=0.; 0088 double s2=20.; 0089 //TH1D* hmass = new TH1D("hmass","", nchan,s1,s2); 0090 //TH1D* hmass0 = new TH1D("hmass0","",nchan,s1,s2); 0091 //TH1D* hmass1 = new TH1D("hmass1","",nchan,s1,s2); 0092 //TH1D* hmass2 = new TH1D("hmass2","",nchan,s1,s2); 0093 //TH1D* hmass3 = new TH1D("hmass3","",nchan,s1,s2); 0094 //TH1D* hmass_flatpt = new TH1D("hmass_flatpt","", nchan,s1,s2); 0095 //TH1D* hmass0_flatpt = new TH1D("hmass0_flatpt","",nchan,s1,s2); 0096 //TH1D* hmass1_flatpt = new TH1D("hmass1_flatpt","",nchan,s1,s2); 0097 //TH1D* hmass2_flatpt = new TH1D("hmass2_flatpt","",nchan,s1,s2); 0098 //TH1D* hmass3_flatpt = new TH1D("hmass3_flatpt","",nchan,s1,s2); 0099 TH1D* hmult1 = new TH1D("hmult1","",20,-0.5,19.5); 0100 TH1D* hmult2 = new TH1D("hmult2","",20,-0.5,19.5); 0101 TH1D* hmultpair = new TH1D("hmultpair","",20,-0.5,19.5); 0102 //hmass->Sumw2(); 0103 //hmass0->Sumw2(); 0104 //hmass1->Sumw2(); 0105 //hmass2->Sumw2(); 0106 //hmass3->Sumw2(); 0107 //hmass_flatpt->Sumw2(); 0108 //hmass0_flatpt->Sumw2(); 0109 //hmass1_flatpt->Sumw2(); 0110 //hmass2_flatpt->Sumw2(); 0111 //hmass3_flatpt->Sumw2(); 0112 0113 const int nbins = 15; 0114 double binlim[nbins+1]; 0115 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);} 0116 //binlim[9] = 10.; 0117 //binlim[10] = 15.; 0118 0119 TH1D* hhmass[nbins+1]; 0120 //TH1D* hhmass_flatpt[nbins+1]; 0121 char hhname[999]; 0122 for(int i=0; i<nbins+1; i++) { 0123 sprintf(hhname,"hhmass_%d",i); 0124 hhmass[i] = new TH1D(hhname,"",nchan,s1,s2); 0125 // sprintf(hhname,"hhmass_flatpt_%d",i); 0126 // hhmass_flatpt[i] = new TH1D(hhname,"",nchan,s1,s2); 0127 (hhmass[i])->Sumw2(); 0128 // (hhmass_flatpt[i])->Sumw2(); 0129 } 0130 0131 // Generate events 0132 0133 for(int ii=0; ii<nstat; ii++) { 0134 for(int i=0; i<stat; i++) { 0135 0136 if(i%10000000==0) cout << "processing event # " << i << endl; 0137 0138 // first introduce 25% fluctuations in multiplicity 0139 float n1smeared = rnd->Gaus(norm2minus,norm2minus*0.25); 0140 float n2smeared = rnd->Gaus(norm2plus,norm2plus*0.25); 0141 if(n1smeared<=0. || n2smeared<=0.) continue; 0142 // float n1smeared = norm2minus; 0143 // float n2smeared = norm2plus; 0144 // then calculate multiplicity for fake electrons and positrons 0145 int n1 = rnd->Poisson(n1smeared); 0146 int n2 = rnd->Poisson(n2smeared); 0147 hmult1->Fill(float(n1)); 0148 hmult2->Fill(float(n2)); 0149 hmultpair->Fill(float(n1*n2)); 0150 if(n1<1 || n2<1) continue; 0151 0152 // loop over pairs 0153 for(int i1=0; i1<n1; i1++) { 0154 0155 /* 0156 double pt1 = rnd->Uniform(start,stop); // flat pT 0157 double eta1 = rnd->Uniform(-1.0,1.0); 0158 double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.); 0159 double phi1 = rnd->Uniform(-3.141592654,3.141592654); 0160 double px1 = pt1*cos(phi1); 0161 double py1 = pt1*sin(phi1); 0162 double pp1 = pt1/sin(theta1); 0163 double pz1 = pp1*cos(theta1); 0164 double ee1 = pp1; 0165 TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1); 0166 */ 0167 double pt11 = minusptdistrAuAu->GetRandom(); // correct pt used for cross-check of histogram normalization 0168 double eta11 = rnd->Uniform(-1.0,1.0); 0169 double theta11 = 3.141592654/2. + (exp(2.*eta11)-1.)/(exp(2.*eta11)+1.); 0170 double phi11 = rnd->Uniform(-3.141592654,3.141592654); 0171 double px11 = pt11*cos(phi11); 0172 double py11 = pt11*sin(phi11); 0173 double pp11 = pt11/sin(theta11); 0174 double pz11 = pp11*cos(theta11); 0175 double ee11 = pp11; 0176 TLorentzVector vec11 = TLorentzVector(px11, py11, pz11, ee11); 0177 0178 for(int i2=0; i2<n2; i2++) { 0179 /* 0180 double pt2 = rnd->Uniform(start,stop); // flat pT 0181 double eta2 = rnd->Uniform(-1.0,1.0); 0182 double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.); 0183 double phi2 = rnd->Uniform(-3.141592654,3.141592654); 0184 double px2 = pt2*cos(phi2); 0185 double py2 = pt2*sin(phi2); 0186 double pp2 = pt2/sin(theta2); 0187 double pz2 = pp2*cos(theta2); 0188 double ee2 = pp2; 0189 TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2); 0190 */ 0191 double pt22 = plusptdistrAuAu->GetRandom(); // correct pt used for cross-check of histogram normalization 0192 double eta22 = rnd->Uniform(-1.0,1.0); 0193 double theta22 = 3.141592654/2. + (exp(2.*eta22)-1.)/(exp(2.*eta22)+1.); 0194 double phi22 = rnd->Uniform(-3.141592654,3.141592654); 0195 double px22 = pt22*cos(phi22); 0196 double py22 = pt22*sin(phi22); 0197 double pp22 = pt22/sin(theta22); 0198 double pz22 = pp22*cos(theta22); 0199 double ee22 = pp22; 0200 TLorentzVector vec22 = TLorentzVector(px22, py22, pz22, ee22); 0201 0202 /* 0203 // flat pT 0204 if(pt1>ptcut && pt2>ptcut) { 0205 TLorentzVector upsilon = vec1 + vec2; // flat pT 0206 double weight1 = minusptdistrAuAu->Eval(pt1)/norm2minus; // weight normalized to 1. 0207 double weight2 = plusptdistrAuAu->Eval(pt2)/norm2plus; 0208 (hhmass_flatpt[nbins])->Fill(upsilon.M(), weight1*weight2); 0209 //if(upsilon.Pt()>0.0 && upsilon.Pt()<=2.0) { hmass0_flatpt->Fill(upsilon.M(), weight1*weight2); } 0210 //if(upsilon.Pt()>2.0 && upsilon.Pt()<=4.0) { hmass1_flatpt->Fill(upsilon.M(), weight1*weight2); } 0211 //if(upsilon.Pt()>4.0 && upsilon.Pt()<=6.0) { hmass2_flatpt->Fill(upsilon.M(), weight1*weight2); } 0212 //if(upsilon.Pt()>6.0 && upsilon.Pt()<=10.0) { hmass3_flatpt->Fill(upsilon.M(), weight1*weight2); } 0213 int mybin = int(upsilon.Pt()); 0214 //if(mybin==9) mybin = 8; 0215 //if(mybin>=11 && mybin<=14) mybin = 9; 0216 if(mybin>=0 && mybin<nbins) { (hhmass_flatpt[mybin])->Fill(upsilon.M(), weight1*weight2); } 0217 } 0218 */ 0219 // // for proper normalization fill histograms with correct pt 0220 if(pt11>ptcut && pt22>ptcut) { 0221 TLorentzVector upsilon2 = vec11 + vec22; 0222 (hhmass[nbins])->Fill(upsilon2.M()); 0223 //if(upsilon2.Pt()>0.0 && upsilon2.Pt()<=2.0) { hmass0->Fill(upsilon2.M()); } 0224 //if(upsilon2.Pt()>2.0 && upsilon2.Pt()<=4.0) { hmass1->Fill(upsilon2.M()); } 0225 //if(upsilon2.Pt()>4.0 && upsilon2.Pt()<=6.0) { hmass2->Fill(upsilon2.M()); } 0226 //if(upsilon2.Pt()>6.0 && upsilon2.Pt()<=10.0) { hmass3->Fill(upsilon2.M()); } 0227 int mybin = int(upsilon2.Pt()); 0228 //if(mybin==9) mybin = 8; 0229 //if(mybin>=11 && mybin<=14) mybin = 9; 0230 if(mybin>=0 && mybin<nbins) { (hhmass[mybin])->Fill(upsilon2.M()); } 0231 } 0232 0233 }} // end loop over pairs 0234 0235 } // end loop over events 0236 } // end second loop 0237 0238 /* 0239 // normalize flat-pt histogram so that number of entries is the same as integral 0240 0241 for(int i=0; i<nbins+1; i++) { 0242 float nn = (hhmass[i])->Integral(1,nchan); 0243 float scerr = 1./sqrt(nn); 0244 float nn_flatpt = (hhmass_flatpt[i])->Integral(1,nchan); 0245 cout << " Scale Factor = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << " " << hhmass[i]->GetEntries() << endl; 0246 (hhmass_flatpt[i])->Scale(nn/nn_flatpt); 0247 } 0248 0249 float nn = hmass->Integral(1,nchan); 0250 float scerr = 1./sqrt(nn); 0251 float nn_flatpt = hmass_flatpt->Integral(1,nchan); 0252 cout << "Scale factor = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << " " << hmass->GetEntries() << endl; 0253 hmass_flatpt->Scale(nn/nn_flatpt); 0254 0255 nn = hmass0->Integral(1,nchan); 0256 float scerr0 = 1./sqrt(nn); 0257 nn_flatpt = hmass0_flatpt->Integral(1,nchan); 0258 cout << "Scale factor 0 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl; 0259 hmass0_flatpt->Scale(nn/nn_flatpt); 0260 0261 nn = hmass1->Integral(1,nchan); 0262 float scerr1 = 1./sqrt(nn); 0263 nn_flatpt = hmass1_flatpt->Integral(1,nchan); 0264 cout << "Scale factor 1 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl; 0265 hmass1_flatpt->Scale(nn/nn_flatpt); 0266 0267 nn = hmass2->Integral(1,nchan); 0268 float scerr2 = 1./sqrt(nn); 0269 nn_flatpt = hmass2_flatpt->Integral(1,nchan); 0270 cout << "Scale factor 2 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl; 0271 hmass2_flatpt->Scale(nn/nn_flatpt); 0272 0273 nn = hmass3->Integral(1,nchan); 0274 float scerr3 = 1./sqrt(nn); 0275 nn_flatpt = hmass3_flatpt->Integral(1,nchan); 0276 cout << "Scale factor 3 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl; 0277 hmass3_flatpt->Scale(nn/nn_flatpt); 0278 0279 0280 // scale to 10B events 0281 float evtscale = 10000000000./float(stat)/float(nstat); 0282 hmass_flatpt ->Scale(evtscale); 0283 hmass0_flatpt->Scale(evtscale); 0284 hmass1_flatpt->Scale(evtscale); 0285 hmass2_flatpt->Scale(evtscale); 0286 hmass3_flatpt->Scale(evtscale); 0287 for(int i=0; i<nbins; i++) { (hhmass_flatpt[i])->Scale(evtscale); } 0288 0289 // count background and errors 0290 double sumbg=0.; double sumbg0=0.; double sumbg1=0.; double sumbg2=0.; double sumbg3=0.; 0291 double ersumbg=0.; double ersumbg0=0.; double ersumbg1=0.; double ersumbg2=0.; double ersumbg3=0.; 0292 float u1start = 9.10; 0293 float u1stop = 9.60; 0294 int fbin = hmass->FindBin(u1start + 0.001); 0295 int lbin = hmass->FindBin(u1stop - 0.001); 0296 cout << "bin range: " << fbin << " - " << lbin << endl; 0297 cout << "inv. mass range: " << u1start << " - " << u1stop << endl; 0298 for(int i=fbin; i<=lbin; i++) { 0299 sumbg += hmass_flatpt->GetBinContent(i); 0300 sumbg0 += hmass0_flatpt->GetBinContent(i); 0301 sumbg1 += hmass1_flatpt->GetBinContent(i); 0302 sumbg2 += hmass2_flatpt->GetBinContent(i); 0303 sumbg3 += hmass3_flatpt->GetBinContent(i); 0304 ersumbg += hmass_flatpt ->GetBinError(i)*hmass_flatpt ->GetBinError(i); 0305 ersumbg0 += hmass0_flatpt->GetBinError(i)*hmass0_flatpt->GetBinError(i); 0306 ersumbg1 += hmass1_flatpt->GetBinError(i)*hmass1_flatpt->GetBinError(i); 0307 ersumbg2 += hmass2_flatpt->GetBinError(i)*hmass2_flatpt->GetBinError(i); 0308 ersumbg3 += hmass3_flatpt->GetBinError(i)*hmass3_flatpt->GetBinError(i); 0309 } 0310 ersumbg = sqrt(ersumbg); 0311 ersumbg0 = sqrt(ersumbg0); 0312 ersumbg1 = sqrt(ersumbg1); 0313 ersumbg2 = sqrt(ersumbg2); 0314 ersumbg3 = sqrt(ersumbg3); 0315 0316 // The main statistical uncertainty comes from histogram normalization 0317 cout << "In 10B 0-10% central events (100B min. bias):" << endl; 0318 cout << " Background(all pT) = " << sumbg << " +- " << ersumbg << "(stat) +- " << sumbg*scerr << "(scale) pairs." << endl; 0319 cout << " Background(0-2) = " << sumbg0 << " +- " << ersumbg0 << "(stat) +- " << sumbg*scerr0 << "(scale) pairs." << endl; 0320 cout << " Background(2-4) = " << sumbg1 << " +- " << ersumbg1 << "(stat) +- " << sumbg*scerr1 << "(scale) pairs." << endl; 0321 cout << " Background(4-6) = " << sumbg2 << " +- " << ersumbg2 << "(stat) +- " << sumbg*scerr2 << "(scale) pairs." << endl; 0322 cout << " Background(6-10) = " << sumbg3 << " +- " << ersumbg3 << "(stat) +- " << sumbg*scerr3 << "(scale) pairs." << endl; 0323 */ 0324 0325 TH1D* hhfakefake[nbins]; 0326 for(int i=0; i<nbins+1; i++) { 0327 sprintf(hhname,"hhfakefake_%d",i); 0328 hhfakefake[i] = (TH1D*)((hhmass[i])->Clone(hhname)); 0329 } 0330 0331 /* 0332 TH1D* hfakefake = (TH1D*)hmass_flatpt->Clone("hfakefake"); 0333 TH1D* hfakefake0 = (TH1D*)hmass0_flatpt->Clone("hfakefake0"); 0334 TH1D* hfakefake1 = (TH1D*)hmass1_flatpt->Clone("hfakefake1"); 0335 TH1D* hfakefake2 = (TH1D*)hmass2_flatpt->Clone("hfakefake2"); 0336 TH1D* hfakefake3 = (TH1D*)hmass3_flatpt->Clone("hfakefake3"); 0337 0338 //TF1* mymass = new TF1("mymass","[0]*pow(x,[3])*pow(1.+x*x/[1]/[1],[2])",2.0,15.0); 0339 TF1* mymass = new TF1("mymass","[0]*pow(1.+x*x/[1]/[1],[2])",2.0,15.0); 0340 mymass->SetParameter(0,5.0e+09); 0341 mymass->SetParameter(1,2.0); 0342 mymass->SetParameter(2,-10.0); 0343 mymass->SetParLimits(2,-50.0,0.0); 0344 //mymass->SetParameter(3,2.0); 0345 0346 double mylim = 0.5; 0347 double x,a,b,erb; 0348 double fitstart = 7.0; 0349 double fitstop = 12.0; 0350 TH1F* hmynorm = (TH1F*)hmass_flatpt->Clone("hmynorm"); 0351 */ 0352 /* 0353 cout << "kuku2" << endl; 0354 for(int j=0; j<nbins; j++) { 0355 cout << "fitting " << j << " " << (hhmass_flatpt[j])->GetEntries() << endl; 0356 (hhmass_flatpt[j])->Fit("mymass","qr","",fitstart,fitstop); 0357 cout << " fit done." << endl; 0358 for(int i=1; i<=nchan; i++) { 0359 cout << " bin " << i << endl; 0360 x = (hhmass_flatpt[j])->GetBinCenter(i); 0361 cout << " x = " << x << endl; 0362 a = ((hhmass_flatpt[j])->GetFunction("mymass"))->Eval(x); 0363 cout << " a = " << a << endl; 0364 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.00001; erb=0.;} 0365 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0366 cout << " bin " << i << " " << b << " +- " << erb << endl; 0367 (hhfakefake[i])->SetBinContent(i,b); 0368 (hhfakefake[i])->SetBinError(i,erb); 0369 cout << " all set." << endl; 0370 } 0371 } 0372 */ 0373 0374 // hmass_flatpt->Fit("mymass","qr","",fitstart,fitstop); 0375 /* 0376 for(int i=hmass_flatpt->FindBin(fitstart+0.001); i<=hmass_flatpt->FindBin(fitstop-0.001); i++) { 0377 double tmpa = hmass_flatpt->GetBinContent(i)/(hmass_flatpt->GetFunction("mymass")->Eval(hmass_flatpt->GetBinCenter(i))); 0378 hmynorm->SetBinContent(i,tmpa); 0379 //double tmpb = tmpa*(hmass_flatpt->GetBinError(i)/hmass_flatpt->GetBinContent(i)); 0380 double tmpb = (hmass_flatpt->GetBinError(i)/hmass_flatpt->GetBinContent(i)); 0381 hmynorm->SetBinError(i,tmpb); 0382 cout << "corr: " << tmpa << " +- " << tmpb << endl; 0383 } 0384 hmynorm->Fit("pol0","r","",fitstart,fitstop); 0385 double mycorr = hmynorm->GetFunction("pol0")->GetParameter(0); 0386 cout << "Correction factor = " << mycorr << endl; 0387 */ 0388 /* 0389 for(int i=1; i<=nchan; i++) { 0390 x = hmass_flatpt->GetBinCenter(i); 0391 a = (hmass_flatpt->GetFunction("mymass"))->Eval(x); 0392 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.00001; erb=0.;} 0393 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0394 hfakefake->SetBinContent(i,b); 0395 hfakefake->SetBinError(i,erb); 0396 } 0397 0398 hmass0_flatpt->Fit("mymass","qr","",fitstart,fitstop); 0399 for(int i=1; i<=nchan; i++) { 0400 x = hmass0_flatpt->GetBinCenter(i); 0401 a = (hmass0_flatpt->GetFunction("mymass"))->Eval(x); 0402 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;} 0403 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0404 hfakefake0->SetBinContent(i,b); 0405 hfakefake0->SetBinError(i,erb); 0406 } 0407 hmass1_flatpt->Fit("mymass","qr","",fitstart,fitstop); 0408 for(int i=1; i<=nchan; i++) { 0409 x = hmass1_flatpt->GetBinCenter(i); 0410 a = (hmass1_flatpt->GetFunction("mymass"))->Eval(x); 0411 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;} 0412 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0413 hfakefake1->SetBinContent(i,b); 0414 hfakefake1->SetBinError(i,erb); 0415 } 0416 hmass2_flatpt->Fit("mymass","qr","",fitstart,fitstop); 0417 for(int i=1; i<=nchan; i++) { 0418 x = hmass2_flatpt->GetBinCenter(i); 0419 a = (hmass2_flatpt->GetFunction("mymass"))->Eval(x); 0420 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;} 0421 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0422 hfakefake2->SetBinContent(i,b); 0423 hfakefake2->SetBinError(i,erb); 0424 } 0425 hmass3_flatpt->Fit("mymass","qr","",fitstart,fitstop); 0426 for(int i=1; i<=nchan; i++) { 0427 x = hmass3_flatpt->GetBinCenter(i); 0428 a = (hmass3_flatpt->GetFunction("mymass"))->Eval(x); 0429 //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;} 0430 b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b); 0431 hfakefake3->SetBinContent(i,b); 0432 hfakefake3->SetBinError(i,erb); 0433 } 0434 */ 0435 0436 0437 fout->Write(); 0438 fout->Close(); 0439 0440 return 0; 0441 } 0442 0443 0444
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|