void plot(const int istation=1, const int logy=1, const char *desc="Aug15-16")
{
  char descrip[100];
  if (logy==1) sprintf(descrip,"%s",desc);
  else sprintf(descrip,"%s_liny",desc);

  char fname[100]; 
  sprintf(fname, "result%s.root", desc);

  TFile *inFile = TFile::Open(fname);

  const int debug = 0;

  // constants for arrays 
  const int kNChan = 32;
  const int kNQR = 26;
  const int kNVRef = 6;
  const int kNACurr = 5;
  const int kNSync = 8;

  TTree *tree = (TTree*) inFile->Get("tree");

  //Declaration of leaves types
  //Declaration of leaves types
  Char_t          filename[200];
  Char_t          datetime[100];
  Char_t          cQR[kNQR];
  Int_t           iChip;
  Int_t           iSta;
  Float_t         fVRef[kNVRef];
  Float_t         fVCurr[kNACurr];
  Float_t         fACurr[kNACurr];
  Float_t         fGain20mV[32];
  Float_t         fGain30mV[32];
  Float_t         fRT20mV[32];
  Float_t         fRT30mV[32];
  Float_t         fPed20mV[32];
  Float_t         fPed30mV[32];
  Float_t         fRms20mV[32];
  Float_t         fRms30mV[32];
  Float_t         fPedAna20mV[32];
  Float_t         fPedAna30mV[32];
  Float_t         fRmsAna20mV[32];
  Float_t         fRmsAna30mV[32];
  Float_t         fX20mVeven;
  Float_t         fX20mVodd;
  Float_t         fX30mVeven;
  Float_t         fX30mVodd;
  Int_t           iSync[kNSync];
  Int_t           iDFT;
  Int_t           iJtag;
  Int_t           iMem;
  Float_t         fRingOsc;
  Int_t           i2cErr;
  Int_t           iErrCode;
  Int_t           iRCode;
  
  // Set branch addresses.
  tree->SetBranchAddress("filename",filename);
  tree->SetBranchAddress("datetime",datetime);
  tree->SetBranchAddress("cQR",&cQR);
  tree->SetBranchAddress("iChip",&iChip);
  tree->SetBranchAddress("iSta",&iSta);
  tree->SetBranchAddress("fVRef",fVRef);
  tree->SetBranchAddress("fVCurr",fVCurr);
  tree->SetBranchAddress("fACurr",fACurr);
  tree->SetBranchAddress("fGain20mV",fGain20mV);
  tree->SetBranchAddress("fGain30mV",fGain30mV);
  tree->SetBranchAddress("fRT20mV",fRT20mV);
  tree->SetBranchAddress("fRT30mV",fRT30mV);
  tree->SetBranchAddress("fPed20mV",fPed20mV);
  tree->SetBranchAddress("fPed30mV",fPed30mV);
  tree->SetBranchAddress("fRms20mV",fRms20mV);
  tree->SetBranchAddress("fRms30mV",fRms30mV);
  tree->SetBranchAddress("fPedAna20mV",fPedAna20mV);
  tree->SetBranchAddress("fPedAna30mV",fPedAna30mV);
  tree->SetBranchAddress("fRmsAna20mV",fRmsAna20mV);
  tree->SetBranchAddress("fRmsAna30mV",fRmsAna30mV);
  tree->SetBranchAddress("fX20mVeven",&fX20mVeven);
  tree->SetBranchAddress("fX20mVodd",&fX20mVodd);
  tree->SetBranchAddress("fX30mVeven",&fX30mVeven);
  tree->SetBranchAddress("fX30mVodd",&fX30mVodd);
  tree->SetBranchAddress("iSync",iSync);
  tree->SetBranchAddress("iDFT",&iDFT);
  tree->SetBranchAddress("iJtag",&iJtag);
  tree->SetBranchAddress("iMem",&iMem);
  tree->SetBranchAddress("fRingOsc",&fRingOsc);
  tree->SetBranchAddress("i2cErr",&i2cErr);
  tree->SetBranchAddress("iErrCode",&iErrCode);
  tree->SetBranchAddress("iRCode",&iRCode);
  
  //     This is the loop skeleton
  //       To read only selected branches, Insert statements like:
  // tree->SetBranchStatus("*",0);  // disable all branches
  // TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname
  Long64_t nentries = tree->GetEntries();

  const int kNHistMax = 50;
  TH1F *h[kNHistMax];
  int nbins[kNHistMax] = { 100 };
  float min[kNHistMax] = { 0 };
  float max[kNHistMax] = { 100 };

  // 
  vector <string> hTitle;
  vector <string> hFill;

  int nHists = 0;
  hTitle.push_back("VDDA (V)"); min[nHists] = 1.2; max[nHists] = 1.3; hFill.push_back("fVRef[0]"); nHists++;
  hTitle.push_back("VREF (V)"); min[nHists] = 1.07; max[nHists] = 1.11; hFill.push_back("fVRef[1]"); nHists++;
  hTitle.push_back("VDD (V)"); min[nHists] = 1.2; max[nHists] = 1.3; hFill.push_back("fVRef[2]"); nHists++;
  hTitle.push_back("V450 (V)"); min[nHists] = 0.45; max[nHists] = 0.55; hFill.push_back("fVRef[3]"); nHists++;
  hTitle.push_back("V600 (V)"); min[nHists] = 0.6; max[nHists] = 0.7; hFill.push_back("fVRef[4]"); nHists++;
  hTitle.push_back("V750 (V)"); min[nHists] = 0.75; max[nHists] = 0.85; hFill.push_back("fVRef[5]"); nHists++;
  hTitle.push_back("V750 - V450 (V)"); min[nHists] = 0.2; max[nHists] = 0.4; hFill.push_back("fVRefDiff53"); nHists++;
  hTitle.push_back("V750 - V600 (V)"); min[nHists] = 0.05; max[nHists] = 0.25; hFill.push_back("fVRefDiff54"); nHists++;

  hTitle.push_back("FE1_CURR (V)"); min[nHists] = 0.25; max[nHists] = 0.55; hFill.push_back("fVCurr[0]"); nHists++;
  hTitle.push_back("FE2_CURR (V)"); min[nHists] = 0.2; max[nHists] = 0.4; hFill.push_back("fVCurr[1]"); nHists++;
  hTitle.push_back("AD_CURR (V)"); min[nHists] = 0.0; max[nHists] = 0.2; hFill.push_back("fVCurr[2]"); nHists++;
  hTitle.push_back("DR_CURR (V)"); min[nHists] = 0.1; max[nHists] = 0.2; hFill.push_back("fVCurr[3]"); nHists++;
  hTitle.push_back("DG_CURR (V)"); min[nHists] = 0.00; max[nHists] = 0.04; hFill.push_back("fVCurr[4]"); nHists++;
  hTitle.push_back("FE1_CURR (mA)"); min[nHists] = 80; max[nHists] = 110; hFill.push_back("fACurr[0]"); nHists++;
  hTitle.push_back("FE2_CURR (mA)"); min[nHists] = 50; max[nHists] = 80; hFill.push_back("fACurr[1]"); nHists++;
  hTitle.push_back("AD_CURR (mA)"); min[nHists] = 0; max[nHists] = 50; hFill.push_back("fACurr[2]"); nHists++;
  hTitle.push_back("DR_CURR (mA)"); min[nHists] = 30; max[nHists] = 60; hFill.push_back("fACurr[3]"); nHists++;
  hTitle.push_back("DG_CURR (mA)"); min[nHists] = 10; max[nHists] = 50; hFill.push_back("fACurr[4]"); nHists++;

  hTitle.push_back("X20mVeven"); min[nHists] = 0; max[nHists] = 0.2; hFill.push_back("fX20mVeven"); nHists++;
  hTitle.push_back("X30mVeven"); min[nHists] = 0; max[nHists] = 0.2; hFill.push_back("fX30mVeven"); nHists++;
  hTitle.push_back("X20mVodd"); min[nHists] = 0; max[nHists] = 0.2; hFill.push_back("fX20mVodd"); nHists++;
  hTitle.push_back("X30mVodd"); min[nHists] = 0; max[nHists] = 0.2; hFill.push_back("fX30mVodd"); nHists++;

  hTitle.push_back("RingOsc"); min[nHists] = 130; max[nHists] = 170; hFill.push_back("fRingOsc"); nHists++;

  const int kNHistGlob = nHists; // this is the number of histograms where there is one entry per chip

  /*  
hTitle.push_back("Gain20mV (mV/fC)"); min[nHists] = 18; max[nHists] = 22; nbins[nHists] = 100*(max[nHists]-min[nHists]); hFill.push_back("fGain20mV"); nHists++;
  hTitle.push_back("Gain30mV (mV/fC)"); min[nHists] = 27; max[nHists] = 33; nbins[nHists] = 100*(max[nHists]-min[nHists]); hFill.push_back("fGain30mV"); nHists++;
  */
hTitle.push_back("Gain20mV (mV/fC)"); min[nHists] = 18; max[nHists] = 22; hFill.push_back("fGain20mV"); nHists++;
  hTitle.push_back("Gain30mV (mV/fC)"); min[nHists] = 27; max[nHists] = 33; hFill.push_back("fGain30mV"); nHists++;

  hTitle.push_back("RiseTime20mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("fRT20mV"); nHists++;
  hTitle.push_back("RiseTime30mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("fRT30mV"); nHists++;
  hTitle.push_back("Ped20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("fPed20mV"); nHists++;
  hTitle.push_back("Ped30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("fPed30mV"); nHists++;
  hTitle.push_back("Rms20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("fRms20mV"); nHists++;
  hTitle.push_back("Rms30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("fRms30mV"); nHists++;
  hTitle.push_back("PedAna20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("fPedAna20mV"); nHists++;
  hTitle.push_back("PedAna30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("fPedAna30mV"); nHists++;
  hTitle.push_back("RmsAna20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("fRmsAna20mV"); nHists++;
  hTitle.push_back("RmsAna30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("fRmsAna30mV"); nHists++;

  const int kNHistChan = nHists - kNHistGlob; // this is the number of histograms where there is one entry per channel
  if (debug>1) cout << " nHists " << nHists << endl;

  // Print Strings stored in Vector
  if (debug>0) { 
    for (int ii=0; ii<hTitle.size(); ii++) {   
      cout << " ii " << ii << " title " << hTitle[ii] << " - fill " << hFill[ii] << "\n"; 
    }
  }

  int j = 0;
  int k = 0;
  char id[100];
  char fillstr[100];

  const double kNSigma = 4;
  FILE* fout = fopen( Form("sta%d_cuts2.txt", istation), "w");
  fprintf(fout, "# Min and Max cut values generated as Mean +- %1.0f*RMS -- Station %d\n", kNSigma, istation);
  fprintf(fout, "%20s %8s %8s %8s %8s\n", "Variable", "Mean", "RMS", "MinCut", "MaxCut");

  // first, create all histos
  for (j = 0; j<nHists; j++) {
    sprintf(id, "h%d", j);
    if (nbins[j]<100) nbins[j]=100; 
    h[j] = new TH1F(id, hTitle[j].c_str(), nbins[j], min[j], max[j]);
  }

  // second, fill all histos
  Long64_t nbytes = 0;
  if (debug>1) cout << " nentries " << nentries << endl;
  for (Long64_t i=0; i<nentries;i++) {
    if (debug>1) cout << " get entry i " << i << endl;
    nbytes += tree->GetEntry(i);
    if (debug>1) cout << " entry i " << i << " iSta " << iSta << endl;

    if (iSta == istation) { // selected station
      if (debug>1) cout << " selected entry i " << i << " iSta " << iSta << endl;
      int nSyncFail = 0;
      int nSyncOK = 0;
      for (k=0; k<kNSync; k++) {
	if ( (200<iSync[k]) && (iSync[k]<300) ) { nSyncOK++; }
	else { nSyncFail++; }
      }
      if (debug>1) cout << " nSyncOK " << nSyncOK << endl;
      if ( nSyncOK == kNSync ) {
	if (debug>1) cout << " ok entry for station " << iSta << endl;
	
	// need to fill histos in the same order as they were declared above
	j=0;
	if (debug>1) cout << " line " << __LINE__ << " j " << j << endl; 
	for (k=0; k<kNVRef; k++) {
	  h[j]->Fill(fVRef[k]); j++;
	}
	h[j]->Fill(fVRef[5] - fVRef[3]); j++;
	h[j]->Fill(fVRef[5] - fVRef[4]); j++;

	if (debug>1) cout << " line " << __LINE__ << " j " << j << endl; 
	for (k=0; k<kNACurr; k++) {
	  h[j]->Fill(fVCurr[k]); j++;
	}
	if (debug>1) cout << " line " << __LINE__ << " j " << j << endl; 
	for (k=0; k<kNACurr; k++) {
	  h[j]->Fill(fACurr[k]); j++;
	}
	if (debug>1) cout << " line " << __LINE__ << " j " << j << endl; 
	
	h[j]->Fill(fX20mVeven); j++;
	h[j]->Fill(fX30mVeven); j++;
	h[j]->Fill(fX20mVodd); j++;
	h[j]->Fill(fX30mVodd); j++;
	
	h[j]->Fill(fRingOsc); j++;
	
	if (j != kNHistGlob) cerr << " line " << __LINE__ << " mismatch! j " << j << " kNHistGlob " << kNHistGlob << endl; 
	
	// all histos with entries per channel
	for (int ichan=0; ichan<32; ichan++) {
	  j = kNHistGlob;
	  if (debug>1) cout << " ichan " << ichan << endl;
	  h[j]->Fill(fGain20mV[ichan]); j++;
	  h[j]->Fill(fGain30mV[ichan]); j++;
	  h[j]->Fill(fRT20mV[ichan]); j++;
	  h[j]->Fill(fRT30mV[ichan]); j++;
	  h[j]->Fill(fPed20mV[ichan]); j++;
	  h[j]->Fill(fPed30mV[ichan]); j++;
	  h[j]->Fill(fRms20mV[ichan]); j++;
	  h[j]->Fill(fRms30mV[ichan]); j++;
	  h[j]->Fill(fPedAna20mV[ichan]); j++;
	  h[j]->Fill(fPedAna30mV[ichan]); j++;
	  h[j]->Fill(fRmsAna20mV[ichan]); j++;
	  h[j]->Fill(fRmsAna30mV[ichan]); j++;
	}
	if (j != nHists) cerr << " line " << __LINE__ << " mismatch! j " << j << " nHists " << nHists << endl; 
      } // OK event (Sync)
    } // station
  }

  // third, create all plots, and check cut info 
  if (debug>1) cout << " create canvas " << endl;
  TCanvas *c1 = new TCanvas("c1","c1");
  if (logy==1) c1->SetLogy(1);
  if (debug>1) cout << " canvas created " << endl;

  for (j = 0; j<nHists; j++) {
    if (debug>1) cout << " line " << __LINE__ << " j " << j << endl; 
    c1->cd(1);
    h[j]->Draw();
    
    c1->Modified();
    sprintf(id, "png/cuts2/%s/sta%d_%s.png", descrip, istation, hFill[j].c_str());
    c1->SaveAs(id);
    sprintf(id, "pdf/cuts2/%s/sta%d_%s.pdf", descrip, istation, hFill[j].c_str());
    c1->SaveAs(id);

    double mean = h[j]->GetMean();
    double rms = h[j]->GetRMS();
    double minCut = mean - kNSigma*rms;
    double maxCut = mean + kNSigma*rms;
    if (minCut<min[j]) minCut=min[j]; 
    if (maxCut>max[j]) maxCut=max[j]; 
    fprintf(fout, "%20s %8.3f %8.3f %8.3f %8.3f\n", hTitle[j].c_str(), mean, rms, minCut, maxCut);
  }

  fclose(fout);

}