cut_L(int nrun = 0)
{
// 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/e06010_det_L_%d.root",nrun));

	TFile* f = new TFile(Form("ROOTfiles/e06010_det_L_Sieve.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')
	{
		gSystem->Exec(Form("cp -vf ./cutfiles/gcut_L_%d.root ./cutfiles/gcut_L_%d.root.old",nrun,nrun));
		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.06,0.06,800,-0.4,0.4);
    TH2F* h2 = new TH2F ("h2","theta_target vs. phi_target", 600, -0.06,0.06,700,-0.07,0.07);

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

// Set y target number
    int nfoil=6, 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; i < nfoil+fmin; i++)
    {
		TCutG* cutg, * tmpcut;
		
        T->Draw("L.tr.vz:L.tr.r_ph>>h1","abs(L.gold.dp)<.05","COLZ");
        c1->Update();
		TPad *cpad = gPad;

		cout<<"Testing "<<Form("fcut_L_%d",i)<<endl;
		f1->cd();
		TCutG* tmpcut=(TCutG*)gROOT->FindObject(Form("fcut_L_%d",i)); //looking for old cut definition
		f->cd();
		if (tmpcut)
		{
			cout<<Form("fcut_L_%d",i)<<" is found, using old one"<<endl;
			tmpcut->Draw("PL");
	// 		delete tmpcut; //delete old cut
		}
		else
		{
		
			cout <<"making cut for foil No."<<i<<", waiting ..."<<endl;
			
			cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
			c1->Update();
			cout <<"done!"<<endl;
			cutg->SetName(Form("fcut_L_%d",i)); //
  // set axises' name
			cutg->SetVarX("L.tr.r_ph");
			cutg->SetVarY("L.tr.vz");

  // 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
			f1->Write();
		}

  // make cuts for each holes
        T->Draw("L.gold.th:L.gold.ph>>h2",Form("fcut_L_%d&&abs(L.gold.dp)<.05",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();
                }
            }

			f1->Write();
        }

        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
}