StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
FitPt2Mass.C
1 // After we test the parameter values and fit function from FitPtMass.C, we fix some parameters and start generating four spin-depdendent pi0 yields.
2 // The default input file is the root file after normalization from normalizemass.C. The default file name is 'allfill.hist.root'.
3 // We are looping over four spin-dependent plot from the root file "hMassPtUU","hMassPtUD","hMassPtDU","hMassPtDD" to generate yield.
4 // The output will be 4*7 fitted histograms and a pi0 yield txt file. By changing the fit ranges three times, we get three yield txts, which will be used by Calall.C to calculate the double spin asymmetry.
5 // Author: Weihong He.
6 
7 const int mx=4;
8 const int mt=1;
9 TFile *fd[mx];
10 TCanvas *c1=0;
11 FILE* fout=fopen("yield.txt","w");assert(fout);
12 FILE* fmass=fopen("mass.txt","w");assert(fout);
13 TF1* f1[7];
14 
15 FitPt2Mass(){
16  fprintf(fout,"ss pt_mean totpi realpi background\n");
17  float fptmin[7]={4.0,5.0,6.0,7.0,8.0,9.0,10.0};
18  float fptmax[7]={5.0,6.0,7.0,8.0,9.0,10.0,25.0};
19  float ffit_x1[7]={0.03,0.03,0.03,0.05,0.05,0.05,0.07};
20  float ffit_x2=0.4;
21  float fpar0[7]={8.0,7.8,7.2,6.4,5.6,4.9,4.9};
22  float fpar1[7]={-8.33,-6.82,-5.15,-4.14,-3.52,-3.0,-2.9};
23  float fpar2[7]={3000,3000,2000,1000,500,400,400};
24  float fpar3[7]={0.1308,0.1356,0.136,0.137,0.138,0.1375,0.1385};
25  float fpar4[7]={0.0208,0.02166,0.0231,0.025,0.027,0.027,0.029};
26  f1[0]=new TF1("ffit","exp([0]-8.262*x)+[1]*exp(-0.5*((x-0.1305)/(0.02073*(1.+6.11*(x-0.1305))))**2)",ffit_x1[0],ffit_x2);
27  f1[1]=new TF1("ffit","exp([0]-6.752*x)+[1]*exp(-0.5*((x-0.1353)/(0.02164*(1.+6.11*(x-0.1353))))**2)",ffit_x1[1],ffit_x2);
28  f1[2]=new TF1("ffit","exp([0]-5.095*x)+[1]*exp(-0.5*((x-0.1359)/(0.02307*(1.+6.11*(x-0.1359))))**2)",ffit_x1[2],ffit_x2);
29  f1[3]=new TF1("ffit","exp([0]-4.041*x)+[1]*exp(-0.5*((x-0.1375)/(0.02498*(1.+6.11*(x-0.1375))))**2)",ffit_x1[3],ffit_x2);
30  f1[4]=new TF1("ffit","exp([0]-3.37*x)+[1]*exp(-0.5*((x-0.1374)/(0.02666*(1.+6.11*(x-0.1374))))**2)",ffit_x1[4],ffit_x2);
31  f1[5]=new TF1("ffit","exp([0]-2.942*x)+[1]*exp(-0.5*((x-0.137)/(0.0266*(1.+6.11*(x-0.137))))**2)",ffit_x1[5],ffit_x2);
32  f1[6]=new TF1("ffit","exp([0]-2.668*x)+[1]*exp(-0.5*((x-0.1392)/(0.02792*(1.+6.11*(x-0.1392))))**2)",ffit_x1[6],ffit_x2);
33  for(int i=0;i<7;i++){
34  if(i<3){
35  const int fnend=int(ffit_x1[i]*100)+1;
36  }
37  else{
38  const int fnend=int(ffit_x1[i]*100);
39  }
40  const int fnstart=int(ffit_x2*100);
41  cout<<"i="<<i<<" nend="<<fnend<<" nstart="<<fnstart<<endl;
42 
43  FitPtMass(i, fptmin[i], fptmax[i], ffit_x1[i], ffit_x2, fpar0[i], fpar1[i], fpar2[i], fnend, fnstart);
44  }
45 }
46 
47 void FitPtMass(int key, float ptmin, float ptmax, float fit_x1, float fit_x2, float par0, float par1, float par2, int nend, int nstart) {
48 
49  int Realpi=0,backpi=0,totpi=0;
50  //int minbin,maxbin;
51  //TH1F* h1= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
52  //TH1F* h2= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
53  //TH1F* hResidual= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
54 
55  //fprintf(fout,"Bin# Entry\n");
56  char *PlotName[mx]={"hMassPtUU","hMassPtUD","hMassPtDU","hMassPtDD"};
57  char *fName[mt]={"allfill.hist"};
58  TString inPath="";
59 
60  gStyle->SetPalette(1,0);
61  int i;
62  for(i=0;i<mt;i++) {
63  TString hFile=inPath+fName[i];
64  hFile+=".root";
65  fd0=new TFile(hFile); assert(fd0->IsOpen());
66  fd[i]=fd0;
67  }
68  //int ploti;
69  int channel_yield[mx][40];
70  int channel_pi_yield[mx][40];
71  int channel_bg_yield[mx][40];
72  for(int cmx=0;cmx<mx;cmx++)
73  {
74  for(int i=0;i<40;i++)
75  {
76  channel_yield[cmx][i]=0;
77  channel_pi_yield[cmx][i]=0;
78  channel_bg_yield[cmx][i]=0;
79  }
80  }
81  for(int ploti=0;ploti<mx;ploti++)
82  {
83  if(c1) delete c1;
84  c1=new TCanvas("c1","c1",500,500);
85  c1->Divide(1,1);
86  c1->cd(1);
87  //gStyle->SetPalette(1,0);
88  // .... inp eve
89  int minbin,maxbin;
90  TH1F* h1= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
91  TH1F* h2= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
92  TH1F* hResidual= new TH1F("hMassAny","diphoton invariant mass",120,0.,1.2 );
93  TString hname=PlotName[ploti];
94  printf("i= %d =%s=\n",ploti,hname.Data());
95  h0=(TH2F*)fd[0]->Get(hname); assert(h0);
96  minbin=h0->GetYaxis()->FindBin(ptmin);
97  maxbin=h0->GetYaxis()->FindBin(ptmax)-1;
98  h1->Add(h0->ProjectionX("htemp",minbin,maxbin,""));
99  int ent=h1->Integral();
100  cout<<"entry="<<ent<<endl;
101  h1->SetEntries(ent);
102  int ccount=0;
103  float sum=0;
104  float tpt=ptmin;
105  for(int i=minbin;i<=maxbin;i++){
106  ccount+=h0->ProjectionX("",i,i,"")->Integral();
107  sum+=(h0->ProjectionX("",i,i,"")->Integral())*tpt;
108  tpt+=(ptmax-ptmin)/(maxbin+1-minbin);
109  }
110  float ptmean=sum/ccount;
111  cout<<"raw="<<int(h1->Integral(11,19))<<" ptmean="<<ptmean<<endl;
112  h1->Draw();
113  h1->SetMinimum(-ent/100);
114  //TF1* f1=new TF1("ffit","exp([0]+par1*x)+[2]*exp(-0.5*((x-par3)/(par4*(1.+6.11*(x-par3))))**2)",fit_x1,fit_x2);
115  f1[key]->SetParameters(par0,par2);
116  h1->Fit(f1[key],"R+","",fit_x1,fit_x2);
117  int nx=h1->GetNbinsX();
118  cout<<"nx="<<nx<<endl;
119  float sumx=0,sum2;
120  for(int ix=1;ix<=40;ix++)
121  {
122  channel_yield[ploti][ix-1]+=h1->GetBinContent(ix);
123 
124  sumx+=h1->GetBinContent(ix);
125  }
126  //cout<<"sumx="<<sumx<<endl;
127  TF1* f2=new TF1("fit2","expo",0.,1.2);
128  //TF1* f2=new TF1("fit2","exp(5.7-6.4*x)",0.,1.2);
129  h2->Add(h1);
130  hResidual->Add(h1);
131  hResidual->Add(f1[key],-1);
132  for (int ix= 1; ix <= nend; ix++) {
133  hResidual->SetBinContent(ix, 0);
134  }
135  for (int ix= nstart; ix <= nx; ix++) {
136  hResidual->SetBinContent(ix, 0);
137  }
138  hResidual->SetLineColor(11);
139  hResidual->Draw("same");
140  f2->SetParameters(f1[key]->GetParameter(0),par1);
141  h2->Add(f2,-1);
142  for (int ix= 1; ix <= 8; ix++) {
143  //if (h2->GetBinContent(ix) < 0) h2->SetBinContent(ix, 0);
144  h2->SetBinContent(ix, 0);
145  }
146  //for(int ix=9;ix<=19;ix++){
147  //sum2+=h2->GetBinContent(ix);
148  //}
149  h2->SetLineColor(kBlue);
150  TF1* f3=new TF1("f3","expo(0)",fit_x1,fit_x2);
151  f3->SetParameters(f1[key]->GetParameter(0),par1);
152  cout<<"p0="<<f1[key]->GetParameter(0)<<endl;
153  //h2->Fit("gaus","","",0.42,0.66);
154  Realpi=h2->Integral(11,19);
155  backpi=h1->Integral(11,19)-h2->Integral(11,19);
156  totpi=h1->Integral(11,19);
157 
158  for(int ix=1;ix<=40;ix++)
159  {
160  channel_pi_yield[ploti][ix-1]+=h2->Integral(ix,ix);
161  //cout<<"channel_pi_yield[ploti][ix-1]="<<channel_pi_yield[ploti][ix-1]<<endl;
162  channel_bg_yield[ploti][ix-1]+=h1->Integral(ix,ix)-h2->Integral(ix,ix);
163  //cout<<"channel_bg_yield[ploti][ix-1]="<<channel_bg_yield[ploti][ix-1]<<endl;
164  }
165  //ttotpi=h1->GetBinContent(11)+h1->GetBinContent(12)+h1->GetBinContent(13)+h1->GetBinContent(14)+h1->GetBinContent(15)+h1->GetBinContent(16)+h1->GetBinContent(17)+h1->GetBinContent(18)+h1->GetBinContent(19);
166  cout<<"total pi="<<totpi<<" Real Pi0 integral="<<Realpi<<" backpi="<<backpi<<endl;
167  //cout<<"check tot pi="<<ttotpi<<endl;
168  fprintf(fout,"%d %6f %d %d %d\n",ploti,ptmean,totpi,Realpi,backpi);
169  h2->Draw("same");
170  f3->SetLineColor(6);
171  f3->Draw("same");
172  //h2->Draw("error");
173  l1=new TLine(0.135,0.,0.135,7000);
174  l1->SetLineColor(kRed);
175  l1->Draw();
176  l2=new TLine(0.18,0.,0.18,7000);
177  l2->SetLineColor(kGreen);
178  l2->Draw();
179  l3=new TLine(0.1,0.,0.1,7000);
180  l3->SetLineColor(kGreen);
181  l3->Draw();
182  //h2->SetStats(kFALSE);
183  TString hMassAny="mass";
184  hMassAny+=key;
185  hMassAny+=ploti;
186  c1->Print(hMassAny+".gif");
187 
188  }//end of mx
189  if(key==5)
190  {
191  for(int j=0;j<40;j++)
192  {
193  int ySame=channel_yield[0][j]+channel_yield[3][j];
194  int yAnti=channel_yield[1][j]+channel_yield[2][j];
195  float ypi=channel_pi_yield[0][j]+channel_pi_yield[1][j]+channel_pi_yield[2][j]+channel_pi_yield[3][j];
196  if(ypi<0.0) ypi=0.0;
197  float ybg=channel_bg_yield[0][j]+channel_bg_yield[1][j]+channel_bg_yield[2][j]+channel_bg_yield[3][j];
198  float ratio=0.0;
199  if(ybg==0)
200  {
201  ratio=0;
202  }
203  else{
204  ratio=ypi/ybg;
205  }
206  fprintf(fmass,"%d %d %d %f\n",j+1,ySame,yAnti,ratio);
207  }
208  }
209  return;
210 
211 }
212 
Generated by   doxygen 1.8.5