//
// Pibero Djawwotho <pibero@indiana.edu>
// Indiana University
// 24 April 2007
//

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <iostream>

#include <cstring>




#include <cstdlib>
#include <ctime>
#include <cerrno>

#ifndef __ROOT__
#include "/usr/src/linux-2.4/include/asm/msr.h"
#else
#include "/usr/src/kernels/2.6.9-42.0.10.EL-smp-i686/include/asm-i386/msr.h"
#endif

class TH1F;
#include "TH1F.h"
#include "TFile.h"

#include "L2Result.h"
#include "Upsilon.hh"

const char* Upsilon::fVersion = "2007.04.24";
const float Upsilon::CPU_MHZ = 1603.680;
const char* Upsilon::fHeader = "runNumber:vpdVertex:zdcVertex:bbcVertex:id1:id2:eta1:eta2:phi1:phi2:mult1:mult2:etrig1:etrig2:ecluster1:ecluster2:invMass:cosTheta:time:evSeen:evAccept:numberOfHotTowers:L2candidates\n";
//mUps = new TNtuple("ups","ups", fHeader);
const int Upsilon::fSizeOfHeader = strlen(fHeader);

//mCandidates = new TH1F("candidates","candidates",1000000,0.0,1000000.);

void Upsilon::init(int runNumber, int* userInt, float* userFloat)
{
  printf("Upsilon: Entering init(run=%d)\n", runNumber);
 //mCandidates = new TH1F("candidates","candidates",1000000,0.0,1000000.);
 //f = new TFile ("file.root","RECREATE");
  //
  // Check if finish was called after last run
  //
  if (fTuple) finish();

  //
  // Set run number, float and int run control parameters
  //
  fRunNumber = runNumber;

  fMinL0ClusterEnergy = userFloat[0];
  fMinL2ClusterEnergy = userFloat[1];
  fMinInvMass         = userFloat[2];
  fMaxInvMass         = userFloat[3];
  fMaxCosTheta        = userFloat[4];

  fL0SeedThreshold = userInt[0];
  fL2SeedThreshold = userInt[1];
  fHotThreshold    = userInt[2];
  fUseVertex       = userInt[3];
  fMaxHotTowers    = userInt[4];

  printf("Upsilon: run number, userInt and userFloat loaded\n");

  mTest = 10;// this is for test ddas

  cout<<" This is the INIT in SIMULATOR:"<<endl;
  //
  // fL0SeedThreshold is always greather than fL2SeedThreshold
  //
  if (fL0SeedThreshold < fL2SeedThreshold) {
    int temp = fL0SeedThreshold;
    fL0SeedThreshold = fL2SeedThreshold;
    fL2SeedThreshold = temp;
    printf("Upsilon: Swapped L0 and L2 thresholds\n");
  }

  //
  // Reset event counters
  //
  fEventsSeen = 0;
  fEventsAccepted = 0;
  printf("Upsilon: Reset event counters\n");

  //
  // Clear BEMC status tables
  //
  memset(bemcTower, 0, sizeof(bemcTower));
  memset(bemcPatch, 0, sizeof(bemcPatch));
  memset(daqIdFromSoftId, 0, sizeof(daqIdFromSoftId));
  printf("Upsilon: Cleared BEMC status tables\n");

  //
  // Dynamic masking of hot towers
  //
  memset(fHotLog, 0, sizeof(fHotLog));
  memset(fHits, 0, sizeof(fHits));
  memset(fHot, 0, sizeof(fHot));
  fNumberOfHotTowers = 0;
  printf("Upsilon: Cleared data structures for dynamic masking of hot towers\n");

  //
  // Log init
  //
  char pathname[1024];
  sprintf(pathname, "data/run%d.l2ups.log", fRunNumber);
  fLog = fopen(pathname, "w");
  if (!fLog) {
    printf("Can't open %s: %s\n", pathname, strerror(errno));
    fLog = stdout;
  }
  printf("Upsilon: Log file opened\n");

  // Start of run summary
  fStartTime = time(0);

  fprintf(fLog, "Version %s\n", fVersion);
  fprintf(fLog, "Compiled on %s %s\n", __DATE__, __TIME__);
  fprintf(fLog, "gcc version %s\n", __VERSION__);
  fprintf(fLog, "Source %s\n", __FILE__);
  fprintf(fLog, "Run Number: %d\n", fRunNumber);
  fprintf(fLog, "Start Time: %s", ctime(&fStartTime));
  fprintf(fLog, "userInt[0]: L0 seed threshold: %d\n", fL0SeedThreshold);
  fprintf(fLog, "userInt[1]: L2 seed threshold: %d\n", fL2SeedThreshold);
  fprintf(fLog, "userInt[2]: Hot tower threshold: %d\n", fHotThreshold);
  fprintf(fLog, "userInt[3]: Use vertex: %d\n", fUseVertex);
  fprintf(fLog, "userInt[4]: Maximum number of hot towers: %d\n", fMaxHotTowers);
  fprintf(fLog, "userFloat[0]: L0 cluster energy threshold (GeV): %f\n", fMinL0ClusterEnergy);
  fprintf(fLog, "userFloat[1]: L2 cluster energy threshold (GeV): %f\n", fMinL2ClusterEnergy);
  fprintf(fLog, "userFloat[2]: Minimum invariant mass (GeV): %f\n", fMinInvMass);
  fprintf(fLog, "userFloat[3]: Maximum invariant mass (GeV): %f\n", fMaxInvMass);
  fprintf(fLog, "userFloat[4]: Maximum cos(theta): %f\n", fMaxCosTheta);

  printf("Upsilon: Version %s\n", fVersion);
  printf("Upsilon: Compiled on %s %s\n", __DATE__, __TIME__);
  printf("Upsilon: gcc version %s\n", __VERSION__);
  printf("Upsilon: Source %s\n", __FILE__);
  printf("Upsilon: Run Number: %d\n", fRunNumber);
  printf("Upsilon: Start Time: %s", ctime(&fStartTime));
  printf("Upsilon: userInt[0]: L0 seed threshold: %d\n", fL0SeedThreshold);
  printf("Upsilon: userInt[1]: L2 seed threshold: %d\n", fL2SeedThreshold);
  printf("Upsilon: userInt[2]: Hot tower threshold: %d\n", fHotThreshold);
  printf("Upsilon: userInt[3]: Use vertex: %d\n", fUseVertex);
  printf("Upsilon: userInt[4]: Maximum number of hot towers: %d\n", fMaxHotTowers);
  printf("Upsilon: userFloat[0]: L0 cluster energy threshold (GeV): %f\n", fMinL0ClusterEnergy);
  printf("Upsilon: userFloat[1]: L2 cluster energy threshold (GeV): %f\n", fMinL2ClusterEnergy);
  printf("Upsilon: userFloat[2]: Minimum invariant mass (GeV): %f\n", fMinInvMass);
  printf("Upsilon: userFloat[3]: Maximum invariant mass (GeV): %f\n", fMaxInvMass);
  printf("Upsilon: userFloat[4]: Maximum cos(theta): %f\n", fMaxCosTheta);

  //
  // Convert thresholds from 6-bit to 12-bit ADC
  //
  fL0SeedThreshold <<= BIT_SHIFT;
  fL2SeedThreshold <<= BIT_SHIFT;
  printf("Upsilon: Bit shifted thresholds for use with 12-bit raw ADCs\n");

  //
  // Warning: These tables must be loaded in this order!!!
  //
  readBemcTower();
  readBemcStatus();
  readBemcKillTable();
  printf("Upsilon: Loaded BEMC status tables\n");

  //
  // Allocate memory for ntuple and write header
  //
  fMaxEvents = 25000;
  fTuple = (char*)malloc(fSizeOfHeader+fMaxEvents*sizeof(Eventups));
  if (!fTuple) {
    fprintf(fLog, "Can't allocate memory for ntuple: %s\n", strerror(errno));
    printf("Can't allocate memory for ntuple: %s\n", strerror(errno));
  }
  else {
    memcpy(fTuple, fHeader, fSizeOfHeader);
    printf("Upsilon: Memory allocated for %d events\n", fMaxEvents);
  }

  //
  // Flush log file to disk
  //
  if (fflush(fLog)) {
    fprintf(fLog, "Can't flush %s: %s", pathname, strerror(errno));
    printf("Can't flush %s: %s", pathname, strerror(errno));    
  }

  printf("Upsilon: Leaving init(run=%d)\n", runNumber);
}

