{

char tdchar[48],omchar[48],tdcchar[48];


#include <iomanip.h>
#include <fstream.h>

//char tdcchar[24],tdchar[24],chchar[24];
ofstream out_file("om.dat",ios::out);

hh1 = new TH1F("Runs","Runs",700.,0.0,700.);
moni->Draw("run >> Runs");
Stat_t m = hh1->GetMinimum();
cerr << "Max:"<<m<<'\n';

TCanvas *c1 = new TCanvas("c1", "c1",10,10,500,500); 

sprintf(tdchar,"noiseeventrate >> noiserate");
h1 = new TH1F("noiserate","NoiseRate",1000.,50.0,2000.);
moni->Draw(tdchar); 

//
// Fit Gauss two the two peaks to find global 
// mean noise rate 
//

TF1 *gfit = new TF1("gfit","[0]*exp(-0.5*((x-[1])/[2])**2.)");

gfit->SetParameters(1000,350,35.);
gfit->SetLineColor(3);
h1->Fit("gfit","L","",280,500);
Float_t ampl_low = gfit->GetParameter(0);
Float_t mean_low = gfit->GetParameter(1);
Float_t sigma_low = gfit->GetParameter(2);
Float_t chi2_low = gfit->GetChisquare();
 
// Fit High Noise Sample.
gfit->SetParameters(1500,1450.,150.);
gfit->SetLineColor(3);

h1->Fit("gfit","L","",1100,1800);
Float_t ampl_high = gfit->GetParameter(0);
Float_t mean_high = gfit->GetParameter(1);
Float_t sigma_high = gfit->GetParameter(2);
Float_t chi2_high = gfit->GetChisquare();

// end global means



Stat_t NoiseRate[303][701];
Stat_t TDCRate[303][701];
Float_t s = 5.;
Float_t s2 = 10.;
Int_t  status_array[303][701];
Int_t  status_array2[303][701];

Int_t bad_om[303];
Int_t bad_array[303];
Int_t bad_run[701];

Float_t Ratio[303];


// 5 sigma map
hxy     = new TH2F("hxy","OM status: 5 sig",700,1.,700.,302,1.,302.);

// 10 sigma map 
hxy2     = new TH2F("hxy2","OM status: 10 sig ",700,1.,700.,302,1.,302.);

// TDC over Noise Ratio MAP
hxy3     = new TH2F("hxy3","TDC over Noise",700,1.,700.,302,1.,302.);


// loop over OMs

h1 = new TH1F("noiserate","Rate",5000.,0.0,5000.);
h2 = new TH1F("TDCrate","tdc",10000.,0.0,10000.);

out_file << "OM STATUS LIST: " << '\n';

 for(Int_t i=1;i<= 302;i++){

  sprintf(tdchar,"noiseeventrate[%1d] >> noiserate",i);
  moni->Draw(tdchar); 
  Stat_t meanyr = h1->GetMean();
  Stat_t sigmayr = h1 ->GetRMS();

  Int_t count = 0;
  Int_t help = 0;
  bad_om[i] = 0;

 
  for(Int_t r=1;r<=700;r++){

    // loop over runs

    h1.Reset();
    h2.Reset();

    status_array[i][r] = 0.;
    bad_array[i] = 0;
    bad_run[r] = 0;

    //HOW DO I GET THE NOISEVENTRATE FOR RUN X for OM Y ?

 
    sprintf(tdcchar,"tdc1stedgeeventrate[%1d] >> TDCrate",i);
    sprintf(tdchar,"noiseeventrate[%1d] >> noiserate",i);
    sprintf(omchar,"run==%1d",r);

    // cerr <<tdchar;
    // cerr <<chchar;

    moni->Draw(tdchar,omchar); 
    Stat_t l = h1->GetMean();
    NoiseRate[i][r] = l;
    
    moni->Draw(tdcchar,omchar);
    Stat_t t = h2->GetMean();
    TDCRate[i][r] = t;




    if (l != 0) Float_t ratio = (TDCRate[i][r] / NoiseRate[i][r]) ;
  

    Float_t mean = mean_high;
    Float_t sigma = sigma_high;


    if (i < 86) {
      mean = mean_low;
      sigma = sigma_low;
    }

    Int_t cond1 = 0;
    Int_t cond2 = 0;
    Int_t cond4 = 0;


    if ((NoiseRate[i][r] >= mean +s*sigma)||(NoiseRate[i][r] <= mean-s*sigma) || (NoiseRate[i][r] >= meanyr +s*sigmayr)||(NoiseRate[i][r] <= meanyr-s*sigmayr)) cond1 = 1;

    if ((NoiseRate[i][r] >= mean +s2*sigma)||(NoiseRate[i][r] <= mean-s2*sigma) || (NoiseRate[i][r] >= meanyr +s2*sigmayr)||(NoiseRate[i][r] <= meanyr-s2*sigmayr)) cond2 = 1;

    if (ratio < 1.0) cond4 = 1;



    // if noise >< 5 sig (global/individual) or noise >< 3 sig in to consecutive runs:

    if (cond1 == 1) 
      { 	hxy->Fill(r,i,1.);
 
      status_array[i][r] = 1.;         
 
      if (cond2 == 1) 
	hxy2->Fill(r,i,1.); 
      else 
	hxy2->Fill(r,i,3.);

      //        bad run counter
      count = count+1;

      if (help != i) {
	out_file << '\n' << "OM #" << " " << i << " " << '\n'; 
	out_file << r << " " << "(" << NoiseRate[i][r] << ")" << " ";
	if (cond2 == 1) {
         out_file << "*" << " ";
	 status_array2[i][r] = 1.; 
	}
	help = i;
      }
      else  
	{
	  out_file << r << " " << "(" << NoiseRate[i][r] << ")" << " "; 
	  
	  if (cond2 == 1) {
	    out_file << "*" << " ";
            status_array2[i][r] = 1.; 
	  }
	}
	//}
	//    else
	//    hxy->Fill(r,i,3.);
      }
    else {
      hxy->Fill(r,i,3.);
      hxy2->Fill(r,i,3.);
    }
  

    if (cond4 == 1) {
      out_file << "T" << " ";
      hxy3->Fill(r,i,1.);
    } 
    else
      hxy3->Fill(r,i,3.);


    cerr << "." << '\n';
  }

  // Ratio of Runs where OM is bad:

  Ratio[i] = count/700.; 

  if (Ratio[i] > 0.5) {
    bad_array[i] = 999;
    bad_om[i] = i;
    out_file << status_array[i] << " "; 
  }

 }

out_file <<'\n' << '\n' << "List of Bad OMs:"  <<'\n';

for(Int_t i=1;i<= 302;i++){
  if (bad_om[i] != 0) out_file << bad_om[i] << " ";
}

out_file <<'\n' << '\n' << "List of BAD RUNs: "  <<'\n';

for (Int_t r = 1; r<=700;r++){
  
  Int_t counter2 = 0;
  Int_t counter3 = 0;  

  for (Int_t i =1; i <= 302; i++){
    counter2 = counter2 + status_array[i][r];
    counter3 = counter3 + status_array2[i][r];

 }

  if ((counter2>151) || (counter3>60)) { 
    bad_run[r] = r;
    out_file << bad_run[r] << " ";
  }

}

out_file <<'\n' << '\n' << "List of 10 sig OMs:"  <<'\n';  

for (Int_t r = 1; r<=700;r++){
  for (Int_t i =1; i <= 302; i++){
    if ((status_array2[i][r] == 1) && (bad_run[r] != r) && (bad_om[i] !=i))
     out_file << i << "(" << r<< ")" << " ";
  }
}



out_file.close();   

TFile g("OMstatus_n_5_10_TDC.root","recreate");
hxy ->Write("hxy");
hxy2 ->Write("hxy2");
hxy3 ->Write("hxy3");
g.Close();

}