////////////////////////////////////////////////////////////////////
//Xavier Portell from original code of Kordas Kostas              //
//Code modified february 20th 2004                                //
//                                                                //
//Offline validation per each run using only Poisson distribution //
//with fluctuating mean (tuned)                                   //
//                                                                //
//High occupancy plots (1 GeV) to look for hot channels           //
//Low occupancy plots (0.5 GeV) to look for dead channels         //
//                                                                //
//Last modified in 6/1/04:                                        //
// - Introduced correction for special towers                     //
// - Output in HTML format                                        //
////////////////////////////////////////////////////////////////////
// Update 10/27/04                                                //
//This program only runs for looking into empty rings             //
////////////////////////////////////////////////////////////////////

#include <iostream.h>
#include <string>
#include <iomanip.h>
#include <fstream.h>
#include <cstdio>
#include <cstdlib>
#include <math.h>

#include <vector>

#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TFile.h>
#include <TPostScript.h>
#include <TGaxis.h>
#include <TLegend.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TChain.h>
#include <TGraphErrors.h>

#include "TGraph.h"
#include "TMath.h"
#include "TRandom.h"

using namespace std;

/* Including other functions */
int isEMChimney (int ieta, int iphi);
int isHADChimney (int ieta, int iphi);
int isDoubledPhiBins(int ieta);
int isKnownProblem(int ihist, int ieta, int iphi, int run_number);
float SigmaRatio(float observedCounts, float  meanCounts, float rmsOnMean);
/* Functions for EM/HAD special towers */
float isEMCorrFactor(int ieta);
float isHADCorrFactor(int ieta);
int isEMSpecialTow(int ieta,int iphi);
int isHADSpecialTow(int ieta,int iphi);
/* Function to know if put CCAL or PCAL */
int TypeOfBit(int emhad,int ieta); //Receives if EM or HAD and the eta.
                                   //Returns 1 for CCAL bit or 2 for PCAL bit

//void GaussCalc(float data_inside,float &Corrdatainside,float &CorrFactor);
int PoissonCalc(float mu,int Mode);
double logfactorial(int x, float mu);

/* Declarations for the GUI and main */
extern void InitGui();
VoidFuncPtr_t initfuncs[]={InitGui, 0};

TROOT root ("GUI", "GUI test environment", initfuncs);

/*********************************************************************/