bool Upsilon::run(TrgDataType* trgData, unsigned short* bemcData)
{
  if (!trgData || !bemcData) return false;

  //
  // We have a valid event
  //
  unsigned int timeStart;
  unsigned int timeStop;

  //
  // Start timer here
  //
  rdtscl(timeStart);
  ++fEventsSeen;

  //
  // Save the trigger data and bemc data pointers for use in member functions
  //
  this->trgData = trgData;
  this->bemcData = bemcData;

  //
  // Set z-vertex
  //
  fVpdVertex = vpdVertex();
  fZdcVertex = zdcVertex();
  fBbcVertex = bbcVertex();

  myevent.VpdVertex=999;
  myevent.ZdcVertex=999;
  myevent.BbcVertex=999;
  myevent.all_tower_list.clear();
  myevent.matched_track_list.clear();
  fill_all_towers();
  myevent.VpdVertex=vpdVertex();
  myevent.ZdcVertex=zdcVertex();
  myevent.BbcVertex=bbcVertex();

  int numberOfTowersL0 = 0;
  int numberOfTowersL2 = 0;
  //int numberOfL2Candidates = 0; //this is to test no of candidates for QM08
  int L0_HT[4800];
  int L2_HT[4800];

  findSeedTowers(L0_HT, L2_HT, numberOfTowersL0, numberOfTowersL2);

  bool accept = false;

  int id1 = -1;//this is a DAQ ID
  int id2 = -1;//this is a DAQ ID
  float cosTheta = 0;
  float invMass = 0;
  BemcCluster cluster1;
  BemcCluster cluster2;

  fNumberOfL2Candidates=0;
  mL2Candidates.clear();
//  cout << "At beginning of candidate loop, container size = " << mL2Candidates.size() << endl;

  //cout<<"ID1==="<<id1<<endl;
//  cout<<" No of L0 Towers: "<< numberOfTowersL0 << endl;	
//  cout<<" No of L2 Towers: "<< numberOfTowersL2 << endl;

  for (int i = 0; i < numberOfTowersL0; ++i) {
       
    
    id1 = L0_HT[i];
    
//    cout<<" L0 tower SOFT ID :"<<bemcTower[id1].softId<<endl;     
   
    calcCluster(id1, cluster1);

//     if(bemcTower[id1].softId==662)cout<<"Cluster1 position x:y:z:E ="<<cluster1.x<<" "<<cluster1.y<<" "<<cluster1.z<<" "<<cluster1.E<<endl;

    if (cluster1.E < fMinL0ClusterEnergy || cluster1.E == 0) continue;
    for (int j = 0; j < numberOfTowersL2; ++j) {
      id2 = L2_HT[j];
      if (id1 == id2) continue;
      calcCluster(id2, cluster2);

 //     cout<<" L2 tower Soft ID :"<< bemcTower[id2].softId <<endl;

 //    if(bemcTower[id2].softId == 2695)cout<<"Cluster2 position x:y:z:E ="<<cluster2.x<<" "<<cluster2.y<<" "<<cluster2.z<<" "<<cluster2.E<<endl;

      if (cluster2.E < fMinL2ClusterEnergy || cluster2.E == 0) continue;
      cosTheta = (cluster1.x * cluster2.x + cluster1.y * cluster2.y + cluster1.z * cluster2.z);
      invMass = sqrt(2 * cluster1.E * cluster2.E * (1 - cosTheta));
      if ((cosTheta < fMaxCosTheta) && (fMinInvMass < invMass && invMass < fMaxInvMass)) {
	accept = true;
        //++fEventsAccepted;
	if (fNumberOfL2Candidates == 0) ++fEventsAccepted;//for qm08
        ++fNumberOfL2Candidates;// for qm08
       // mCandidates->Fill(numberOfL2Candidates);
         //cout<<numberOfL2Candidates<<endl;
	//goto fini;
        pair<int,int> l2Pair(bemcTower[id1].softId,bemcTower[id2].softId);
 //       cout << "first, second towerID from simulator:"<<l2Pair.first << ", " << l2Pair.second << endl;
        mL2Candidates.push_back(l2Pair);
//	cout<<"Soft id1 :"<<bemcTower[id1].softId << ", E1:" << cluster1.E <<endl;
//	cout<<"Soft id2 :"<<bemcTower[id2].softId << ", E2:" << cluster2.E <<endl;
//	cout <<"pair inv m: " << invMass << ", cosTheta: " <<cosTheta << endl;
      }
    }
    //cout<<"id1 :"<<id1<<endl;
    //cout<<"id2 :"<<id2<<endl;
  }

// cout<<"upsilon candidates: "<<fNumberOfL2Candidates<<endl;
// cout<<"id1 outside loop:"<<id1<<endl;
// cout<<"id2 outside loop:"<<id2<<endl;
 
// fini:

  L2Result result;
  memset(&result, 0, sizeof(L2Result));

  if (accept) result.tag |= 1;

  result.numberOfTowersL0 = numberOfTowersL0;
  result.numberOfTowersL2 = numberOfTowersL2;

  if (id1 >= 0) result.triggerTowerL0 = bemcTower[id1].softId;
  if (id2 >= 0) result.triggerTowerL2 = bemcTower[id2].softId;

//  cout<<" Within fini: triggerTowerL0 "<<  result.triggerTowerL0 <<endl;
//  cout<<" Within fini: triggerTowerL2 "<<  result.triggerTowerL2 <<endl;

  result.cosTheta = cosTheta;
  result.invMass = (unsigned short)(invMass * 1e3);
  result.eventsSeen = fEventsSeen;
  result.eventsAccepted = fEventsAccepted;
  //result.numberOfL2Candidates = numberOfL2Candidates; 
  if (id1 >= 0) result.energyOfClusterL0 = (unsigned short)(cluster1.E * 1e3);
  if (id2 >= 0) result.energyOfClusterL2 = (unsigned short)(cluster2.E * 1e3);

  result.vertexZ = fVpdVertex;

  result.numberOfHotTowers = fNumberOfHotTowers;

  //
  // Stop timer here
  //
  rdtscl(timeStop);
  if (timeStop > timeStart) result.time = (unsigned short)((timeStop - timeStart) / CPU_MHZ);

  //
  // Save L2 result into trigger data
  //
 // memcpy(trgData->TrgSum.L2Result+L2RESULTS_2007_OFFSET_UPS, &result, sizeof(L2Result));
  memcpy(trgData->TrgSum.L2Result+L2RESULTS_2008_OFFSET_UPS, &result, sizeof(L2Result));

  if (accept) {
    //
    // Write to nutple
    //
    Eventups event = {
      fRunNumber,
      fVpdVertex,
      fZdcVertex,
      fBbcVertex,
      bemcTower[id1].softId,
      bemcTower[id2].softId,
      bemcTower[id1].eta,
      bemcTower[id2].eta,
      bemcTower[id1].phi,
      bemcTower[id2].phi,
      numberOfTowersL0,
      numberOfTowersL2,
      calcEnergy(id1),
      calcEnergy(id2),
      cluster1.E,
      cluster2.E,
      invMass,
      cosTheta,
      result.time,
      fEventsSeen,
      fEventsAccepted,
      fNumberOfHotTowers,
      fNumberOfL2Candidates
    };

    if (fEventsAccepted > fMaxEvents) {
      fMaxEvents *= 2;
      fTuple = (char*)realloc(fTuple, fSizeOfHeader+fMaxEvents*sizeof(Eventups));
    }

    memcpy(fTuple+fSizeOfHeader+(fEventsAccepted-1)*sizeof(Eventups), &event, sizeof(Eventups));
  }

  return accept;
}

