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