int main(int argc, char *argv[])
{

  /* GUI Initialization */
  //TApplication TheApp("App", 0, 0);
  
  /* To work without plotting into the screen */
  TRint TheApp("App",0,0);
  gROOT->Reset();
  
  gStyle->SetPadBorderMode(0); 
  gStyle->SetOptFit(111111);
  //Int_t type=111;//portrait ps
  
  // *************************************************************************
  // SET the Thresholds and the Details
  // *************************************************************************
  
  // ATTENTION!!!  Root histogram bin numbers start from 1:  1, 2, 3, 4, .... 52
  // but the actual values on the Histo axes 
  // (in acordance to CDF-convention) start from 0:          0, 1, 2, 3, .... 51
  const int nbins_eta = 52; // # of eta_slices in the histograms (with the miniplug)
  const int min_eta_bin =  4;  // 8;
  const int max_eta_bin = 47;  // 43;
  
  // # of phi_slices in the histograms: 
  // even in the central, were there are 24 wedges, the number of phi bins is 48 
  // by giving two consecutive phi bins the same occupnacy (not half).
  const int nbins_phi = 48;
  
  //__________ sigma cuts for finding COLD/DEAD/HOT towers ____________
  // One-sided (as in Poisson tests), 6 sigma equivalent probability:
  Double_t prob_cut_flag = 1.00 * TMath::Power(10.0, -9); // 6 sigma
  // float prob_cut_flag = 2.85 * TMath::Power(10.0, -7); // 5 sigma 
  // float prob_cut_flag_dead = 3.15 * TMath::Power(10.0, -5); // 4 sigma
  

  int histoentries;  
  int min_rms = 2.0; // if the occupancy is very close to Zero, 
  // the Mean and RMS are close to zero, and then only the 0 occupancy
  // bins are left to participate in the recalculated Mean and RMS... 
  // and then we'll end up with a Mean = RMS = 0 ==> BAD!!!
  // so, instead of doing the Poisson thing right, set a min_rms to work 
  // with.
 
  //Minimum entries  
  int min_entries = 1;  // do not check histograms with entries  <= min_entries 

  //Value of the fluctuation of the mean
  //-------------------------------------
  //This value need to be tuned
  //For the moment take fluctuation times the RMS of the truncated distribution
  //without correcting for truncation
  //  float fluctuation = 1.0;

  //Number of excluded towers for each side (needs to be *2 when 48 towers)
  //The dead channels need to be cut also
  int NumTruncTowers=3;

  //_______ Define the arrays and histograms holding the information _______
  //------------------------------------------------------------------------
  float eta_phi_bin[nbins_eta][nbins_phi]; //Take values
  float eta_phi_bin_ordered[nbins_eta][nbins_phi]; //Ordered values
  float eta_slice_mean[nbins_eta]; //mean value for each eta slice
  float eta_slice_rms[nbins_eta];  //rms
  float eta_slice_mean_trunc[nbins_eta]; //mean value for truncated
  float eta_slice_rms_trunc[nbins_eta]; //rms corrected for truncated

  float eta_slice_mean_err[nbins_eta]; //error on the mean value for each eta slice
  float eta_slice_rms_err[nbins_eta]; //error on rms 
  float eta_slice_mean_trunc_err[nbins_eta];// error on mean for truncated
  float eta_slice_rms_trunc_err[nbins_eta]; //error on rms for corrected

  float MaxTowerOccup[nbins_eta]; //Maximum occupancy
  float MinTowerOccup[nbins_eta]; //Minimum occupancy

  //  int eta_phi_bin_ncold[nbins_eta][nbins_phi];
  // int eta_phi_bin_ndead[nbins_eta][nbins_phi];
  // int eta_phi_bin_nhot[nbins_eta][nbins_phi];
  // int eta_phi_bin_nfixable[nbins_eta][nbins_phi];

  //  int eta_phi_bin_ncold_runs[nbins_eta][nbins_phi];
  //  int eta_phi_bin_ndead_runs[nbins_eta][nbins_phi];
  //  int eta_phi_bin_nhot_runs[nbins_eta][nbins_phi];
  //  int eta_phi_bin_nfixable_runs[nbins_eta][nbins_phi];

  int em_chimney = 0; // false;
  int had_chimney = 0; // false;
  int ndead[nbins_eta]; //Number of dead histos per eta


  int ColdEta[52],HotEta[52];
  int ColdPhi[52],HotPhi[52];
  int ColdProb[52],HotProb[52];
  int ncold;
  int nhot;

  float CorrFactor; 

  int pmt0,pmt1; //To store the values of PMTs
    
  if(argc!=3) {
      cout<<"ERROR: Program need to be call with an argument:\n OfflineCalo <run_number> <filename>\n";
      return(0);
    }
  string run_number=argv[1];

  cout<<"ANALYSING RUN "<<run_number<<endl;
 
  //To calculate the semi-interval of a gaussian for a 87.5% of data
  //  float HalfInterval=90;// to have a dummy init
  //float Corrdatainside=0.875167; //percentage of data inside the gauss curve
  // float CorrFactor=1.32399;
  ////GaussCalc(0.875,Corrdatainside,CorrFactor);
  //cout<<Corrdatainside<<" and "<<CorrFactor<<endl;
  // exit(0);

  //Depending on the luminosity of the histo (high/low), looking for hot/dead towers
  int FindMode; 
  
  //Root staff
  gROOT->Reset();
  string filename(argv[2]);
  cout << " Reading from file " << filename << endl;
  //string filename="gmbs.root";
  string dir=":/DQMValidationModule/Tower/";
  string directory;
  string HistoName[]={"TowOccEM","TowOccHad","TowOccEMHigh","TowOccHadHigh"};

  int MaxOccup,MinOccup;
 
  //Opening root file
  char *num=new char[2];
  
  TH2F* histo[4]; //pointer to the read histos
  TH2F* PMT0histo[4]; //pointer to PMT0 histogram
  TH2F* PMT1histo[4]; //pointer to PMT1 histogram
 
  int bad_data=0; //Modify for more histos


  // ============================= STEP 1 ================================  
  // Loop over the histograms, get info into arrays, 
  // and calulate initial Mean, RMS for each eta_slice
  // with ALL THE BINS participating in the calculation
  // =====================================================================  
  //

  // Prepare outputs
  string aux="Calo.txt";
  
  ofstream listout;
  int flag=0; //Title not written

  ofstream listout2;
  string TagTow; //Variable to buffer the output of tagged towers
  int flag2=0; //Flag to not open file more than one time
  int CCAL=1; //Flag to test if CCAL has been tagged (1 means not tagged)
  int PCAL=1; //Flag to test if PCAL has been tagged (1 means not tagged)

  //Vector for emptyness check
  vector<int> EmptyRing_emcold;
  vector<int> EmptyRing_hadcold;
  vector<int> EmptyRing_emhot;
  vector<int> EmptyRing_hadhot;

  //Start loop
  for(int hloop=0;hloop<4;hloop++) {
    //Read from file 
    TFile f(filename.c_str());
    
    directory=filename;
    directory.append(dir);
    f.cd(directory.c_str());
    histo[hloop] = (TH2F*)gDirectory->Get(HistoName[hloop].c_str());
    histo[hloop]->SetDirectory(gROOT); 
    //Read PMT histograms
    string PMTdir=filename;
    PMTdir.append(":/TowerValidationModule/occurrence/");
    f.cd(PMTdir.c_str());
    if(hloop==0 || hloop==2){ //EM  
      PMT0histo[hloop] = (TH2F*)gDirectory->Get("Tower Occ. em pmt0");
      PMT0histo[hloop]->SetDirectory(gROOT); 
      PMT1histo[hloop] = (TH2F*)gDirectory->Get("Tower Occ. em pmt1");
      PMT1histo[hloop]->SetDirectory(gROOT); 
    }else{ //Had
      PMT0histo[hloop] = (TH2F*)gDirectory->Get("Tower Occ. ha pmt0");
      PMT0histo[hloop]->SetDirectory(gROOT); 
      PMT1histo[hloop] = (TH2F*)gDirectory->Get("Tower Occ. ha pmt1");
      PMT1histo[hloop]->SetDirectory(gROOT); 
      //cout<<"file "<<PMTdir.c_str()<<endl;
    }
    f.Close();
    
    //If looking for dead towers -> 1
    //If looking for hot towers -> 2
    FindMode=(hloop<2? 1 : 2); //Depending on low and high luminosity
    
    bad_data = 0; // false = "histogram is good enough to work with"
    if( !histo[hloop] ) {
      bad_data = 1; // Can not get the pointer to the histogram
    } else {
      histoentries = histo[hloop]->GetEntries();
      if ( histoentries <= min_entries ) {
        bad_data = 2; // too few entries to work with
      }
    }
    
    if ( bad_data != 0 ){
      if ( bad_data == 1 )
	cout <<"Histogram not accessible"<<endl;
      if ( bad_data == 2 )
	cout<<"Histogram with too few statistics"<< endl;
    } else {
      //Histogram is OK so go on...
    
      for (int ieta = min_eta_bin; ieta != (max_eta_bin + 1); ++ieta) {
	int doubled_phi_bins;
	ndead[ieta]=0;
	for (int iphi = 0; iphi != nbins_phi; ++iphi) { 
		  
	  // number of entries in each bin
	  // the info for eta_tower=0 in bin#1 (bin#'s start from 1 in Root)
	  // same is true for the phi_wedges...
	  eta_phi_bin[ieta][iphi]=histo[hloop]->GetBinContent(ieta+1 , iphi+1); 
	  
	  /******************************************************************/
	  // Deal with known problems first (e.g., chimney)
	  // to avoid reporting known problematic towers all the time...
	  //
	  if(hloop%2==0)
	    em_chimney  =  isEMChimney(ieta, iphi);
	  else
	    had_chimney = isHADChimney(ieta, iphi);
	    
	  //___ Flag the artificially Doubled number of Wedges in the Histogram
	  doubled_phi_bins = isDoubledPhiBins(ieta);
	  
	  // In the Eta bins that we have 48 different entries in Phi
	  // (i.e., for  8 <= eta <= 13   and   38 <= eta <= 43 ), 
	  // the occupancies have been doubled in the histograms
	  // (see CalorMods/src/TowerValidationModule.cc, 
	  //      lines 241, 296-297 and 335-336 in version 4.8.4)
	  // To correct for this, because I only care about the raw 
	  // occupancies here,  I divide these entries by 2 (two) to get 
	  // back the true occupnacy of these towers.
	  //
	  if ( !doubled_phi_bins ) eta_phi_bin[ieta][iphi] /= 2.0;
	      
	  // Count number of ZERO occupancy channels.
	  // If too many, then histogram is useless...
	  //
	  if ( !em_chimney && !had_chimney && 
	       eta_phi_bin[ieta][iphi] == 0 ) {
	    ndead[ieta]++;
	  }

	  // New at 5/28/04
	  // Look if tower is an special case to be corrected for hot/cold
	  // If so, the correction factor is returned; otherwise a 1 is returned
	  //  if(hloop%2==0)
	  //   CorrFactor=isEMSpecialTow(ieta, iphi);
	  // else
	  // CorrFactor=isHADSpecialTow(ieta, iphi);
	  //eta_phi_bin[ieta][iphi] /= CorrFactor;      

	  /*****************************************************************/
	  //_____________ Ordered array ______________
	  // Make a new array which is going to have the entries ordered in
	  // lowest to highest occupancy
	  // Make sure that in this ordered array the chimney bins get 
	  // the lowest occupancy
	  //
	  
	  eta_phi_bin_ordered[ieta][iphi] = eta_phi_bin[ieta][iphi];
	  if ( em_chimney || had_chimney )  {
	    eta_phi_bin_ordered[ieta][iphi] = -1;
	  }
	  //
	  // order the array as it gets filled: loop over all the elements up
	  // to now and order them from lowest occupancy to highest occupancy
	  //
	  for (int ia = 0; ia != iphi; ++ia) {
	    for ( int ib = ia+1; ib != iphi+1 ; ++ib ) {
	      if ( eta_phi_bin_ordered[ieta][ib] < 
		       eta_phi_bin_ordered[ieta][ia] ) {
		// swap the entries: lower occupancies go to lower bin numbers
		    int entry1 = eta_phi_bin_ordered[ieta][ib];
		    int entry2 = eta_phi_bin_ordered[ieta][ia];
		    eta_phi_bin_ordered[ieta][ia] = entry1;
		    eta_phi_bin_ordered[ieta][ib] = entry2;
	      }
	    }
	  } //_________ end of ordering the array for now __________
	} //_________ PHI for-loop ____________
	

	/***********************************************************/
	//_______ Statistics for this eta_slice: Mean and RMS
	//For cut and not cut towers (highest and lowest occupancy + ndead)
	//To make the plots (not truncated data), do not consider chimneys

	float nwedges = 0.0; //counter of meaningful wedges for each eta-slice
	float sum_eta = 0.0;  // sum of bin contents in each eta-slice 
	float sum2_eta = 0.0;  // sum of squares      "  "        "
	float nwedges_trunc = 0.0; //the same as before but for truncated data
	float sum_eta_trunc = 0.0;  // "
	float sum2_eta_trunc = 0.0;  // "
	int MinTower, MaxTower; //To get extreme towers
	int chimney=0;

	if(!doubled_phi_bins) {
	  MinTower=NumTruncTowers*2-1; //Case 48 towers we neglect double
	  MaxTower=nbins_phi-NumTruncTowers*2;
	  for(int iphi=0;iphi<nbins_phi;iphi++) {
	    if(eta_phi_bin_ordered[ieta][iphi]>=0) { //not tagged as bad
	      chimney=(hloop%2==0? isEMChimney(ieta,iphi) : isHADChimney(ieta,iphi));
	      if(!chimney){
		sum_eta=sum_eta+eta_phi_bin_ordered[ieta][iphi];
		sum2_eta = sum2_eta +
		  eta_phi_bin_ordered[ieta][iphi]*eta_phi_bin_ordered[ieta][iphi];
		++nwedges;
	      }
	      if(iphi>=MinTower && iphi<MaxTower 
		 && eta_phi_bin_ordered[ieta][iphi]>0) {//not dead tower
		sum_eta_trunc=sum_eta_trunc+eta_phi_bin_ordered[ieta][iphi];
		sum2_eta_trunc = sum2_eta_trunc +
		  eta_phi_bin_ordered[ieta][iphi]*eta_phi_bin_ordered[ieta][iphi];
		++nwedges_trunc;
	      }
	    }
	  } //end for loop
	} else {
	  MinTower=NumTruncTowers-1; //Case 24 towers
	  MaxTower=nbins_phi-NumTruncTowers;
	  for(int iphi=0;iphi<nbins_phi;iphi++) {
	    if ( iphi%2 == 0 ) {	
	      if(eta_phi_bin_ordered[ieta][iphi]>=0) { //not tagged as bad
		chimney=(hloop%2==0? isEMChimney(ieta,iphi) : isHADChimney(ieta,iphi));
		if(!chimney){
		  sum_eta=sum_eta+eta_phi_bin_ordered[ieta][iphi];
		  sum2_eta = sum2_eta +
		    eta_phi_bin_ordered[ieta][iphi]*eta_phi_bin_ordered[ieta][iphi];
		  ++nwedges;
		}  
		if(iphi>=MinTower && iphi<MaxTower
		   && eta_phi_bin_ordered[ieta][iphi]>0) {//not dead tower
		  sum_eta_trunc=sum_eta_trunc+eta_phi_bin_ordered[ieta][iphi];
		  sum2_eta_trunc = sum2_eta_trunc +
		    eta_phi_bin_ordered[ieta][iphi]*eta_phi_bin_ordered[ieta][iphi];
		  ++nwedges_trunc;
		}
	      }
	    }
	  } //end for loop
	}

	///////////////////////////////////////////////////////////////////
	//Looking for empty rings (Modified 10/27/04
	//
	//Check if the ring is empty (only if looking for cold towers!!!)
	if(FindMode==1 && (sum_eta==0 || sum_eta_trunc==0)){
	  //mark CCAL/PCAL
	  int bit=TypeOfBit(hloop,ieta); //to know if CCAL or PCAL
	  if(bit==1)
	    CCAL=0;
	  else
	    PCAL=0;
	  //Check had,em
	    if(hloop%2==0){ //EM
	      EmptyRing_emcold.push_back(ieta);
	    }else{ //HAD
	      EmptyRing_hadcold.push_back(ieta);
	    }

	}//End if empty ring

	///////////////////////////////////////////////////////////////////
	//For raw data
	if ( nwedges == 0.0 ) {
	  nwedges = 1.0;
	  sum_eta = 0.5;
	  sum2_eta = sum_eta * sum_eta;
	}
	
	float mean=sum_eta/nwedges;
	float rms2 = (sum2_eta + nwedges * mean * mean - 2.0 * mean * sum_eta) / nwedges;
	eta_slice_mean[ieta] = mean; 

	//cout<<"mean: "<<eta_slice_mean[ieta]<<" occup0: "<<eta_phi_bin[ieta][0]<<" occup44 "<<eta_phi_bin[ieta][44]<<" ieta: "<<ieta<<endl;

	if ( sqrt(rms2) < min_rms ) {
	  eta_slice_rms[ieta] = (sqrt(mean)>min_rms) ? sqrt(mean) : min_rms ;
	} else {
	  eta_slice_rms[ieta] = sqrt(rms2);
	}
	
	//Errors on mean and RMS
	eta_slice_mean_err[ieta] = eta_slice_rms[ieta]/sqrt(nwedges-1.);
	eta_slice_rms_err[ieta] = eta_slice_rms[ieta]/sqrt(2.*(nwedges-1.));

	////////////////////////////////////////////////////////////
	//For truncated data
	if ( nwedges_trunc == 0.0 ) {
	  nwedges_trunc = 1.0;
	  sum_eta_trunc = 0.5;
	  sum2_eta_trunc = sum_eta_trunc * sum_eta_trunc;;
	}
	
	mean=sum_eta_trunc/nwedges_trunc;
	rms2 = (sum2_eta_trunc + nwedges_trunc * mean * mean - 
		2.0 * mean * sum_eta_trunc) / nwedges_trunc;
	eta_slice_mean_trunc[ieta] = mean; 
	// cout<<"mean: "<<mean<<" nwedges_trunc: "<<nwedges_trunc<<endl;
	
	if ( sqrt(rms2) < min_rms ) {
	  eta_slice_rms_trunc[ieta] = (sqrt(mean)>min_rms) ? sqrt(mean) : min_rms ;
	} else {
	  eta_slice_rms_trunc[ieta] = sqrt(rms2);
	}
	//eta_slice_rms_trunc[ieta] = eta_slice_rms_trunc[ieta]*CorrFactor; //corrected RMS
	
	//Errors on mean and RMS
	eta_slice_mean_trunc_err[ieta] = eta_slice_rms_trunc[ieta]/sqrt(nwedges_trunc-1.);
	eta_slice_rms_trunc_err[ieta] = eta_slice_rms_trunc[ieta]/sqrt(2.*(nwedges_trunc-1.));
	
	
	/***********************************************************/
	//_______ Looking for the highest and lowest occupancy tower
	//
	MaxTowerOccup[ieta]=eta_phi_bin_ordered[ieta][nbins_phi-1];
	for(int i=0;i<nbins_phi;i++) {
	  if(eta_phi_bin_ordered[ieta][i]>=0) { //not tagged as bad
	    MinTowerOccup[ieta]=eta_phi_bin_ordered[ieta][i];
	    break;
	  }
	}
      } // eta loop
    } // end-if reasonable histogram
    
    //cout<<"Loop for calculation finished.\n";
        
    ///////////////////////////////////////////////////////////////////////////////////////

    // ============================= STEP 2 ================================  
    // Look for dead or hot channels (depending on the histo)
    // Calculate the deviation from the Poisson mean fluctuated 
    // =====================================================================  
    //
    ncold=0;
    nhot=0;
    int chimney;
    if(FindMode==1) { //Looking for cold towers
      float mu2,occup_mu2,occup_mu2_2,miscalib;
      for (int ieta = min_eta_bin; ieta < (max_eta_bin + 1); ++ieta) {
	int doubled_phi_bins = isDoubledPhiBins(ieta);

	//Calculate corrections and occupancy limits
	if(hloop%2==0){ //EM
	  miscalib= 0.04; //EM (4%)
	  mu2=eta_slice_mean_trunc[ieta]-3.*miscalib*eta_slice_mean_trunc[ieta];
	  occup_mu2=PoissonCalc(mu2,FindMode);

	  //Correct the mean value for the case of special towers
	  CorrFactor=isEMCorrFactor(ieta);
	  if(CorrFactor != 1.)
	    occup_mu2_2=PoissonCalc(CorrFactor*mu2,FindMode);
	} else { //HAD
	  miscalib= 0.03; //HAD (3%)
	  mu2=eta_slice_mean_trunc[ieta]-3.*miscalib*eta_slice_mean_trunc[ieta];
	  occup_mu2=PoissonCalc(mu2,FindMode);
	  
	  //Correct the mean value for the case of special towers
	  CorrFactor=isHADCorrFactor(ieta);
	  if(CorrFactor != 1.)
	    occup_mu2_2=PoissonCalc(CorrFactor*mu2,FindMode);
	}
	
	//Start loop for phi to look for cold towers
	for(int iphi=0;iphi<nbins_phi;iphi++) {
	  if(!doubled_phi_bins) {
	    if(hloop%2==0) { //EM case
	      if(!isEMChimney(ieta,iphi)) {
		if(!isEMSpecialTow(ieta,iphi)) { //it is a normal EM tower
		  if(eta_phi_bin[ieta][iphi] < occup_mu2) {
		    ColdEta[ncold]=ieta;
		    ColdPhi[ncold]=iphi;
		    ColdProb[ncold]=occup_mu2;
		    ++ncold;
		  } 
		} else {// it is a EM special tower
		  if(eta_phi_bin[ieta][iphi] < occup_mu2_2) {
		    ColdEta[ncold]=ieta;
		    ColdPhi[ncold]=iphi;
		    ColdProb[ncold]=occup_mu2_2;
		    ++ncold;
		  }
		}
	      }
	    } else { //HAD case
	      if(!isHADChimney(ieta,iphi)) {
		if(!isHADSpecialTow(ieta,iphi)) { //it is a normal HAD tower
		  if(eta_phi_bin[ieta][iphi] < occup_mu2) {
		    ColdEta[ncold]=ieta;
		    ColdPhi[ncold]=iphi;
		    ColdProb[ncold]=occup_mu2;
		    ++ncold;
		  } 
		} else { //it is a special HAD tower
		  if(eta_phi_bin[ieta][iphi] < occup_mu2_2) {
		    ColdEta[ncold]=ieta;
		    ColdPhi[ncold]=iphi;
		    ColdProb[ncold]=occup_mu2_2;
		    ++ncold;
		  }
		}
	      }
	    } //End EM/HAD
	    
	  } else { //phi bins doubled so only work when iphi%2==0
	    if ( iphi%2 == 0 ) {
	      if(hloop%2==0) { //EM case
		if(!isEMChimney(ieta,iphi)) {
		  if(!isEMSpecialTow(ieta,iphi)) { //it is a normal EM tower
		    if(eta_phi_bin[ieta][iphi] < occup_mu2) {
		      ColdEta[ncold]=ieta;
		      ColdPhi[ncold]=iphi;
		      ColdProb[ncold]=occup_mu2;
		      ++ncold;
		    } 
		  } else {// it is a EM special tower
		    if(eta_phi_bin[ieta][iphi] < occup_mu2_2) {
		      ColdEta[ncold]=ieta;
		      ColdPhi[ncold]=iphi;
		      ColdProb[ncold]=occup_mu2_2;
		      ++ncold;
		    }
		  }
		}
	      } else { //HAD case
		if(!isHADChimney(ieta,iphi)) {
		  if(!isHADSpecialTow(ieta,iphi)) { //it is a normal HAD tower
		    if(eta_phi_bin[ieta][iphi] < occup_mu2) {
		      ColdEta[ncold]=ieta;
		      ColdPhi[ncold]=iphi;
		      ColdProb[ncold]=occup_mu2;
		      ++ncold;
		    } 
		  } else { //it is a special HAD tower
		    if(eta_phi_bin[ieta][iphi] < occup_mu2_2) {
		      ColdEta[ncold]=ieta;
		      ColdPhi[ncold]=iphi;
		      ColdProb[ncold]=occup_mu2_2;
		    ++ncold;
		    }
		  }
		}
	      } //End EM/HAD
	    }
	  } //End simple/double bins
	}//End loop over phi
      } //End loop over eta
    } else { //Now hot towers
      float mu2,occup_mu2,occup_mu2_2,miscalib;
      for (int ieta = min_eta_bin; ieta < (max_eta_bin + 1); ++ieta) {
	int doubled_phi_bins = isDoubledPhiBins(ieta);
	
	//Calculate corrections and occupancy limits
	if(hloop%2==0){ //EM
	  miscalib= 0.04; //EM (4%)
	  mu2=eta_slice_mean_trunc[ieta]+3.*miscalib*eta_slice_mean_trunc[ieta];
	  occup_mu2=PoissonCalc(mu2,FindMode);
	  
	  //Correct the mean value for the case of special towers
	  CorrFactor=isEMCorrFactor(ieta);
	  if(CorrFactor != 1.)
	    occup_mu2_2=PoissonCalc(CorrFactor*mu2,FindMode);
	} else { //HAD
	  miscalib= 0.03; //HAD (3%)
	  mu2=eta_slice_mean_trunc[ieta]+3.*miscalib*eta_slice_mean_trunc[ieta];
	  occup_mu2=PoissonCalc(mu2,FindMode);
          //cout<<"occup_mu2= "<<occup_mu2<<" mean= "<<eta_slice_mean_trunc[ieta]<<" mu2= "<<mu2<<" eta= "<<ieta<<endl;

	  //Correct the mean value for the case of special towers
	  CorrFactor=isHADCorrFactor(ieta);
	  if(CorrFactor != 1.)
	    occup_mu2_2=PoissonCalc(CorrFactor*mu2,FindMode);
	}
	
	//Start loop for phi to look for hot towers
	for(int iphi=0;iphi<nbins_phi;iphi++) {
	  if(!doubled_phi_bins) {
	    if(hloop%2==0) { //EM case
	      if(!isEMChimney(ieta,iphi)) {
		if(!isEMSpecialTow(ieta,iphi)) { //it is a normal EM tower
		  if(eta_phi_bin[ieta][iphi] > occup_mu2) {
		    HotEta[nhot]=ieta;
		    HotPhi[nhot]=iphi;
		    HotProb[nhot]=occup_mu2;
		    ++nhot;
		  } 
		} else {// it is a EM special tower
		  if(eta_phi_bin[ieta][iphi] > occup_mu2_2) {
		    HotEta[nhot]=ieta;
		    HotPhi[nhot]=iphi;
		    HotProb[nhot]=occup_mu2_2;
		    ++nhot;
		  }
		}
	      }
	    } else { //HAD case
	      if(!isHADChimney(ieta,iphi)) {
		if(!isHADSpecialTow(ieta,iphi)) { //it is a normal HAD tower
		  if(eta_phi_bin[ieta][iphi] > occup_mu2) {
		    HotEta[nhot]=ieta;
		    HotPhi[nhot]=iphi;
		    HotProb[nhot]=occup_mu2;
		    ++nhot;
		  } 
		} else { //it is a special HAD tower
		  if(eta_phi_bin[ieta][iphi] > occup_mu2_2) {
		    HotEta[nhot]=ieta;
		    HotPhi[nhot]=iphi;
		    HotProb[nhot]=occup_mu2_2;
		    ++nhot;
		  }
		}
	      }
	    } //End EM/HAD
	    
	  } else { //phi bins doubled so only work when iphi%2==0
	    if ( iphi%2 == 0 ) {
	      if(hloop%2==0) { //EM case
		if(!isEMChimney(ieta,iphi)) {
		  if(!isEMSpecialTow(ieta,iphi)) { //it is a normal EM tower
		    if(eta_phi_bin[ieta][iphi] > occup_mu2) {
		      HotEta[nhot]=ieta;
		      HotPhi[nhot]=iphi;
		      HotProb[nhot]=occup_mu2;
		      ++nhot;
		    } 
		  } else {// it is a EM special tower
		    if(eta_phi_bin[ieta][iphi] > occup_mu2_2) {
		      HotEta[nhot]=ieta;
		      HotPhi[nhot]=iphi;
		      HotProb[nhot]=occup_mu2_2;
		      ++nhot;
		    }
		  }
		}
	      } else { //HAD case
		if(!isHADChimney(ieta,iphi)) {
		  if(!isHADSpecialTow(ieta,iphi)) { //it is a normal HAD tower
		    if(eta_phi_bin[ieta][iphi] > occup_mu2) {
		      HotEta[nhot]=ieta;
		      HotPhi[nhot]=iphi;
		      HotProb[nhot]=occup_mu2;
		      ++nhot;
		    } 
		  } else { //it is a special HAD tower
		    if(eta_phi_bin[ieta][iphi] > occup_mu2_2) {
		      HotEta[nhot]=ieta;
		      HotPhi[nhot]=iphi;
		      HotProb[nhot]=occup_mu2_2;
		      ++nhot;
		    }
		  }
		}
	      } //End EM/HAD
	    }
	  } //End simple/double bins
	}//End loop over phi
      } //End loop over eta
    } //End hot towers
    
    
    // ============================= STEP 3 ================================  
    // Writing the outputs with the PMT values
    // =====================================================================  
    //
    int big,skip;
    string Format;
    char *ToString=new char[8];

    //Write the title if it is necessary and has not been done
    if(flag==0 && (nhot+ncold>0 || CCAL==0 || PCAL==0)) {
      listout.open(aux.c_str());
      listout<<"-------------"<<endl;
      listout<<"Run: "<<run_number.c_str()<<endl;
      listout<<"-------------"<<endl;
      flag=1;
    }

    if(flag2==0) { //Only done at the beginning
      TagTow="<td><font size=-1>"; 
      flag2=1;
     // CCAL=PCAL=1; //Not tagged
    }

    //Prepare for reading PMT values from plots
    if( !PMT0histo[hloop] ) { // Can not get the pointer to the histogram
      TagTow.append("Histogram PMT0 ");
      sprintf(ToString,"%d",hloop);
      TagTow.append(ToString);
      TagTow.append(" is not accessible!!!</td>");
      CCAL=PCAL=0;
      exit(1);
    }
    if( !PMT1histo[hloop] ) {
      TagTow.append("Histogram PMT1 ");
      sprintf(ToString,"%d",hloop);
      TagTow.append(ToString);
      TagTow.append(" is not accessible!!!</td>");
      CCAL=PCAL=0;
      exit(1);
    }
    
    int meanaux; //temporal variable
    switch(hloop) {
    case 0:
      //empty rings
      if(EmptyRing_emcold.size()!=0){
	TagTow.append("EMCOLD EMPTY RINGS: ");
	listout<<"\nEMPTY RINGS: ";
	char *ToString=new char[7];
	for(vector<int>::size_type pos=0;pos!=EmptyRing_emcold.size();pos++){
	  listout<<EmptyRing_emcold[pos]<<" ";
	  sprintf(ToString,"%d",EmptyRing_emcold[pos]);
	  TagTow.append(ToString);
	  TagTow.append(" ");
	}
	listout<<"\n";      
	TagTow.append("<br>");
      }
      //towers
      if(ncold!=0) {
	//One file
	listout<<"\nEM cold towers (occ. threshold of 0.5 GeV)"<<endl;
	listout<<"------------------------------------------"<<endl;
	listout<<endl;
	listout<<setiosflags(ios::right);
	listout.precision(0);

	listout<<setw(7)<<"Eta/Phi"<<setw(6)<<"mu"<<setw(8)<<"mu-12%"<<setw(11)<<"sqrt(mu*)"<<setw(6)<<"RMS"<<setw(8)<<"Occup."<<setw(9)<<"PoissLim"<<endl;
	
	//Second output file
	TagTow.append("EMCOLD ");

	for(register int i=0;i<ncold;i++) {
	  // FindPMTs(PMT0histo[hloop],PMT1histo[hloop],ColdEta[i],ColdPhi[i],pmt0,pmt1);
	  // cout<<"ColdEta= "<<ColdEta[i]<<" ColdPhi= "<<ColdPhi[i]<<endl;
	  //cout<<PMT0histo[hloop]->GetEntries()<<endl;
	  pmt0=PMT0histo[hloop]->GetBinContent(ColdEta[i]+1 , ColdPhi[i]+1);
	  // cout<<"pmt0 "<<pmt0<<endl;
	  pmt1=PMT1histo[hloop]->GetBinContent(ColdEta[i]+1 , ColdPhi[i]+1);
	  // cout<<"pmt1 "<<pmt1<<endl;
	  // cout<<"here4"<<endl;
	  meanaux=int(eta_slice_mean_trunc[ColdEta[i]]-0.12*eta_slice_mean_trunc[ColdEta[i]]);
	  listout<<setw(4)<<ColdEta[i]<<"/"<<setw(2)<<ColdPhi[i]
		 <<setw(6)<<(int)eta_slice_mean_trunc[ColdEta[i]]
		 <<setw(8)<<meanaux
		 <<setw(11)<<(int)sqrt(meanaux)
		 <<setw(6)<<(int)eta_slice_rms[ColdEta[i]]
		 <<setw(8)<<(int)eta_phi_bin[ColdEta[i]][ColdPhi[i]]
		 <<setw(9)<<ColdProb[i]
		 <<endl;

	  sprintf(ToString,"%d",ColdEta[i]); 
	  Format=ToString; Format.append("/");
	  sprintf(ToString,"%d",ColdPhi[i]);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt0);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt1);
	  Format.append(ToString);

	  if(i!=0 && i%4==0){
	    TagTow.append("<br>.................... ");
	    TagTow.append(Format);
	  }else{
	    TagTow.append("  ");
	    TagTow.append(Format);
	  }
	  
	  //Check CCAL and PCAL bits
	  int bit=TypeOfBit(hloop,ColdEta[i]);
	  if(bit==1) //CCAL
	    CCAL=0;
	  else //PCAL
	    PCAL=0;
	}
	TagTow.append("<br>");
      }
      break;
    case 1:
      //Empty rings 
      if(EmptyRing_hadcold.size()!=0){
	TagTow.append("HADCOLD EMPTY RINGS: ");
	listout<<"\nEMPTY RINGS: ";
	char *ToString=new char[7];
	for(vector<int>::size_type pos=0;pos!=EmptyRing_hadcold.size();pos++){
	  listout<<EmptyRing_hadcold[pos]<<" ";
	  sprintf(ToString,"%d",EmptyRing_hadcold[pos]);
	  TagTow.append(ToString);
	  TagTow.append(" ");
	}
	listout<<"\n";      
	TagTow.append("<br>");
      }
      //Towers
      if(ncold!=0) {
	listout<<"\nHAD cold towers (occ. threshold of 0.5 GeV)"<<endl;
	listout<<"-------------------------------------------"<<endl;
	listout<<endl;
	listout<<setiosflags(ios::right);
	listout.precision(0);
	
	listout<<setw(7)<<"Eta/Phi"<<setw(6)<<"mu"<<setw(8)<<"mu-9%"<<setw(11)<<"sqrt(mu*)"<<setw(6)<<"RMS"<<setw(8)<<"Occup."<<setw(9)<<"PoissLim"<<endl;
	
	//Second output file
	TagTow.append("HADCOLD ");

	for(register int i=0;i<ncold;i++) {
	  pmt0=PMT0histo[hloop]->GetBinContent(ColdEta[i]+1 , ColdPhi[i]+1);
	  pmt1=PMT1histo[hloop]->GetBinContent(ColdEta[i]+1 , ColdPhi[i]+1);
	  meanaux=int(eta_slice_mean_trunc[ColdEta[i]]-0.09*eta_slice_mean_trunc[ColdEta[i]]);

	  listout<<setw(4)<<ColdEta[i]<<"/"<<setw(2)<<ColdPhi[i]
		 <<setw(6)<<(int)eta_slice_mean_trunc[ColdEta[i]]
		 <<setw(8)<<meanaux
		 <<setw(11)<<(int)sqrt(meanaux)
		 <<setw(6)<<(int)eta_slice_rms[ColdEta[i]]
		 <<setw(8)<<(int)eta_phi_bin[ColdEta[i]][ColdPhi[i]]
		 <<setw(9)<<ColdProb[i]
		 <<endl;

	  sprintf(ToString,"%d",ColdEta[i]); 
	  Format=ToString; Format.append("/");
	  sprintf(ToString,"%d",ColdPhi[i]);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt0);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt1);
	  Format.append(ToString);
	  
	  if(i!=0 && i%4==0){
	    TagTow.append("<br>.................... ");
	    TagTow.append(Format);
	  }else{
	    TagTow.append("  ");
	    TagTow.append(Format);
	  }

	  //Check CCAL and PCAL bits
	  int bit=TypeOfBit(hloop,ColdEta[i]);
	  if(bit==1) //CCAL
	    CCAL=0;
	  else //PCAL
	    PCAL=0;
	}
	TagTow.append("<br>");
      }
      break;
    case 2:
      if(nhot!=0) { 
	listout<<"\nEM hot towers (occ. threshold of 1.0 GeV)"<<endl;
	listout<<"-----------------------------------------"<<endl;
	listout<<endl;
	listout<<setiosflags(ios::right);
	listout.precision(0);
	
	listout<<setw(7)<<"Eta/Phi"<<setw(6)<<"mu"<<setw(8)<<"mu+12%"<<setw(11)<<"sqrt(mu*)"<<setw(6)<<"RMS"<<setw(8)<<"Occup."<<setw(9)<<"PoissLim"<<endl;
	
	//Second output file
	TagTow.append("EMHOT ");

	for(register int i=0;i<nhot;i++) {
	  pmt0=PMT0histo[hloop]->GetBinContent(HotEta[i]+1 , HotPhi[i]+1);
	  pmt1=PMT1histo[hloop]->GetBinContent(HotEta[i]+1 , HotPhi[i]+1);
	  meanaux=int(eta_slice_mean_trunc[HotEta[i]]+0.12*eta_slice_mean_trunc[HotEta[i]]);

	  listout<<setw(4)<<HotEta[i]<<"/"<<setw(2)<<HotPhi[i]
		 <<setw(6)<<(int)eta_slice_mean_trunc[HotEta[i]]
		 <<setw(8)<<meanaux
		 <<setw(11)<<(int)sqrt(meanaux)
		 <<setw(6)<<(int)eta_slice_rms[HotEta[i]]
		 <<setw(8)<<(int)eta_phi_bin[HotEta[i]][HotPhi[i]]
		 <<setw(9)<<HotProb[i]
		 <<endl;

	  sprintf(ToString,"%d",HotEta[i]); 
	  Format=ToString; Format.append("/");
	  sprintf(ToString,"%d",HotPhi[i]);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt0);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt1);
	  Format.append(ToString);

	 if(i!=0 && i%4==0){
	    TagTow.append("<br>.................... ");
	    TagTow.append(Format);
	  }else{
	    TagTow.append("  ");
	    TagTow.append(Format);
	  }

	 //Check CCAL and PCAL bits
	 int bit=TypeOfBit(hloop,HotEta[i]);
	 if(bit==1) //CCAL
	   CCAL=0;
	 else //PCAL
	   PCAL=0;
	}
	TagTow.append("<br>");
      }
      break;
    case 3:
      if(nhot!=0) { 
	listout<<"\nHAD hot towers (occ. threshold of 1.0 GeV)"<<endl;
	listout<<"-------------------------------------------"<<endl;
	listout<<endl;
	listout<<setiosflags(ios::right);
	listout.precision(0);
	
	listout<<setw(7)<<"Eta/Phi"<<setw(6)<<"mu"<<setw(8)<<"mu+9%"<<setw(11)<<"sqrt(mu*)"<<setw(6)<<"RMS"<<setw(8)<<"Occup."<<setw(9)<<"PoissLim"<<endl;
	
	//Second output file
	TagTow.append("HADHOT ");
	
	for(register int i=0;i<nhot;i++) {
	  pmt0=PMT0histo[hloop]->GetBinContent(HotEta[i]+1 , HotPhi[i]+1);
	  pmt1=PMT1histo[hloop]->GetBinContent(HotEta[i]+1 , HotPhi[i]+1);
	  meanaux=int(eta_slice_mean_trunc[HotEta[i]]+0.09*eta_slice_mean_trunc[HotEta[i]]);
	  listout<<setw(4)<<HotEta[i]<<"/"<<setw(2)<<HotPhi[i]
		 <<setw(6)<<(int)eta_slice_mean_trunc[HotEta[i]]
		 <<setw(8)<<meanaux
		 <<setw(11)<<(int)sqrt(meanaux)
		 <<setw(6)<<(int)eta_slice_rms[HotEta[i]]
		 <<setw(8)<<(int)eta_phi_bin[HotEta[i]][HotPhi[i]]
		 <<setw(9)<<HotProb[i]
		 <<endl;

	  sprintf(ToString,"%d",HotEta[i]); 
	  Format=ToString; Format.append("/");
	  sprintf(ToString,"%d",HotPhi[i]);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt0);
	  Format.append(ToString); Format.append("/");
	  sprintf(ToString,"%d",pmt1);
	  Format.append(ToString);
	  
	  if(i!=0 && i%4==0){
	    TagTow.append("<br>.................... ");
	    TagTow.append(Format);
	  }else{
	    TagTow.append("  ");
	    TagTow.append(Format);
	  }

	  //Check CCAL and PCAL bits
	  int bit=TypeOfBit(hloop,HotEta[i]);
	  if(bit==1) //CCAL
	    CCAL=0;
	  else //PCAL
	    PCAL=0;
	}
	TagTow.append("<br>");
      }
      break;
    default: 
      cout<<"ERROR!"<<endl;
      exit(0);
    }    
  }//End hloop loop
  
  //Write the output to a file
  listout2.open("Calo.out");

  if(CCAL==0 && PCAL==0)
    listout2<<"<tr><td><font color=ORANGE> "<<run_number<<"</font> </td><td> -1 </td><td> -1 </td><td> -1 </td>"<<TagTow.c_str()<<"</font></td>";
  if(CCAL==0 && PCAL==1) //tagged only CCAL
    listout2<<"<tr><td><font color=ORANGE> "<<run_number<<"</font> </td><td> -1 </td><td> 1 </td><td> -1 </td>"<<TagTow.c_str()<<"</font></td>";
  if(CCAL==1 && PCAL==0) //tagged only PCAL
    listout2<<"<tr><td><font color=ORANGE> "<<run_number<<"</font> </td><td> 1 </td><td> -1 </td><td> -1 </td>"<<TagTow.c_str()<<"</font></td>";
  if(CCAL==1 && PCAL==1)
    listout2<<"<tr><td><font color=LIME> "<<run_number<<" </font></td><td> 1 </td><td> 1 </td><td> 1 </td><td> --- </td>";

  if(nhot+ncold>0 || CCAL==0 || PCAL==0)
    listout2<<"<td> <a href=\"Calo_"<<run_number.c_str()<<".txt\">Details</a></td>";
  listout2<< "</tr>"<< std::endl;

  listout2.close();
}