void Upsilon::finish()
{

// TFile f("file.root","RECREATE");
  //f->Write();
  //f->Close();

  printf("Upsilon: Entering finish(run=%d)\n", fRunNumber);

  // Write ntuple to file
  if (fLog && fTuple) {
    // Open bin file
    char pathname[1024];
    sprintf(pathname, "data/run%d.l2ups.bin", fRunNumber);
    int fd = creat(pathname, 0644 | O_SYNC);
    if (fd == -1) {
      fprintf(fLog, "Can't open %s: %s\n", pathname, strerror(errno));
    }
    else {
      // Write ntuple to bin file
      int count = fSizeOfHeader+fEventsAccepted*sizeof(Eventups);
      int nWritten = write(fd, fTuple, count);
      if (nWritten == -1)
	fprintf(fLog, "Can't write to %s: %s\n", pathname, strerror(errno));
      else if (nWritten != count)
	fprintf(fLog, "Wrote only %d out of %d bytes to %s\n", nWritten, count, pathname);
      else
	fprintf(fLog, "%s written successfully\n", pathname);

      if (close(fd) == -1)
	fprintf(fLog, "Can't close %s\n", pathname);
    }

    // Free ntuple memory
    free(fTuple);
    fTuple = 0;

    printf("Upsilon: %d events written\n", fEventsAccepted);

    // End of run summary
    fStopTime = time(0);
    fprintf(fLog, "Stop Time: %s", ctime(&fStopTime));
    fprintf(fLog, "Events seen: %d\n", fEventsSeen);
    fprintf(fLog, "Events accepted: %d\n", fEventsAccepted);
    //fprintf(fLog, "upsilon candidates:%d\n",fnumberOfL2Candidates);
    printf("Upsilon: Stop Time: %s", ctime(&fStopTime));
    printf("Upsilon: Events seen: %d\n", fEventsSeen);
    printf("Upsilon: Events accepted: %d\n", fEventsAccepted);
    //printf("upsilon:candidates passed L2:%d\n", fnumberOfL2Candidates);
    // Dynamically masked hot towers
    fprintf(fLog, fHotLog);
    printf(fHotLog);
    fprintf(fLog, "Number of hot towers dynamically masked: %d\n", fNumberOfHotTowers);
    printf("Upsilon: Number of hot towers dynamically masked: %d\n", fNumberOfHotTowers);

    // Hot towers not masked
    findHotTowers();

    if (fLog != stdout && fclose(fLog))
      printf("Can't close log file: %s\n", strerror(errno));

    fLog = 0;

    printf("Upsilon: Log file closed\n");
  }

  printf("Upsilon: Leaving finish(run=%d)\n", fRunNumber);
}

