Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:18 

0001 #include <iostream>
0002 #include <fstream>
0003 #include <iomanip>
0004 #include <string>
0005 #include <list>
0006 
0007 #include <math.h>
0008 #include <vector>
0009 
0010 #include "Rtypes.h"
0011 #include "TROOT.h"
0012 #include "TRint.h"
0013 #include "TObject.h"
0014 #include "TFile.h"
0015 #include "TTree.h"
0016 #include "TH1F.h"
0017 #include "TCanvas.h"
0018 #include "TApplication.h"
0019 #include "TRefArray.h"
0020 #include "TStyle.h"
0021 #include "TGraph.h"
0022 
0023 //const int modelsITEP=5, modelsBNL=5, nEnergy0=8, nEnergy1=7, nEnergy2=10, nEnergy3=8;
0024 //const int modelsITEP=2, modelsBNL=2, nEnergy0=8, nEnergy1=7, nEnergy2=10, nEnergy3=8;
0025 const int modelsITEP=6, modelsBNL=5, nEnergy0=8, nEnergy1=7, nEnergy2=10, nEnergy3=8;
0026 //const int modelsITEP=4, modelsBNL=4, nEnergy0=8, nEnergy1=7, nEnergy2=10, nEnergy3=8;
0027 //const int modelsITEP=5, modelsBNL=5, nEnergy0=8, nEnergy1=7, nEnergy2=10, nEnergy3=8;
0028 std::string ModelsI[6]      = {"LEP", "FTF", "Bertini", "Binary", "QGSC", "FTFP"};
0029 std::string ModelFilesI[6]  = {"LEP", "FTF", "Bertini", "Binary", "QGSC", "FTFP"};
0030 std::string ModelNamesI[6]  = {"LEP", "FTF-Binary", "Bertini", "Binary", "QGS-CHIPS", "FTF-Preco"};
0031 std::string ModelsB[6]      = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "FTFP"};
0032 std::string ModelFilesB[6]  = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "FTFP"};
0033 std::string ModelNamesB[6]  = {"LEP", "FTF-Binary", "Bertini", "QGS-Preco", "QGS-CHIPS", "FTF-Preco"};
0034 //std::string ModelsB[6]      = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0035 //std::string ModelFilesB[6]  = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0036 //std::string ModelNamesB[6]  = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0037 //std::string ModelsI[5]      = {"LEP", "FTF", "Bertini", "Binary", "QGSC"};
0038 //std::string ModelFilesI[5]  = {"LEP", "FTF", "Bertini", "Binary", "QGSC"};
0039 //std::string ModelNamesI[5]  = {"LEP", "FTF", "Bertini", "Binary", "QGSC"};
0040 //std::string ModelsI[6]      = {"Bertini1", "Bertini2", "Bertini3", "Bertini4", "Bertini5", "Bertini6"};
0041 //std::string ModelFilesI[6]  = {"Bertini1", "Bertini2", "Bertini3", "Bertini4", "Bertini5", "Bertini6"};
0042 //std::string ModelNamesI[6]  = {"Bertini (9.2.b01)", "Bertini (9.1.ref08)", "Bertini (+ Coulomb)", "Bertini (9.2)", "Bertini (9.2.ref02)", "Bertini (9.2.p01)"};
0043 //std::string ModelsI[2]      = {"Bertini", "Bertini"};
0044 //std::string ModelFilesI[2]  = {"OldBertini", "NewBertini"};
0045 //std::string ModelNamesI[2]  = {"Bertini (Old)", "Bertini (New)"};
0046 //std::string ModelsI[6]      = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0047 //std::string ModelFilesI[6]  = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0048 //std::string ModelNamesI[6]  = {"LEP", "FTF", "Bertini", "QGSP", "QGSC", "RPG"};
0049 //std::string ModelsI[4]      = {"LEP", "RPG", "FTF", "Bertini"};
0050 //std::string ModelFilesI[4]  = {"LEP", "RPG", "FTF", "Bertini"};
0051 //std::string ModelNamesI[4]  = {"LEP", "RPG", "FTF", "Bertini"};
0052 //std::string ModelsI[5]      = {"LEP", "Binary", "FTFP", "QGSC", "Bertini"};
0053 //std::string ModelFilesI[5]  = {"LEP", "Binary", "FTFP", "QGSC", "Bertini"};
0054 //std::string ModelNamesI[5]  = {"LEP", "Binary", "FTFP", "QGSC", "Bertini"};
0055 //std::string ModelsI[5]      = {"LEP", "QGSP", "FTFP", "QGSC", "Bertini"};
0056 //std::string ModelFilesI[5]  = {"LEP", "QGSP", "FTFP", "QGSC", "Bertini"};
0057 //std::string ModelNamesI[5]  = {"LEP", "QGSP", "FTFP", "QGSC", "Bertini"};
0058 int         colModel[6]    = {1, 2, 6, 4, 7, 11};
0059 int         symbModel[6]   = {24, 29, 25, 27, 26, 23};
0060 int         stylModel[6]   = {1, 2, 3, 4, 5, 6};
0061 double      keproton[4]    = {0.09, 0.15, 0.19, 0.23};
0062 double      keneutron[4]   = {0.07, 0.11, 0.15, 0.17};
0063 double      energyScan0[8] = {5.7, 6.2, 6.5, 7.0, 7.5, 8.2, 8.5, 9.0};
0064 double      energyScan1[7] = {6.0, 6.5, 7.0, 7.5, 8.2, 8.5, 9.0};
0065 double      energyScan2[10]= {1.0, 2.0, 3.0,  5.0, 6.0, 6.5,
0066                   7.0, 7.5, 8.25, 9.0};
0067 bool        debug=true;
0068 
0069 void plotData(char element[2], char ene[6], char angle[6], 
0070           char beam[8]="proton", char particle[8]="neutron", int save=0) {
0071 
0072   char file1[50], file2[50];
0073   sprintf (file1, "itep/%s/%s/%s%sGeV%sdeg.dat",  beam, particle, element, ene, angle);
0074   sprintf (file2, "itep/%s/%s/%s%sGeV%sdeg.dat2", beam, particle, element, ene, angle);
0075   cout << file1 << "  " << file2 << "\n";
0076   ifstream infile1, infile2;
0077   infile1.open(file1);
0078   infile2.open(file2);
0079   
0080   Int_t   q1, i=0;
0081   Float_t m1, r1, x1[30], y1[30], stater1[30], syser1[30];
0082   infile1 >> m1 >> r1 >> q1;
0083   for (i=0; i<q1; i++) infile1 >> x1[i] >> y1[i] >> stater1[i] >> syser1[i];
0084 
0085   Int_t   q2, n=0;
0086   Float_t m2, r2, x2[30], y2[30], stater2[30], syser2[30];
0087   infile2 >> m2 >> r2 >> q2;
0088   for (i=0; i<q2; i++) infile2 >> x2[i] >> y2[i] >> stater2[i] >> syser2[i];
0089 
0090   Float_t x[30], chi[30], dif[30], edif[30], chi2=0.0, diff=0., dify=0. ;
0091   for (i=0; i<q1; i++) {
0092     for (int j=0; j<q2; j++) {
0093       double dx = ((x1[i]-x2[j]) >= 0 ? (x1[i]-x2[j]): -(x1[i]-x2[j]));
0094       if (dx < 0.0001) {
0095     double d1 = stater1[i]/y1[i];
0096     double d2 = stater2[j]/y2[j];
0097     double dd = sqrt(stater1[i]*stater1[i] + stater2[j]*stater2[j]);
0098     x[n]    = x1[i];
0099     dif[n]  = 200.*(y1[i]-y2[j])/(y1[i]+y2[j]);
0100     double da = (dif[n] > 0 ? dif[n]: -dif[n]);
0101     edif[n] = da*sqrt(d1*d1+d2*d2);
0102     chi[n]  = pow(((y1[i]-y2[j])/dd),2);
0103     diff   += da;
0104     chi2   += chi[n];
0105     dify   += dif[n];
0106     n++;
0107     break;
0108       }
0109     }
0110   }
0111 
0112   for (i=0; i<n; i++)
0113     std::cout << "Data " << i << " E " << x[i] << " Difference " << dif[i] << " +- " << edif[i] << " Chi2 " << chi[i] << "\n";
0114 
0115   dify /= n;
0116   std::cout << "Chi-Square = " << chi2 << "/" << n << " Mean Difference = " << diff/n << "\n";
0117 
0118   setStyle(); gStyle->SetOptStat(1111);
0119 
0120   char name[30], title[60];
0121   sprintf (name, "%s%sGeV%sdeg", element, ene, angle);
0122   sprintf (title, "%s from p+%s at %s GeV (#theta = %s^{o})", particle, element, ene, angle);
0123   TCanvas* c1  = new TCanvas("c1",name,400,300); c1->SetLeftMargin(0.15);
0124   TGraph*  gr1 = new TGraphErrors(q1,x1,y1,0,stater1);
0125   gr1->SetTitle(""); gr1->SetMarkerColor(4);  // blue
0126   gr1->SetMarkerStyle(22);  gr1->SetMarkerSize(1.4); gr1->Draw("ALP");
0127   gr1->GetXaxis()->SetTitle("Energy (GeV)"); 
0128   gr1->GetYaxis()->SetTitle("E#frac{d^{3}#sigma}{dp^{3}} (mb/GeV^{2})"); 
0129 
0130   TGraph* gr2 = new TGraphErrors(q2,x2,y2,0,stater2);
0131   gr2->SetMarkerColor(2);  // red
0132   gr2->SetMarkerStyle(23);  gr2->SetMarkerSize(1.4); gr2->Draw("LP");
0133 
0134   TLegend *leg1 = new TLegend(0.55,0.80,0.90,0.90);
0135   leg1->SetHeader(title); leg1->SetFillColor(0);
0136   leg1->SetTextSize(0.04);
0137   leg1->Draw();
0138 
0139   sprintf (name, "Chi%s%sGeV%sdeg", element, ene, angle);
0140   TCanvas *c2  = new TCanvas("c2",name,400,300); c2->SetLeftMargin(0.15);
0141   TGraph  *gr3 = new TGraph(n,x,chi);
0142   gr3->SetTitle("");  gr3->SetMarkerStyle(22);  gr3->SetMarkerColor(4);
0143   gr3->GetXaxis()->SetTitle("Energy (GeV)"); gr3->GetYaxis()->SetTitle("#chi^{2}");  
0144   gr3->Draw("ALP"); gr3->SetMarkerSize(1.4);
0145   leg1->Draw();
0146 
0147   TGraph*  gr4 = new TGraphErrors(n,x,dif,0,edif);
0148   gr4->SetTitle("");  gr4->SetMarkerStyle(20); gr4->SetMarkerSize(1.25);
0149   gr4->SetMarkerColor(6);
0150   gr4->GetYaxis()->SetRangeUser(-25.,25.);
0151   gr4->GetXaxis()->SetTitle("Energy (GeV)"); 
0152   gr4->GetYaxis()->SetTitle("Difference (%)");
0153   double xmin = gr4->GetXaxis()->GetXmin();
0154   double xmax = gr4->GetXaxis()->GetXmax();
0155   if (debug) std::cout << " Xmin " << xmin << " " << xmax << "\n";
0156   sprintf (name, "Diff%s%sGeV%sdeg", element, ene, angle);
0157   TCanvas *c3  = new TCanvas("c3",name,400,300); c3->SetLeftMargin(0.15);
0158   TLine   *line = new TLine(xmin,dify,xmax,dify);
0159   line->SetLineStyle(2); line->SetLineColor(4);
0160   gr4->Draw("AP"); line->Draw();
0161   leg1->Draw();
0162 
0163   if (save > 0) {
0164     char fname[60];
0165     sprintf (fname, "%s%sGeV%sdeg_1.eps",  element, ene, angle);
0166     c1->SaveAs(fname);
0167     sprintf (fname, "%s%sGeV%sdeg_2.eps",  element, ene, angle);
0168     c2->SaveAs(fname);
0169     sprintf (fname, "%s%sGeV%sdeg_3.eps",  element, ene, angle);
0170     c3->SaveAs(fname);
0171   }
0172 }
0173 
0174 void plotKEx(char ene[6], char angle[6], int first=0, int logy=0, int save=0, 
0175          double ymin=-1., char particle[8]="proton", char beam[8]="proton",
0176          bool ratio='false', int leg1=1, int leg2=1, char dir[20]=".", 
0177          char dird[40]=".", char mark=' ') {
0178 
0179   setStyle();
0180   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0181 
0182   int leg=leg1; if (leg2 == 0) leg=leg2;
0183   char markf[4]=" ";
0184   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0185   if (mark == 'y') sprintf(markf, "(a)");
0186   if (ratio)
0187     plotKERatio("C", ene,angle,first,logy,ymin,particle,beam,leg1,dir,dird,markf);
0188   else
0189     plotKE("C", ene,angle,first,logy,ymin,particle,beam, leg1, dir,dird,markf);
0190   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0191   if (mark == 'y') sprintf(markf, "(b)");
0192   if (ratio)
0193     plotKERatio("Cu",ene,angle,first,logy,ymin,particle,beam,leg2,dir,dird,markf);
0194   else
0195     plotKE("Cu",ene,angle,first,logy,ymin,particle,beam,leg2,dir,dird,markf);
0196   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0197   if (mark == 'y') sprintf(markf, "(c)");
0198   if (ratio)
0199     plotKERatio("Pb",ene,angle,first,logy,ymin,particle,beam,leg, dir,dird,markf);
0200   else
0201     plotKE("Pb",ene,angle,first,logy,ymin,particle,beam,leg, dir,dird,markf);
0202   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0203   if (mark == 'y') sprintf(markf, "(d)");
0204   if (ratio)
0205     plotKERatio("U", ene,angle,first,logy,ymin,particle,beam,leg, dir,dird,markf);
0206   else
0207     plotKE("U", ene,angle,first,logy,ymin,particle,beam,leg, dir,dird,markf);
0208 
0209   char anglx[6], fname[60];
0210   int nx = 0;
0211   for (int i=0; i<6; i++) {
0212     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0213   }
0214   if (save != 0) {
0215     std::string tag=".gif";
0216     if (ratio) {
0217       if (save > 0) tag = "R.eps";
0218       else          tag = "R.gif";
0219     } else {
0220       if (save > 0) tag = ".eps";
0221     }
0222     sprintf (fname, "%sCCuPbUto%sat%sGeV%sdeg%s", beam, particle, ene, anglx, tag.c_str());
0223     myc->SaveAs(fname);
0224   }
0225 }
0226 
0227 
0228 void plotKEn(char ene[6], int first=0, int logy=0, int save=0, double ymin=-1.,
0229          char beam[8]="proton", bool ratio=false, int leg1=1, int leg2=1,
0230          char dir[20]=".", char dird[40]=".", char mark=' ') {
0231 
0232   setStyle();  
0233   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0234 
0235   int leg=leg1; if (leg2 == 0) leg=leg2;
0236   char markf[4]=" ";
0237   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0238   if (mark == 'y') sprintf(markf, "(a)");
0239   if (ratio)
0240     plotKERatio("C","1.4","119.0",first,logy,ymin,"neutron",beam,leg1,dir,dird,markf);
0241   else
0242     plotKE("C","1.4","119.0",first,logy,ymin,"neutron",beam,leg1,dir,dird,markf);
0243   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0244   if (mark == 'y') sprintf(markf, "(b)");
0245   if (ratio)
0246     plotKERatio("C",ene,   "119.0",first,logy,ymin,"neutron",beam,leg2,dir,dird,markf);
0247   else
0248     plotKE("C",ene,   "119.0",first,logy,ymin,"neutron",beam,leg2,dir,dird,markf);
0249   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0250   if (mark == 'y') sprintf(markf, "(c)");
0251   if (ratio)
0252     plotKERatio("U","1.4","119.0",first,logy,ymin,"neutron",beam,leg,dir,dird,markf);
0253   else
0254     plotKE("U","1.4","119.0",first,logy,ymin,"neutron",beam,leg,dir,dird,markf);
0255   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0256   if (mark == 'y') sprintf(markf, "(d)");
0257   if (ratio)
0258     plotKERatio("U",ene,   "119.0",first,logy,ymin,"neutron",beam,leg,dir,dird,markf);
0259   else
0260     plotKE("U",ene,   "119.0",first,logy,ymin,"neutron",beam,leg,dir,dird,markf);
0261   
0262   char fname[60];
0263   if (save != 0) {
0264     std::string tag=".gif";
0265     if (ratio) {
0266       if (save > 0) tag = "R.eps";
0267       else          tag = "R.gif";
0268     } else {
0269       if (save > 0) tag = ".eps";
0270     }
0271     sprintf (fname, "%sCUtoneutron_1%s", beam, tag.c_str());
0272     myc->SaveAs(fname);
0273   }
0274 
0275 }
0276 
0277 void plotKEp(char element[2], char ene[6], int first=0, int logy=0, int save=0,
0278          double ymin=-1., char beam[8]="proton", bool ratio=false, 
0279          int leg1=1, int leg2=1, char dir[20]=".", char dird[40]=".", 
0280          char mark=' ') {
0281 
0282   setStyle();  
0283   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0284 
0285   int leg=leg1; if (leg2 == 0) leg=leg2;
0286   char markf[4]=" ";
0287   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0288   if (mark == 'y') sprintf(markf, "(a)");
0289   if (ratio)
0290     plotKERatio(element,"1.4"," 59.1",first,logy,ymin,"proton",beam,leg1,dir,dird,markf);
0291   else
0292     plotKE(element,"1.4"," 59.1",first,logy,ymin,"proton",beam,leg1,dir,dird,markf);
0293   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0294   if (mark == 'y') sprintf(markf, "(b)");
0295   if (ratio)
0296     plotKERatio(element,ene,   " 59.1",first,logy,ymin,"proton",beam,leg2,dir,dird,markf);
0297   else
0298     plotKE(element,ene,   " 59.1",first,logy,ymin,"proton",beam,leg2,dir,dird,markf);
0299   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0300   if (mark == 'y') sprintf(markf, "(c)");
0301   if (ratio)
0302     plotKERatio(element,"1.4","119.0",first,logy,ymin,"proton",beam,leg,dir,dird,markf);
0303   else
0304     plotKE(element,"1.4","119.0",first,logy,ymin,"proton",beam,leg,dir,dird,markf);
0305   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0306   if (mark == 'y') sprintf(markf, "(d)");
0307   if (ratio)
0308     plotKERatio(element,ene,   "119.0",first,logy,ymin,"proton",beam,leg,dir,dird,markf);
0309   else
0310     plotKE(element,ene,   "119.0",first,logy,ymin,"proton",beam,leg,dir,dird,markf);
0311 
0312   char fname[60];
0313   if (save != 0) {
0314     std::string tag=".gif";
0315     if (ratio) {
0316       if (save > 0) tag = "R.eps";
0317       else          tag = "R.gif";
0318     } else {
0319       if (save > 0) tag = ".eps";
0320     }
0321     sprintf (fname, "%s%stoproton_1%s", beam, element, tag.c_str());
0322     myc->SaveAs(fname);
0323   }
0324 
0325 }
0326 
0327 void plotKE4(char element[2], char ene[6], int first=0, int logy=0, int save=0,
0328          double ymin=-1., char particle[8]="proton", char beam[8]="proton",
0329          bool ratio=false, int leg1=1, int leg2=1, char dir[20]=".", 
0330          char dird[40]=".", char mark=' ') {
0331 
0332   setStyle();  
0333   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0334 
0335   int leg=leg1; if (leg2 == 0) leg=leg2;
0336   char markf[4]=" ";
0337   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0338   if (mark == 'y') sprintf(markf, "(a)");
0339   if (ratio)
0340     plotKERatio(element,ene," 59.1",first,logy,ymin,particle,beam,leg1,dir,dird,markf);
0341   else
0342     plotKE(element,ene," 59.1",first,logy,ymin,particle,beam,leg1,dir,dird,markf);
0343   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0344   if (mark == 'y') sprintf(markf, "(b)");
0345   if (ratio)
0346     plotKERatio(element,ene," 89.0",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
0347   else
0348     plotKE(element,ene," 89.0",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
0349   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0350   if (mark == 'y') sprintf(markf, "(c)");
0351   if (ratio)
0352     plotKERatio(element,ene,"119.0",first,logy,ymin,particle,beam,leg,dir,dird,markf);
0353   else
0354     plotKE(element,ene,"119.0",first,logy,ymin,particle,beam,leg,dir,dird,markf);
0355   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0356   if (mark == 'y') sprintf(markf, "(d)");
0357   if (ratio)
0358     plotKERatio(element,ene,"159.6",first,logy,ymin,particle,beam,leg,dir,dird,markf);
0359   else
0360     plotKE(element,ene,"159.6",first,logy,ymin,particle,beam,leg,dir,dird,markf);
0361 
0362   char fname[60];
0363   if (save != 0) {
0364     std::string tag=".gif";
0365     if (ratio) {
0366       if (save > 0) tag = "R.eps";
0367       else          tag = "R.gif";
0368     } else {
0369       if (save > 0) tag = ".eps";
0370     }
0371     sprintf (fname, "%s%sto%sat%sGeV_1%s", beam, element, particle, ene, tag.c_str());
0372     myc->SaveAs(fname);
0373   }
0374 
0375 }
0376 
0377 void plotKE1(char element[2], char ene[6], char angle[6], int first=0, 
0378          int logy=0, int save=0, double ymin=-1, char particle[8]="proton",
0379          char beam[8]="proton", bool ratio=false, int legend=1, 
0380          char dir[20]=".", char dird[40]=".", char markf[4]=" ") {
0381 
0382   setStyle();
0383   TCanvas *myc = new TCanvas("myc","",800,600); myc->SetLeftMargin(0.15);
0384   if (logy != 0) gPad->SetLogy(1);
0385   if (ratio) 
0386     plotKERatio(element,ene,angle,first,logy,ymin,particle,beam,legend,dir,dird,markf);
0387   else
0388     plotKE(element,ene,angle,first,logy,ymin,particle,beam,legend,dir,dird,markf);
0389 
0390   char anglx[6], fname[120];
0391   int nx = 0;
0392   for (int i=0; i<6; i++) {
0393     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0394   }
0395   if (save != 0) {
0396     std::string tag=".gif";
0397     if (ratio) {
0398       if (save > 0) tag = "R.eps";
0399       else          tag = "R.gif";
0400     } else {
0401       if (save > 0) tag = ".eps";
0402     }
0403     sprintf (fname, "%s%sto%sat%sGeV%sdeg%s", beam, element, particle, ene, anglx, tag.c_str());
0404     myc->SaveAs(fname);
0405   }
0406 }
0407 
0408 void plotKE(char element[2], char ene[6], char angle[6], int first=0, 
0409         int logy=0, double ymin=-1, char particle[8]="proton",
0410             char beam[8]="proton", int legend=1, char dir[20]=".", 
0411         char dird[40]=".", char markf[4]=" ") {
0412 
0413   char fname[120], list[40], hname[60], titlx[50];
0414   TH1F *hi[6];
0415   int i=0, icol=1;
0416   sprintf (titlx, "Kinetic Energy of %s (GeV)", particle);
0417   double  ymx0=1, ymi0=100., xlow=0.06, xhigh=0.26;
0418   if (particle == "neutron") {xlow= 0.0; xhigh=0.20;}
0419   for (i=0; i<modelsITEP; i++) {
0420     sprintf (list, "%s", ModelFilesI[i].c_str()); 
0421     sprintf (fname, "%s/%s/%s/%s%s%sGeV_1.root", dir, beam, particle, element, list, ene);
0422     sprintf (list, "%s", ModelsI[i].c_str()); icol = colModel[i];
0423     sprintf (hname, "KE0%s%s%sGeV%s", element, list, ene, angle);
0424     TFile *file = new TFile(fname);
0425     hi[i] = (TH1F*) file->Get(hname);
0426     std::cout << "Get " << hname << " from " << fname <<" as " << hi[i] <<"\n";
0427     int nx = hi[i]->GetNbinsX();
0428     for (int k=1; k <= nx; k++) {
0429       double xx = hi[i]->GetBinCenter(k);
0430       double yy = hi[i]->GetBinContent(k);
0431       if (xx > xlow && xx < xhigh) {
0432     if (yy > ymx0) ymx0 = yy;
0433     if (yy < ymi0 && yy > 0) ymi0 = yy;
0434       }
0435     }
0436     hi[i]->GetXaxis()->SetRangeUser(xlow, xhigh); hi[i]->SetTitle("");
0437     hi[i]->GetXaxis()->SetTitle(titlx);
0438     hi[i]->SetLineStyle(1);  hi[i]->SetLineWidth(2); hi[i]->SetLineColor(icol);
0439     //    file->Close();
0440   }
0441 
0442   char anglx[6];
0443   int nx = 0;
0444   for (i=0; i<6; i++) {
0445     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0446   }
0447   sprintf (fname, "%s/itep/%s/%s/%s%sGeV%sdeg.dat", dird, beam, particle, element, ene, anglx);
0448   std::cout << "Reads data from file " << fname << "\n";
0449   ifstream infile;
0450   infile.open(fname);
0451   
0452   int     q1;
0453   float   m1, r1, x1[30], y1[30], stater1[30], syser1[30];
0454   infile >> m1 >> r1 >> q1;
0455   for (i=0; i<q1; i++) {
0456     infile >> x1[i] >> y1[i] >> stater1[i] >> syser1[i];
0457     syser1[i] *= y1[i];
0458     double err = sqrt(syser1[i]*syser1[i]+stater1[i]*stater1[i]);
0459     stater1[i] = err;
0460     if (y1[i]+stater1[i] > ymx0) ymx0 = y1[i]+stater1[i];    
0461     if (y1[i]-stater1[i] < ymi0 && y1[i]-stater1[i] > 0) ymi0=y1[i]-stater1[i];
0462     if (debug) std::cout << i << " " << x1[i] << " " << y1[i] << " " << stater1[i] << "\n";
0463   }
0464   TGraph*  gr1 = new TGraphErrors(q1,x1,y1,0,stater1);
0465   gr1->SetMarkerColor(4);  gr1->SetMarkerStyle(22);
0466   gr1->SetMarkerSize(1.6);
0467 
0468   if (logy == 0) {ymx0 *= 1.5; ymi0 *= 0.8;}
0469   else           {ymx0 *=10.0; ymi0 *= 0.2; }
0470   if (ymin > 0) ymi0 = ymin;
0471   for (i = 0; i<modelsITEP; i++) {
0472     if (debug) std::cout << "Model " << i << " " << hi[i] << " " << ymi0 << " " << ymx0 << "\n";
0473     hi[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
0474   }
0475 
0476   hi[first]->GetYaxis()->SetTitleOffset(1.6);
0477   hi[first]->Draw();
0478   for (i=0; i<modelsITEP; i++) {
0479     if (i != first) hi[i]->Draw("same");
0480   }
0481   gr1->Draw("p");
0482 
0483   TLegend *leg1;
0484   if (legend < 0) {
0485     leg1 = new TLegend(0.60,0.55,0.90,0.90);
0486   } else {
0487     if (markf == " ") leg1 = new TLegend(0.42,0.55,0.90,0.90);
0488     else              leg1 = new TLegend(0.38,0.70,0.90,0.90);
0489   }
0490   for (i=0; i<modelsITEP; i++) {
0491     sprintf (list, "%s", ModelNamesI[i].c_str()); 
0492     leg1->AddEntry(hi[i],list,"F");
0493   }
0494   char header[120], beamx[8], partx[2];
0495   if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
0496   else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
0497   else                        sprintf (beamx, "p");
0498   if      (particle == "neutron") sprintf (partx, "n");
0499   else                            sprintf (partx, "p");
0500   if (legend < 0) {
0501     sprintf (header,"%s+A #rightarrow %s+X", beamx, partx);
0502   } else {
0503     if (markf == " ") 
0504       sprintf (header,"%s+%s #rightarrow %s+X at %s GeV (#theta = %s^{o})", beamx, element, partx, ene, angle);
0505     else 
0506       sprintf (header,"%s %s+%s #rightarrow %s+X at %s GeV (#theta = %s^{o})", markf, beamx, element, partx, ene, angle);
0507   }
0508   leg1->SetHeader(header); leg1->SetFillColor(0);
0509   leg1->SetTextSize(0.04);
0510   if (legend != 0) leg1->Draw("same");
0511   if (debug) std::cout << "End\n";
0512 }
0513 
0514 void plotKERatio(char element[2], char ene[6], char angle[6], int first=0, 
0515          int logy=0, double ymin=-1, char particle[8]="proton",
0516          char beam[8]="proton", int legend=1, char dir[20]=".", 
0517          char dird[40]=".", char markf[4]=" ") {
0518 
0519   // First open the data file
0520   char anglx[6], fname[120];
0521   int nx = 0;
0522   for (int i=0; i<6; i++) {
0523     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0524   }
0525   sprintf (fname, "%s/itep/%s/%s/%s%sGeV%sdeg.dat", dird, beam, particle, element, ene, anglx);
0526   if (debug) std::cout << "Reads data from file " << fname << "\n";
0527   ifstream infile;
0528   infile.open(fname);
0529 
0530   // Read contents of the data file
0531   int     q1;
0532   float   m1, r1, x1[30], y1[30], er1[30], staterr, syserr;
0533   infile >> m1 >> r1 >> q1;
0534   for (i=0; i<q1; i++) {
0535     infile >> x1[i] >> y1[i] >> staterr >> syserr;
0536     syserr *= y1[i];
0537     er1[i]  = sqrt(syserr*syserr+staterr*staterr);
0538     if (debug) std::cout << i << " " << x1[i] << " " << y1[i] << " " << er1[i] << "\n";
0539   }
0540 
0541   char list[40], hname[60], titlx[50];
0542   TGraphErrors *gr[6];
0543   int icol=1, ityp=20;
0544   sprintf (titlx, "Kinetic Energy of %s (GeV)", particle);
0545   double  ymx0=0.1, ymi0=100., xlow=0.06, xhigh=0.26;
0546   if (particle == "neutron") {xlow= 0.0; xhigh=0.20;}
0547   for (int i=0; i<modelsITEP; i++) {
0548     icol = colModel[i]; ityp = symbModel[i];
0549     sprintf (list, "%s", ModelFilesI[i].c_str()); 
0550     sprintf (fname, "%s/%s/%s/%s%s%sGeV_1.root", dir, beam, particle, element, list, ene);
0551     sprintf (list, "%s", ModelsI[i].c_str()); 
0552     sprintf (hname, "KE0%s%s%sGeV%s", element, list, ene, angle);
0553 
0554     TFile *file = new TFile(fname);
0555     TH1F *hi = (TH1F*) file->Get(hname);
0556     if (debug) std::cout << "Get " << hname << " from " << fname <<" as " << hi <<"\n";
0557             
0558     if (hi != 0 && q1 > 0) {
0559       float xx[30], dx[30], rat[30], drt[30];
0560       int   nx = hi->GetNbinsX();
0561       int   np = 0;
0562       if (debug) std::cout << "Start with " << nx << " bins\n";
0563       for (int k=1; k <= nx; k++) {
0564     double xx1 = hi->GetBinLowEdge(k);
0565     double xx2 = hi->GetBinWidth(k);
0566     for (int j=0; j<q1; j++) {
0567       if (xx1 < x1[j] && xx1+xx2 > x1[j]) {
0568         double yy = hi->GetBinContent(k);
0569         xx[np]    = x1[j];
0570         dx[np]    = 0;
0571         rat[np]   = yy/y1[j];
0572         drt[np]   = er1[j]*rat[j]/y1[j];
0573         if (xx[np] > xlow && xx[np] < xhigh) {
0574           if (rat[np]+drt[np] > ymx0) ymx0 = rat[np]+drt[np];
0575           if (rat[np]-drt[np] < ymi0) ymi0 = rat[np]-drt[np];
0576         }
0577         if (debug) std::cout << np << "/" << j << "/" << k << " x " << xx[np] << " (" << xx1 << ":" << xx1+xx2 << ")" << " y " << yy << "/" << y1[j] << " = " << rat[np] << " +- " << drt[np] << "\n";
0578         np++;
0579         break;
0580       }
0581     }
0582       }
0583       gr[i] = new TGraphErrors(np, xx, rat, dx, drt);
0584       gr[i]->GetXaxis()->SetRangeUser(xlow, xhigh); gr[i]->SetTitle("");
0585       gr[i]->GetXaxis()->SetTitle(titlx);
0586       gr[i]->GetYaxis()->SetTitle("MC/Data");
0587       gr[i]->SetLineStyle(stylModel[i]); gr[i]->SetLineWidth(2); 
0588       gr[i]->SetLineColor(icol);         gr[i]->SetMarkerColor(icol); 
0589       gr[i]->SetMarkerStyle(ityp);       gr[i]->SetMarkerSize(1.0); 
0590     } else {
0591       gr[i] = 0;
0592     }
0593     file->Close();
0594   }
0595 
0596   if (logy == 0) {ymx0 *= 1.5; ymi0 *= 0.8;}
0597   else           {ymx0 *=10.0; ymi0 *= 0.2; }
0598   if (ymin > 0)   ymi0 = ymin;
0599   for (i = 0; i<modelsITEP; i++) {
0600     if (debug) std::cout << "Model " << i << " " << gr[i] << " " << ymi0 << " " << ymx0 << "\n";
0601     if (gr[i] != 0) gr[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
0602   }
0603 
0604   gr[first]->GetYaxis()->SetTitleOffset(1.6);
0605   gr[first]->Draw("APl");
0606   for (i=0; i<modelsITEP; i++) {
0607     if (i != first && gr[i] != 0) gr[i]->Draw("Pl");
0608   }
0609 
0610   TLegend *leg1;
0611   if (legend < 0) {
0612     leg1 = new TLegend(0.60,0.55,0.90,0.90);
0613   } else {
0614     if (markf == " ") leg1 = new TLegend(0.42,0.55,0.90,0.90);
0615     else              leg1 = new TLegend(0.38,0.70,0.90,0.90);
0616   }
0617   for (i=0; i<modelsITEP; i++) {
0618     if (gr[i] != 0) {
0619       sprintf (list, "%s", ModelNamesI[i].c_str()); 
0620       leg1->AddEntry(gr[i],list,"lP");
0621     }
0622   }
0623   char header[120], beamx[8], partx[2];
0624   if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
0625   else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
0626   else                        sprintf (beamx, "p");
0627   if      (particle == "neutron") sprintf (partx, "n");
0628   else                            sprintf (partx, "p");
0629   if (legend < 0) {
0630     sprintf (header,"%s+A #rightarrow %s+X", beamx, partx);
0631   } else {
0632     if (markf == " ") 
0633       sprintf (header,"%s+%s #rightarrow %s+X at %s GeV (#theta = %s^{o})", beamx, element, partx, ene, angle);
0634     else 
0635       sprintf (header,"%s %s+%s #rightarrow %s+X at %s GeV (#theta = %s^{o})", markf, beamx, element, partx, ene, angle);
0636   }
0637   leg1->SetHeader(header); leg1->SetFillColor(0);
0638   leg1->SetTextSize(0.04);
0639   if (legend != 0) leg1->Draw("same");
0640 
0641   xx[0]=xlow; xx[1]=xhigh; rat[0]=rat[1]=1.0;
0642   TGraph *gr0 = new TGraph(2, xx, rat);
0643   gr0->GetXaxis()->SetRangeUser(xlow, xhigh); gr0->SetTitle("");
0644   gr0->SetLineStyle(1);   gr0->SetLineWidth(1.4); 
0645   gr0->SetLineColor(1);   gr0->SetMarkerColor(1); 
0646   gr0->SetMarkerStyle(20);gr0->SetMarkerSize(1.6);
0647   gr0->Draw("l");
0648 }
0649 
0650 void plotCT4(char element[2], char ene[6], int first=0, int scan=1, int logy=0,
0651          int save=0, char particle[8]="proton", char beam[8]="proton", 
0652          int leg1=1, int leg2=1, char dir[20]=".", char dird[40]=".") {
0653 
0654   setStyle();
0655   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0656 
0657   for (int i=0; i<4; i++) {
0658     myc->cd(i+1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0659     double ke = keproton[i];
0660     if (particle == "neutron") ke = keneutron[i];
0661     if (i == 0) 
0662       plotCT(element, ene,ke, first, scan, logy, particle,beam,leg1,dir,dird); 
0663     else
0664       plotCT(element, ene,ke, first, scan, logy, particle,beam,leg2,dir,dird); 
0665   }
0666 
0667   char fname[60];
0668   if (save != 0) {
0669     if (save > 0) sprintf (fname, "%s%sto%sat%sGeV_2.eps", beam, element, particle, ene);
0670     else          sprintf (fname, "%s%sto%sat%sGeV_2.gif", beam, element, particle, ene);
0671     myc->SaveAs(fname);
0672   }
0673 }
0674  
0675 void plotCT1(char element[2], char ene[6], double ke, int first=0, int scan=1,
0676          int logy=0, int save=0, char particle[8]="proton", 
0677          char beam[8]="proton", int legend=1, char dir[20]=".",
0678          char dird[40]=".") {
0679 
0680   setStyle();
0681   TCanvas *myc = new TCanvas("myc","",800,600); myc->SetLeftMargin(0.15);
0682   if (logy != 0) gPad->SetLogy(1);
0683   plotCT(element, ene,ke, first, scan, logy, particle,beam, legend, dir,dird);
0684 
0685   char fname[60];
0686   if (save != 0) {
0687     if (save > 0) sprintf (fname, "%s%sto%sat%sGeV%4.2fGeV.eps", beam, element, particle, ene, ke);
0688     else          sprintf (fname, "%s%sto%sat%sGeV%4.2fGeV.gif", beam, element, particle, ene, ke);
0689     myc->SaveAs(fname);
0690   }
0691 }
0692 
0693 void plotCT(char element[2], char ene[6], double ke, int first=0, int scan=1,
0694         int logy=0, char particle[8]="proton", char beam[8]="proton", 
0695         int legend=1, char dir[20]=".", char dird[40]=".") {
0696 
0697   static double pi  = 3.1415926;
0698   static double deg = pi/180.; 
0699   if (debug) std::cout << "Scan " << scan;
0700   std::vector<double> angles = angleScan(scan);
0701   int    nn = (int)(angles.size());
0702   if (debug) std::cout << " gives " << nn << " angles\n";
0703 
0704   char fname[120], list[40], hname[60];
0705   TH1F *hi[6];
0706   int i=0, icol=1;
0707   double  ymx0=1, ymi0=100., xlow=-1.0, xhigh=1.0;
0708   for (i=0; i<modelsITEP; i++) {
0709     sprintf (list, "%s", ModelFilesI[i].c_str()); icol = colModel[i];
0710     sprintf (fname, "%s/%s/%s/%s%s%sGeV_1.root", dir, beam, particle, element, list, ene);
0711     sprintf (list, "%s", ModelsI[i].c_str()); icol = colModel[i];
0712     sprintf (hname, "CT0%s%s%sGeV%4.2f", element, list, ene, ke);
0713     TFile *file = new TFile(fname);
0714     hi[i] = (TH1F*) file->Get(hname);
0715     if (debug) std::cout << "Get " << hname << " from " << fname <<" as " << hi[i] <<"\n";
0716     int nx = hi[i]->GetNbinsX();
0717     for (int k=1; k <= nx; k++) {
0718       double xx = hi[i]->GetBinCenter(k);
0719       double yy = hi[i]->GetBinContent(k);
0720       if (xx > xlow && xx < xhigh) {
0721     if (yy > ymx0)           ymx0 = yy;
0722     if (yy < ymi0 && yy > 0) ymi0 = yy;
0723       }
0724     }
0725     if (debug) std::cout << "Y limit " << ymi0 << " " << ymx0 << " after " << i;
0726     hi[i]->GetXaxis()->SetRangeUser(xlow, xhigh); hi[i]->SetTitle("");
0727     hi[i]->SetLineStyle(1);  hi[i]->SetLineWidth(2); hi[i]->SetLineColor(icol);
0728     //    file->Close();
0729   }
0730 
0731   int     q1, kk0=0;
0732   float   m1, r1, x1[30], y1[30], stater1[30], syser1[30];
0733 
0734   for (int kk=0; kk<nn; kk++) {
0735     char angle[6], anglx[6];
0736     sprintf (angle, "%5.1f", angles[kk]);
0737     int nx = 0;
0738     for (i=0; i<6; i++) {
0739       if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0740     }
0741     sprintf (fname, "%s/itep/%s/%s/%s%sGeV%sdeg.dat", dird, beam, particle, element, ene, anglx);
0742     ifstream infile;
0743     infile.open(fname);
0744   
0745     infile >> m1 >> r1 >> q1;
0746     for (i=0; i<q1; i++) {
0747       float xx1, yy1, stater, syser;
0748       infile >> xx1 >> yy1 >> stater >> syser;
0749       if (xx1 > ke-0.001 && xx1 < ke+0.001) {
0750     x1[kk0] = cos(deg*angles[kk]);
0751     y1[kk0] = yy1; stater1[kk0] = stater; syser1[kk0] = syser;
0752     syser *= yy1;
0753     double err = sqrt(syser*syser+stater*stater);
0754     stater1[kk0] = err; 
0755     if (y1[kk0]+stater1[kk0] > ymx0) ymx0 = y1[kk0]+stater1[kk0];
0756     if (y1[kk0]-stater1[kk0] < ymi0 && y1[kk0]-stater1[kk0] > 0) ymi0 = y1[kk0]-stater1[kk0];
0757     kk0++;
0758       }
0759     }
0760     infile.close();
0761     if (debug) std::cout << kk << " File " << fname << " X " << x1[kk] << " Y " << y1[kk] << " DY " << stater1[kk] << "\n";
0762   }
0763 
0764   TGraph*  gr1 = new TGraphErrors(kk0,x1,y1,0,stater1);
0765   gr1->SetMarkerColor(4);  gr1->SetMarkerStyle(22);
0766   gr1->SetMarkerSize(1.6);
0767 
0768   if (logy == 0) {ymx0 *= 1.5; ymi0 *= 0.8;}
0769   else           {ymx0 *=10.0; ymi0 *= 0.2; }
0770   for (i = 0; i<modelsITEP; i++)
0771     hi[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
0772   
0773   hi[first]->GetYaxis()->SetTitleOffset(1.6);
0774   hi[first]->Draw();
0775   for (i=0; i<modelsITEP; i++) {
0776     if (i != first)  hi[i]->Draw("same");
0777   }
0778   gr1->Draw("p");
0779 
0780   TLegend *leg1;
0781   if (legend == 1) leg1 = new TLegend(0.15,0.70,0.62,0.90);
0782   else             leg1 = new TLegend(0.15,0.55,0.62,0.90);
0783   for (i=0; i<modelsITEP; i++) {
0784     sprintf (list, "%s", ModelNamesI[i].c_str());
0785     leg1->AddEntry(hi[i],list,"F");
0786   }
0787   char header[80], beamx[8], partx[2];
0788   if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
0789   else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
0790   else                        sprintf (beamx, "p");
0791   if      (particle == "neutron") sprintf (partx, "n");
0792   else                            sprintf (partx, "p");
0793   sprintf (header, "%s+%s #rightarrow %s+X at %s GeV (%4.2f GeV)", beamx, element, partx, ene, ke);
0794   leg1->SetHeader(header); leg1->SetFillColor(0);
0795   leg1->SetTextSize(0.04);
0796   if (legend != 0) leg1->Draw();
0797 
0798 }
0799 
0800 void plotBE4(char element[2], int logy=0, int scan=1, int save=0, 
0801          char particle[8]="proton", char beam[8]="proton", 
0802          char dir[20]=".", char dird[40]=".") {
0803 
0804   setStyle();
0805   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0806 
0807   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0808   plotBE(element, " 59.1", 0.11, logy, scan, particle, beam, dir, dird);
0809   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0810   plotBE(element, " 59.1", 0.21, logy, scan, particle, beam, dir, dird);
0811   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0812   plotBE(element, "119.0", 0.11, logy, scan, particle, beam, dir, dird);
0813   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
0814   plotBE(element, "119.0", 0.21, logy, scan, particle, beam, dir, dird);
0815 
0816   char fname[60];
0817   if (save != 0) {
0818     if (save > 0) sprintf (fname, "%s%sto%s_1.eps", beam, element, particle);
0819     else          sprintf (fname, "%s%sto%s_1.gif", beam, element, particle);
0820     myc->SaveAs(fname);
0821   }
0822 }
0823 
0824 void plotBE1(char element[2], char angle[6], double ke, int logy=0, int scan=1,
0825          int save=0, char particle[8]="proton", char beam[8]="proton", 
0826          char dir[20]=".", char dird[40]=".") {
0827 
0828   setStyle();
0829   TCanvas *myc = new TCanvas("myc","",800,600); myc->SetLeftMargin(0.15);
0830   if (logy != 0) gPad->SetLogy(1);
0831   plotBE(element, angle, ke, logy, scan, particle, beam, dir, dird);
0832 
0833   char anglx[6], fname[60];
0834   int i=0, nx=0;
0835   for (i=0; i<6; i++) {
0836     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0837   }
0838   if (save != 0) {
0839     if (save>0) sprintf (fname, "%s%sto%sat%sdeg%4.2fGeV.eps", beam, element, particle, anglx, ke);
0840     else        sprintf (fname, "%s%sto%sat%sdeg%4.2fGeV.gif", beam, element, particle, anglx, ke);
0841     myc->SaveAs(fname);
0842   }
0843 }
0844 
0845 void plotBE(char element[2], char angle[6], double ke, int logy=0, int scan=1,
0846         char particle[8]="proton", char beam[8]="proton", 
0847         char dir[20]=".", char dird[40]=".") {
0848 
0849   double ene[15];
0850   int    nene=0;
0851   if (scan == 0) {
0852     nene = nEnergy0;
0853     for (int i=0; i<nene; i++) ene[i] = energyScan0[i];
0854   } else if (scan <= 1) {
0855     nene = nEnergy1;
0856     for (int i=0; i<nene; i++) ene[i] = energyScan1[i];
0857   } else if (scan == 2) {
0858     nene = nEnergy2;
0859     for (int i=0; i<nene; i++) ene[i] = energyScan2[i];
0860   } else {
0861     nene = nEnergy3;
0862     for (int i=0; i<nene; i++) ene[i] = energyScan2[i];
0863   }
0864  
0865   char anglx[6];
0866   int i=0, nx=0;
0867   for (i=0; i<6; i++) {
0868     if (angle[i] != ' ') { anglx[nx] = angle[i]; nx++;}
0869   }
0870 
0871   TGraph *gr[4];
0872   char fname[120], list[40], hname[60];
0873   int j=0, icol=1, ityp=20;
0874   double  ymx0=1, ymi0=10000., xmi=5.0, xmx=10.0;
0875   if (scan > 1) { 
0876     xmi = 0.5;
0877     xmx = 9.5;
0878   }
0879   for (i=0; i<modelsITEP; i++) {
0880     icol = colModel[i]; ityp = symbModel[i];
0881     double yt[15];
0882     for (j=0; j<nene; j++) {
0883       sprintf (list, "%s", ModelFilesI[i].c_str()); 
0884       sprintf (fname, "%s/%s/%s/%s%s%3.1fGeV_1.root", dir, beam, particle, element, list, ene[j]);
0885       sprintf (list, "%s", ModelsI[i].c_str()); 
0886       sprintf (hname, "KE0%s%s%3.1fGeV%s", element, list, ene[j], angle);
0887       TFile *file = new TFile(fname);
0888       TH1F *hi = (TH1F*) file->Get(hname);
0889       if (hi == 0) {
0890     sprintf (hname, "KE0%s%s%4.2fGeV%s", element, list, ene[j], angle);
0891     hi = (TH1F*) file->Get(hname);
0892       }
0893       if (debug) std::cout << "Get " << hname << " from " << fname <<" as " << hi <<"\n";
0894       int    nk=0, nx = hi->GetNbinsX();
0895       double yy0=0;
0896       for (int k=1; k <= nx; k++) {
0897     double xx0 = hi->GetBinCenter(k);
0898     if (xx0 > ke-0.01 && xx0 < ke+0.01) {
0899       yy0 += hi->GetBinContent(k);
0900       nk++;
0901     }
0902       }
0903       if (nk > 0 )                 yy0 /= nk;
0904       if (yy0 > ymx0)              ymx0 = yy0;
0905       if (yy0 < ymi0 && yy0 > 0.1) ymi0 = yy0;
0906       if (debug) std::cout << hname << " # " << nk << " Y " << yy0 << " min " << ymi0 << " max " << ymx0 << "\n";
0907       yt[j] = yy0;
0908       file->Close();
0909     }
0910     gr[i] = new TGraph(nene, ene, yt); gr[i]->SetMarkerSize(1.2);
0911     gr[i]->SetTitle(list); gr[i]->SetLineColor(icol); 
0912     gr[i]->SetLineStyle(i+1); gr[i]->SetLineWidth(2);
0913     gr[i]->SetMarkerColor(icol);  gr[i]->SetMarkerStyle(ityp); 
0914     gr[i]->GetXaxis()->SetTitle("Beam Energy (GeV)");
0915     if (debug) {
0916       std::cout << "Graph " << i << " with " << nene << " points\n";
0917       for (j=0; j<nene; j++) std::cout << j << " x " << ene[j] << " y " << yt[j] << "\n";
0918     }
0919   }
0920 
0921   double ye[15], dy[15];
0922   for (j=0; j<nene; j++) {
0923     sprintf (fname, "%s/itep/%s/%s/%s%3.1fGeV%sdeg.dat", dird, beam, particle, element, ene[j], anglx);
0924     if (debug) std::cout << "Reads data from file " << fname << "\n";
0925     ifstream infile;
0926     infile.open(fname);
0927   
0928     int     q1;
0929     float   m1, r1, xx, yy, stater, syser;
0930     infile >> m1 >> r1 >> q1;
0931     for (i=0; i<q1; i++) {
0932       infile >> xx >> yy >> stater >> syser;
0933       if (xx > ke-0.01 && xx < ke+0.01) {
0934     ye[j] = yy;
0935     syser *= yy;
0936     double err = sqrt(syser*syser+stater*stater);
0937     dy[j] = err;
0938       }
0939     }
0940     infile.close();
0941     if (ye[j]+dy[j] > ymx0) ymx0 = ye[j]+dy[j];
0942     if (ye[j]-dy[j] < ymi0 && ye[j]-dy[j] > 0) ymi0 = ye[j]-dy[j];
0943   }
0944   if (debug) {
0945     std::cout << "Graph Data with " << nene << " points\n";
0946     for (j=0; j<nene; j++) std::cout << j << " x " << ene[j] << " y " << ye[j] << " +- " << dy[j] << "\n";
0947   }
0948   TGraph*  gr1 = new TGraphErrors(nene,ene,ye,0,dy);
0949   gr1->SetMarkerColor(1);  gr1->SetMarkerStyle(22);
0950   gr1->SetMarkerSize(1.6);
0951 
0952   if (logy == 0) {
0953     ymx0 *= 1.8; ymi0 *= 0.8;
0954   } else {
0955     ymx0 *= 50.0; ymi0 *= 0.2;
0956     if (scan > 1) ymx0 *= 4;
0957   }
0958   for (i = 0; i<modelsITEP; i++) {
0959     gr[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
0960     gr[i]->GetXaxis()->SetRangeUser(xmi,xmx);
0961   }
0962   gr1->GetXaxis()->SetRangeUser(xmi,xmx);
0963   gr1->GetYaxis()->SetRangeUser(ymi0,ymx0);
0964   gr1->GetXaxis()->SetTitle("Energy (GeV)"); 
0965   gr1->GetYaxis()->SetTitle("E#frac{d^{3}#sigma}{dp^{3}} (mb/GeV^{2})"); 
0966 
0967   gr1->GetYaxis()->SetTitleOffset(1.6); gr1->SetTitle("");
0968   gr1->Draw("ap");
0969   for (i=0; i<modelsITEP; i++)
0970     gr[i]->Draw("lp");
0971   
0972   TLegend *leg1 = new TLegend(0.35,0.60,0.90,0.90);
0973   for (i=0; i<modelsITEP; i++) {
0974     sprintf (list, "%s", ModelNamesI[i].c_str());
0975     leg1->AddEntry(gr[i],list,"LP");
0976   }
0977   char header[80], beamx[8], partx[2];
0978   if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
0979   else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
0980   else                        sprintf (beamx, "p");
0981   if      (particle == "neutron") sprintf (partx, "n");
0982   else                            sprintf (partx, "p");
0983   sprintf (header, "%s+%s #rightarrow %s+X at (KE = %3.1f GeV, #theta = %s^{o})", beamx, element, partx, ke, angle);
0984   leg1->SetHeader(header); leg1->SetFillColor(0);
0985   leg1->SetTextSize(0.04);
0986   leg1->Draw();
0987 }
0988  
0989 void plotMT4(char ene[6], int first=0, int logy=0, int save=0, double ymin=-1,
0990          char particle[8]="piplus", char beam[8]="proton",bool ratio=false,
0991          int leg1=1, int leg2=1, char dir[20]=".", char dird[40]=".", 
0992          char mark=' ') {
0993 
0994   setStyle();
0995   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
0996 
0997   char markf[4]=" ";
0998   int leg=leg1; if (leg2 == 0) leg=leg2;
0999   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1000   if (mark == 'y') sprintf(markf, "(a)");
1001   if (ratio)
1002     plotMTRatio("Be",ene,"1.10", first,logy,ymin,particle,beam,leg1,dir,dird,markf);
1003   else
1004     plotMT("Be",ene,"1.10", first,logy,ymin,particle,beam,leg1,dir,dird,markf);
1005   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1006   if (mark == 'y') sprintf(markf, "(b)");
1007   if (ratio)
1008     plotMTRatio("Be",ene,"2.30", first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1009   else
1010     plotMT("Be",ene,"2.30", first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1011   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1012   if (mark == 'y') sprintf(markf, "(c)");
1013   if (ratio)
1014     plotMTRatio("Au",ene,"1.10", first,logy,ymin,particle,beam,leg,dir,dird,markf);
1015   else
1016     plotMT("Au",ene,"1.10", first,logy,ymin,particle,beam,leg,dir,dird,markf);
1017   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1018   if (mark == 'y') sprintf(markf, "(d)");
1019   if (ratio)
1020     plotMTRatio("Au",ene,"2.30", first,logy,ymin,particle,beam,leg,dir,dird,markf);
1021   else
1022     plotMT("Au",ene,"2.30", first,logy,ymin,particle,beam,leg,dir,dird,markf);
1023 
1024   char fname[60];
1025   if (save != 0) {
1026     std::string tag=".gif";
1027     if (ratio) {
1028       if (save > 0) tag = "R.eps";
1029       else          tag = "R.gif";
1030     } else {
1031       if (save > 0) tag = ".eps";
1032     }
1033     sprintf (fname, "%sBeAuto%sat%sGeV%s", beam, particle, ene, tag.c_str());
1034     myc->SaveAs(fname);
1035   }
1036 }
1037  
1038 void plotMT4(char element[2], char ene[6], int first=0, int logy=0, int save=0,
1039          double ymin=-1, char particle[8]="piplus", char beam[8]="proton", 
1040          bool ratio=false, int leg1=1, int leg2=1, char dir[20]=".", 
1041          char dird[40]=".", char mark=' ') {
1042 
1043   setStyle();
1044   TCanvas *myc = new TCanvas("myc","",800,600); myc->Divide(2,2);
1045 
1046   char markf[4]=" ";
1047   int leg=leg1; if (leg2 == 0) leg=leg2;
1048   myc->cd(1); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1049   if (mark == 'y') sprintf(markf, "(a)");
1050   if (ratio)
1051     plotMTRatio(element,ene,"1.10",first,logy,ymin,particle,beam,leg1,dir,dird,markf);
1052   else
1053     plotMT(element,ene,"1.10",first,logy,ymin,particle,beam,leg1,dir,dird,markf);
1054   myc->cd(2); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1055   if (mark == 'y') sprintf(markf, "(b)");
1056   if (ratio)
1057     plotMTRatio(element,ene,"1.50",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1058   else
1059     plotMT(element,ene,"1.50",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1060   myc->cd(3); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1061   if (mark == 'y') sprintf(markf, "(c)");
1062   if (ratio)
1063     plotMTRatio(element,ene,"1.90",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1064   else
1065     plotMT(element,ene,"1.90",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1066   myc->cd(4); if (logy != 0) gPad->SetLogy(1); gPad->SetLeftMargin(0.15);
1067   if (mark == 'y') sprintf(markf, "(d)");
1068   if (ratio)
1069     plotMTRatio(element,ene,"2.30",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1070   else
1071     plotMT(element,ene,"2.30",first,logy,ymin,particle,beam,leg2,dir,dird,markf);
1072 
1073   char fname[60];
1074   if (save != 0) {
1075     std::string tag=".gif";
1076     if (ratio) {
1077       if (save > 0) tag = "R.eps";
1078       else          tag = "R.gif";
1079     } else {
1080       if (save > 0) tag = ".eps";
1081     }
1082     sprintf (fname, "%s%sto%sat%sGeV%s", beam, element, particle, ene, tag.c_str());
1083     myc->SaveAs(fname);
1084   }
1085 }
1086  
1087 void plotMT1(char element[2], char ene[6], char rapid[6], int first=0, 
1088          int logy=0, int save=0, double ymin=-1, char particle[8]="piplus",
1089          char beam[8]="proton", bool ratio=false, int legend=1, 
1090          char dir[20]=".", char dird[40]=".", char markf[4]=" ") {
1091 
1092   setStyle();
1093   TCanvas *myc = new TCanvas("myc","",800,600); myc->SetLeftMargin(0.15);
1094   if (logy != 0) gPad->SetLogy(1);
1095   if (ratio)
1096     plotMTRatio(element,ene,rapid,first,logy,ymin,particle,beam,legend,dir,dird,markf);
1097   else
1098     plotMT(element,ene,rapid,first,logy,ymin,particle,beam,legend,dir,dird,markf);
1099 
1100   char fname[60];
1101   if (save != 0) {
1102     std::string tag=".gif";
1103     if (ratio) {
1104       if (save > 0) tag = "R.eps";
1105       else          tag = "R.gif";
1106     } else {
1107       if (save > 0) tag = ".eps";
1108     }
1109     sprintf (fname, "%s%sto%sat%sGeVY%s%s", beam, element, particle, ene, rapid, tag.c_str());
1110     myc->SaveAs(fname);
1111   }
1112 }
1113 
1114 void plotMT(char element[2], char ene[6], char rapid[6], int first=0, 
1115         int logy=0, double ymin=-1, char particle[8]="piplus", 
1116         char beam[8]="proton", int legend=0, char dir[20]=".",
1117         char dird[40]=".", char markf[4]=" ") {
1118 
1119   char fname[120], list[40], hname[60], titlx[50], sym[6];
1120   TH1F *hi[6];
1121   int i=0, icol=1;
1122   if      (particle=="piminus") sprintf(sym, "#pi^{-}");
1123   else if (particle=="piplus")  sprintf(sym, "#pi^{+}");
1124   else if (particle=="kminus")  sprintf(sym, "K^{-}");
1125   else if (particle=="kplus")   sprintf(sym, "K^{+}");
1126   else                          sprintf(sym, "p");
1127   sprintf (titlx, "Reduced m_{T} (GeV)");
1128   double  ymx0=1, ymi0=100., xlow=0.1, xhigh=1.6;
1129   for (i=0; i<modelsBNL; i++) {
1130     sprintf (list, "%s", ModelFilesB[i].c_str()); 
1131     sprintf (fname, "%s/%s/%s/%s%s%sGeV_1.root", dir, beam, particle, element, list, ene);
1132     sprintf (list, "%s", ModelsB[i].c_str()); icol = colModel[i];
1133     sprintf (hname, "KE0%s%s%sGeVy%s", element, list, ene, rapid);
1134     TFile *file = new TFile(fname);
1135     hi[i] = (TH1F*) file->Get(hname);
1136     std::cout << "Get " << hname << " from " << fname <<" as " << hi[i] <<"\n";
1137     int nx = hi[i]->GetNbinsX();
1138     for (int k=1; k <= nx; k++) {
1139       double xx = hi[i]->GetBinCenter(k);
1140       double yy = hi[i]->GetBinContent(k);
1141       if (xx > xlow && xx < xhigh) {
1142     if (yy > ymx0) ymx0 = yy;
1143     if (yy < ymi0 && yy > 0) ymi0 = yy;
1144       }
1145     }
1146     if (debug) std::cout << "ylimit " << ymi0 << ":" << ymx0 << "\n";
1147     hi[i]->GetXaxis()->SetRangeUser(xlow, xhigh); hi[i]->SetTitle("");
1148     hi[i]->GetXaxis()->SetTitle(titlx);
1149     hi[i]->SetLineStyle(1);  hi[i]->SetLineWidth(2); hi[i]->SetLineColor(icol);
1150     //    file->Close();
1151   }
1152 
1153   sprintf (fname, "%s/bnl802/%s/%s/%s%sGeVRap%s.dat", dird, beam, particle, element, ene, rapid);
1154   std::cout << "Reads data from file " << fname << "\n";
1155   ifstream infile;
1156   infile.open(fname);
1157   int     q1;
1158   float   ym1, ym2, sys, x1[50], y1[50], stater1[50], syser1[50];
1159   infile >> q1 >> ym1 >> ym2 >> sys;
1160   for (i=0; i<q1; i++) {
1161     infile >> x1[i] >> y1[i] >> stater1[i];
1162     syser1[i] = sys*y1[i];
1163     double err = sqrt(syser1[i]*syser1[i]+stater1[i]*stater1[i]);
1164     stater1[i] = err;
1165     if (y1[i]+stater1[i] > ymx0) ymx0 = y1[i]+stater1[i];    
1166     if (y1[i]-stater1[i] < ymi0 && y1[i]-stater1[i] > 0) ymi0=y1[i]-stater1[i];
1167     if (debug) std::cout << i << " " << x1[i] << " " << y1[i] << " " << stater1[i] << "\n";
1168   }
1169   TGraph*  gr1 = new TGraphErrors(q1,x1,y1,0,stater1);
1170   gr1->SetMarkerColor(4);  gr1->SetMarkerStyle(22);
1171   gr1->SetMarkerSize(1.6);
1172 
1173   if (logy == 0) {ymx0 *= 1.5; ymi0 *= 0.8;}
1174   else           {ymx0 *=10.0; ymi0 *= 0.2; }
1175   if (ymin > 0) ymi0 = ymin;
1176   for (i = 0; i<modelsBNL; i++) {
1177     if (debug) std::cout << "Model " << i << " " << hi[i] << " " << ymi0 << " " << ymx0 << "\n";
1178     hi[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
1179   }
1180 
1181   hi[first]->GetYaxis()->SetTitleOffset(1.1);
1182   hi[first]->Draw();
1183   for (i=0; i<modelsBNL; i++) {
1184     if (i != first) hi[i]->Draw("same");
1185   }
1186   gr1->Draw("p");
1187 
1188   TLegend *leg1;
1189   if (legend < 0) {
1190     leg1 = new TLegend(0.50,0.55,0.90,0.90);
1191   } else {
1192     if (markf == " " ) leg1 = new TLegend(0.42,0.70,0.90,0.90);
1193     else               leg1 = new TLegend(0.38,0.70,0.90,0.90);
1194   }
1195   for (i=0; i<modelsBNL; i++) {
1196     sprintf (list, "%s", ModelNamesB[i].c_str()); 
1197     leg1->AddEntry(hi[i],list,"F");
1198   }
1199   char header[120], beamx[8], partx[2];
1200   if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
1201   else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
1202   else                        sprintf (beamx, "p");
1203   if (legend < 0) {
1204     sprintf (header,"%s+%s #rightarrow %s+X at %s GeV", beamx, element, sym, ene);
1205   } else {
1206     if (markf == " ")
1207       sprintf (header,"%s+%s #rightarrow %s+X at %s GeV (y = %s)", beamx, element, sym, ene, rapid);
1208     else
1209       sprintf (header,"%s %s+%s #rightarrow %s+X at %s GeV (y = %s)", markf, beamx, element, sym, ene, rapid);
1210   }
1211   leg1->SetHeader(header); leg1->SetFillColor(0);
1212   leg1->SetTextSize(0.04);
1213   if (legend != 0) leg1->Draw("same");
1214 
1215   if (debug) std::cout << "End\n";
1216 }
1217 
1218 void plotMTRatio(char element[2], char ene[6], char rapid[6], int first=0, 
1219          int logy=0, double ymin=-1, char particle[8]="piplus", 
1220          char beam[8]="proton", int legend=0, char dir[20]=".",
1221          char dird[40]=".", char markf[4]=" ") {
1222 
1223   char titlx[50], sym[6];
1224   int  i=0, icol=1, ityp=20;
1225   if      (particle=="piminus") sprintf(sym, "#pi^{-}");
1226   else if (particle=="piplus")  sprintf(sym, "#pi^{+}");
1227   else if (particle=="kminus")  sprintf(sym, "K^{-}");
1228   else if (particle=="kplus")   sprintf(sym, "K^{+}");
1229   else                          sprintf(sym, "p");
1230   sprintf (titlx, "Reduced m_{T} (GeV)");
1231 
1232   //Read in the data files
1233   char fname[120];
1234   sprintf (fname, "%s/bnl802/%s/%s/%s%sGeVRap%s.dat", dird, beam, particle, element, ene, rapid);
1235   if (debug) std::cout << "Reads data from file " << fname << "\n";
1236   ifstream infile;
1237   infile.open(fname);
1238   int     q1;
1239   float   ym1, ym2, sys, x1[50], y1[50], er1[50], staterr, syserr;
1240   infile >> q1 >> ym1 >> ym2 >> sys;
1241   for (i=0; i<q1; i++) {
1242     infile >> x1[i] >> y1[i] >> staterr;
1243     syserr = sys*y1[i];
1244     er1[i] = sqrt(syserr*syserr+staterr*staterr);
1245     if (debug) std::cout << i << " " << x1[i] << " " << y1[i] << " " << er1[i] << "\n";
1246   }
1247 
1248   char          list[40], hname[60];
1249   TGraphErrors *gr[6];
1250   double        ymx0=1, ymi0=100., xlow=0.1, xhigh=1.6;
1251   for (i=0; i<modelsBNL; i++) {
1252     icol = colModel[i]; ityp = symbModel[i];
1253     sprintf (list, "%s", ModelFilesB[i].c_str()); 
1254     sprintf (fname, "%s/%s/%s/%s%s%sGeV_1.root", dir, beam, particle, element, list, ene);
1255     sprintf (list, "%s", ModelsB[i].c_str());
1256     sprintf (hname, "KE0%s%s%sGeVy%s", element, list, ene, rapid);
1257 
1258     TFile *file = new TFile(fname);
1259     TH1F  *hi   = (TH1F*) file->Get(hname);
1260     if (debug) std::cout << "Get " << hname << " from " << fname <<" as " << hi <<"\n";
1261 
1262     if (hi != 0 && q1 > 0) {
1263       float xx[50], dx[50], rat[50], drt[50];
1264       int   nx = hi->GetNbinsX();
1265       int   np = 0;
1266       if (debug) std::cout << "Start with " << nx << " bins\n";
1267       for (int k=1; k <= nx; k++) {
1268     double xx1 = hi->GetBinLowEdge(k);
1269     double xx2 = hi->GetBinWidth(k);
1270     for (int j=0; j<q1; j++) {
1271       if (xx1 < x1[j] && xx1+xx2 > x1[j]) {
1272         double yy = hi->GetBinContent(k);
1273         xx[np]    = x1[j];
1274         dx[np]    = 0;
1275         rat[np]   = yy/y1[j];
1276         drt[np]   = er1[j]*rat[j]/y1[j];
1277         if (xx[np] > xlow && xx[np] < xhigh) {
1278           if (rat[np]+drt[np] > ymx0) ymx0 = rat[np]+drt[np];
1279           if (rat[np]-drt[np] < ymi0) ymi0 = rat[np]-drt[np];
1280         }
1281         if (debug) std::cout << np << "/" << j << "/" << k << " x " << xx[np] << " (" << xx1 << ":" << xx1+xx2 << ")" << " y " << yy << "/" << y1[j] << " = " << rat[np] << " +- " << drt[np] << "\n";
1282         np++;
1283         break;
1284       }
1285     }
1286       }
1287       gr[i] = new TGraphErrors(np, xx, rat, dx, drt);
1288       gr[i]->GetXaxis()->SetRangeUser(xlow, xhigh); gr[i]->SetTitle("");
1289       gr[i]->GetXaxis()->SetTitle(titlx);
1290       gr[i]->GetYaxis()->SetTitle("MC/Data");
1291       gr[i]->SetLineStyle(stylModel[i]); gr[i]->SetLineWidth(2); 
1292       gr[i]->SetLineColor(icol);         gr[i]->SetMarkerColor(icol); 
1293       gr[i]->SetMarkerStyle(ityp);       gr[i]->SetMarkerSize(1.0); 
1294     } else {
1295       gr[i] = 0;
1296     }
1297     file->Close();
1298   }
1299 
1300   if (logy == 0) {ymx0 *= 1.5; ymi0 *= 0.8;}
1301   else           {ymx0 *=10.0; ymi0 *= 0.2; }
1302   if (ymin > 0) ymi0 = ymin;
1303   for (i = 0; i<modelsBNL; i++) {
1304     if (gr[i] != 0) {
1305       if (debug) std::cout << "Model " << i << " " << gr[i] << " " << ymi0 << " " << ymx0 << "\n";
1306       gr[i]->GetYaxis()->SetRangeUser(ymi0,ymx0);
1307     }
1308   }
1309 
1310   if (gr[first] > 0) {
1311     gr[first]->GetYaxis()->SetTitleOffset(1.1);
1312     gr[first]->Draw("APl");
1313     for (i=0; i<modelsBNL; i++) {
1314       if (i != first && gr[i] != 0) gr[i]->Draw("Pl");
1315     }
1316 
1317     TLegend *leg1;
1318     if (legend < 0) {
1319       leg1 = new TLegend(0.50,0.55,0.90,0.90);
1320     } else {
1321       if (markf == " " ) leg1 = new TLegend(0.42,0.70,0.90,0.90);
1322       else               leg1 = new TLegend(0.38,0.70,0.90,0.90);
1323     }
1324     for (i=0; i<modelsBNL; i++) {
1325       if (gr[i] != 0) {
1326     sprintf (list, "%s", ModelNamesB[i].c_str()); 
1327     leg1->AddEntry(gr[i],list,"lP");
1328       }
1329     }
1330     char header[120], beamx[8], partx[2];
1331     if      (beam == "piplus")  sprintf (beamx, "#pi^{+}");
1332     else if (beam == "piminus") sprintf (beamx, "#pi^{-}");
1333     else                        sprintf (beamx, "p");
1334     if (legend < 0) {
1335       sprintf (header,"%s+%s #rightarrow %s+X at %s GeV", beamx, element, sym, ene);
1336     } else {
1337       if (markf == " ")
1338     sprintf (header,"%s+%s #rightarrow %s+X at %s GeV (y = %s)", beamx, element, sym, ene, rapid);
1339       else
1340     sprintf (header,"%s %s+%s #rightarrow %s+X at %s GeV (y = %s)", markf, beamx, element, sym, ene, rapid);
1341     }
1342     leg1->SetHeader(header); leg1->SetFillColor(0);
1343     leg1->SetTextSize(0.04);
1344     if (legend != 0) leg1->Draw("same");
1345 
1346     xx[0]=xlow; xx[1]=xhigh; rat[0]=rat[1]=1.0;
1347     TGraph *gr0 = new TGraph(2, xx, rat);
1348     gr0->GetXaxis()->SetRangeUser(xlow, xhigh); gr0->SetTitle("");
1349     gr0->SetLineStyle(1);   gr0->SetLineWidth(1.4); 
1350     gr0->SetLineColor(1);   gr0->SetMarkerColor(1); 
1351     gr0->SetMarkerStyle(20);gr0->SetMarkerSize(1.6);
1352     gr0->Draw("l");
1353   }
1354 }
1355 
1356 void setStyle() {
1357 
1358   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1359   gStyle->SetPadColor(kWhite);    gStyle->SetFrameBorderMode(0);
1360   gStyle->SetFrameBorderSize(1);  gStyle->SetFrameFillColor(0);
1361   gStyle->SetFrameFillStyle(0);   gStyle->SetFrameLineColor(1);
1362   gStyle->SetFrameLineStyle(1);   gStyle->SetFrameLineWidth(1);
1363   gStyle->SetTitleOffset(1.6,"Y");  gStyle->SetOptStat(0);
1364   gStyle->SetLegendBorderSize(1);
1365 
1366 }
1367 
1368 std::vector<double> angleScan(int scan) {
1369 
1370   std::vector<double> tmp;
1371   if (scan <= 1) {
1372     tmp.push_back(59.1);
1373     tmp.push_back(89.0);
1374     tmp.push_back(119.0);
1375     tmp.push_back(159.6);
1376   } else {
1377     tmp.push_back(10.1);
1378     tmp.push_back(15.0);
1379     tmp.push_back(19.8);
1380     tmp.push_back(24.8);
1381     tmp.push_back(29.5);
1382     tmp.push_back(34.6);
1383     tmp.push_back(39.6);
1384     tmp.push_back(44.3);
1385     tmp.push_back(49.3);
1386     tmp.push_back(54.2);
1387     tmp.push_back(59.1);
1388     tmp.push_back(64.1);
1389     tmp.push_back(69.1);
1390     tmp.push_back(74.1);
1391     tmp.push_back(79.1);
1392     tmp.push_back(84.1);
1393     tmp.push_back(89.0);
1394     tmp.push_back(98.9);
1395     tmp.push_back(108.9);
1396     tmp.push_back(119.0);
1397     tmp.push_back(129.1);
1398     tmp.push_back(139.1);
1399     tmp.push_back(149.3);
1400     tmp.push_back(159.6);
1401     tmp.push_back(161.4);
1402     tmp.push_back(165.5);
1403     tmp.push_back(169.5);
1404     tmp.push_back(173.5);
1405     tmp.push_back(177.0);
1406   }
1407   if (debug) {
1408     std::cout << "Scan " << tmp.size() << " angular regions:\n";
1409     for (unsigned int i=0; i<tmp.size(); i++) {
1410       std::cout << tmp[i];
1411       if (i == tmp.size()-1) std::cout << " degrees\n";
1412       else                   std::cout << ", ";
1413     }
1414   }
1415   return tmp;
1416 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!