////////////////////////////////////////////
////////////////////////////////////////////
//Functions (TowerMonitoring_utilities.C)
////////////////////////////////////////////
////////////////////////////////////////////
int isEMChimney (int ieta, int iphi) {

  int wedgeNo = (int) iphi/2;

  //
  //_______ EM chimmey: BOTH PMT's missing at eta = 34 and 35, phi = 5
  // The one at eta = 33 is recording double, so it's hot by design
  // and I flag him as "Chimney tower" as well, because his
  // occupancy is not typical.
	  
  int em_chimney = 0; // false;
  if ((ieta == 33 || ieta == 34 || ieta == 35) && 
       (wedgeNo == 5)) {
    em_chimney = 1; // true;
  }

  return em_chimney;

} // isEMChimney()



int isHADChimney (int ieta, int iphi) {

  int wedgeNo = (int) iphi/2;
  
  //_______ HAD chimney
  // (in reality, only PMT1 is missing, at eta = 30,31,32 && phi = 5
  // and actualy there is a PMT at eta=32, phi=5 but it records little)

  int had_chimney = 0; // false;
  if ((ieta == 30 || ieta == 31 || ieta == 32) && 
       (wedgeNo == 5)) {
    had_chimney = 1; // true;
  }

  return had_chimney;

} // isHADChimney()