void Upsilon::calcCluster(int daqId, BemcCluster& cluster)
{
  memset(&cluster, 0, sizeof(BemcCluster));

  //
  // Calculate the energies of all neighbors
  //
  float energy = calcEnergy(daqId);
  float energies[8];

  memset(energies, 0, sizeof(energies));

  for (int i = 0; i < bemcTower[daqId].numberOfNeighbors; ++i) {
    int daqId2 = bemcTower[daqId].neighbor[i];
    if (!towerUnmasked(daqId2) || bemcData[daqId2] < bemcTower[daqId2].pedestal) {
      energies[i] = 0;
    }
    else {
      energies[i] = calcEnergy(daqId2);
      if (energy < energies[i]) return;
    }
  }

  //
  // Sort the high tower neighbors in order of descending energy.
  // Use an insertion sort algorithm which is fairly efficient
  // for small arrays.
  //
  for (int i = 1; i < bemcTower[daqId].numberOfNeighbors; ++i) {
    int j = i;
    float indexE = energies[j];
    int indexId = bemcTower[daqId].neighbor[j];
    while (j > 0 && indexE > energies[j-1]) {
      energies[j] = energies[j-1];
      bemcTower[daqId].neighbor[j] = bemcTower[daqId].neighbor[j-1];
      --j;
    }
    energies[j] = indexE;
    bemcTower[daqId].neighbor[j] = indexId;
  }

  //
  // Calculate energy-weighted centroid of cluster
  //
  cluster.E = energy;
  cluster.x = bemcTower[daqId].x * energy;
  cluster.y = bemcTower[daqId].y * energy;
  cluster.z = bemcTower[daqId].z * energy;

  for (int i = 0; i < NUMBER_OF_TOWERS_PER_CLUSTER - 1; ++i) {
    int daqId2 = bemcTower[daqId].neighbor[i];
    cluster.E += energies[i];
    cluster.x += energies[i] * bemcTower[daqId2].x;
    cluster.y += energies[i] * bemcTower[daqId2].y;
    cluster.z += energies[i] * bemcTower[daqId2].z;
  }

  cluster.x /= cluster.E;
  cluster.y /= cluster.E;
  cluster.z /= cluster.E;

  //
  // Correct z-vertex
  //
  if (fUseVertex) cluster.z -= fVpdVertex;

  //
  // Normalize (x, y, z) vector to unity
  //
  float norm = sqrt(cluster.x * cluster.x + cluster.y * cluster.y + cluster.z * cluster.z);

  if (norm) {
    cluster.x /= norm;
    cluster.y /= norm;
    cluster.z /= norm;
  }
}

