void Main_layers_graphics(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 sigdy = 0.2;
   // Double_t rstraw=2;
   #define npts 1;
   Double_t phimin=20.3;
   Double_t phimax=phimin+1;
   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
   // kgeom=2 is axial geometry with close-packing

   Int_t kgeom=2;

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

   //
   // 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);

   // 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,5,0,5);

   TH1F *chisq1 = new TH1F("chisq1","Chi2 of fit to true hits",100,0,5);
   TH1F *chisqN = new TH1F("chisqN","Chi2 of fit to all combinations",100,0,5);
   TH1F *prob1 = new TH1F("prob1","Prob of fit to true hits",100,0,1);
   TH1F *probN = new TH1F("probN","Prob of fit to all combinations",100,0,1);

   // 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)) {
              printf ("j=%d, phi=%f\n",j,phi);
          }
          else {
           printf ("*** Main_layers_graphics: plot_layers failed!\n" ); 
          }
   }


   // plot histogram

   c0->cd(1);
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);
   chisq1->Draw();
   c0->cd(2);
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);
   chisqN->Draw();

   c0->cd(3);
   c0_3->SetBorderMode(0);
   c0_3->SetFillColor(0);
   prob1->Draw();
   c0->cd(4);
   c0_4->SetLogy();
   c0_4->SetBorderMode(0);
   c0_4->SetFillColor(0);
   probN->Draw();
   sum = probN->Integral();
   sprintf (string,"sum=%.1f trials\n",sum);
   t1 = new TLatex(0.5,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

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


   // save histogram to file

   sprintf(string,"Main_layers_graphics_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();*/
   out->Close();



}
// 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.

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

// gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
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 *prob1 = (TH1F *) gROOT->FindObject("prob1");
TH1F *probN = (TH1F *) gROOT->FindObject("probN");
// 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];

   // 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 kpi=3.14159;
   Double_t phi=angle*kpi/180;
   // Double_t xoffset = 1;
   // Double_t scale=1.2;
   Double_t xoffset = -6;
   Double_t scale=0.8;

   Double_t ymin=-4*rstraw;

   // Convert to actual cylindrical geometry with circle of integral number of tubes
   // ntubes will be the number of tubes in the first layer.

   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 = shift + Dphi/2;
         }
         elseif (j==2 || j==4) {
            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);

      printf ("angle=%f, r1=%f, Dphi=%f, ntubes=%d, arc=%f, arc1=%f, phidif=%f, shift=%f\n",angle,r1,Dphi,ntubes,arc,arc1,phidif,shift);

      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;
         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;

      // the next statement allows for geometrical inefficiencies to be included in the study.

      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();
   }
   printf ("wire_hitsx1=%f, wire_hitsx2=%f, wire_hitsx3=%f, wire_hitsx4=%f\n",wire_hitsx[0],wire_hitsx[1],wire_hitsx[2],wire_hitsx[3]);
   printf ("true_hitsy1=%f, true_hitsy2=%f, true_hitsy3=%f, true_hitsy4=%f\n",true_hitsy[0],true_hitsy[1],true_hitsy[2],true_hitsy[3]);
   printf ("errorsy1=%f, errorsy2=%f, errorsy3=%f, errorsy4=%f\n",errorsy[0],errorsy[1],errorsy[2],errorsy[3]);
   printf ("errorsy5=%f, errorsy6=%f\n",errorsy[4],errorsy[5]);

   // fit the true hits

   Double_t nsolutions = 0;
   TGraphErrors *combination = new TGraphErrors(klayers,true_hitsx,true_hitsy,errorsx,errorsy);
   TF1 *track = new TF1("track","pol1",rmin-2*rstraw,rmin+klayers*4*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();
   Double_t prob=track->GetProb();
   Int_t dof=track->GetNDF();
   chi2=chi2/dof;
   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");
   // delete combination;
   // delete track;

   // fill histogram

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

   // plot fits to wires only


   TGraphErrors *combination = new TGraphErrors(klayers,wire_hitsx,wire_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","R");
   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();
   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 ("Wire hits x=%f %f %f %f, y=%f %f %f %f, offset=%f +/-%f, slope=%f +/- %f, chi2/dof=%f, dof=%d\n",
           wire_hitsx[0],wire_hitsx[1],wire_hitsx[2],wire_hitsx[3],
           wire_hitsy[0],wire_hitsy[1],wire_hitsy[2],wire_hitsy[3],offset,offset_error,slope,slope_error,chi2,dof);
   track->SetLineColor(2);
   // track->Draw("same");
   // delete combination;
   // delete track;
   
   sprintf (string,"#Phi=%.1f deg\n",phi*180/kpi);
   t1 = new TLatex(0.05,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   sprintf (string,"rstraw=%.1f cm\n",rstraw);
   t1 = new TLatex(0.05,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   sprintf (string,"sigdy=%.2f cm\n",sigdy);
   t1 = new TLatex(0.05,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   sprintf (string,"kgeom=%d\n",kgeom);
   t1 = new TLatex(0.05,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();


   // loop over all possible hit combinations and compute Chi2

   Int_t comb=0;
   Double_t nsolutions = 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];
                     }
                     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, x=%f %f %f %f, y=%f %f %f %f, offset=%f +/-%f, slope=%f +/- %f, prob=%f, chi2/dof=%f, dof=%f\n",comb,select_hitsx[0],select_hitsx[1],select_hitsx[2],select_hitsx[3],
           select_hitsy[0],select_hitsy[1],select_hitsy[2],select_hitsy[3],offset,offset_error,slope,slope_error,prob,chi2,dof);
                    track->SetLineColor(1);
                    // track->SetLineWidth(0);
                    if (prob > prob_cut && comb != 1) track->Draw("same");
                    chisqN->Fill(chi2);
                    probN->Fill(prob);
                    if (prob > prob_cut) nsolutions++;
                    // delete combination;
                    // delete track;
                  }
                }
              }
           }
       }
   }
   sol_comb->Fill(angle,nsolutions);

   
   sprintf(string,"plot_layers_c1_phi%.0f_l%d_g%d.eps",angle,klayers,kgeom);
   printf("Display file=%s\n",string);
   c1->SaveAs(string);
   sprintf(string,"plot_layers_c1_phi%0.f_l%d_g%d.gif",angle,klayers,kgeom);
   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;
}