int isDoubledPhiBins(int ieta) {
 
  //_____Flag the artificially double number of wedges in the histogram
  // e.g., in the CEM we have 24 wedges in the detector, 
  //       but the histogram has 48 phi-bins 
 
  int doubled_phi_bins = 0; // false;

  if ( (  4 <= ieta && ieta <=  7 ) || 
       ( 14 <= ieta && ieta <= 37 ) ||
       ( 44 <= ieta && ieta <= 47 )
       ) doubled_phi_bins = 1; // true;

  return doubled_phi_bins;

} // isDoubledPhiBins()


float SigmaRatio(float observedCounts, float  meanCounts, float rmsOnMean) {

  // observedCountrs: number of counts observed
  // meanCounts:      average number of counts expected
  // rmsOnMean:       uncertainty on the meanCounts: rmsCounts/sqrt(N)


  float ratio = 1.0;
  if ( meanCounts != 0.0 ) { 
    ratio = observedCounts /  meanCounts;
  }

  //
  // the uncertainty on the ratio takes into account the 
  // uncertainties of both values forming the ratio
  //
  float rel1 = 0.0;
  if ( observedCounts != 0.0 ) {
    rel1 = sqrt(observedCounts) / observedCounts;
  }
  float rel2 = 0.0;
  if ( meanCounts != 0.0 ) {
    rel2 = rmsOnMean / meanCounts;
  }
  float temp2 = ratio * ratio * (rel1 * rel1 + rel2 * rel2);

  float sigmaRatio = sqrt(temp2); 

  return sigmaRatio;

} // SigmaRatio()


