StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
makeUA1Ratio.C
1 #include "commonmacro/common.h"
2 #include "common/Name.cc"
3 #include "commonmacro/histutil.h"
4 
5 void makeUA1Ratio(
6  const char* inName=
7  "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
8  const char* psDir="ps",
9  int cut = 1,
10  const char* outDir="./",
11  const char* more = "west",
12  float extraValue = 10)
13 {
14  cout << "--------------------------" << endl;
15  cout << "in name=" << inName << endl
16  << "ps dir=" << psDir << endl
17  << "cut=" << cut << endl;
18  cout << "--------------------------" << endl;
19 
20  inRoot = new TFile(inName);
21 
22  if(!inRoot){
23  cout << "cannot find the infile" << endl;
24  return;
25  }
26  gStyle->SetOptStat(0);gStyle->SetPadTickX(1);gStyle->SetPadTickY(1);
27  TCanvas c1("c1","c1",500,400);
28  //-------------------------------------------------------
29  const int nbin=2;
30  TGraphAsymmErrors* gSpec[2]; //plus+minus
31  TGraphAsymmErrors* gSpecPM[2][2]; // plus,minus
32  TGraphAsymmErrors* gUA1, *gRatio;
33  // get stuff
34 
35  for(int ibin=0;ibin<nbin;ibin++){
36  setName(name,"gSpecCorrected",ibin);
37  cout << name << endl;
38  gSpec[ibin]=(TGraphAsymmErrors*)inRoot.Get(name);
39  cout << gSpec[ibin]->GetName() << endl;
40  if(!gSpec[ibin]) return;
41  for(int ipm=0;ipm<2;ipm++){
42  setName(name,"gSpecCorrected",ibin,sPM[ipm].Data());
43  cout << name << endl;
44  gSpecPM[ibin][ipm]=(TGraphAsymmErrors*)inRoot.Get(name);
45  if(!gSpecPM[ibin][ipm]) return;
46  }
47  }
48  // ratios
49  for(int ibin=0;ibin<2; ibin++){
50  cout << "setting title"<< endl;
51  sprintf(title,"%s / UA1 (bin %d)(cut %d)",
52  gSpec[ibin]->GetName(),ibin,cut);
53  Divide(&c1,1,1,title,inName);
54  cout << "making" << endl;
55  gUA1=makeUA1(gSpec[ibin]);
56  gRatio=makeRatio(gSpec[ibin],gUA1);
57  gRatio->Draw("ap"); sprintf(title,"%sOverUA1",gSpec[ibin]->GetName());
58  gRatio->SetMaximum(1.2); gRatio->SetMinimum(.2);
59  Print(&c1,psDir,title);
60  for(int ipm=0;ipm<2;ipm++){
61  // sName=gSpec[ibin][ipm]->GetName(); // charge sign
62  // sName.Replace(0,sName.First("."),"");
63  //sName.Replace(sName.First("."),sName.Length(),"");
64  sprintf(title,"%s / UA1 (bin %d)(cut %d)",
65  sName.Data(),gSpecPM[ibin][ipm]->GetName(),ibin,cut);
66  Divide(&c1,1,1,title,inName);
67  //gUA1=makeUA1(gSpec[ibin][ipm]);
68  gRatio=makeRatio(gSpecPM[ibin][ipm],gUA1);
69  gRatio->Draw("ap");
70  gRatio->SetMaximum(1.2); gRatio->SetMinimum(.2);
71 
72  sprintf(title,"%sOverUA1",gSpecPM[ibin][ipm]->GetName());
73  Print(&c1,psDir,title);
74 
75  }
76  }
77 
78 
79 
80 
81 }
82 
83 
84 TGraphAsymmErrors*
85 makeUA1(const TGraphAsymmErrors* graph)
86 {
87 
88  //
89  // find the bins
90  //
91  const Int_t nBin = graph->GetN();
92  Double_t* xValues = graph->GetX();
93  Double_t* yValues = graph->GetY();
94  Double_t* errXLow = graph->GetEXlow();
95  Double_t* errXHigh = graph->GetEXhigh();
96 
97  //
98  // ratios
99  //
100  /*
101  1.500000 0.6146703 6.9558643E-02 0.1131642
102  2.000000 0.5883511 8.5044047E-03 1.4454641E-02
103  2.500000 0.5589531 1.7606174E-03 3.1498480E-03
104  3.000000 0.5283710 4.8784356E-04 9.2329737E-04
105  3.500000 0.5029035 1.6314255E-04 3.2440131E-04
106  4.000000 0.4790372 6.2226289E-05 1.2989865E-04
107  4.500000 0.4574749 2.6327352E-05 5.7549289E-05
108  5.000000 0.4378719 1.2056932E-05 2.7535298E-05
109  5.500000 0.4186996 5.9010908E-06 1.4093855E-05
110  6.000000 0.4019553 3.0456495E-06 7.5770845E-06
111  */
112 
113  sprintf(name,"h1tmp%s",graph->GetName());
114  TH1D* h1 = new TH1D(name,name,9,1.5,6.0);
115 
116  h1->SetBinContent(1,.61467);
117  h1->SetBinContent(2,.58835);
118  h1->SetBinContent(3,.55895);
119  h1->SetBinContent(4,.52837);
120  h1->SetBinContent(5,.50290);
121  h1->SetBinContent(6,.47903);
122  h1->SetBinContent(7,.45747);
123  h1->SetBinContent(8,.43787);
124  h1->SetBinContent(9,.40195);
125 
126 
127  // fixed
128  Double_t A = 266.6;
129  Double_t p0 = 1.895;
130  Double_t n = 12.98;
131 
132 
133  // old
134  /*
135  double A=286;
136  double p0=1.8;
137  double n=12.14;
138  */
139  Double_t TAA = 26;
140 
141  char func[200], mean[200];
142 
143  sprintf(func,"2*pi*%.0f*%.0f*(1+x/%.2f)^(-%.2f)",A,TAA,p0,n);
144 
145  sprintf(mean,"2*pi*x*%.0f*%.0f*(1+x/%.2f)^(-%.2f)",A,TAA,p0,n);
146 
147 
148  TF1* fUA1 = new TF1("fUA1",func,lowPt,highPt);
149  TF1* fMean = new TF1("fUA1Mean",mean,lowPt,highPt);
150 
151  gUA1 = new TGraphAsymmErrors;
152 
153  for(Int_t i=0; i<nBin; i++){
154  //
155  // find the x value
156  //
157  Double_t lowEdge = xValues[i]-errXLow[i];
158  Double_t highEdge = xValues[i]+errXHigh[i];
159 
160 
161  Double_t num = fMean->Integral(lowEdge,highEdge);
162  Double_t denom = fUA1->Integral(lowEdge,highEdge);
163 
164  //
165  // just use star's x
166  //
167 
168  Double_t UA1ErrLow = fabs(lowEdge-xValues[i]);
169  Double_t UA1ErrHigh= fabs(highEdge-xValues[i]);
170 
171  //Double_t y = fUA1->Eval(xValues[i]);
172  double y = denom/(highEdge-lowEdge);
173 
174  // ratio
175 
176  //Int_t iBin = h1->GetXaxis()->FindBin(xValues[i]);
177  //Double_t fraction = h1->GetBinContent(iBin);
178 
179  // cout << "ratio : " << h1->GetXaxis()->GetBinLowEdge(iBin) << " "
180  // << h1->GetXaxis()->GetBinUpEdge(iBin) << " = " << fraction << endl;
181 
182 
183  //Double_t yCorrected = fraction*y;
184  double yCorrected = y;
185 
186  cout << "UA1 uncorrected " << y << ", corrected y " << yCorrected
187  << " graph y " << yValues[i] << endl;
188 
189  gUA1->SetPoint(i,xValues[i],yCorrected);
190  gUA1->SetPointError(i,UA1ErrLow,UA1ErrHigh,0,0);
191  }
192 
193  gUA1->SetName("gUA1Spectrum");
194  gUA1->SetTitle("gUA1Spectrum");
195 
196  return gUA1;
197 
198 }
199 
200 
201 TGraphAsymmErrors*
202 makeRatio(TGraphAsymmErrors* gStar, TGraphAsymmErrors* gUA1)
203 {
204  cout << "********************************" << endl;
205  cout << gStar->GetName() << endl;
206  const Int_t nBin = gStar->GetN();
207 
208  Double_t* starY = gStar->GetY();
209  Double_t* starX = gStar->GetX();
210  Double_t* starXErrLow = gStar->GetEXlow();
211  Double_t* starXErrHigh = gStar->GetEXhigh();
212  Double_t* starYErrHigh = gStar->GetEYhigh();
213 
214  Double_t* ua1Y = gUA1->GetY();
215  Double_t* ua1YErrHigh = gUA1->GetEYhigh();
216 
217  TGraphAsymmErrors* gRatio = new TGraphAsymmErrors;
218 
219  for(Int_t i=0; i<nBin; i++){
220 
221  Double_t ratio = starY[i]/ua1Y[i];
222 
223  Double_t yError =
224  TMath::Sqrt(ratio*ratio*(starYErrHigh[i]*starYErrHigh[i]/(starY[i]*starY[i])));
225  // ua1YErrHigh[i]*ua1YErrHigh[i]/(ua1Y[i]*ua1Y[i])));
226 
227  cout << "star x : " << starX[i] << endl;
228  cout << "star y : " << starY[i] << " ua1 Y : " << ua1Y[i]
229  << " ratio : " << ratio << endl
230  << "error : " << yError << endl;
231 
232  gRatio->SetPoint(i,starX[i],ratio);
233  gRatio->SetPointError(i,starXErrLow[i],starXErrHigh[i],yError,yError);
234  }
235 
236  TString sName = gStar->GetName();
237  sName.Append("UA1Ratio");
238 
239  gRatio->SetName(sName.Data());
240  gRatio->SetTitle(sName.Data());
241 
242  return gRatio;
243 
244 }
Generated by   doxygen 1.8.5