{
// Script for showing the optics figures and generating cuts on Left Arm
// according to paw kumac from Nilanga
// by Yi Qiang, 4/28/2004

int nrun = 1467;
cout << "Run number?"<<endl;
cin >> nrun;

// Read root files
TFile* f = new TFile(Form("./rootfiles/optics_%d.root",nrun));

// Select Output format: root files or C files
char OTPT = 'root'; // options are 'root' or 'C'

// Create output root file if needed
if (OTPT == 'root')
{
  TFile *f1 = new TFile(Form("./cutfiles/gcut_L_%d.root",nrun),"UPDATE");
  f.cd(); // goback to f directory
}

// Define Canvas
TCanvas* c1 = new TCanvas("c1","Optics",600,800);
gStyle->SetPalette(1,0);

// Define Histograms
TH2F* h1 = new TH2F ("h1","ReactZ vs. Phi_rotate", 600, -0.03,0.03,800,-0.4,0.4);
TH2F* h2 = new TH2F ("h2","theta_target vs. phi_target", 600, -0.03,0.03,700,-0.07,0.07);

// Draw ReactZ vs. Phi_rotate
T->Draw("rpl.z:L.tr.r_ph>>h1","abs(L.gold.dp)<.05","COLZ"); // need finer delta cut later
c1->Update();

// Set y target number
int nfoil=0, fmin =1;
cout << "How many foils?"<<endl;
cin >> nfoil;
if (nfoil<=0) return ();
cout << "foil number begins with=?"<<endl;
cin >> fmin;


// output ps files
//c1->Print(Form("./cutfiles/gcut_L_%d.ps[",nrun));

// Choose the foil you want to make cut)
for (int i = fmin, CutG* cutg, CutG* tmpcut; i < nfoil+fmin; i++)
{
  T->Draw("rpl.z:L.tr.r_ph>>h1","abs(L.gold.dp)<.05","COLZ"); 
  c1->Update();  
  cout <<"making cut for foil No."<<i<<", waiting ..."<<endl;
  cutg=(TCutG*)gPad->WaitPrimitive("CUTG","CutG"); // making cut, store to CUTG
  c1->Update();
  cout <<"done!"<<endl;
  TCutG* tmpcut=(TCutG*)gROOT->FindObject(Form("fcut_L_%d",i)); //looking for old cut definition
  delete tmpcut; //delete old cut
  cutg->SetName(Form("fcut_L_%d",i)); //

  // set axises' name
  cutg->SetVarX("L.tr.r_ph");
  cutg->SetVarY("rpl.z");

  // output cut to disk
  if (OTPT == 'root'){
    f1.cd();
    cutg->Write("",TObject::kOverwrite); // Overwrite old cut
    f.cd();
  }
  if (OTPT == 'C'){
    ofstream out(Form("./cutfiles/%d/fcut_L_%d.C",nrun,i));
    cutg->SavePrimitive(out,"");
    out.close();
  }

  // output ps file
  if (i==fmin)  c1->Print(Form("./cutfiles/gcut_L_%d.ps(",nrun));
  else c1->Print(Form("./cutfiles/gcut_L_%d.ps",nrun)); // insert page into ps file


  // make cuts for each holes
  T->Draw("L.gold.th:L.gold.ph>>h2",Form("fcut_L_%d&&abs(L.gold.dp)<.05&&D.evtypebits&0x08",i),"COLZ");
  c1->Update();
  
  int ncol = 0, cmin=1;
  cout << "How many columns of holes do you see?" <<endl;
  cin >> ncol;
  if (ncol<=0) continue;
  if (ncol!=7){
    cout << "column number starts with=?" <<endl;
    cin >> cmin;
  }

  // make cuts from the first available column
  for (int j = cmin; j < ncol+cmin; j++) //column number j
    {
      
      int nhol=0, hmin=1;
      cout << "How many holes in this No."<<j<<" column?" <<endl;
      cin >> nhol;
      if (nhol <= 0) continue;
      if (nhol!=7){
	cout << "hole number starts with=?" <<endl;
	cin >> hmin;
      }

      //make cuts for each holes from begin number
      for (int k=hmin; k < (nhol+hmin); k++)
	{
          // define cuts
	  cout <<"making cut for foil="<<i<<", column="<<j<<", No."<<k<<" hole: waiting ..."<<endl;
	  cutg=(TCutG*)gPad->WaitPrimitive("CUTG","CutG");
	  c1->Update();
	  cout <<"done!"<<endl;
	  TCutG* tmpcut=(TCutG*)gROOT->FindObject(Form("hcut_L_%d_%d_%d",i,j,k));
	  delete tmpcut;
	  cutg->SetName(Form("hcut_L_%d_%d_%d",i,j,k));
 
	  // set axises
	  cutg->SetVarX("L.gold.ph");
	  cutg->SetVarY("L.gold.th");

	  // output
	  if (OTPT == 'root'){
	    f1.cd();
	    cutg->Write("",TObject::kOverwrite);
	    f.cd();
	  }
	  if (OTPT == 'C'){
	    ofstream out(Form("./cutfiles/%d/hcut_L_%d_%d_%d.C",nrun,i,j,k));
	    cutg->SavePrimitive(out,"");
	    out.close();
	  }
	}
    }

  if(i==fmin+nfoil-1) c1->Print(Form("./cutfiles/gcut_L_%d.ps)",nrun));
  else c1->Print(Form("./cutfiles/gcut_L_%d.ps",nrun)); // insert page into ps file
}
//gROOT->Reset();
f1.Close();
//c1->Print(Form("./cutfiles/gcut_L_%d.ps]",nrun)); // close ps file
}