int PoissonCalc(float mu,int Mode)
{
  int k=0;
  double sum=0.;
  if(mu>0)
    double logmu=TMath::Log(mu);
  else return 0;

  if(mu>500) //For long calculations, start with a shortcut (not enough resolution)
    k=int(mu-10*sqrt(mu));
  if(Mode==1) { //Look for dead towers
    while(sum < TMath::Power(10,-9)) {
      sum=sum + TMath::Exp(logfactorial(k,mu)-mu);
      ++k;
      //  cout<<k<<" "<<mu<<" "<<sum<<endl;
      // exit(0);
    }
  } 
  else { //Look for dead towers
    while(sum < 1.-TMath::Power(10,-9)) {
      sum=sum + TMath::Exp(logfactorial(k,mu)-mu);
      ++k;
    }
  }
  //  cout<<"k= "<<k<<endl;
  return k-1;
}


double logfactorial(int x, float mu) {
  double sum=0.;
  if(x==0)
    return sum;
  for(int i=1;i<=x;i++)
    sum=sum+TMath::Log(mu)-TMath::Log(i);
  //  cout<<sum<<endl;
  return sum;
}


// void FindPMTs(TH2F *h0,TH2F *h1,int eta,int phi,float& pmt0,float& pmt1);
// {
//   // number of entries in each bin
//   // the info for eta_tower=0 in bin#1 (bin#'s start from 1 in Root)
//   // same is true for the phi_wedges...
//   pmt0=histo[hloop]->GetBinContent(eta+1 , phi+1); 
  
