void cutStudy(const int istation=1, const char *desc="Aug15-16")
{
  char fname[100]; 
  sprintf(fname, "result%s.root", desc);

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

  const int debug = 1;

  const int kNChips = 100000; // max number of chips so far, for any station
  const int kChipOffset = 100000;

  // 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 = 100;
  TH1F *h[kNHistMax];
  TH1F *hCutLow[kNHistMax];
  TH1F *hCutHigh[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

  // for the channel histos, keep a min and max version for each variable
  // first the min versions
  hTitle.push_back("minGain20mV (mV/fC)"); min[nHists] = 18; max[nHists] = 22; nbins[nHists] = 100; hFill.push_back("minGain20mV"); nHists++;
  hTitle.push_back("minGain30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 60; nbins[nHists] = 600; hFill.push_back("minGain30mV"); nHists++;

  hTitle.push_back("minRiseTime20mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("minRT20mV"); nHists++;
  hTitle.push_back("minRiseTime30mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("minRT30mV"); nHists++;
  hTitle.push_back("minPed20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 1023.5; nbins[nHists] = 1024; hFill.push_back("minPed20mV"); nHists++;
  hTitle.push_back("minPed30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 1023.5; nbins[nHists] = 1024; hFill.push_back("minPed30mV"); nHists++;
  hTitle.push_back("minRms20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("minRms20mV"); nHists++;
  hTitle.push_back("minRms30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("minRms30mV"); nHists++;
  hTitle.push_back("minPedAna20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("minPedAna20mV"); nHists++;
  hTitle.push_back("minPedAna30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("minPedAna30mV"); nHists++;
  hTitle.push_back("minRmsAna20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 3; hFill.push_back("minRmsAna20mV"); nHists++;
  hTitle.push_back("minRmsAna30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 3; hFill.push_back("minRmsAna30mV"); nHists++;

  const int kNHistChanMin = nHists - kNHistGlob; // this is the number of channel histograms - Min

  // then the max versions
  hTitle.push_back("maxGain20mV (mV/fC)"); min[nHists] = 18; max[nHists] = 22; nbins[nHists] = 100; hFill.push_back("maxGain20mV"); nHists++;
  hTitle.push_back("maxGain30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 60; nbins[nHists] = 600; hFill.push_back("maxGain30mV"); nHists++;

  hTitle.push_back("maxRiseTime20mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("maxRT20mV"); nHists++;
  hTitle.push_back("maxRiseTime30mV (ns)"); min[nHists] = 150; max[nHists] = 190; hFill.push_back("maxRT30mV"); nHists++;
  hTitle.push_back("maxPed20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 1023.5; nbins[nHists] = 1024; hFill.push_back("maxPed20mV"); nHists++;
  hTitle.push_back("maxPed30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 1023.5; nbins[nHists] = 1024; hFill.push_back("maxPed30mV"); nHists++;
  hTitle.push_back("maxRms20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("maxRms20mV"); nHists++;
  hTitle.push_back("maxRms30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 5; hFill.push_back("maxRms30mV"); nHists++;
  hTitle.push_back("maxPedAna20mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("maxPedAna20mV"); nHists++;
  hTitle.push_back("maxPedAna30mV (mV/fC)"); min[nHists] = -0.5; max[nHists] = 150.5; nbins[nHists] = 151; hFill.push_back("maxPedAna30mV"); nHists++;
  hTitle.push_back("maxRmsAna20mV (mV/fC)"); min[nHists] = 0; max[nHists] = 3; hFill.push_back("maxRmsAna20mV"); nHists++;
  hTitle.push_back("maxRmsAna30mV (mV/fC)"); min[nHists] = 0; max[nHists] = 3; hFill.push_back("maxRmsAna30mV"); nHists++;

  const int kNHistChan = nHists - kNHistGlob; // this is the number of channel histogram
  const int kNHistChanMax = kNHistChan - kNHistChanMin; // this is the number of channel histograms - Max

  if (debug>0) {
    cout << " nHists " << nHists << endl;
    cout << " kNHistGlob " << kNHistGlob << endl;
    cout << " kNHistChanMin " << kNHistChanMin
	 << " kNHistChanMax " << kNHistChanMax << endl;
  }
  float chanMin[kNHistChanMin] = {0};
  float chanMax[kNHistChanMax] = {0};

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

  // 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]);
    sprintf(id, "hCutLow%d", j);
    hCutLow[j] = new TH1F(id, hTitle[j].c_str(), nbins[j], min[j], max[j]);
    sprintf(id, "hCutHigh%d", j);
    hCutHigh[j] = new TH1F(id, hTitle[j].c_str(), nbins[j], min[j], max[j]);
  }

  // second, fill all histos
  int minChip = kNChips;
  int maxChip = 0;
  int iEntry[kNChips] = {0};
  int locChip = 0;  

  // we do a loop over entries first, to get the latest OK entry
  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;

	locChip = iChip - kChipOffset;  
	if ( (locChip>=0) && (locChip < kNChips) ) { // right range of chips

	  if (minChip>locChip) minChip = locChip;
	  if (maxChip<locChip) maxChip = locChip;    

	  iEntry[locChip] = i; // keep track of latest entry for each chip
	}

      } // sync
    } // station
  } // entries
  
  // then we loop over the last OK entry per chip						 
  for (int ichip = minChip; ichip<=maxChip; ichip++) {
    if (iEntry[ichip] > 0) {
      nbytes += tree->GetEntry(iEntry[ichip]);
      locChip = iChip - kChipOffset;  

      // 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; 
      
      for (k=0; k<kNHistChanMin; k++) chanMin[k] = max[k+kNHistGlob];
      for (k=0; k<kNHistChanMax; k++) chanMax[k] = min[k+kNHistGlob+kNHistChanMin];
      
      // find all the min and max values
      for (int ichan=0; ichan<32; ichan++) {
	k = 0;
	if (debug>1) cout << " ichan " << ichan << endl;
	if (chanMin[k] > fGain20mV[ichan]) chanMin[k] = fGain20mV[ichan];
	if (chanMax[k] < fGain20mV[ichan]) chanMax[k] = fGain20mV[ichan];
	k++;
	if (chanMin[k] > fGain30mV[ichan]) chanMin[k] = fGain30mV[ichan];
	if (chanMax[k] < fGain30mV[ichan]) chanMax[k] = fGain30mV[ichan];
	k++;
	if (chanMin[k] > fRT20mV[ichan]) chanMin[k] = fRT20mV[ichan];
	if (chanMax[k] < fRT20mV[ichan]) chanMax[k] = fRT20mV[ichan];
	k++;
	if (chanMin[k] > fRT30mV[ichan]) chanMin[k] = fRT30mV[ichan];
	if (chanMax[k] < fRT30mV[ichan]) chanMax[k] = fRT30mV[ichan];
	k++;
	if (chanMin[k] > fPed20mV[ichan]) chanMin[k] = fPed20mV[ichan];
	if (chanMax[k] < fPed20mV[ichan]) chanMax[k] = fPed20mV[ichan];
	k++;
	if (chanMin[k] > fPed30mV[ichan]) chanMin[k] = fPed30mV[ichan];
	if (chanMax[k] < fPed30mV[ichan]) chanMax[k] = fPed30mV[ichan];
	k++;
	if (chanMin[k] > fRms20mV[ichan]) chanMin[k] = fRms20mV[ichan];
	if (chanMax[k] < fRms20mV[ichan]) chanMax[k] = fRms20mV[ichan];
	k++;
	if (chanMin[k] > fRms30mV[ichan]) chanMin[k] = fRms30mV[ichan];
	if (chanMax[k] < fRms30mV[ichan]) chanMax[k] = fRms30mV[ichan];
	k++;
	
	if (chanMin[k] > fPedAna20mV[ichan]) chanMin[k] = fPedAna20mV[ichan];
	if (chanMax[k] < fPedAna20mV[ichan]) chanMax[k] = fPedAna20mV[ichan];
	k++;
	if (chanMin[k] > fPedAna30mV[ichan]) chanMin[k] = fPedAna30mV[ichan];
	if (chanMax[k] < fPedAna30mV[ichan]) chanMax[k] = fPedAna30mV[ichan];
	k++;
	if (chanMin[k] > fRmsAna20mV[ichan]) chanMin[k] = fRmsAna20mV[ichan];
	if (chanMax[k] < fRmsAna20mV[ichan]) chanMax[k] = fRmsAna20mV[ichan];
	k++;
	if (chanMin[k] > fRmsAna30mV[ichan]) chanMin[k] = fRmsAna30mV[ichan];
	if (chanMax[k] < fRmsAna30mV[ichan]) chanMax[k] = fRmsAna30mV[ichan];
      }

      // fill histos
      j = kNHistGlob;
      for (k=0; k<kNHistChanMin; k++) {
	if (debug>1) cout << " min-hist " << j+k << endl;
	h[j+k]->Fill(chanMin[k]);
      }
      j = kNHistGlob + kNHistChanMin;
      for (k=0; k<kNHistChanMax; k++) {
	if (debug>1) cout << " max-hist " << j+k << endl;
	h[j+k]->Fill(chanMax[k]);
      }

    } // ok entry
  } // chip

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

  gStyle->SetOptStat(0);

  int ibin = 1;
  float fraction = 1;
  float norm = 1;
  float edge = 100;
  const float kLimit = 0.995;
  char lowText[100];
  char highText[100];
  TText t(.5,.5,"");
  t.SetTextColor(kRed);
  //    tLow->SetTextSize(40);

  for (j = 0; j<kNHistGlob; j++) {
    if (debug>1) cout << " line " << __LINE__ << " j " << j << " " << hTitle[j].c_str() << endl; 
    
    norm = 100.0/h[j]->GetEntries();

    edge = kLimit * norm * h[j]->Integral(0, nbins[j]);
    for (ibin=0; ibin<=nbins[j]; ibin++) {
      fraction = norm * h[j]->Integral(ibin, nbins[j]);
      if (fraction < edge) {
	sprintf(lowText, "x=%.1f val %4.2f", h[j]->GetBinLowEdge(ibin), fraction);
	if (debug>0) cout << " lowText " << lowText << endl;
	edge = -1; // only print the first drop
      }
      hCutLow[j]->SetBinContent(ibin, fraction);
    }

    edge = kLimit * norm * h[j]->Integral(0, nbins[j]);
    for (ibin=nbins[j]; ibin>=0; ibin--) {
      fraction = norm * h[j]->Integral(0, ibin);
      if (fraction < edge) {
	sprintf(highText, "x=%.1f val %4.2f", h[j]->GetBinLowEdge(ibin) + h[j]->GetBinWidth(ibin), fraction);
	if (debug>0) cout << " highText " << highText << endl;
	edge = -1; // only print the first drop
      }
      hCutHigh[j]->SetBinContent(ibin, fraction);
    }

    c1->cd(1);
    hCutLow[j]->Draw();
    t.DrawTextNDC(0.6, 0.8, lowText);
    c1->Modified();
    sprintf(id, "pdf/cutStudy/%s/sta%d_%sLow.pdf", desc, istation, hFill[j].c_str());
    c1->SaveAs(id);

    hCutHigh[j]->Draw();
    t.DrawTextNDC(0.2, 0.8, highText);

    c1->Modified();
    sprintf(id, "pdf/cutStudy/%s/sta%d_%sHigh.pdf", desc, istation, hFill[j].c_str());
    c1->SaveAs(id);

  } // j

  j = kNHistGlob;
  for (k=0; k<kNHistChanMin; k++) {
    norm = 100.0/h[j+k]->GetEntries();

    edge = kLimit * norm * h[j]->Integral(0, nbins[j]);
    for (ibin=0; ibin<=nbins[j]; ibin++) {
      fraction = norm * h[j+k]->Integral(ibin, nbins[j]);
      if (fraction < edge) {
	sprintf(lowText, "x=%.1f val %4.2f", h[j+k]->GetBinLowEdge(ibin), fraction);
	if (debug>0) cout << " lowText " << lowText << endl;
	edge = -1; // only print the first drop
      }
      hCutLow[j+k]->SetBinContent(ibin, fraction);
    }

    c1->cd(1);
    hCutLow[j+k]->Draw();
    t.DrawTextNDC(0.6, 0.8, lowText);
    c1->Modified();
    sprintf(id, "pdf/cutStudy/%s/sta%d_%sLow.pdf", desc, istation, hFill[j+k].c_str());
    c1->SaveAs(id);

  }

  j = kNHistGlob + kNHistChanMin;
  for (k=0; k<kNHistChanMax; k++) {
    norm = 100.0/h[j+k]->GetEntries();

    edge = kLimit * norm * h[j+k]->Integral(0, nbins[j+k]);
    for (ibin=nbins[j+k]; ibin>=0; ibin--) {
      fraction = norm * h[j+k]->Integral(0, ibin);
      if (fraction < edge) {
	sprintf(highText, "x=%.1f val %4.2f", h[j+k]->GetBinLowEdge(ibin) + h[j]->GetBinWidth(ibin), fraction);
	if (debug>0) cout << " highText " << highText << endl;
	edge = -1; // only print the first drop
      }
      hCutHigh[j+k]->SetBinContent(ibin, fraction);
    }

    c1->cd(1);
    hCutHigh[j+k]->Draw();
    t.DrawTextNDC(0.2, 0.8, highText);
    c1->Modified();
    sprintf(id, "pdf/cutStudy/%s/sta%d_%sHigh.pdf", desc, istation, hFill[j+k].c_str());
    c1->SaveAs(id);
  }

}