void Upsilon::readBemcTower()
{
  const char* pathname = "bemcTower.dat";
  FILE* fp = fopen(pathname, "r");
  if (!fp) {
    fprintf(fLog, "Can't open %s: %s\n", pathname, strerror(errno));
    fprintf(fLog, "Can't run without %s\n", pathname);
    return;
  }
  fprintf(fLog, "Reading %s\n", pathname);
  char line[1024];
  int lineNumber = 0;
  while (fgets(line, 1024, fp)) {
    ++lineNumber;
    if (line[0] == '#') continue;
    int daqId, softId, crate, crateSeq;
    float gain, x, y, z, eta, phi;
    int numberOfNeighbors;
    int neighbors[8];
    sscanf(line, "%d %d %d %d %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d\n", &daqId, &softId, &crate, &crateSeq, &gain, &x, &y, &z, &eta, &phi, &numberOfNeighbors, &neighbors[0], &neighbors[1], &neighbors[2], &neighbors[3], &neighbors[4], &neighbors[5], &neighbors[6], &neighbors[7]);
    if (!checkDaqId(daqId))
      fprintf(fLog, "Bad daqId=% in %s:%d\n", daqId, pathname, lineNumber);
    if (!checkSoftId(softId))
      fprintf(fLog, "Bad softId=%d in %s:%d\n", softId, pathname, lineNumber);
    if (!checkCrate(crate))
      fprintf(fLog, "Bad crate=%d in %s:%d\n", crate, pathname, lineNumber);
    if (!checkCrateSeq(crateSeq))
      fprintf(fLog, "Bad crateSeq=%d in %s:%d\n", crateSeq, pathname, lineNumber);
    if (!checkGain(gain))
      fprintf(fLog, "Bad gain=%f in %s:%d\n", gain, pathname, lineNumber);
    if (!checkXYZ(x, y, z))
      fprintf(fLog, "Bad x=%f, y=%f, z=%f in %s:%d\n", x, y, z, pathname, lineNumber);
    if (!checkEtaPhi(eta, phi))
      fprintf(fLog, "Bad eta=%f, phi=%f in %s:%d\n", eta, phi, pathname, lineNumber);
    if (!checkNumberOfNeighbors(numberOfNeighbors))
      fprintf(fLog, "Bad number of tower neighbors=%d in %s:%d\n", numberOfNeighbors, pathname, lineNumber);
    bemcTower[daqId].softId = softId;
    bemcTower[daqId].crate = crate;
    bemcTower[daqId].crateSeq = crateSeq;
    bemcTower[daqId].gain = gain;
    bemcTower[daqId].x = x;
    bemcTower[daqId].y = y;
    bemcTower[daqId].z = z;
    bemcTower[daqId].eta = eta;
    bemcTower[daqId].phi = phi;
    bemcTower[daqId].numberOfNeighbors = numberOfNeighbors;
    daqIdFromSoftId[softId] = daqId;
    for (int i = 0; i < numberOfNeighbors; ++i) {
      if (!checkDaqId(neighbors[i]))
	fprintf(fLog, "Bad daqIdOfNeighbor=%d in %s:%d\n", neighbors[i], pathname, lineNumber);
      bemcTower[daqId].neighbor[i] = neighbors[i];
    }
  }
  if (fclose(fp))
    fprintf(fLog, "Can't close %s: %s\n", pathname, strerror(errno));
  fprintf(fLog, "Done reading %s\n", pathname);
  if (fflush(fLog))
    fprintf(fLog, "Can't flush log file: %s\n", strerror(errno));
}