//   /******************************************************************/
//   // Deal with known problems first (e.g., chimney)
//   // to avoid reporting known problematic towers all the time...
//   //
//   em_chimney  =  isEMChimney(ieta, iphi);
//   had_chimney = isHADChimney(ieta, iphi);
  
//   //___ Flag the artificially Doubled number of Wedges in the Histogram
//   doubled_phi_bins = isDoubledPhiBins(ieta);
  
// }

float isEMCorrFactor(int ieta)
{
  float CorrFact=1.;
  if (ieta == 33)
    CorrFact=1.71;
  return CorrFact;
}


float isHADCorrFactor(int ieta)
{
  float CorrFact=1.;
  if (ieta == 33)
    CorrFact=0.275;
  if (ieta == 34)
    CorrFact=1.25;
  if (ieta == 35)
    CorrFact=1.54;
  return CorrFact;
}

int isEMSpecialTow(int ieta,int iphi)
{
 int wedgeNo= (int)iphi/2;
 if (ieta == 33 && wedgeNo == 5)
   return 1;
 return 0;
}

int isHADSpecialTow(int ieta,int iphi)
{
 int wedgeNo= (int)iphi/2;
 if (ieta == 33 && wedgeNo == 5)
   return 1;
 if (ieta == 34 && wedgeNo == 5)
   return 1;
 if (ieta == 35 && wedgeNo == 5)
   return 1;
 return 0;
}
  
