StRoot: StHighptPool/macros/makeUA1Ratio.C Source File

00001 #include "commonmacro/common.h"
00002 #include "common/Name.cc"
00003 #include "commonmacro/histutil.h"
00004 
00005 void makeUA1Ratio(
00006                   const char* inName=
00007                   "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
00008                   const char* psDir="ps",
00009                   int cut = 1,
00010                   const char* outDir="./",
00011                   const char* more = "west",
00012                   float extraValue = 10)
00013 {
00014   cout << "--------------------------" << endl;
00015   cout << "in name=" << inName << endl
00016        << "ps dir=" << psDir << endl
00017        << "cut=" << cut << endl;
00018   cout << "--------------------------" << endl;
00019  
00020   inRoot = new TFile(inName);
00021   
00022   if(!inRoot){
00023     cout << "cannot find the infile" << endl;
00024     return;
00025   }
00026   gStyle->SetOptStat(0);gStyle->SetPadTickX(1);gStyle->SetPadTickY(1);
00027   TCanvas c1("c1","c1",500,400);
00028   //-------------------------------------------------------
00029   const int nbin=2;
00030   TGraphAsymmErrors* gSpec[2]; //plus+minus
00031   TGraphAsymmErrors* gSpecPM[2][2]; // plus,minus
00032   TGraphAsymmErrors* gUA1, *gRatio;
00033   // get stuff
00034   
00035   for(int ibin=0;ibin<nbin;ibin++){
00036     setName(name,"gSpecCorrected",ibin);
00037     cout << name << endl;
00038     gSpec[ibin]=(TGraphAsymmErrors*)inRoot.Get(name); 
00039     cout << gSpec[ibin]->GetName() << endl;
00040     if(!gSpec[ibin]) return;
00041     for(int ipm=0;ipm<2;ipm++){
00042       setName(name,"gSpecCorrected",ibin,sPM[ipm].Data());
00043       cout << name << endl;
00044       gSpecPM[ibin][ipm]=(TGraphAsymmErrors*)inRoot.Get(name); 
00045       if(!gSpecPM[ibin][ipm]) return;
00046     }
00047   }
00048   // ratios
00049   for(int ibin=0;ibin<2; ibin++){
00050     cout << "setting title"<< endl;
00051     sprintf(title,"%s / UA1 (bin %d)(cut %d)",
00052             gSpec[ibin]->GetName(),ibin,cut);
00053     Divide(&c1,1,1,title,inName);
00054     cout << "making" << endl;
00055     gUA1=makeUA1(gSpec[ibin]);
00056     gRatio=makeRatio(gSpec[ibin],gUA1);
00057     gRatio->Draw("ap"); sprintf(title,"%sOverUA1",gSpec[ibin]->GetName());
00058     gRatio->SetMaximum(1.2); gRatio->SetMinimum(.2);
00059     Print(&c1,psDir,title);
00060     for(int ipm=0;ipm<2;ipm++){
00061       //     sName=gSpec[ibin][ipm]->GetName(); // charge sign
00062       // sName.Replace(0,sName.First("."),""); 
00063       //sName.Replace(sName.First("."),sName.Length(),"");
00064       sprintf(title,"%s / UA1 (bin %d)(cut %d)",
00065               sName.Data(),gSpecPM[ibin][ipm]->GetName(),ibin,cut);
00066       Divide(&c1,1,1,title,inName);
00067       //gUA1=makeUA1(gSpec[ibin][ipm]);
00068       gRatio=makeRatio(gSpecPM[ibin][ipm],gUA1);
00069       gRatio->Draw("ap");
00070       gRatio->SetMaximum(1.2); gRatio->SetMinimum(.2);
00071  
00072       sprintf(title,"%sOverUA1",gSpecPM[ibin][ipm]->GetName());
00073       Print(&c1,psDir,title);
00074 
00075     }
00076   }
00077     
00078 
00079 
00080 
00081 }
00082 
00083 
00084 TGraphAsymmErrors*
00085 makeUA1(const TGraphAsymmErrors* graph)
00086 {
00087   
00088   //
00089   // find the bins
00090   //
00091   const Int_t nBin = graph->GetN();
00092   Double_t* xValues = graph->GetX();
00093   Double_t* yValues = graph->GetY();
00094   Double_t* errXLow = graph->GetEXlow();
00095   Double_t* errXHigh = graph->GetEXhigh();
00096 
00097   // 
00098   // ratios
00099   //
00100   /*
00101    1.500000      0.6146703      6.9558643E-02  0.1131642    
00102    2.000000      0.5883511      8.5044047E-03  1.4454641E-02
00103    2.500000      0.5589531      1.7606174E-03  3.1498480E-03
00104    3.000000      0.5283710      4.8784356E-04  9.2329737E-04
00105    3.500000      0.5029035      1.6314255E-04  3.2440131E-04
00106    4.000000      0.4790372      6.2226289E-05  1.2989865E-04
00107    4.500000      0.4574749      2.6327352E-05  5.7549289E-05
00108    5.000000      0.4378719      1.2056932E-05  2.7535298E-05
00109    5.500000      0.4186996      5.9010908E-06  1.4093855E-05
00110    6.000000      0.4019553      3.0456495E-06  7.5770845E-06
00111   */
00112 
00113   sprintf(name,"h1tmp%s",graph->GetName());
00114   TH1D* h1 = new TH1D(name,name,9,1.5,6.0);
00115 
00116   h1->SetBinContent(1,.61467);
00117   h1->SetBinContent(2,.58835);
00118   h1->SetBinContent(3,.55895);
00119   h1->SetBinContent(4,.52837);
00120   h1->SetBinContent(5,.50290);
00121   h1->SetBinContent(6,.47903);
00122   h1->SetBinContent(7,.45747);
00123   h1->SetBinContent(8,.43787);
00124   h1->SetBinContent(9,.40195);
00125 
00126   
00127   // fixed
00128   Double_t A  = 266.6;
00129   Double_t p0 = 1.895;
00130   Double_t n  = 12.98;
00131   
00132 
00133   // old
00134   /*
00135   double A=286;
00136   double p0=1.8;
00137   double n=12.14;
00138   */
00139   Double_t TAA = 26;
00140 
00141   char func[200], mean[200];
00142   
00143   sprintf(func,"2*pi*%.0f*%.0f*(1+x/%.2f)^(-%.2f)",A,TAA,p0,n);
00144 
00145   sprintf(mean,"2*pi*x*%.0f*%.0f*(1+x/%.2f)^(-%.2f)",A,TAA,p0,n);
00146 
00147 
00148   TF1* fUA1 = new TF1("fUA1",func,lowPt,highPt);
00149   TF1* fMean = new TF1("fUA1Mean",mean,lowPt,highPt);
00150 
00151   gUA1 = new TGraphAsymmErrors;
00152 
00153   for(Int_t i=0; i<nBin; i++){
00154     //
00155     // find the x value
00156     //
00157     Double_t lowEdge = xValues[i]-errXLow[i];
00158     Double_t highEdge = xValues[i]+errXHigh[i];
00159 
00160 
00161     Double_t num = fMean->Integral(lowEdge,highEdge);
00162     Double_t denom = fUA1->Integral(lowEdge,highEdge);
00163 
00164     //
00165     // just use star's x
00166     //
00167 
00168     Double_t UA1ErrLow = fabs(lowEdge-xValues[i]);
00169     Double_t UA1ErrHigh= fabs(highEdge-xValues[i]);
00170 
00171     //Double_t y = fUA1->Eval(xValues[i]);
00172     double y = denom/(highEdge-lowEdge);
00173 
00174     // ratio
00175 
00176     //Int_t iBin = h1->GetXaxis()->FindBin(xValues[i]);
00177     //Double_t fraction = h1->GetBinContent(iBin);
00178     
00179     //    cout << "ratio : " << h1->GetXaxis()->GetBinLowEdge(iBin) << " " 
00180     // << h1->GetXaxis()->GetBinUpEdge(iBin) << " = " << fraction << endl;
00181 
00182     
00183     //Double_t yCorrected = fraction*y;
00184     double yCorrected = y;
00185 
00186     cout << "UA1 uncorrected " << y << ", corrected y " << yCorrected  
00187          << " graph y " << yValues[i] << endl;
00188 
00189     gUA1->SetPoint(i,xValues[i],yCorrected);
00190     gUA1->SetPointError(i,UA1ErrLow,UA1ErrHigh,0,0);
00191   }    
00192   
00193   gUA1->SetName("gUA1Spectrum");
00194   gUA1->SetTitle("gUA1Spectrum");
00195 
00196   return gUA1;
00197 
00198 }
00199 
00200 
00201 TGraphAsymmErrors*
00202 makeRatio(TGraphAsymmErrors* gStar, TGraphAsymmErrors* gUA1)
00203 {
00204   cout << "********************************" << endl;
00205   cout << gStar->GetName() << endl;
00206   const Int_t nBin = gStar->GetN();
00207 
00208   Double_t* starY = gStar->GetY();
00209   Double_t* starX = gStar->GetX();
00210   Double_t* starXErrLow = gStar->GetEXlow();
00211   Double_t* starXErrHigh = gStar->GetEXhigh();
00212   Double_t* starYErrHigh = gStar->GetEYhigh();
00213 
00214   Double_t* ua1Y  = gUA1->GetY();
00215   Double_t* ua1YErrHigh = gUA1->GetEYhigh();
00216 
00217   TGraphAsymmErrors* gRatio = new TGraphAsymmErrors;
00218 
00219   for(Int_t i=0; i<nBin; i++){
00220     
00221     Double_t ratio = starY[i]/ua1Y[i];
00222     
00223     Double_t yError = 
00224       TMath::Sqrt(ratio*ratio*(starYErrHigh[i]*starYErrHigh[i]/(starY[i]*starY[i])));
00225     //                         ua1YErrHigh[i]*ua1YErrHigh[i]/(ua1Y[i]*ua1Y[i])));
00226     
00227     cout << "star x : " << starX[i] << endl;
00228     cout << "star y : " << starY[i] << " ua1 Y : " << ua1Y[i]
00229          << " ratio  : " << ratio << endl
00230          << "error : " << yError << endl;
00231     
00232     gRatio->SetPoint(i,starX[i],ratio);
00233     gRatio->SetPointError(i,starXErrLow[i],starXErrHigh[i],yError,yError);
00234   }
00235 
00236   TString sName = gStar->GetName();
00237   sName.Append("UA1Ratio");
00238   
00239   gRatio->SetName(sName.Data());
00240   gRatio->SetTitle(sName.Data());
00241 
00242   return gRatio;
00243 
00244 }