void Main_layers(void)
{
// Main_test inputing arguments to a script
//

#include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kTRUE);
gStyle->SetOptStat(kFALSE);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetOptStat(111111);
// gStyle->SetOptFit(1111);
//

   Int_t j,jj;
   char string[256];

   Double_t rmin=10;
   Int_t krows=20;
   Int_t klayers=6;
   Double_t sigdy = 0.02;
   Double_t rstraw=0.8;
   // Double_t rstraw=2;
   #define npts 10000;
   // #define npts 1000;
   Double_t phimin=0;
   Double_t phimax=360;
   Double_t prob_cut=0.01;
   // Reff: Fraction of radius for efficient operation. 0.97-> Thresh=1KeV/4.16 KeV for 1.6 cm
   Double_t Reff=0.97; 

   // use kgeom to tag geometry
   // kgeom=0 is original, unshifted axial geometry
   // kgeom=1 is axial geometry with shifted layers

   Int_t kgeom=2;


   //
   // define canvas 

   TCanvas *c0 = new TCanvas("c0","Histogram of Chi2 of fits",200,10,700,700);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);
   // c0->SetGridx();
   // c0->SetGridy();
   c0->Divide(2,2);

   TCanvas *c2 = new TCanvas("c2","Good Solutions vs angle",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   // c2->SetGridx();
   // c2->SetGridy();
   c2->Divide(1,2);

   TCanvas *c3 = new TCanvas("c3","Extended range Chi2 plots",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   // c3->SetGridx();
   // c3->SetGridy();
   c3->Divide(2,2);

   TCanvas *c4 = new TCanvas("c4","Chi2 of parameters",200,10,700,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   // c4->SetGridx();
   // c4->SetGridy();
   c4->Divide(1,2);

   TCanvas *c5 = new TCanvas("c5","Chi2 of parameters",200,10,700,700);
   c5->SetBorderMode(0);
   c5->SetFillColor(0);
   // c5->SetGridx();
   // c5->SetGridy();
   c5->Divide(2,2);

   // Define Chi2 histogram 

   TH2F *sol_true = new TH2F("sol_true","Good solutions with true hits vs angle",360,0,360,5,0,5);
   TH2F *sol_comb = new TH2F("sol_comb","Good solutions with all combinations vs angle",360,0,360,20,0,20);

   Int_t nbins=100;
   TH1F *chisq1 = new TH1F("chisq1","Chi2 of fit to true hits",nbins,0,5);
   TH1F *chisqN = new TH1F("chisqN","Chi2 of fit to all combinations",nbins,0,5);
   TH1F *Schisq1 = new TH1F("Schisq1","SChi2 of fit to true hits",10*nbins,0,100);
   TH1F *SchisqN = new TH1F("SchisqN","SChi2 of fit to all combinations",10*nbins,0,100);
   TH1F *intchi1 = new TH1F("intchi1","chisq1 bins integral",10*nbins,0,100);
   TH1F *intchiN = new TH1F("intchiN","chisqN bins integral",10*nbins,0,100);
   TH1F *prob1 = new TH1F("prob1","Prob of fit to true hits",nbins,0,1);
   TH1F *probN = new TH1F("probN","Prob of fit to all combinations",nbins,0,1);
   TH1F *solN = new TH1F("solN","Good solutions (p>cut)",20,0,20);
   TH1F *goodN = new TH1F("goodN","Good solutions (p>cut && chi2par>5)",20,0,20);

   TH2F *parchi2 = new TH2F("parchi2","Parameter chi2 vs angle (p>cut)",360,0,360,20,0,20);
   TH2F *trkchi2 = new TH2F("trkchi2","Track chi2 vs angle",360,0,360,nbins,0,5);
   TH2F *pchi2chi2 = new TH2F("pchi2chi2","Parameter chi2 vs track chi2",20,0,20,20,0,20);
   TH1F *pchi2 = new TH1F("pchi2","Parameter chi2 fit (p>cut)",20,0,20);

   sprintf(string,"Starting Test Main_layers_l%d_g%d.root",klayers,kgeom);
   printf ("filename=%s\n",string);

   // Loop over phi.

   for (j=0;j<npts;j++) {
      Double_t phi = phimin + (phimax-phimin)*j/npts;
          if (! plot_layers(kgeom,rmin,rstraw,krows,klayers,phi,sigdy,Reff)) {
              if ((j/100)*100 == j) printf ("j=%d, phi=%f\n",j,phi);
          }
          else {
           printf ("*** Main_layers: plot_layers failed!\n" ); 
          }
   }

   // get integral of chi2 distributions.

   Double_t sum = 0;
   for (j=1;j<=10*nbins;j++) {
       sum += Schisq1->GetBinContent(j);
       intchi1->SetBinContent(j,sum);;
       // printf("j=%d, sum=%f\n",j,sum);
       }
   Double_t sum = 0;
   for (j=1;j<=10*nbins;j++) {
       sum += SchisqN->GetBinContent(j);
       intchiN->SetBinContent(j,sum);
       // printf("j=%d, sum=%f\n",j,sum);
       }

   Int_t jcut=prob_cut*nbins;
   Double_t sum = 0;
   Double_t sumcut = 0;
   Double_t sumtrue = 0;
   Double_t sumcut_true = 0;
   for (j=1;j<=nbins;j++) {
       sum += probN->GetBinContent(j);
       sumtrue += prob1->GetBinContent(j);
       if (j > jcut) {
           sumcut+ = probN->GetBinContent(j);
           sumcut_true+ = prob1->GetBinContent(j);
       }
   }
   printf("sum=%f, sumcut=%f\n",sum,sumcut);

   // plot histogram

   c0->cd(1);
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);
   Double_t ymax = chisq1->GetMaximum();
   chisq1->GetYaxis()->SetRangeUser(0,1.2*ymax);
   chisq1->Draw();
   sprintf (string,"Rmin=%.1f cm \n",rmin);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"klayers=%d\n",klayers);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sigdy=%.3f cm\n",sigdy);
   t1 = new TLatex(0.5,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm \n",rstraw);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"Reff=%.2f\n",Reff);
   t1 = new TLatex(0.5,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c0->cd(2);
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);
   Double_t ymax = chisqN->GetMaximum();
   chisqN->GetYaxis()->SetRangeUser(0,1.2*ymax);
   chisqN->Draw();

   c0->cd(3);
   c0_3->SetBorderMode(0);
   c0_3->SetFillColor(0);
   Double_t ymax = prob1->GetMaximum();
   prob1->GetYaxis()->SetRangeUser(0,1.2*ymax);
   prob1->Draw();
   sprintf (string,"sum=%.1f trials \n",sumtrue);
   t1 = new TLatex(0.4,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sum=%.1f, p>%.3f \n",sumcut_true,prob_cut);
   t1 = new TLatex(0.4,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c0->cd(4);
   c0_4->SetLogy();
   c0_4->SetBorderMode(0);
   c0_4->SetFillColor(0);
   probN->Draw();
   sprintf (string,"sum=%.1f trials \n",sum);
   t1 = new TLatex(0.4,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sum=%.1f, p>%.3f \n",sumcut,prob_cut);
   t1 = new TLatex(0.4,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c2->cd(1);
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   sol_true->Draw("colz");


   sprintf (string,"Rmin=%.1f cm \n",rmin);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"klayers=%d\n",klayers);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sigdy=%.3f cm\n",sigdy);
   t1 = new TLatex(0.5,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm \n",rstraw);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"Reff=%.2f\n",Reff);
   t1 = new TLatex(0.5,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   sol_comb->Draw("colz");

   c3->cd(1);
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);
   // Schisq1->GetXaxis()->SetRangeUser(0,10);
   Double_t ymax = Schisq1->GetMaximum();
   Schisq1->GetYaxis()->SetRangeUser(0,1.2*ymax);
   Schisq1->Draw();

   sprintf (string,"Rmin=%.1f cm \n",rmin);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"klayers=%d\n",klayers);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sigdy=%.3f cm\n",sigdy);
   t1 = new TLatex(0.5,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm \n",rstraw);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"Reff=%.2f\n",Reff);
   t1 = new TLatex(0.5,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c3->cd(2);
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   // SchisqN->GetXaxis()->SetRangeUser(0,10);
   Double_t ymax = SchisqN->GetMaximum();
   SchisqN->GetYaxis()->SetRangeUser(0,1.2*ymax);
   SchisqN->Draw();

   c3->cd(3);
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);
   Double_t ymax = intchi1->GetMaximum();
   // intchi1->GetXaxis()->SetRangeUser(1,10*npts);
   intchi1->GetYaxis()->SetRangeUser(0,1.2*ymax);
   c3_3->SetLogx();
   sum = Schisq1->Integral();
   intchi1->Scale(1/sum);
   intchi1->Draw();
   sprintf (string,"sum=%.1f trials\n",sum);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c3->cd(4);
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   c3_4->SetLogx();
   Double_t ymax = intchiN->GetMaximum();
   // intchiN->GetXaxis()->SetRangeUser(1,10*npts);
   intchiN->GetYaxis()->SetRangeUser(0,1.2*ymax);
   sum = SchisqN->Integral();
   intchiN->Scale(1/sum);
   intchiN->Draw();
   sprintf (string,"sum=%.1f trials \n",sum);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c4->cd(1);
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);
   // parchi2->GetXaxis()->SetRangeUser(0,10);
   Double_t ymax = Schisq1->GetMaximum();
   parchi2->GetYaxis()->SetRangeUser(0,1.2*ymax);
   parchi2->Draw("colz");

   sprintf (string,"Rmin=%.1f cm \n",rmin);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"klayers=%d\n",klayers);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sigdy=%.3f cm\n",sigdy);
   t1 = new TLatex(0.5,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm \n",rstraw);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"Reff=%.2f\n",Reff);
   t1 = new TLatex(0.5,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c4->cd(2);
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   // sol_comb->GetXaxis()->SetRangeUser(0,10);
   c4_2->SetLogz();
   Double_t ymax = chisqN->GetMaximum();
   trkchi2->GetYaxis()->SetRangeUser(0,1.2*ymax);
   trkchi2->Draw("colz");


   c5->cd(1);
   c5_1->SetBorderMode(0);
   c5_1->SetFillColor(0);
   // pchi2->GetXaxis()->SetRangeUser(0,10);
   Double_t ymax = pchi2->GetMaximum();
   c5_1->SetLogy();
   // pchi2->GetYaxis()->SetRangeUser(0,1.2*ymax);
   pchi2->Draw();

   sprintf (string,"Rmin=%.1f cm \n",rmin);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"klayers=%d\n",klayers);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"sigdy=%.3f cm\n",sigdy);
   t1 = new TLatex(0.5,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm \n",rstraw);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf (string,"Reff=%.2f\n",Reff);
   t1 = new TLatex(0.5,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c5->cd(2);
   c5_2->SetBorderMode(0);
   c5_2->SetFillColor(0);
   // pchi2chi2->GetXaxis()->SetRangeUser(0,10);
   Double_t ymax = pchi2chi2->GetMaximum();
   c5_2->SetLogz();
   // pchi2chi2->GetYaxis()->SetRangeUser(0,1.2*ymax);
   pchi2chi2->Draw("colz");
   c5->cd(3);
   c5_3->SetLogy();
   c5_3->SetBorderMode(0);
   c5_3->SetFillColor(0);
   Double_t sol_mean = solN->GetMean();
   Double_t sol_meanerr = solN->GetMeanError();
   solN->Draw();
   sprintf (string,"mean=%.3f #pm %.3f\n",sol_mean,sol_meanerr);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   Double_t sol_sum = solN->GetSum();
   sprintf (string,"sum=%.2f\n",sol_sum);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c5->cd(4);
   c5_4->SetLogy();
   c5_4->SetBorderMode(0);
   c5_4->SetFillColor(0);
   Double_t sol_mean = goodN->GetMean();
   Double_t sol_meanerr = goodN->GetMeanError();
   goodN->Draw();
   sprintf (string,"mean=%.3f #pm %.3f\n",sol_mean,sol_meanerr);
   t1 = new TLatex(0.5,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   Double_t sol_sum = goodN->GetSum();
   sprintf (string,"sum=%.2f\n",sol_sum);
   t1 = new TLatex(0.5,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   

   // save histogram to file

   sprintf(string,"Main_layers_l%d_g%d.root",klayers,kgeom);
   printf ("filename=%s\n",string);
   TFile *out = new TFile(string,"recreate");
   chisq1->Write();
   chisqN->Write();
   intchi1->Write();
   intchiN->Write();
   Schisq1->Write();
   SchisqN->Write();
   prob1->Write();
   probN->Write();
   sol_true->Write();
   sol_comb->Write();
   solN->Write();
   goodN->Write();
   parchi2->Write();
   trkchi2->Write();
   pchi2->Write();
   out->Close();

   /*sprintf(string,"Main_layers_c0_l%d_g%d.eps",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_c0_l%d_g%d.gif",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_c2_l%d_g%d.eps",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_c2_l%d_g%d.gif",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_c3_l%d_g%d.eps",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_c3_l%d_g%d.gif",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_c4_l%d_g%d.eps",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_c4_l%d_g%d.gif",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_c5_l%d_g%d.eps",klayers,kgeom);
   c5->SaveAs(string);
   sprintf(string,"Main_layers_c5_l%d_g%d.gif",klayers,kgeom);
   c5->SaveAs(string);*/

   /*sprintf(string,"Main_layers_eff2_c0_l%d_g%d.eps",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c0_l%d_g%d.gif",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c2_l%d_g%d.eps",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c2_l%d_g%d.gif",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c3_l%d_g%d.eps",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c3_l%d_g%d.gif",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c4_l%d_g%d.eps",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c4_l%d_g%d.gif",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c5_l%d_g%d.eps",klayers,kgeom);
   c5->SaveAs(string);
   sprintf(string,"Main_layers_eff2_c5_l%d_g%d.gif",klayers,kgeom);
   c5->SaveAs(string);*/

   sprintf(string,"Main_layers_eff3_c0_l%d_g%d.eps",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c0_l%d_g%d.gif",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c2_l%d_g%d.eps",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c2_l%d_g%d.gif",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c3_l%d_g%d.eps",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c3_l%d_g%d.gif",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c4_l%d_g%d.eps",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c4_l%d_g%d.gif",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c5_l%d_g%d.eps",klayers,kgeom);
   c5->SaveAs(string);
   sprintf(string,"Main_layers_eff3_c5_l%d_g%d.gif",klayers,kgeom);
   c5->SaveAs(string);

   /*sprintf(string,"Main_layers_test_c0_l%d_g%d.eps",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_test_c0_l%d_g%d.gif",klayers,kgeom);
   c0->SaveAs(string);
   sprintf(string,"Main_layers_test_c2_l%d_g%d.eps",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_test_c2_l%d_g%d.gif",klayers,kgeom);
   c2->SaveAs(string);
   sprintf(string,"Main_layers_test_c3_l%d_g%d.eps",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_test_c3_l%d_g%d.gif",klayers,kgeom);
   c3->SaveAs(string);
   sprintf(string,"Main_layers_test_c4_l%d_g%d.eps",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_test_c4_l%d_g%d.gif",klayers,kgeom);
   c4->SaveAs(string);
   sprintf(string,"Main_layers_test_c5_l%d_g%d.eps",klayers,kgeom);
   c5->SaveAs(string);
   sprintf(string,"Main_layers_test_c5_l%d_g%d.gif",klayers,kgeom);
   c5->SaveAs(string);*/



}
// Int_t plot_layers(Int_t kgeom, Double_t rmin=10, Double_t rstraw=2, Int_t krows=20, Int_t klayers=4, Double_t angle=39, Double_t sigdy=0.02, Double_t Reff=0.97)
Int_t plot_layers(Int_t kgeom, Double_t rmin, Double_t rstraw, Int_t krows, Int_t klayers, Double_t angle, Double_t sigdy, Double_t Reff)
{
// Plot hits in four consecutive rows of CDC layers with fitted tracks.
// This version has no graphics (commented out) to avoid memory issues.

#include <TMath.h>
#include <TRandom.h>

// gROOT->Reset();
TH2F *sol_true = (TH2F *) gROOT->FindObject("sol_true");
TH2F *sol_comb = (TH2F *) gROOT->FindObject("sol_comb");
TH1F *chisq1 = (TH1F *) gROOT->FindObject("chisq1");
TH1F *chisqN = (TH1F *) gROOT->FindObject("chisqN");
TH1F *Schisq1 = (TH1F *) gROOT->FindObject("Schisq1");
TH1F *SchisqN = (TH1F *) gROOT->FindObject("SchisqN");
TH1F *prob1 = (TH1F *) gROOT->FindObject("prob1");
TH1F *probN = (TH1F *) gROOT->FindObject("probN");
TH1F *solN = (TH1F *) gROOT->FindObject("solN");
TH1F *goodN = (TH1F *) gROOT->FindObject("goodN");
TH2F *parchi2 = (TH2F *) gROOT->FindObject("parchi2");
TH2F *trkchi2 = (TH2F *) gROOT->FindObject("trkchi2");
TH2F *pchi2chi2 = (TH2F *) gROOT->FindObject("pchi2chi2");
TH1F *pchi2 = (TH1F *) gROOT->FindObject("pchi2");
// printf(" plot_layers- All histograms found in memory\n");
// gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kFALSE);
// gStyle->SetPadRightMargin(0.15);
// gStyle->SetPadLeftMargin(0.2);
// gStyle->SetPadBottomMargin(0.15);
// gStyle->SetOptFit(1111);
//

   char string[256];
   Int_t j,jj,k,kk,kkk,kkkk;
   Double_t prob_cut=0.01;

   Double_t *wire_hitsx =new Double_t [klayers] ;
   Double_t *wire_hitsy =new Double_t [klayers];
   Double_t *true_hitsx =new Double_t [klayers];
   Double_t *true_hitsy =new Double_t [klayers];
   Double_t *amb_hitsx =new Double_t [klayers];
   Double_t *amb_hitsy =new Double_t [klayers];
   Double_t *select_hitsx =new Double_t [klayers];
   Double_t *select_hitsy =new Double_t [klayers];
   Double_t *errorsx =new Double_t [klayers];
   Double_t *errorsy =new Double_t [klayers];
   
   // #define klayersmax 7;
   // Double_t shifts[klayersmax]={0.018265,-0.011701,0.006756,-0.004651,0.002856,-0.002335,0.001727}

   // fill error histograms

   for (j=0;j<klayers;j++) {
      errorsx[j] = 0.0001;
      errorsy[j] = sigdy;
   }

   //
   // define canvas 

   /*TCanvas *c1 = new TCanvas("c1","Plot hits in four CDC layers",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetGridx();
   c1->SetGridy();
   // gPad->SetFixedAspectRatio(kTRUE);*/

   // Double_t rstraw=0.8;
   // Double_t rmin=10;
   Double_t rmax=55;
   Double_t kpi=3.14159;
   // Int_t krows=20;
   // Int_t klayers=4;
   Double_t phi=angle*kpi/180;
   // Double_t sigdy = 0.02;
   Double_t xoffset = 1;
   Double_t scale=2;

   Double_t ymin=-4*rstraw;

   // Convert to actual cylindrical geometry

   Double_t Circ=2*kpi*rmin;
   Int_t ntubes = ceil(Circ/(2*rstraw));
   // printf ("rmin=%f, ntubes=%d\n",rmin,ntubes);
   Double_t r1=2*rstraw*ntubes/(2*kpi);
   Double_t Dphi = 2*kpi/ntubes;
   Double_t ymin = 2*r1*sin(Dphi);
   Double_t xmin = r1*cos(Dphi);

   /* gPad->Range(0-xoffset,-rmin*scale,2*rmin*scale-xoffset,rmin*scale);
   TEllipse *circle = new TEllipse(xmin,ymin,rstraw);
   circle->Draw();
   TMarker *origin = new TMarker(0,0,2);
   origin->Draw();
   TLine *horizon = new TLine(0-xoffset,0,2*rmin*scale-xoffset,0);
   horizon->SetLineStyle(2);
   horizon->Draw();*/

   TRandom1 *rand = new TRandom1();
   // rand->SetSeed();


   // draw row

   Double_t shift=0;
   Int_t jfirst = krows/2 -1;
   Double_t phidif0=0;
   Double_t phidif=0;
   for (j=0;j<klayers;j++) {
      if (kgeom < 2) {

         // normal configuration

         ntubes = ntubes + 7;
         Double_t r1=2*rstraw*ntubes/(2*kpi);
         Double_t Dphi = 2*kpi/ntubes;

         // find location of second tube for reference  
         Double_t arc=phi*r1;
         Double_t arc1 = 2*rstraw*floor(arc/(2*rstraw));
         Double_t phidif0 = (arc-arc1)/r1;

      }
      elseif (kgeom == 2) {

         // close packed configuration

         if (j == 0) {
            ntubes = ntubes + 7;
            Double_t r1=2*rstraw*ntubes/(2*kpi);
            Double_t Dphi = 2*kpi/ntubes;
            Double_t arc=phi*r1;
            Double_t arc1 = 2*rstraw*floor(arc/(2*rstraw));
            Double_t phidif0 = (arc-arc1)/r1;
            Double_t shift = 0.;
         }
         elseif (j == 1 || j ==3 || j==5) {
            r1 = r1*cos(Dphi/2) + sqrt(3)*rstraw;
            // Double_t shift = Dphi/2;
            Double_t shift = shift + Dphi/2;
         }
         elseif (j==2) {
            Double_t Circ=2*kpi*r1;
            Int_t ntubes = ceil(Circ/(2*rstraw)) + 7;
            Double_t r1=2*rstraw*ntubes/(2*kpi);
            Double_t Dphi = 2*kpi/ntubes;
            shift = 0;
            Double_t arc=phi*r1;
            Double_t arc1 = 2*rstraw*floor(arc/(2*rstraw));
            Double_t phidif0 = (arc-arc1)/r1;

         // calculate additional the shift for "shifted geometry"

         switch(j%8){
           case 0: shift =       0.0;    break;
           case 1: shift = +Dphi/2.0;    break;
           case 2: shift = -Dphi/3.0;    break;
           case 3: shift = +Dphi/5.0;    break;
           case 4: shift = -Dphi/7.0;    break;
           case 5: shift = +Dphi/11.0;    break;
           case 6: shift = -Dphi/13.0;    break;
           case 7: shift = +Dphi/17.0;    break;
         }
         }
         // printf ("\n j=%d, ntubes=%d, r1=%f, Dphi=%f, phidif0=%f, shift=%f\n",j,ntubes,r1,Dphi,phidif0,shift);
      }   

      // printf ("layer=%d, ntubes=%d, radius=%f cm, Dphi=%f rad\n",j,ntubes,r1,Dphi);


      if (kgeom == 1) {

         // calculate the shift for "shifted geometry"

         switch(j%8){
           case 0: shift =       0.0;    break;
           case 1: shift = +Dphi/2.0;    break;
           case 2: shift = -Dphi/3.0;    break;
           case 3: shift = +Dphi/5.0;    break;
           case 4: shift = -Dphi/7.0;    break;
           case 5: shift = +Dphi/11.0;    break;
           case 6: shift = -Dphi/13.0;    break;
           case 7: shift = +Dphi/17.0;    break;
         }
      }

      // add shift to layer

      if (kgeom==1 || kgeom==2) {
         phidif = phidif0 + shift;
      }
      else {
         phidif = phidif0;
      }

      // place hit on track
      Double_t dy = rand->Gaus(0,sigdy);
      // TMarker *hit = new TMarker(r1,dy,5);
      true_hitsx[j] = r1;
      true_hitsy[j] = dy;
      // hit->SetMarkerColor(4);
      // hit->Draw();
      // printf ("r1=%f, dy=%f\n",r1,dy);

      for (jj=0;jj<krows;jj++) {
         Double_t x = r1*cos((jj-jfirst)*Dphi-phidif);
         Double_t y = r1*sin((jj-jfirst)*Dphi-phidif);
         Double_t test = sqrt(x*x+y*y);
         // printf ("rmin=%f, ntubes=%d, r1=%f, Dphi=%f, x=%f, y=%f, test=%f\n",rmin,ntubes,r1,Dphi,x,y,test);
         // circle->DrawEllipse(x,y,rstraw,0,0,360,0);
         // TMarker *dot = new TMarker(x,y,7);
         // dot->Draw();
      }
      x = r1*cos(Dphi-phidif);
      // compute dt 
      y = r1*sin(Dphi-phidif);
      if (y > rstraw+dy) {
         x = r1*cos(-phidif);
         y = r1*sin(-2*phidif) - dy;
         wire_hitsy[j] = r1*sin(-1*phidif);
      }
      else {
         y = r1*sin(2*(Dphi-phidif))-dy;
         // y = r1*sin(-Dphi-phidif);
         wire_hitsy[j] = r1*sin(1*(Dphi-phidif));
      }
      // printf ("j=%d, dy=%f, y=%f\n",j,dy,y);
      // TMarker *hit2 = new TMarker(x,y,5);
      amb_hitsx[j] = x;
      amb_hitsy[j] = y;

      // this line allows for inefficiencies at the edges of the straws.

      if (fabs(y-wire_hitsy[j]) > rstraw*Reff) errorsy[j] = 1.5*rstraw;
      wire_hitsx[j] = x;
      // hit2->SetMarkerColor(2);
      // hit2->Draw();
      // TMarker *hit3 = new TMarker(wire_hitsx[j],wire_hitsy[j],4);
      // hit3->SetMarkerColor(1);
      // hit3->Draw();
   }

   // fit the true hits

   Int_t comb=0;
   Double_t nsolutions = 0;

   TGraphErrors *combination = new TGraphErrors(klayers,true_hitsx,true_hitsy,errorsx,errorsy);

   // remove cells that miss hits

   for (Int_t m=klayers-1;m>-1;m--) {
       if (errorsy[m] > 1.1*sigdy) {
           combination->RemovePoint(m);
       }
   }
   TF1 *track = new TF1("track","pol1",rmin-2*rstraw,rmin+klayers*4*rstraw);
   combination->Fit("track","RQ");
   Double_t true_offset=track->GetParameter(0);
   Double_t true_slope=track->GetParameter(1);
   Double_t true_offset_error=track->GetParError(0);
   Double_t true_slope_error=track->GetParError(1);
   Double_t chi2=track->GetChisquare();
   Double_t prob=track->GetProb();
   Int_t dof=track->GetNDF();
   if (dof > 0) {
      chi2=chi2/dof;
   }
   else {
      chi2=0;
      if(combination->GetN() > 1) prob=1;
   }
   // printf ("   True ntps=%d, chi2=%f, prob=%f\n",combination->GetN(),chi2,prob);
   /* printf ("True hits x=%f %f %f %f, y=%f %f %f %f, offset=%f +/-%f, slope=%f +/- %f, chi2/dof=%f, dof=%d\n",
             true_hitsx[0],true_hitsx[1],true_hitsx[2],true_hitsx[3],
             true_hitsy[0],true_hitsy[1],true_hitsy[2],true_hitsy[3],offset,offset_error,slope,slope_error,chi2,dof);*/
   // track->SetLineColor(4);
   // track->Draw("same");

   // fill histogram

   prob1->Fill(prob);
   chisq1->Fill(chi2);
   Schisq1->Fill(chi2);
   if (prob > prob_cut) {
      nsolutions = 1;
   }
   else {
      nsolutions = 0;
   }
   sol_true->Fill(angle,nsolutions);
   delete combination;
   delete track;


   // loop over all possible hit combinations and compute Chi2

   Int_t comb=0;
   Double_t nsolutions = 0;
   Double_t nsolutionspar = 0;
   for (j=0;j<2;j++) {
       for (jj=0;jj<2;jj++) {
            for (k=0;k<2;k++) {
                 for (kk=0;kk<2;kk++) {
                    for (kkk=0;kkk<2;kkk++) {
                    for (kkkk=0;kkkk<2;kkkk++) {
                     if (j == 0) {
                       select_hitsx[0] = true_hitsx[0];
                       select_hitsy[0] = true_hitsy[0];
                       if (errorsy[0] > 1.1*sigdy) {
                          j++;
                       }
                     }
                     else {
                       select_hitsx[0] = amb_hitsx[0];
                       select_hitsy[0] = amb_hitsy[0];
                     }
                     if (jj == 0) {
                       select_hitsx[1] = true_hitsx[1];
                       select_hitsy[1] = true_hitsy[1];
                       if (errorsy[1] > 1.1*sigdy) {
                          jj++;
                       }
                     }
                     else {
                       select_hitsx[1] = amb_hitsx[1];
                       select_hitsy[1] = amb_hitsy[1];
                     }
                     if (k == 0) {
                       select_hitsx[2] = true_hitsx[2];
                       select_hitsy[2] = true_hitsy[2];
                       if (errorsy[2] > 1.1*sigdy) {
                          k++;
                       }
                     }
                     else {
                       select_hitsx[2] = amb_hitsx[2];
                       select_hitsy[2] = amb_hitsy[2];
                     }

                     /*   NOTE: for different number of layers inner loops must be commented/uncommented by hand. */

                     if (kk == 0) {
                       select_hitsx[3] = true_hitsx[3];
                       select_hitsy[3] = true_hitsy[3];
                       if (errorsy[3] > 1.1*sigdy) {
                          kk++;
                       }
                     }
                     else {
                       select_hitsx[3] = amb_hitsx[3];
                       select_hitsy[3] = amb_hitsy[3];
                     }
                     if (kkk == 0) {
                       select_hitsx[4] = true_hitsx[4];
                       select_hitsy[4] = true_hitsy[4];
                       if (errorsy[4] > 1.1*sigdy) {
                          kkk++;
                       }
                     }
                     else {
                       select_hitsx[4] = amb_hitsx[4];
                       select_hitsy[4] = amb_hitsy[4];
                     }
                     if (kkkk == 0) {
                       select_hitsx[5] = true_hitsx[5];
                       select_hitsy[5] = true_hitsy[5];
                       if (errorsy[5] > 1.1*sigdy) {
                          kkkk++;
                       }
                     }
                     else {
                       select_hitsx[5] = amb_hitsx[5];
                       select_hitsy[5] = amb_hitsy[5];
                     }
                    comb++;
                    TGraphErrors *combination = new TGraphErrors(klayers,select_hitsx,select_hitsy,errorsx,errorsy);

                    // remove cells that miss hits

                    for (Int_t m=klayers-1;m>-1;m--) {
                       if (errorsy[m] > 1.1*sigdy) {
                          combination->RemovePoint(m);
                       }
                    }
                    TF1 *track = new TF1("track","pol1",xmin-rstraw,xmin+klayers*3*rstraw);
                    combination->Fit("track","RQ");
                    Double_t offset=track->GetParameter(0);
                    Double_t slope=track->GetParameter(1);
                    Double_t offset_error=track->GetParError(0);
                    Double_t slope_error=track->GetParError(1);
                    Double_t chi2=track->GetChisquare();
                    Int_t dof=track->GetNDF();
                    Double_t prob=track->GetProb();
                    if (dof > 0) {
                       chi2=chi2/dof;
                    }
                    else {
                       chi2=0;
                       if(combination->GetN() > 1) prob=1;
                   }
  // printf ("comb=%d, klayers=%d, x=%f %f %f %f %f %f, y=%f %f %f %f %f %f\n",comb,klayers,
  // select_hitsx[0],select_hitsx[1],select_hitsx[2],select_hitsx[3],select_hitsx[4],select_hitsx[5],
  // select_hitsy[0],select_hitsy[1],select_hitsy[2],select_hitsy[3],select_hitsy[4],select_hitsy[5]);
  // printf ("offset=%f +/-%f, slope=%f +/- %f, chi2/dof=%f, dof=%d\n",offset,offset_error,slope,slope_error,chi2,dof);
                    // track->SetLineColor(1);
                    // track->SetLineWidth(0);
                    // track->Draw("same");
                    // track->DrawCopy();
                    chisqN->Fill(chi2);
                    probN->Fill(prob);
                    // if (prob > prob_cut/10) SchisqN->Fill(chi2);
                    SchisqN->Fill(chi2);

                    // compute a chi2 of parameters
                    if (offset_error> 0 && slope_error> 0) {
                        Double_t chi2par = (true_offset-offset)*(true_offset-offset)/(offset_error*offset_error) +
                                       (true_slope-slope)*(true_slope-slope)/(slope_error*slope_error); 
                    }
                    else {
                        Double_t chi2par = 0;
                    }
                    chi2par = chi2par/2;
                    pchi2chi2->Fill(chi2,chi2par);
                    trkchi2->Fill(angle,chi2);
                    if (prob > prob_cut) {
                        nsolutions++;
                        parchi2->Fill(angle,chi2par);
                        pchi2->Fill(chi2par);
                        if (chi2par > 5) {
                           nsolutionspar++;
                           // printf ("klayers=%d, comb=%d, angle=%f, chi2=%f, dof=%d, chi2par=%f\n",klayers,comb,angle,chi2,dof,chi2par);
                        }
                    }
                    delete combination;
                    delete track;
                   }
                 }
               }
           }
       }
   }
   solN->Fill(nsolutions);
   goodN->Fill(nsolutionspar);
   sol_comb->Fill(angle,nsolutions);
   // printf ("angle=%f, comb=%d, nsolutions=%d\n",angle,comb,nsolutions);

   
   /*sprintf (string,"#Phi=%.1f\n",phi*180/kpi);
   t1 = new TLatex(0.1,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   
   sprintf(string,"plot_layers.eps");
   c1->SaveAs(string);
   sprintf(string,"plot_layers.gif");
   c1->SaveAs(string);
   c1->Close();*/


   delete[] wire_hitsx;
   delete[] wire_hitsy;
   delete[] true_hitsx;
   delete[] true_hitsy;
   delete[] amb_hitsx;
   delete[] amb_hitsy;
   delete[] select_hitsx;
   delete[] select_hitsy;
   delete[] errorsx;
   delete[] errorsy;

   return 0;
}