int TypeOfBit(int emhad,int ieta)
{
  if(emhad%2==0){ //EM
    if(ieta>=16 && ieta<=35)
      return 1;
    else
      return 2;
  } else { //HAD
    if(ieta>=14 && ieta<=37)
      return 1;
    else
      return 2;
  }
}

//Not in use for the moment!!!
int isKnownProblem(int ihist, int ieta, int iphi, int run_number) {

  int wedge = (int) ( iphi/2 ); // 15 degrees wedge number: 0 - 23
	  
  int known_problem = 0; // false;
  if ( 
      //
      //________CENTRAL and WALL calorimeters___________
      //
      //____ CHA/23W/1 : CDF eta/wedge  = 24/23 
      (ihist == 1 && ieta == 24 &&  wedge == 23 && run_number <= 156487)
      //____ CHA/1E/2  : CDF eta/wedge = 28/1
      || (ihist == 1 && ieta == 28 && wedge == 1 && run_number <= 156487)
      //____ CHA/(6,7,11,16)W)/0 : CDF eta/phi = 25/6,7,11,16
      || (ihist == 1 && ieta == 25 &&
	  (wedge == 6 || wedge == 7 || wedge == 11 || wedge == 16) &&
	  run_number <= 156487
	  )
      //____ WHA/1E/9    : CDF eta/wedge = 35/1
      || ( (ihist == 1 && ieta == 35 && wedge == 1) &&
	   (152400 <= run_number && run_number <= 154211)
	   ) 
      //____ WHA/12E/10  : CDF eta/wedge = 36/12
      || ( (ihist == 1 && ieta == 36 && wedge == 12) &&
	   (151443 <= run_number && run_number <= 152382)
	   )
      //____ CEM/19W/7 (pmt0) : CDF eta/wedge =  18/19
      || ( (ihist == 0 && ieta == 18 && wedge == 19) &&
	   (147697 <= run_number && run_number <= 147736)
	   )
      //____ WHA/0W/10 (ptm0): CDF eta/wedge = 15/0  UNPLUGGED(DEAD)
      || ( (ihist == 1 && ieta == 15 && wedge == 0) &&
	   (164724 <= run_number && run_number <= 169000)
	   )

      //	      
      //_______ PLUG calorimeter________________
      // 
      //____ PEM/22E/13-1 :  CDF eta/wedge = 39/45--> 7.5deg wedges 
      ||( (ihist == 0 && ieta == 39 && iphi == 45) &&
	  (0 <= run_number && run_number <= 146313)
	  )
      //____ PHA/12W/11-1 :  CDF eta/wedge = 14/25
      || ( (ihist == 1 && ieta == 14 && iphi == 25) &&
	   (147814 <= run_number && run_number <= 148786)
	   )
      //____ PHA/16W/16-0 :  CDF eta/wedge =  9/32
      || ( (ihist == 1 && ieta ==  9 && iphi == 32) &&
	   (run_number == 147815)
	   )
      //____ PHA/ 2E/12-1 :  CDF eta/wedge = 38/5
      || ( (ihist == 1 && ieta == 38 && iphi == 5) &&
	   (145020 <= run_number && run_number <= 149730)
	   )
      //____ PHA/ 5E/16-1 :  CDF eta/wedge = 42/11
      || (ihist == 1 && ieta == 42 && iphi == 11 && run_number <= 156487)
      //____ PHA/10W/15-1 :  CDF eta/wedge = 10/21
      || (ihist == 1 && ieta == 10 && iphi == 21 && run_number <= 156487)
      ) {
    known_problem = 1; // true;
  } // end-if known problems

  return known_problem;

} // isKnownProblem()