void Upsilon::readBemcStatus()
{
  //
  // Open global bemcStatus.txt if available, otherwise revert to local copy
  //
  const char* pathname = "/home/startrg/trg/cfg/Tier1/DSM_LUT/bemcStatus.txt";
  FILE* fp = fopen(pathname, "r");
  if (!fp) {
    fprintf(fLog, "Can't open %s: %s\n", pathname, strerror(errno));
    pathname = "bemcStatus.txt";
    fp = fopen(pathname, "r");
    if (!fp) {
      fprintf(fLog, "Can't open %s: %s\n", pathname, strerror(errno));
      fprintf(fLog, "Can't run without %s: %s\n", pathname, strerror(errno));
      return;
    }
  }
  fprintf(fLog, "Reading %s\n", pathname);

  //
  // Fill bemcTower and bemcPatch
  //
  char line[1024];
  int lineNumber = 0;
  while (fgets(line, 1024, fp)) {
    ++lineNumber;
    if (line[0] == '#')
      continue;
    else if (strncmp(line, "SoftId", 6) == 0) {
      int softId, crate, crateSeq, towerUnmasked, patchUnmaskedInHT, patchUnmaskedInSum, patch;
      float pedestal;
      sscanf(line, "SoftId %d %d %d %d %d %d %f %d\n", &softId, &crate, &crateSeq, &towerUnmasked, &patchUnmaskedInHT, &patchUnmaskedInSum, &pedestal, &patch);
      int daqId = daqIdFromSoftId[softId];
      if (!checkSoftId(softId))
	fprintf(fLog, "Bad softId=%d in %s:%d\n", softId, pathname, lineNumber);
      if (!checkDaqId(daqId))
	fprintf(fLog, "Bad daqId=%d from softId=%d in %s:%d\n", daqId, softId, pathname, lineNumber);
      if (!checkCrate(crate))
	fprintf(fLog, "Bad crate=%d in %s:%d\n", crate, pathname, lineNumber);
      if (!checkCrateSeq(crateSeq))
	fprintf(fLog, "Bad crateSeq=%d in %s:%d\n", crateSeq, pathname, lineNumber);
      if (!checkPatch(patch))
	fprintf(fLog, "Bad patch=%d in %s:%d\n", patch, pathname, lineNumber);
      bemcTower[daqId].crate = crate;
      bemcTower[daqId].crateSeq = crateSeq;
      bemcTower[daqId].towerUnmasked = towerUnmasked;
      bemcTower[daqId].patchUnmaskedInHT = patchUnmaskedInHT;
      bemcTower[daqId].patchUnmaskedInSum = patchUnmaskedInSum;
      bemcTower[daqId].pedestal = pedestal;
      bemcTower[daqId].patch = patch;
    }
    else if (strncmp(line, "TriggerPedestalShift", 20) == 0) {
      sscanf(line, "TriggerPedestalShift %d\n", &pedestalShift);
    }
    else if (strncmp(line, "triggerPatch", 12) == 0) {
      int patch, crate, cratePatch, unmaskedInHT, unmaskedInSum, bitConvMode, formula, parameters[6];
      sscanf(line, "triggerPatch %d %d %d %d %d %d %d %d %d %d %d %d %d", &patch, &crate, &cratePatch, &unmaskedInHT, &unmaskedInSum, &bitConvMode, &formula, &parameters[0], &parameters[1], &parameters[2], &parameters[3], &parameters[4], &parameters[5]);
      if (!checkPatch(patch))
	fprintf(fLog, "Bad patch=%d in %s:%d\n", patch, pathname, lineNumber);
      if (!checkCrate(crate))
	fprintf(fLog, "Bad crate=%d in %s:%d\n", crate, pathname, lineNumber);
      if (!checkCratePatch(cratePatch))
	fprintf(fLog, "Bad cratePatch=%d in %s:%d\n", cratePatch, pathname, lineNumber);
      bemcPatch[patch].crate = crate;
      bemcPatch[patch].cratePatch = cratePatch;
      bemcPatch[patch].unmaskedInHT = unmaskedInHT;
      bemcPatch[patch].unmaskedInSum = unmaskedInSum;
      bemcPatch[patch].bitConvMode = bitConvMode;
      bemcPatch[patch].formula = formula;
      memcpy(bemcPatch[patch].parameters, parameters, sizeof(parameters));
    }
  }
  if (fclose(fp))
    fprintf(fLog, "Can't close %s: %s\n", pathname, strerror(errno));
  fprintf(fLog, "Done reading %s\n", pathname);
  if (fflush(fLog))
    fprintf(fLog, "Can't flush log file: %s\n", strerror(errno));

  //
  // Save a copy of bemcStatus.txt
  //
  char pathname2[1024];
  sprintf(pathname2, "data/run%d.l2ups.bemcStatus.txt", fRunNumber);
  char command[1024];
  sprintf(command, "/bin/cp %s %s", pathname, pathname2);
  if (system(command) == -1)
    fprintf(fLog, "Can't copy %s to %s\n", pathname, pathname2);
}

