Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fakee_invmass.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fakee_invmass.cc
1 #include "TFile.h"
2 #include "TF1.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TNtuple.h"
6 #include "TTree.h"
7 #include "TRandom3.h"
8 #include "TMCParticle.h"
9 #include "TLorentzVector.h"
10 
11 #include <cstdlib>
12 #include <cstdio>
13 #include <iostream>
14 #include <iomanip>
15 #include <fstream>
16 #include <vector>
17 
18 using namespace std;
19 
20 int fakee_invmass(double, string, long int);
21 
22 // This program will generate invariant mass distribution from fake
23 // electrons (mis-identified charged hadrons) in central AuAu collisions
24 // and scales it up to 5B central events (50B min. bias)
25 // You can choose 0-10% (centr=0) or 10-20% (centr=1) central events
26 // You can choose three different pion rejection functions,
27 // or use constant (Pt-independent) pion rejection factor
28 
29 int main(int argc, char* argv[]) {
30  double eideff = 0.7;
31  string ofname="fakee.root";
32  long int stat=500000000;
33  if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl;
34 // if(argc==2) {eideff = atoi(argv[1]);}
35 // if(argc==3) {eideff = atoi(argv[1]); ofname=argv[2]; stat=atoi(argv[3]);}
36  return fakee_invmass(eideff, ofname, stat);
37 }
38 
39 
40 int fakee_invmass(double eideff, string ofname, long int stat) {
41 
42 int nstat = 20;
43 
44 // This is dN/dpT for fake electrons from all sources (pions, kaons,protons, anti-protons)
45 //float start = 1.;
46 double start = 2.;
47 double stop = 20.0;
48 double ptcut=2.0;
49 string plusdistributionAuAu;
50 string minusdistributionAuAu;
51 
52 if(eideff==0.9) {
53 cout << "eID efficiency = " << eideff << endl;
54 // eid efficiency 90%
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))";
57 
58 }
59 else if(eideff==0.7) {
60 cout << "eID efficiency = " << eideff << endl;
61 // eid efficnecy 70%
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))";
64 }
65 else if(eideff==1.0) { // constant pion rejection factor = 90
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.";
69 }
70 else { cerr << "ERROR: eid efficiency must be 0.9, 0.7 or 1.0 (constant rejection factor = 90)" << endl; }
71 
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;
75 
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;
79 
80 
81 
82 TRandom* rnd = new TRandom3();
83 rnd->SetSeed(0);
84 
85 TFile* fout = new TFile(ofname.c_str(),"RECREATE");
86 int nchan = 400;
87 double s1=0.;
88 double s2=20.;
89 //TH1D* hmass = new TH1D("hmass","", nchan,s1,s2);
90 //TH1D* hmass0 = new TH1D("hmass0","",nchan,s1,s2);
91 //TH1D* hmass1 = new TH1D("hmass1","",nchan,s1,s2);
92 //TH1D* hmass2 = new TH1D("hmass2","",nchan,s1,s2);
93 //TH1D* hmass3 = new TH1D("hmass3","",nchan,s1,s2);
94 //TH1D* hmass_flatpt = new TH1D("hmass_flatpt","", nchan,s1,s2);
95 //TH1D* hmass0_flatpt = new TH1D("hmass0_flatpt","",nchan,s1,s2);
96 //TH1D* hmass1_flatpt = new TH1D("hmass1_flatpt","",nchan,s1,s2);
97 //TH1D* hmass2_flatpt = new TH1D("hmass2_flatpt","",nchan,s1,s2);
98 //TH1D* hmass3_flatpt = new TH1D("hmass3_flatpt","",nchan,s1,s2);
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);
102 //hmass->Sumw2();
103 //hmass0->Sumw2();
104 //hmass1->Sumw2();
105 //hmass2->Sumw2();
106 //hmass3->Sumw2();
107 //hmass_flatpt->Sumw2();
108 //hmass0_flatpt->Sumw2();
109 //hmass1_flatpt->Sumw2();
110 //hmass2_flatpt->Sumw2();
111 //hmass3_flatpt->Sumw2();
112 
113 const int nbins = 15;
114 double binlim[nbins+1];
115 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);}
116 //binlim[9] = 10.;
117 //binlim[10] = 15.;
118 
119 TH1D* hhmass[nbins+1];
120 //TH1D* hhmass_flatpt[nbins+1];
121 char hhname[999];
122 for(int i=0; i<nbins+1; i++) {
123  sprintf(hhname,"hhmass_%d",i);
124  hhmass[i] = new TH1D(hhname,"",nchan,s1,s2);
125 // sprintf(hhname,"hhmass_flatpt_%d",i);
126 // hhmass_flatpt[i] = new TH1D(hhname,"",nchan,s1,s2);
127  (hhmass[i])->Sumw2();
128 // (hhmass_flatpt[i])->Sumw2();
129 }
130 
131 // Generate events
132 
133 for(int ii=0; ii<nstat; ii++) {
134 for(int i=0; i<stat; i++) {
135 
136  if(i%10000000==0) cout << "processing event # " << i << endl;
137 
138 // first introduce 25% fluctuations in multiplicity
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;
142 // float n1smeared = norm2minus;
143 // float n2smeared = norm2plus;
144 // then calculate multiplicity for fake electrons and positrons
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;
151 
152  // loop over pairs
153  for(int i1=0; i1<n1; i1++) {
154 
155 /*
156  double pt1 = rnd->Uniform(start,stop); // flat pT
157  double eta1 = rnd->Uniform(-1.0,1.0);
158  double theta1 = 3.141592654/2. + (exp(2.*eta1)-1.)/(exp(2.*eta1)+1.);
159  double phi1 = rnd->Uniform(-3.141592654,3.141592654);
160  double px1 = pt1*cos(phi1);
161  double py1 = pt1*sin(phi1);
162  double pp1 = pt1/sin(theta1);
163  double pz1 = pp1*cos(theta1);
164  double ee1 = pp1;
165  TLorentzVector vec1 = TLorentzVector(px1, py1, pz1, ee1);
166 */
167  double pt11 = minusptdistrAuAu->GetRandom(); // correct pt used for cross-check of histogram normalization
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);
175  double ee11 = pp11;
176  TLorentzVector vec11 = TLorentzVector(px11, py11, pz11, ee11);
177 
178  for(int i2=0; i2<n2; i2++) {
179 /*
180  double pt2 = rnd->Uniform(start,stop); // flat pT
181  double eta2 = rnd->Uniform(-1.0,1.0);
182  double theta2 = 3.141592654/2. + (exp(2.*eta2)-1.)/(exp(2.*eta2)+1.);
183  double phi2 = rnd->Uniform(-3.141592654,3.141592654);
184  double px2 = pt2*cos(phi2);
185  double py2 = pt2*sin(phi2);
186  double pp2 = pt2/sin(theta2);
187  double pz2 = pp2*cos(theta2);
188  double ee2 = pp2;
189  TLorentzVector vec2 = TLorentzVector(px2, py2, pz2, ee2);
190 */
191  double pt22 = plusptdistrAuAu->GetRandom(); // correct pt used for cross-check of histogram normalization
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);
199  double ee22 = pp22;
200  TLorentzVector vec22 = TLorentzVector(px22, py22, pz22, ee22);
201 
202 /*
203  // flat pT
204  if(pt1>ptcut && pt2>ptcut) {
205  TLorentzVector upsilon = vec1 + vec2; // flat pT
206  double weight1 = minusptdistrAuAu->Eval(pt1)/norm2minus; // weight normalized to 1.
207  double weight2 = plusptdistrAuAu->Eval(pt2)/norm2plus;
208  (hhmass_flatpt[nbins])->Fill(upsilon.M(), weight1*weight2);
209  //if(upsilon.Pt()>0.0 && upsilon.Pt()<=2.0) { hmass0_flatpt->Fill(upsilon.M(), weight1*weight2); }
210  //if(upsilon.Pt()>2.0 && upsilon.Pt()<=4.0) { hmass1_flatpt->Fill(upsilon.M(), weight1*weight2); }
211  //if(upsilon.Pt()>4.0 && upsilon.Pt()<=6.0) { hmass2_flatpt->Fill(upsilon.M(), weight1*weight2); }
212  //if(upsilon.Pt()>6.0 && upsilon.Pt()<=10.0) { hmass3_flatpt->Fill(upsilon.M(), weight1*weight2); }
213  int mybin = int(upsilon.Pt());
214  //if(mybin==9) mybin = 8;
215  //if(mybin>=11 && mybin<=14) mybin = 9;
216  if(mybin>=0 && mybin<nbins) { (hhmass_flatpt[mybin])->Fill(upsilon.M(), weight1*weight2); }
217  }
218 */
219 // // for proper normalization fill histograms with correct pt
220  if(pt11>ptcut && pt22>ptcut) {
221  TLorentzVector upsilon2 = vec11 + vec22;
222  (hhmass[nbins])->Fill(upsilon2.M());
223  //if(upsilon2.Pt()>0.0 && upsilon2.Pt()<=2.0) { hmass0->Fill(upsilon2.M()); }
224  //if(upsilon2.Pt()>2.0 && upsilon2.Pt()<=4.0) { hmass1->Fill(upsilon2.M()); }
225  //if(upsilon2.Pt()>4.0 && upsilon2.Pt()<=6.0) { hmass2->Fill(upsilon2.M()); }
226  //if(upsilon2.Pt()>6.0 && upsilon2.Pt()<=10.0) { hmass3->Fill(upsilon2.M()); }
227  int mybin = int(upsilon2.Pt());
228  //if(mybin==9) mybin = 8;
229  //if(mybin>=11 && mybin<=14) mybin = 9;
230  if(mybin>=0 && mybin<nbins) { (hhmass[mybin])->Fill(upsilon2.M()); }
231  }
232 
233  }} // end loop over pairs
234 
235 } // end loop over events
236 } // end second loop
237 
238 /*
239 // normalize flat-pt histogram so that number of entries is the same as integral
240 
241 for(int i=0; i<nbins+1; i++) {
242  float nn = (hhmass[i])->Integral(1,nchan);
243  float scerr = 1./sqrt(nn);
244  float nn_flatpt = (hhmass_flatpt[i])->Integral(1,nchan);
245  cout << " Scale Factor = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << " " << hhmass[i]->GetEntries() << endl;
246  (hhmass_flatpt[i])->Scale(nn/nn_flatpt);
247 }
248 
249  float nn = hmass->Integral(1,nchan);
250  float scerr = 1./sqrt(nn);
251  float nn_flatpt = hmass_flatpt->Integral(1,nchan);
252  cout << "Scale factor = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << " " << hmass->GetEntries() << endl;
253  hmass_flatpt->Scale(nn/nn_flatpt);
254 
255  nn = hmass0->Integral(1,nchan);
256  float scerr0 = 1./sqrt(nn);
257  nn_flatpt = hmass0_flatpt->Integral(1,nchan);
258  cout << "Scale factor 0 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl;
259  hmass0_flatpt->Scale(nn/nn_flatpt);
260 
261  nn = hmass1->Integral(1,nchan);
262  float scerr1 = 1./sqrt(nn);
263  nn_flatpt = hmass1_flatpt->Integral(1,nchan);
264  cout << "Scale factor 1 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl;
265  hmass1_flatpt->Scale(nn/nn_flatpt);
266 
267  nn = hmass2->Integral(1,nchan);
268  float scerr2 = 1./sqrt(nn);
269  nn_flatpt = hmass2_flatpt->Integral(1,nchan);
270  cout << "Scale factor 2 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl;
271  hmass2_flatpt->Scale(nn/nn_flatpt);
272 
273  nn = hmass3->Integral(1,nchan);
274  float scerr3 = 1./sqrt(nn);
275  nn_flatpt = hmass3_flatpt->Integral(1,nchan);
276  cout << "Scale factor 3 = " << nn << " / " << nn_flatpt << " = " << nn/nn_flatpt << endl;
277  hmass3_flatpt->Scale(nn/nn_flatpt);
278 
279 
280 // scale to 10B events
281  float evtscale = 10000000000./float(stat)/float(nstat);
282  hmass_flatpt ->Scale(evtscale);
283  hmass0_flatpt->Scale(evtscale);
284  hmass1_flatpt->Scale(evtscale);
285  hmass2_flatpt->Scale(evtscale);
286  hmass3_flatpt->Scale(evtscale);
287  for(int i=0; i<nbins; i++) { (hhmass_flatpt[i])->Scale(evtscale); }
288 
289 // count background and errors
290  double sumbg=0.; double sumbg0=0.; double sumbg1=0.; double sumbg2=0.; double sumbg3=0.;
291  double ersumbg=0.; double ersumbg0=0.; double ersumbg1=0.; double ersumbg2=0.; double ersumbg3=0.;
292  float u1start = 9.10;
293  float u1stop = 9.60;
294  int fbin = hmass->FindBin(u1start + 0.001);
295  int lbin = hmass->FindBin(u1stop - 0.001);
296  cout << "bin range: " << fbin << " - " << lbin << endl;
297  cout << "inv. mass range: " << u1start << " - " << u1stop << endl;
298  for(int i=fbin; i<=lbin; i++) {
299  sumbg += hmass_flatpt->GetBinContent(i);
300  sumbg0 += hmass0_flatpt->GetBinContent(i);
301  sumbg1 += hmass1_flatpt->GetBinContent(i);
302  sumbg2 += hmass2_flatpt->GetBinContent(i);
303  sumbg3 += hmass3_flatpt->GetBinContent(i);
304  ersumbg += hmass_flatpt ->GetBinError(i)*hmass_flatpt ->GetBinError(i);
305  ersumbg0 += hmass0_flatpt->GetBinError(i)*hmass0_flatpt->GetBinError(i);
306  ersumbg1 += hmass1_flatpt->GetBinError(i)*hmass1_flatpt->GetBinError(i);
307  ersumbg2 += hmass2_flatpt->GetBinError(i)*hmass2_flatpt->GetBinError(i);
308  ersumbg3 += hmass3_flatpt->GetBinError(i)*hmass3_flatpt->GetBinError(i);
309  }
310  ersumbg = sqrt(ersumbg);
311  ersumbg0 = sqrt(ersumbg0);
312  ersumbg1 = sqrt(ersumbg1);
313  ersumbg2 = sqrt(ersumbg2);
314  ersumbg3 = sqrt(ersumbg3);
315 
316 // The main statistical uncertainty comes from histogram normalization
317  cout << "In 10B 0-10% central events (100B min. bias):" << endl;
318  cout << " Background(all pT) = " << sumbg << " +- " << ersumbg << "(stat) +- " << sumbg*scerr << "(scale) pairs." << endl;
319  cout << " Background(0-2) = " << sumbg0 << " +- " << ersumbg0 << "(stat) +- " << sumbg*scerr0 << "(scale) pairs." << endl;
320  cout << " Background(2-4) = " << sumbg1 << " +- " << ersumbg1 << "(stat) +- " << sumbg*scerr1 << "(scale) pairs." << endl;
321  cout << " Background(4-6) = " << sumbg2 << " +- " << ersumbg2 << "(stat) +- " << sumbg*scerr2 << "(scale) pairs." << endl;
322  cout << " Background(6-10) = " << sumbg3 << " +- " << ersumbg3 << "(stat) +- " << sumbg*scerr3 << "(scale) pairs." << endl;
323 */
324 
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));
329 }
330 
331 /*
332  TH1D* hfakefake = (TH1D*)hmass_flatpt->Clone("hfakefake");
333  TH1D* hfakefake0 = (TH1D*)hmass0_flatpt->Clone("hfakefake0");
334  TH1D* hfakefake1 = (TH1D*)hmass1_flatpt->Clone("hfakefake1");
335  TH1D* hfakefake2 = (TH1D*)hmass2_flatpt->Clone("hfakefake2");
336  TH1D* hfakefake3 = (TH1D*)hmass3_flatpt->Clone("hfakefake3");
337 
338  //TF1* mymass = new TF1("mymass","[0]*pow(x,[3])*pow(1.+x*x/[1]/[1],[2])",2.0,15.0);
339  TF1* mymass = new TF1("mymass","[0]*pow(1.+x*x/[1]/[1],[2])",2.0,15.0);
340  mymass->SetParameter(0,5.0e+09);
341  mymass->SetParameter(1,2.0);
342  mymass->SetParameter(2,-10.0);
343  mymass->SetParLimits(2,-50.0,0.0);
344  //mymass->SetParameter(3,2.0);
345 
346  double mylim = 0.5;
347  double x,a,b,erb;
348  double fitstart = 7.0;
349  double fitstop = 12.0;
350  TH1F* hmynorm = (TH1F*)hmass_flatpt->Clone("hmynorm");
351 */
352 /*
353 cout << "kuku2" << endl;
354 for(int j=0; j<nbins; j++) {
355  cout << "fitting " << j << " " << (hhmass_flatpt[j])->GetEntries() << endl;
356  (hhmass_flatpt[j])->Fit("mymass","qr","",fitstart,fitstop);
357  cout << " fit done." << endl;
358  for(int i=1; i<=nchan; i++) {
359  cout << " bin " << i << endl;
360  x = (hhmass_flatpt[j])->GetBinCenter(i);
361  cout << " x = " << x << endl;
362  a = ((hhmass_flatpt[j])->GetFunction("mymass"))->Eval(x);
363  cout << " a = " << a << endl;
364  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.00001; erb=0.;}
365  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
366  cout << " bin " << i << " " << b << " +- " << erb << endl;
367  (hhfakefake[i])->SetBinContent(i,b);
368  (hhfakefake[i])->SetBinError(i,erb);
369  cout << " all set." << endl;
370  }
371 }
372 */
373 
374 // hmass_flatpt->Fit("mymass","qr","",fitstart,fitstop);
375 /*
376  for(int i=hmass_flatpt->FindBin(fitstart+0.001); i<=hmass_flatpt->FindBin(fitstop-0.001); i++) {
377  double tmpa = hmass_flatpt->GetBinContent(i)/(hmass_flatpt->GetFunction("mymass")->Eval(hmass_flatpt->GetBinCenter(i)));
378  hmynorm->SetBinContent(i,tmpa);
379  //double tmpb = tmpa*(hmass_flatpt->GetBinError(i)/hmass_flatpt->GetBinContent(i));
380  double tmpb = (hmass_flatpt->GetBinError(i)/hmass_flatpt->GetBinContent(i));
381  hmynorm->SetBinError(i,tmpb);
382  cout << "corr: " << tmpa << " +- " << tmpb << endl;
383  }
384  hmynorm->Fit("pol0","r","",fitstart,fitstop);
385  double mycorr = hmynorm->GetFunction("pol0")->GetParameter(0);
386  cout << "Correction factor = " << mycorr << endl;
387 */
388 /*
389  for(int i=1; i<=nchan; i++) {
390  x = hmass_flatpt->GetBinCenter(i);
391  a = (hmass_flatpt->GetFunction("mymass"))->Eval(x);
392  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.00001; erb=0.;}
393  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
394  hfakefake->SetBinContent(i,b);
395  hfakefake->SetBinError(i,erb);
396  }
397 
398  hmass0_flatpt->Fit("mymass","qr","",fitstart,fitstop);
399  for(int i=1; i<=nchan; i++) {
400  x = hmass0_flatpt->GetBinCenter(i);
401  a = (hmass0_flatpt->GetFunction("mymass"))->Eval(x);
402  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;}
403  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
404  hfakefake0->SetBinContent(i,b);
405  hfakefake0->SetBinError(i,erb);
406  }
407  hmass1_flatpt->Fit("mymass","qr","",fitstart,fitstop);
408  for(int i=1; i<=nchan; i++) {
409  x = hmass1_flatpt->GetBinCenter(i);
410  a = (hmass1_flatpt->GetFunction("mymass"))->Eval(x);
411  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;}
412  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
413  hfakefake1->SetBinContent(i,b);
414  hfakefake1->SetBinError(i,erb);
415  }
416  hmass2_flatpt->Fit("mymass","qr","",fitstart,fitstop);
417  for(int i=1; i<=nchan; i++) {
418  x = hmass2_flatpt->GetBinCenter(i);
419  a = (hmass2_flatpt->GetFunction("mymass"))->Eval(x);
420  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;}
421  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
422  hfakefake2->SetBinContent(i,b);
423  hfakefake2->SetBinError(i,erb);
424  }
425  hmass3_flatpt->Fit("mymass","qr","",fitstart,fitstop);
426  for(int i=1; i<=nchan; i++) {
427  x = hmass3_flatpt->GetBinCenter(i);
428  a = (hmass3_flatpt->GetFunction("mymass"))->Eval(x);
429  //if(a>mylim) { b = int(rnd->Gaus(a,sqrt(a))+0.5); erb = sqrt(b); } else {b=0.; erb=0.;}
430  b = int(rnd->Poisson(a)+0.5); if(b<0.){b=0.;} erb = sqrt(b);
431  hfakefake3->SetBinContent(i,b);
432  hfakefake3->SetBinError(i,erb);
433  }
434 */
435 
436 
437  fout->Write();
438  fout->Close();
439 
440  return 0;
441 }
442 
443 
444