void Upsilon::readBemcKillTable()
{
  const char* pathname = "bemcKill.dat";
  FILE* fp = fopen(pathname, "r");
  if (!fp) return;
  fprintf(fLog, "Reading %s\n", pathname);
  char line[1024];
  while (fgets(line, 1024, fp)) {
    if (line[0] == '#') continue;
    int daqId;
    int softId;
    int kill;
    sscanf(line, "%d %d %d\n", &daqId, &softId, &kill);
    if (softId != bemcTower[daqId].softId) {
      fprintf(fLog, "softId=%d from %s does not match softId=%d from bemcTower.dat for daqId=%d\n",
	      softId, pathname, bemcTower[daqId].softId, daqId);
    }
    else {
      bemcTower[daqId].towerUnmasked = !kill;
    }
  }
 if (fclose(fp)) {
    fprintf(fLog, "Can't close %s: %s\n", pathname, strerror(errno));
  }
  fprintf(fLog, "Done reading %s\n", pathname);
  if (fflush(fLog))
    fprintf(fLog, "Can't flush log file: %s\n", strerror(errno));
}

void Upsilon::findSeedTowers(int* L0_HT, int* L2_HT, int& numberOfTowersL0, int& numberOfTowersL2)
{
  for (int daqId = 0; daqId < NUMBER_OF_TOWERS; ++daqId) {
    if (towerUnmasked(daqId)) {
      float adc = bemcData[daqId] - bemcTower[daqId].pedestal;
      if (adc >= fL2SeedThreshold) {
	if (fHotThreshold && hotTower(daqId)) continue;
	L2_HT[numberOfTowersL2++] = daqId;
	if (adc >= fL0SeedThreshold) {
	  L0_HT[numberOfTowersL0++] = daqId;
	}
      }
    }
  }
}

void Upsilon::fill_all_towers()
{
	  for (int daqId = 0; daqId < NUMBER_OF_TOWERS; ++daqId) {
			if (towerUnmasked(daqId)) {
				  float adc = bemcData[daqId] - bemcTower[daqId].pedestal;
				  float energy = calcEnergy(daqId);
				  if (energy >0.0) {
					  struct_BemcTower tower;
					  tower.x=bemcTower[daqId].x;
					  tower.y=bemcTower[daqId].y;
					  tower.z=bemcTower[daqId].z;
					  tower.E=energy;
					  tower.adc=adc;
//					  tower.softid=daqId;
					  tower.softid=bemcTower[daqId].softId;
					  tower.eta=bemcTower[daqId].eta;
					  tower.phi=bemcTower[daqId].phi;
					  tower.numberOfNeighbors=bemcTower[daqId].numberOfNeighbors;
					  for (int i = 0; i < tower.numberOfNeighbors; ++i) {
						  tower.neighbor[i]= bemcTower[daqId].neighbor[i];
					  }
					  myevent.all_tower_list.push_back(tower);
				   }
			}
	  }
}

//
// Dynamic masking of hot towers
//
bool Upsilon::hotTower(int daqId)
{
  ++fHits[daqId];
  // Switch tower on and off
  if (fHits[daqId] > fHotThreshold) {
    // Turned on and has a lot of hits
    if ((float)fHits[daqId] / fEventsSeen > (float)fHotThreshold / NUMBER_OF_TOWERS) {
      // Has too many hits
      if (!fHot[daqId]) {
	// Don't switch off more than fMaxHotTowers
	if (fNumberOfHotTowers >= fMaxHotTowers) return false;
	// Is not yet bad
	fHot[daqId] = true;
	++fNumberOfHotTowers;
        // Log as turned bad with event number
	sprintf(fHotLog+strlen(fHotLog), "Tower %d turned bad with event %d\n", bemcTower[daqId].softId, fEventsSeen);
      }
    }
    else {
      // Has reasonable amount of hits
      if (fHot[daqId]) {
	// Is now good again but was previously bad
	fHot[daqId] = false;
	--fNumberOfHotTowers;
        // Log as turned good with event number
	sprintf(fHotLog+strlen(fHotLog), "Tower %d turned good with event %d\n", bemcTower[daqId].softId, fEventsSeen);
      }
    }
  }
  return fHot[daqId];
}

void Upsilon::findHotTowers()
{
  float mean = 0;
  float sigma = 0;
  const float SPREAD = 3;

  moment(NUMBER_OF_TOWERS, fHits, mean, sigma);

  fprintf(fLog, "List of hot towers (frequency = %f sigma above mean):\n", SPREAD);
  fprintf(fLog, "mean  : %f\n", mean );
  fprintf(fLog, "sigma : %f\n", sigma);
  fprintf(fLog, "+---------------------------------+\n");
  fprintf(fLog, "|  softId  |  nHits  |  nSigma    |\n");
  fprintf(fLog, "+---------------------------------+\n");
  for (int daqId = 0; daqId < NUMBER_OF_TOWERS; ++daqId) {
    float nSigma = sigma ? (fHits[daqId] - mean) / sigma : 0;
    if (nSigma > SPREAD) {
      fprintf(fLog, "| %8d | %7d | %10f |\n", bemcTower[daqId].softId, fHits[daqId], nSigma);
    }
  }
  fprintf(fLog, "+---------------------------------+\n");
  if (fflush(fLog))
    fprintf(fLog, "Can't flush log file: %s\n", strerror(errno));
}

//
// Calculate average (mean) and standard deviation (sigma) using
// the corrected two-pass algorithm to reduce roundoff errors.
//
// W.H. Press et al., Numerical Recipes in C, 2nd Ed, 1992, section 14.1, page 613
//
// http://www.nrbook.com/b/bookcpdf.php
//
void Upsilon::moment(int n, int* data, float& mean, float& sigma)
{
  mean = 0;
  sigma = 0;
  for (int i = 0; i < n; ++i) mean += data[i];
  mean /= n;
  float err = 0;
  for (int i = 0; i < n; ++i) {
    float s = data[i] - mean;
    err += s;
    sigma += s * s;
  }
  sigma = sqrt((sigma - err * err / n) / (n - 1));
}
