// Author: Andrew M. Rogers <mailto, NSCL 07/01/2006
//* Copyright (C) 2006-2008 Andrew M. Rogers
// Many thanks to Mark Wallace and Mike Famiano for their knowledge
// and any contributions to this code.  

#include <fstream>

#include <TCanvas.h>
#include <TClonesArray.h>
#include <TFriendElement.h>
#include <THtml.h>
#include <TObject.h>
#include <TStopwatch.h>
#include <TTree.h>

#include "THiRA.h"
  //#include "TCrange.h"


ClassImp(THiRA);


////////////////////////////////////////////////////////////////////////////////
/* BEGIN_HTML
<p> The THiRA class inherits from TObject.  This Class provides the structure for the 
HiRA device.  Raw event data is filled into the class members.                                   
</p>

<h3><a name="phil">I. Philosophy</a></h3>

<h3><a name="unpack">II. Unpacking</a></h3>

<h3><a name="coord">III. Coordinate System</a></h3>

<h3><a name="vars">IV. Important Data Members</a></h3>

<h4><a name="order"> III.a Readout Order</a></h4>
<p>While using the ASIC electronics for digitization of signals generated in the HiRA 
Si detectors it has become apparent that corrections based on readout order must be applied.  
The following data members and function methods are important in dealing with these, as we 
have observed so far, constant shifts of channels.
<h5><a name="orderMembers"> Data Members</a></h5>
<pre>
    (THiRASiDet) fCbSequence:    
    (THiRASiDet) fChSequence:
    (THiRASiDet) fHitSequence:
    (THiRASiDet) fHitSequence:
                 fHits[]:
</pre>
<h5><a name="orderMeth"> Methods</a></h5>
</p>

<h3><a name="rndm">IV. Random Number Generator</a></h3>
<p>A few notes should be made regarding the radomization of integer type data to floating 
point types.
<h5><a name="rClass"></h5>

</p>

<h3><a name="examples">V. Examples</a></h3>

<h3>Related References</h3>
<h4>[1] Wallace, M.S., Famiano, M.A., van Goethem, M.J., Rogers, A.M., Lynch, W.G., Clifford, J., Delaunay, F., Lee, J., Labostov, S., Mocko, M., Morris, L., Moroni, A., Nett, B., Oostdyk, D.J., Krishnasamy, R., Tsang, M.B., de Souza, R.T., Hudan, S., Sobotka, L.G., Charity, R.J., Elson, J., and Engel, G.L.: The high resolution array (HiRA) for rare isotope beam experiments, Nuclear Inst. and Methods in Physics Research, A, Elsevier, 2007. <a href="http://dx.doi.org/10.1016/j.nima.2007.08.248">[Full Text]</a></h4>
END_HTML */
////////////////////////////////////////////////////////////////////////////////


TClonesArray *THiRA::fgPixel = 0;

//______________________________________________________________________________
THiRA::THiRA(const THiRA &hira) : TObject(hira)
{
  // -- Copy constructor.
  ((THiRA&)hira).Copy(*this);
}

//______________________________________________________________________________
void THiRA::Copy(TObject &hira) const
{
  // -- Copy this method.
  
  TObject::Copy((TObject&)hira);
  ((THiRA&)hira).fName  = fName;
  ((THiRA&)hira).fTitle = fTitle;
  ((THiRA&)hira).fAnalysisState = fAnalysisState;
  // fTarget.Copy(((THiRA&)hira).fTarget);
  // ((THiRA&)hira).fTarget.parent = &((THiRA&)hira);
  ((THiRA&)hira).chipMapFile       = chipMapFile;
  ((THiRA&)hira).csiMapFile        = csiMapFile;
  ((THiRA&)hira).csiMapPath        = csiMapPath;
  ((THiRA&)hira).fReadoutShiftPath = fReadoutShiftPath;
  for(Int_t i=0; i<21; i++){for(Int_t j=0; j<32; j++){
      ((THiRA&)hira).adcTeleMap[i][j] = adcTeleMap[i][j];
      ((THiRA&)hira).CsIMap[i][j]     = CsIMap[i][j];}}
  ((THiRA&)hira).adcList  = adcList;
  ((THiRA&)hira).fAddRndm = fAddRndm;
  ((THiRA&)hira).fCorrectReadout = fCorrectReadout;
  ((THiRA&)hira).SiRawThresh = SiRawThresh;
  ((THiRA&)hira).CsIRawThresh = CsIRawThresh;
  ((THiRA&)hira).fEulerPhi = fEulerPhi;
  ((THiRA&)hira).fEulerTheta = fEulerTheta;
  ((THiRA&)hira).fEulerPsi = fEulerPsi;
  //  asic.Copy(((THiRA&)hira).asic);
  //  ((THiRA&)hira).asic.parent = &((THiRA&)hira);
  teleTList.Copy(((THiRA&)hira).teleTList);
  for(Int_t i=0; i<THIRA_NTELES; i++) {
    tele[i].Copy(((THiRA&)hira).tele[i]); 
    ((THiRA&)hira).tele[i].parent = &((THiRA&)hira);
    ((THiRA&)hira).teleTList.Add(&((THiRA&)hira).tele[i]);}

  
  for(Int_t i=0; i<THIRA_NTELES; i++) ((THiRA&)hira).fTele[i] = &((THiRA&)hira).tele[i];

  ((THiRA&)hira).fgPixel = new TClonesArray("THiRAPixel",10);
  ((THiRA&)hira).fPixel  = ((THiRA&)hira).fgPixel;

  ((THiRA&)hira).Clear();
}


//______________________________________________________________________________
void THiRA::Calculate()
{
  // -- Calculates all THiRA quantities.
  // This method calls all calculate methods of any member classes of THiRA such
  // as the Calculate() method of TFront, TTele, etc.  Any calculations related
  // to this class, THiRA, are also done here.


}


//______________________________________________________________________________
Int_t THiRA::Calibrate(Option_t* prefix,Option_t* path,Option_t* treeName,Option_t* option)
{
  // -- Calibrates THiRA components.
  // Calibrates all elements in the THiRA class such as Si and CsI slopes and
  // offsets.  Note that a calibration file must be loaded for any members that
  // are to be calibrated.

  if(!fChain){printf("TTree or TChain has not been initialized!!!\n"); return -1;}

  TStopwatch stopWatch;
  stopWatch.Start();

  Int_t  nTrees    = 1;
  TTree *inifChain = fChain;
  if(strcmp(fChain->ClassName(),"TChain")==0) nTrees = ((TChain*)fChain)->GetNtrees();

  // Call the copy constructor.
  THiRA *newHiRA = new THiRA(*this);   
  for(Int_t i=0; i<THIRA_NTELES; i++) newHiRA->tele[i].fCon = &tele[i];

  // Loop over all trees.
  Long64_t curEntry = 0;
  for(Int_t tt=0; tt<nTrees; tt++){
    fChain->LoadTree(curEntry);
    TTree *cTree = fChain->GetTree();
    // Check the tree analysis state.
    //    if(strcmp(cTree->GetUserInfo()->FindObject("Analysis State")->GetTitle(),"0")!=0) {printf("* ERROR: Analysis State!=0\n"); return -2;}
    curEntry = curEntry+cTree->GetEntries()+1;
    InitTree(cTree);

    // Set the output filename.
    TString fileName;
    TNamed *runNumber = (TNamed*)cTree->GetUserInfo()->FindObject("Run Number");
    TString fPath(path);
    TString fPrefix(prefix);
    fileName = runNumber->GetTitle();
    if(prefix == "") fileName = fPath+"run"+fileName+"Cal.root";
    else             fileName = fPath+fPrefix+fileName+"Cal.root";
    
    // File options.
    if(option=="RECREATE"){ cout << "Creating file: " << fileName << endl;
    }else {option="UPDATE"; cout << "Updating file: " << fileName << endl;}

    // Make a new file for our calibrated TTree.
    TString newTName((Char_t*)treeName);
    TString oldTName(cTree->GetName());
    TFile   newFile(fileName,option);
    TTree  *newTree = new TTree("Cal","Calibrated "+oldTName+" Tree",99);

    // Check for existing TTree.
    TNamed aState("Analysis State","1");
    if(newFile.FindObject(oldTName+"Cal")!=0){printf("Tree already exists!\n"); return -2;}
    else if(strcmp(newTName,"")==0){
      newTName = oldTName+"Cal";
      newTree->SetName(newTName);
      newTree->GetUserInfo()->Add(&aState);
    }else{
      newTree->SetName(newTName+"Cal");
      newTree->GetUserInfo()->Add(&aState);
    }

    // Copy UserInfo
    for(Int_t j=0; j<cTree->GetUserInfo()->GetEntries(); j++){
      TNamed *ntmp = (TNamed*)cTree->GetUserInfo()->At(j);
      TNamed *ctmp = (TNamed*)ntmp->Clone();
      if(strcmp(ntmp->GetName(),"Analysis State")!=0) newTree->GetUserInfo()->Add(ctmp);
    }

    // Branches.
    newHiRA->InitBranch(newTree,"Cal");
    newHiRA->InitTree(newTree);

    // Calibrate
    Long64_t nEntries = cTree->GetEntries();
    for(Long64_t ientry=0; ientry<nEntries; ientry++){
      Double_t prog = ((Double_t)ientry/nEntries)*100;
      if(ientry%100==0){printf("Entry: %15lli \t%4.2lf %%\r",ientry,prog); cout.flush();}
      for(Int_t ii=0; ii<THIRA_NTELES; ii++) newHiRA->tele[ii].Calibrate(ientry);
      newTree->Fill();
    }
    newTree->AutoSave();
    newFile.Close();
    InitTree(inifChain);
  }
  stopWatch.Stop();
  stopWatch.Print();

  return nTrees;
}


//______________________________________________________________________________
void THiRA::Clear(Option_t *option)
{
  // -- Clears data members.
  //

  fNwords  = 0;
  fRawMult = 0;
  memset(fHits,'\0',sizeof(fHits));

  fPixel->Clear();

  for (Int_t i=0; i<THIRA_NTELES; i++) tele[i].Clear();
}


//______________________________________________________________________________
void THiRA::ClearErr()
{
  // -- Clear all data members related to counting errors.
  //

  for(Int_t i=0; i<THIRA_NMBS; i++){
    for(Int_t j=0; j<THIRA_MB_SLOTS*2; j++){
      fChipErr[i][j] = 0;
      for(Int_t k=0; k<THIRA_SI_CHANNELS; k++){
	fChanErr[i][j][k] = 0;
      }
    }
  }
}


//______________________________________________________________________________
void THiRA::CreateFolders()
{
  // -- Create HiRA related folders in the current root file.
  //

  f_HiRA        = gROOT->GetRootFolder()->AddFolder("HiRA","The HiRA folder");
  f_pixelMatrix = f_HiRA->AddFolder("Pixel Matrices","Contains the matrices for the pixel positions");
  f_parameters  = f_HiRA->AddFolder("Parameters","Calibration and other HiRA parameters used for this tree");
  f_input       = f_HiRA->AddFolder("Input-Files","Input files used to generate this tree");

  f_HiRA->Write();

  printf("HiRA Folder Creation \t\t\t\t\t\t\t\t [OK]\n"); 
}


//______________________________________________________________________________
Long64_t THiRA::Draw(const Char_t* varexp, const Char_t* selection,
		     Option_t* option,Option_t* drawType)
{
  // -- Draw an expression for multiple telescope.
  // Many times when we wish to draw a histogram for HiRA we want to create it for 
  // all telescopes.  This method allows one to easily accomplish this task.
  //
  // 1. varexp is the normal expression one would use EXCEPT without the tele#. prefix.
  // 2. selection is any condition put on the output.
  // 3. option is the typical Draw option such "colz","surf", etc.
  // 4. drawType must be set to "all" as this is the only option implemented.  This 
  //    option produces a canvas with all telescopes drawn.
  //

  Char_t var[2000],*v1=NULL,*v2=NULL,*v3=NULL;
  //  Char_t var[2000],v1[2000],v2[2000],v3[2000];
  sprintf(var,"%s",varexp);

  Char_t delims[] = ":";
  Char_t *result = NULL;
  Int_t p=0;
  result = strtok(var,delims);
  v1 = result;
  while(result != NULL) {
    p++;
    result = strtok( NULL, delims );
    if(p==1) v2=result;
    if(p==2) v3=result;
  }  

  //  printf( "result is %s %s %s\n",v1,v2,v3 );
  cout << drawType << endl;
  if(strcmp(drawType,"all")==0){
    //    TCanvas *c_HiRA = new TCanvas("c_HiRA","HiRA telescope canvas");
    fHiRACanvas = new TCanvas("c_HiRA","HiRA telescope canvas");
    TCanvas *c_HiRA = fHiRACanvas;

    c_HiRA->Divide(THIRA_MAXCOL,THIRA_MAXROW);
    
    for(Int_t i=0; i<THIRA_NTELES; i++){
      c_HiRA->cd(GetPadIndex(i));
      Char_t expression[900];
      cout << p << endl;
      if(p==1) sprintf(expression,"tele%i.%s",teleID[i],v1);
      if(p==2) sprintf(expression,"tele%i.%s:tele%i.%s",teleID[i],v1,teleID[i],v2);
      if(p==3) sprintf(expression,"tele%i.%s:tele%i.%s:tele%i.%s",teleID[i],v1,teleID[i],v2,teleID[i],v3);
      cout << expression << " " << selection << " " << option << endl;
      fChain->Draw(expression,selection,option);
      c_HiRA->Modified();
      c_HiRA->Update();
    }
    return 1;
  }else if(drawType == ""){
    return 1;
  }

  return -1;
}


//______________________________________________________________________________
Bool_t THiRA::Fill(Short_t energy, Short_t time, Short_t teleNum, Short_t faceNum, 
		   Short_t stripNum, Short_t cbSeq, Short_t chSeq, UInt_t hitSeq)
{
  // -- Fill the telescope branches with data.
  // The following operations occur,
  //     1. The CB and Hit sequences are set.
  //     2. 

  switch(faceNum){
  case 0:
    #if THIRA_OP_DE
    tele[teleNum].DE.fCbSequence  = cbSeq;
    tele[teleNum].DE.fChSequence  = chSeq;
    tele[teleNum].DE.fHitSequence = hitSeq;
    tele[teleNum].DE.ERaw[stripNum] = energy;
    tele[teleNum].DE.TRaw[stripNum] = time;
    if(energy>SiRawThresh) tele[teleNum].DE.rawMult++;
    if(energy>tele[teleNum].DE.Emax){ tele[teleNum].DE.Emax = energy; tele[teleNum].DE.Emaxch = stripNum;}
    #endif
    break;
  case 1:
    tele[teleNum].EF.fCbSequence  = cbSeq;
    tele[teleNum].EF.fChSequence  = chSeq;
    tele[teleNum].EF.fHitSequence = hitSeq;
    tele[teleNum].EF.ERaw[stripNum] = energy;
    tele[teleNum].EF.TRaw[stripNum] = time;
    if(energy>SiRawThresh) {
      tele[teleNum].EF.rawMult++;
      fRawMult++;
    }
    if(energy>tele[teleNum].EF.Emax){ tele[teleNum].EF.Emax = energy; tele[teleNum].EF.Emaxch = stripNum;}
    break;
  case 2:
    tele[teleNum].EB.fCbSequence  = cbSeq;
    tele[teleNum].EB.fChSequence  = chSeq;
    tele[teleNum].EB.fHitSequence = hitSeq;
    tele[teleNum].EB.ERaw[stripNum] = energy;
    tele[teleNum].EB.TRaw[stripNum] = time;
    if(energy>SiRawThresh) tele[teleNum].EB.rawMult++;
    if(energy>tele[teleNum].EB.Emax){ tele[teleNum].EB.Emax = energy; tele[teleNum].EB.Emaxch = stripNum;}
    break;
  }

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::FillCsI(Short_t energy, Short_t teleNum, Short_t csiNum)
{
  // -- Fill the telescope CsI branches with data.
  //

  if(teleNum == -99){ printf("Fill CsI error: \n"); return kFALSE;}

  tele[teleNum].CsI.ERaw[csiNum] = energy;
  if(energy>CsIRawThresh) tele[teleNum].CsI.rawMult++;
  if(energy>tele[teleNum].CsI.Emax) { tele[teleNum].CsI.Emax  = (Double_t)energy; tele[teleNum].CsI.Emaxch  = csiNum;}
  if(energy>tele[teleNum].CsI.E2max){ tele[teleNum].CsI.E2max = (Double_t)energy; tele[teleNum].CsI.E2maxch = csiNum;}
  
  return kTRUE;
}


//______________________________________________________________________________
Int_t THiRA::GetCsIQuadIndex(Int_t xStrip, Int_t yStrip)
{
  // -- Get the index for the quadrant of silicon located behind a pixel.

  if(xStrip<16){
    if(yStrip<16) return 0;
    if(yStrip>15) return 3;
  }else if(xStrip>15){
    if(yStrip<16) return 1;
    if(yStrip>15) return 2;
  }

  return -1;
}


//______________________________________________________________________________
Int_t THiRA::GetPadIndex(Int_t index)
{
  // -- Get the TPad index corresponding to the physical position of a telescope.
  // When drawing histograms for HiRA it is sometimes desireable to create a 
  // canvas containing pads corresponding to each detector.  It is also convenient 
  // to draw these histograms in the actual order of the setup.  This method returns 
  // the index of a TPad defined by the towerMatrix[][] in the rootDefine.h file.  
  // This pad number should be the physical position of the detector in the setup.

  Int_t padIndex = 0;
  Int_t num = teleID[index];

  for(Int_t i=0; i<THIRA_MAXCOL; i++){
    for(Int_t j=0; j<THIRA_MAXROW; j++){
      padIndex++;
      if(towerMatrix[i][j] == num) return padIndex;
    }
  }

  return -1;
}


//______________________________________________________________________________
Int_t THiRA::GenPrime(Option_t* prefix,Option_t* path,Option_t* treeName,Option_t* option)
{
  // -- Generate HiRA primary data tree.
  //

  if(!fChain){printf("TTree or TChain has not been initialized!!!\n"); return -1;}

  TStopwatch stopWatch;
  stopWatch.Start();

  Int_t  nTrees    = 1;
  TTree *inifChain = fChain;
  if(strcmp(fChain->ClassName(),"TChain")==0) nTrees = ((TChain*)fChain)->GetNtrees();

  // Call the copy constructor.
  THiRA *newHiRA = new THiRA(*this);   
  for(Int_t i=0; i<THIRA_NTELES; i++) newHiRA->tele[i].fCon = &tele[i];

  // Loop over all trees.
  //  TNamed *uName=0;
  Long64_t curEntry = 0;
  for(Int_t tt=0; tt<nTrees; tt++){
    fChain->LoadTree(curEntry);
    TTree *cTree = fChain->GetTree();
    // Check the tree analysis state.
    if(strcmp(cTree->GetUserInfo()->FindObject("Analysis State")->GetTitle(),"1")!=0) {printf("* ERROR: Analysis State!=1\n"); return -2;}
    curEntry = curEntry+cTree->GetEntries()+1;
    InitTree(cTree);

    // Set the output filename.
    TString fileName;
    //    TNamed *runNumber = (TNamed*)cTree->GetFriend("E02023")->GetUserInfo()->FindObject("Run Number");
    //    TNamed *runNumber = (TNamed*)(((TFriendElement*)cTree->GetListOfFriends()->At(0))->GetTree()->GetUserInfo()->FindObject("Run Number"));
    TNamed *runNumber = (TNamed*)cTree->GetUserInfo()->FindObject("Run Number");
    TString fPath(path);
    TString fPrefix(prefix);
    fileName = runNumber->GetTitle();
    //    Char_t tmpC[100];
    //    sprintf(tmpC,"%i",371+tt);
    if(prefix == "") fileName = fPath+"run"+fileName+"Pri.root";
    else             fileName = fPath+fPrefix+fileName+"Pri.root";
    
    // File options.
    if(option=="RECREATE"){ cout << "Creating file: " << fileName << endl;
    }else {option="UPDATE"; cout << "Updating file: " << fileName << endl;}

    // Make a new file for our calibrated TTree.
    TString newTName((Char_t*)treeName);
    TString oldTName(cTree->GetName());
    TFile   newFile(fileName,option);
    TTree  *newTree = new TTree("Pri","Primary "+oldTName+" Tree",99);
    // Remove Cal ending to be replaced with Pri ending.
    oldTName.Chop();
    oldTName.Chop();
    oldTName.Chop();

    // Check for existing TTree.
    TNamed aState("Analysis State","2");
    if(newFile.FindObject(oldTName+"Pri")!=0){printf("Tree already exists!\n"); return -2;}
    else if(strcmp(newTName,"")==0){
      newTName = oldTName+"Pri";
      newTree->SetName(newTName);
      newTree->GetUserInfo()->Add(&aState);
    }else{
      newTree->SetName(newTName+"Pri");
      newTree->GetUserInfo()->Add(&aState);
    }

    // Copy UserInfo
    for(Int_t j=0; j<cTree->GetUserInfo()->GetEntries(); j++){
      TNamed *ntmp = (TNamed*)cTree->GetUserInfo()->At(j);
      TNamed *ctmp = (TNamed*)ntmp->Clone();
      if(strcmp(ntmp->GetName(),"Analysis State")!=0) newTree->GetUserInfo()->Add(ctmp);
    }

    //    uName = (TNamed*)newTree->GetUserInfo()->FindObject("Analysis State");
    //    uName->SetNameTitle("Analysis State","2");

    // Branches.
    //    newHiRA->InitBranch(newTree,"Pri");
    //    newHiRA->InitTree(newTree);
    newTree->Branch("fPixel",&fPixel);

    // Generate Primary Data.
    Int_t materials[3];
    materials[0] = fCRange->targetID("Si");
    materials[1] = fCRange->targetID("Mylar");
    materials[2] = fCRange->targetID((Char_t*)fTarget->GetName());
    Int_t    nHits;
    //    TBranch *b_hit;
    newTree->Branch("hira.nHits",&nHits,"nHits/I");
    //    newTree->SetBranchAddress("hira.fPixel",&newHiRA->fPixel);
    //    newTree->SetBranchAddress("hira.nHits",&nHits,&b_hit);
    Long64_t nEntries = cTree->GetEntries();
    for(Long64_t ientry=0; ientry<nEntries; ientry++){
      Double_t prog = ((Double_t)ientry/nEntries)*100;
      if(ientry%100==0){printf("Entry: %15lli \t%4.2lf %%\r",ientry,prog); cout.flush();}

      fPixel->Clear();
      for(Int_t j=0; j<THIRA_NTELES; j++) {
	#if THIRA_OP_DE
	tele[j].DE.b_E->GetEntry(ientry);
	#endif
	tele[j].EF.b_E->GetEntry(ientry);
	tele[j].EB.b_E->GetEntry(ientry);
	tele[j].CsI.b_E->GetEntry(ientry);
	tele[j].CsI.b_Emax->GetEntry(ientry);
      }

      nHits=0;
      UShort_t nCsIHits;
      Int_t    subEMaxCh[3];
      Double_t subEMax[3];
      // Generate hit pixels.
      // Currently no sorting done for high mult events or dE events.
      for(Int_t k=0; k<THIRA_NTELES; k++){
	nCsIHits=0;
	for(Int_t r=0; r<THIRA_CSI_CHANNELS; r++){
	  THiRAPixel aPix;
	  if(tele[k].CsI.E[r]>0){
	    for(Int_t ss=0; ss<3; ss++) subEMaxCh[ss]=99;
	    memset(subEMax,'\0',sizeof(subEMax));
	    nCsIHits++;
	    switch(r){
	    case 0:
	      for(Int_t tt=0;  tt<16; tt++) if(tele[k].EB.E[tt]>subEMax[2]){subEMax[2]=tele[k].EB.E[tt]; subEMaxCh[2]=tt;}
	      for(Int_t tt=0;  tt<16; tt++) if(tele[k].EF.E[tt]>subEMax[1]){subEMax[1]=tele[k].EF.E[tt]; subEMaxCh[1]=tt;}
	      break;
	    case 1:
	      for(Int_t tt=16; tt<32; tt++) if(tele[k].EB.E[tt]>subEMax[2]){subEMax[2]=tele[k].EB.E[tt]; subEMaxCh[2]=tt;}
	      for(Int_t tt=0;  tt<16; tt++) if(tele[k].EF.E[tt]>subEMax[1]){subEMax[1]=tele[k].EF.E[tt]; subEMaxCh[1]=tt;}
	      break;
	    case 2:
	      for(Int_t tt=16; tt<32; tt++) if(tele[k].EB.E[tt]>subEMax[2]){subEMax[2]=tele[k].EB.E[tt]; subEMaxCh[2]=tt;}
	      for(Int_t tt=16; tt<32; tt++) if(tele[k].EF.E[tt]>subEMax[1]){subEMax[1]=tele[k].EF.E[tt]; subEMaxCh[1]=tt;}
	      break;
	    case 3:
	      for(Int_t tt=0;  tt<16; tt++) if(tele[k].EB.E[tt]>subEMax[2]){subEMax[2]=tele[k].EB.E[tt]; subEMaxCh[2]=tt;}
	      for(Int_t tt=16; tt<32; tt++) if(tele[k].EF.E[tt]>subEMax[1]){subEMax[1]=tele[k].EF.E[tt]; subEMaxCh[1]=tt;}
	      break;
	    }
	    if(subEMaxCh[1]!=99 && subEMaxCh[2]!=99){
	      aPix.AddHit(fCRange,materials,&tele[k],subEMaxCh[1],subEMaxCh[2],(UShort_t)r,nHits,nCsIHits);
	      nHits++;
	    }
	  }else{
	    // Do nothing for the moment.
	  }
	}
      }
// 	for(Int_t m=0; m<THIRA_SI_CHANNELS; m++){
// 	  if(tele[k].DE.E[m]>0 || tele[k].EB.E[m]>0){
// 	    for(Int_t n=0; n<THIRA_SI_CHANNELS; n++){
// 	      if(tele[k].EF.E[n]>0){
		
// 		THiRAPixel aPix;
// 		THiRAPixel *bPix;
// 		bPix = aPix.AddHit(fCRange,materials,&tele[k],n,m,nHits);
// 		nHits++;
// 	      }
// 	    }
// 	  }
// 	}
//       }
      newTree->Fill();
    }
    newTree->AutoSave();
    newTree->GetUserInfo()->Clear();
    newFile.Close();
    InitTree(inifChain);
  }
  stopWatch.Stop();
  stopWatch.Print();

  return nTrees;
}

//______________________________________________________________________________
//void THiRA::GenPrime(Option_t *output,Option_t *treeName,Option_t *option)
// {
//   // -- Generate HiRA primary data tree.
//   //

//   TStopwatch stopWatch;
//   stopWatch.Start();

//   // Set the output filename.
//   TString fileName;
//   if(output == ""){
//     TNamed *runNumber = (TNamed*)fChain->GetUserInfo()->FindObject("Run Number");
//     fileName = runNumber->GetTitle();
//     fileName = "run"+fileName+"Pri.root";
//   }else fileName = output;

//   // File options.
//   if(option=="RECREATE"){ cout << "Creating file: " << fileName << endl;
//   }else {option="UPDATE"; cout << "Updating file: " << fileName << endl;}

//   // Make a new file for our prime TTree.
//   TFile   *newFile = new TFile(fileName,option);
//   TTree   *newTree;
//   TString oldTName(fChain->GetName());
//   // Remove Cal ending to be replaced with Pri ending.
//   oldTName.Chop();
//   oldTName.Chop();
//   oldTName.Chop();
  
//   // Check for existing TTree.
//   if(newFile->FindObject(oldTName+"Pri")!=0){newTree = (TTree*)newFile->FindObject(oldTName+"Pri");
//   }else{
//     newTree = new TTree(oldTName+"Pri","Primary "+oldTName+" Tree",99);
//     // We must set the entries of the calibrated tree since we will be filling individual branches.
//     newTree->SetEntries(fChain->GetEntries());
//   }

//   // Add friends.
//   //  newTree->AddFriend(fChain,"E02023Cal",kTRUE);

//   //  TTree *oldTree = fChain;
//   THiRA newHira(*this);
//   newHira.InitBranch(newTree,"Pri");
  
//   //  newHira.InitTree(newTree);

//   // GENERATE PRIMARY TTREE HERE!!!
//   // Initialize the energy-range calculator fCRange.
//   // Calculate the pix TClonesArray.
//   Int_t materials[3];
//   materials[0] = fCRange->targetID("Si");
//   materials[1] = fCRange->targetID("Mylar");
//   materials[2] = fCRange->targetID((Char_t*)fTarget->GetName());

// //   printf("Calulations of energy losses will be done for the following layers,\n");
  
// //   // Initialize the range tabels.
// //   // Proton
// //   fCRange->init_range_table(0, materials[0], 1, 1);
// //   fCRange->init_range_table(0, materials[1], 1, 1);
// //   fCRange->init_range_table(0, materials[2], 1, 1);
// //   // Deuetron
// //   fCRange->init_range_table(1, materials[0], 1, 2);
// //   fCRange->init_range_table(1, materials[1], 1, 2);
// //   fCRange->init_range_table(1, materials[2], 1, 2);
// //   // Triton
// //   fCRange->init_range_table(2, materials[0], 1, 3);
// //   fCRange->init_range_table(2, materials[1], 1, 3);
// //   fCRange->init_range_table(2, materials[2], 1, 3);
// //   // He3
// //   fCRange->init_range_table(3, materials[0], 2, 3);
// //   fCRange->init_range_table(3, materials[1], 2, 3);
// //   fCRange->init_range_table(3, materials[2], 2, 3);
// //   // He4
// //   fCRange->init_range_table(4, materials[0], 2, 4);
// //   fCRange->init_range_table(4, materials[1], 2, 4);
// //   fCRange->init_range_table(4, materials[2], 2, 4);
// //   printf("* Range tables initialized.\n");


//   Int_t nHits;
//   TBranch *b_hit;
//   newTree->Branch("hira.nHits",&nHits,"nHits/I");

//   newTree->SetBranchAddress("hira.fPixel",&newHira.fPixel,&newHira.b_fPixel);
//   newTree->SetBranchAddress("hira.nHits",&nHits,&b_hit);

//   newHira.b_fPixel->Reset();
//   Int_t counter = 0;
//   printf("Entries to be processed:  %15lli\n",newTree->GetEntries());
//   for(Int_t i=0; i<newTree->GetEntries(); i++){
//     newHira.Clear();
//     this->Clear();

//     if(counter == 100){
//       cout << "Processing Buffer:         " << i << "\r";
//       cout.flush();
//       counter=0;
//     }
//     counter++;

//     for(Int_t j=0; j<THIRA_NTELES; j++) {
//       tele[j].EF.b_E->GetEntry(i);
//       tele[j].EB.b_E->GetEntry(i);
//       tele[j].CsI.b_E->GetEntry(i);
//       tele[j].CsI.b_Emax->GetEntry(i);
//     }

//     nHits=0;
//     for(Int_t k=0; k<THIRA_NTELES; k++){
//       for(Int_t m=0; m<THIRA_SI_CHANNELS; m++){
// 	// Generate hit pixels.
// 	if(tele[k].DE.E[m]>0 || tele[k].EB.E[m]>0){
// 	  for(Int_t n=0; n<THIRA_SI_CHANNELS; n++){
// 	    if(tele[k].EF.E[n]>0){
// 	      //	      newHira.tele[k].DE.E[m]  = tele[k].DE.E[n];
// 	      newHira.tele[k].DE.E[m]  = 0;
// 	      newHira.tele[k].EF.E[n]  = tele[k].EF.E[n];
// 	      newHira.tele[k].EB.E[m]  = tele[k].EB.E[n];
// 	      newHira.tele[k].CsI.Emax = tele[k].CsI.Emax;
	      
// 	      THiRAPixel aPix;
// 	      THiRAPixel *bPix;
// 	      bPix = aPix.AddHit(fCRange,materials,&newHira.tele[k],n,m,nHits);
// 	      nHits++;
// 	    }
// 	  }
// 	}
//       }
//     }
//     newHira.b_fPixel->Fill();
//     b_hit->Fill();
//   }
  
//   newTree->AutoSave();
//   newFile->Close();

//   stopWatch.Stop();
//   stopWatch.Print();
  
//   printf("*Generated primary data.\n");
//}


//______________________________________________________________________________
Int_t THiRA::GetTeleIndex(Int_t id)
{
  // -- Gets the telescope index of a telescope with id.
  // There are two main integer numbers that describe a telescope.  The first is 
  // what we call the telescope id.  This is the number that is typically located 
  // on the back of the detector box.  It is the label for the telescope.
  //
  // The second is the telescope index.  This is the postion in the tele[] array.  
  // The purpose of this array is that it allows us to easily loop over all 
  // telescopes using a for() loop.
  // 
  // For example, if we have 16 telescopes then the tele[] array will have 16 
  // elements.  The telescope index then runs from 0-15.  However, we may be using 
  // telescopes 7,0,1,2,3,4,5,6,8,9,10,11,12,16,17,19.  These are not consecutive 
  // numbers that are difficult to loop over.
  //
  // NOTE:  While this is one way to loop over telescopes there is a more object
  // oriented method.  This method involves using the TIter class and iterating 
  // over the elements in the THiRA TList or an additional TObjArray.

  Int_t index = -99;

  //  return fTeleIndexMap[id];
  for(Int_t i=0; i<THIRA_NTELES; i++){
    if(teleID[i] == id) {index = i; break;}
  }
 if(index==-99){ printf("Telescope indexing error:  telescope ID=%i\n",id);}
 return index;
}


//______________________________________________________________________________
void THiRA::InitBranch(TTree *tree, Option_t *stage)
{
  // -- 
  // Initializes or creates all the branches and member branches for the THiRA
  // class on the tree.  This is typically called after instanciating TTree and
  // THiRA.  The branch structure will then be generated on the TTree.

  Char_t *stageOp = (Char_t*)stage;
  Char_t tmpVar1[100],tmpVar2[100];
  sprintf(tmpVar1,"%s.fHits[%i]",GetName(),THIRA_NMBS);
  sprintf(tmpVar2,"fHits[%i]/S",THIRA_NMBS);

  tree->Branch(fName+".fNwords",&fNwords,"fNwords/s");
  tree->Branch(fName+".fRawMult",&fRawMult,"fRawMult/S");
  tree->Branch(fName+".fHits[4]",&fHits,tmpVar2);
  tree->Branch(fName+".fPixel","TClonesArray",&fPixel,32000,99);

  tree->Branch(&teleTList,32000,3);

  #if THIRA_OP_USETIMES==0
  for(Int_t j=0; j<THIRA_NTELES; j++){
    TString nameT(tele[j].GetName());
    tree->SetBranchStatus(nameT+"DE.TRaw[32]",0);
    tree->SetBranchStatus(nameT+"EF.TRaw[32]",0);
    tree->SetBranchStatus(nameT+"EB.TRaw[32]",0);
    tree->SetBranchStatus(nameT+"DE.T[32]",0);
    tree->SetBranchStatus(nameT+"EF.T[32]",0);
    tree->SetBranchStatus(nameT+"EB.T[32]",0);
  }
  #endif

  if(strcmp(stageOp,"Raw")==0){
    printf("* Disabling HiRA calibrated and primary branches.\n");
    printf("  (These branches will be filled during the calibration and primary data generation stages.)\n");
    for(Int_t k=0; k<THIRA_NTELES; k++){
      TString nameD(tele[k].GetName());
      tree->SetBranchStatus(nameD+"DE.E[32]",0);
      tree->SetBranchStatus(nameD+"DE.T[32]",0);
      tree->SetBranchStatus(nameD+"EF.E[32]",0);
      tree->SetBranchStatus(nameD+"EF.T[32]",0);
      tree->SetBranchStatus(nameD+"EB.E[32]",0);
      tree->SetBranchStatus(nameD+"EB.T[32]",0);
      tree->SetBranchStatus(nameD+"CsI.E[4]",0);
      tree->SetBranchStatus(nameD+"CsI.E2[4]",0);

      //      tree->SetBranchStatus(fName+".fPixel*",0);
    }
  }else if(strcmp(stageOp,"Cal")==0){
    //    printf("* Disabling HiRA raw branches.\n");
    for(Int_t k=0; k<THIRA_NTELES; k++){
      TString nameD(tele[k].GetName());
      // Not working correctly.  Disabling after branches are filled.
//       tree->SetBranchStatus(nameD+"DE.ERaw[32]",0);
//       tree->SetBranchStatus(nameD+"DE.TRaw[32]",0);
//       tree->SetBranchStatus(nameD+"EF.ERaw[32]",0);
//       tree->SetBranchStatus(nameD+"EF.TRaw[32]",0);
//       tree->SetBranchStatus(nameD+"EB.ERaw[32]",0);
//       tree->SetBranchStatus(nameD+"EB.TRaw[32]",0);
//       tree->SetBranchStatus(nameD+"CsI.ERaw[4]",0);
    }
  }else if(strcmp(stageOp,"Pri")==0){
    printf("* Enabling HiRA primary branches.\n");
    for(Int_t k=0; k<THIRA_NTELES; k++){
      TString nameD(tele[k].GetName());
      tree->SetBranchStatus(nameD+"DE.ERaw[32]",0);      tree->SetBranchStatus(nameD+"DE.E[32]",0);
      tree->SetBranchStatus(nameD+"DE.TRaw[32]",0);      tree->SetBranchStatus(nameD+"DE.T[32]",0);
      tree->SetBranchStatus(nameD+"EF.ERaw[32]",0);      tree->SetBranchStatus(nameD+"EF.E[32]",0);
      tree->SetBranchStatus(nameD+"EF.TRaw[32]",0);      tree->SetBranchStatus(nameD+"EF.T[32]",0);
      tree->SetBranchStatus(nameD+"EB.ERaw[32]",0);      tree->SetBranchStatus(nameD+"EB.E[32]",0);
      tree->SetBranchStatus(nameD+"EB.TRaw[32]",0);      tree->SetBranchStatus(nameD+"EB.T[32]",0);
      tree->SetBranchStatus(nameD+"CsI.ERaw[4]",0);      tree->SetBranchStatus(nameD+"CsI.E[4]",0); tree->SetBranchStatus(nameD+"CsI.E2[4]",0);
    }
  }
}


//______________________________________________________________________________
Bool_t THiRA::InitClass()
{
  // --
  // Loads the chip/telescope mapping for the THiRA class as well as initializing 
  // class variables for THiRA. 
  
  for(Int_t i=0; i<THIRA_NTELES; i++) fTeleIndexMap[teleID[i]] = i;

  SetNameTitle("hira","HiRA Class");

  TString chain;

  // Load the ASIC chip mapping file.
  printf("* ASIC Chip Map Loading . . . ");
  if(asic.LoadMap()) printf("\t\t\t\t\t\t\t [OK]\n");
  else{              printf("\t\t\t\t\t\t\t [FAILED!!!]\n"); return kFALSE;}

  // Load the CsI ADC mapping file.
  printf("* CsI ADC Map Loading . . .   ");
  if(LoadAdcMap())   printf("\n");
  else{              printf("\t\t\t\t\t\t\t [FAILED!!!]\n"); return kFALSE;}

  // Initialize the individual HiRA telescopes.
  for(Int_t i=0; i<THIRA_NTELES; i++){
    chain = "tele";
    // This is used to set the telescope title.
    chain += teleID[i];                      
    chain += ".";
    tele[i].InitClass(this,chain,i,teleID[i]);
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++){
      tele[i].DE_Eslope[j]  = 1;
      tele[i].EF_Eslope[j]  = 1;
      tele[i].EB_Eslope[j]  = 1;
      tele[i].CsI_Eslope[j] = 1;
      tele[i].DE_Tslope[j]  = 1;
      tele[i].EF_Tslope[j]  = 1;
      tele[i].EB_Tslope[j]  = 1;

      tele[i].DE_Eoffset[j]  = 0;
      tele[i].EF_Eoffset[j]  = 0;
      tele[i].EB_Eoffset[j]  = 0;
      tele[i].CsI_Eoffset[j] = 0;
      tele[i].DE_Toffset[j]  = 0;
      tele[i].EF_Toffset[j]  = 0;
      tele[i].EB_Toffset[j]  = 0;
    }
    
    // Add each telescope to the TList.  This inherits from TObject and is 
    // the class that will be written to the tree.
    teleTList.Add(&tele[i]);
  }

  // Set the motherboard number to the Si detectors.
  Int_t id;
  Int_t face;
  for(Int_t i=0; i<THIRA_NMBS; i++){
    for(Int_t j=0; j<(Int_t)(asic.tower[i].slotList.size()); j++){
      id   = asic.tower[i].slot[asic.tower[i].slotList[j]].telenum;
      face = asic.tower[i].slot[asic.tower[i].slotList[j]].facenum;

      #ifdef THIRA_DEBUG
      printf("%i %i %i %i\n",i,j,id,face);
      #endif
      
      switch (face){
      case 0:
	tele[GetTeleIndex(id)].DE.fMb   = i;
	tele[GetTeleIndex(id)].DE.fSlot = asic.tower[i].slotList[j];
	break;
      case 1:
	tele[GetTeleIndex(id)].EF.fMb   = i;
	tele[GetTeleIndex(id)].EF.fSlot = asic.tower[i].slotList[j];
	break;
      case 2:
	tele[GetTeleIndex(id)].EB.fMb   = i;
	tele[GetTeleIndex(id)].EB.fSlot = asic.tower[i].slotList[j];
	break;
      default:
	printf("Invalid detector face %i\n",face);
	break;
      }
    }
  }

  // Set the CsI crystal postions.
  for(Int_t k=0; k<(Int_t)adcList.size(); k++){
    Int_t slot = adcList[k];
    for(Int_t chan=0; chan<32; chan++){
      for(Int_t ii=0; ii<4; ii++){
	tele[GetTeleIndex(adcTeleMap[slot][chan])].CsI.fPos[ii] = CsIMap[slot][chan];
      }
    }
  }

  fAddRndm        = kTRUE;
  fCorrectReadout = kTRUE;
  fNwords         = 0;

  SiRawThresh  = 0;
  CsIRawThresh = 0;
  fRawMult     = 0;

  // Set the reference frame for HiRA.
  fEulerPhi   = 0.00;
  fEulerTheta = 0.00;
  fEulerPsi   = 0.00;

  fgPixel = new TClonesArray("THiRAPixel",10);
  fPixel  = fgPixel;

  fCRange = new TCrange();

  return kTRUE;
}


//______________________________________________________________________________
void THiRA::InitTree(TTree *tree)
{
  // -- Initializes a tree containing THiRA branches.
  //

  TString bName(fName);
  Char_t name[1000];
  fChain   = tree;
  fCurrent = -1;
  //   if(strstr(fChain->GetName(),"Cal")!=0){
  //     bName = "E02023."+bName;
  //     fChain->GetTree()->GetFriend("E02023")->SetMakeClass(1);
  //   }

  if(!fChain->GetTree()) printf("Tree not loaded or does not exist!!! \n");
  else{
    TNamed *aState = (TNamed*)fChain->GetTree()->GetUserInfo()->FindObject("Analysis State");
    if(strcmp(aState->GetTitle(),"2")!=0){
      fChain->SetBranchAddress(bName+".fNwords",&fNwords,&b_fNwords);
      fChain->SetBranchAddress(bName+".fRawMult",&fRawMult,&b_fRawMult);
      fChain->SetBranchAddress(bName+".fPixel",&fPixel,&b_fPixel);
      sprintf(name,".fHits[%i]",THIRA_NMBS); fChain->SetBranchAddress(bName+name,&fHits,&b_fHits);
    
      // Initialize the telescopes with the tree.
      for (Int_t i=0; i<THIRA_NTELES; i++){
	sprintf(name,"tele%i.",teleID[i]);
	fTele[i] = &tele[i];
	fChain->SetBranchAddress(name,&fTele[i]);
	tele[i].InitTree(tree);
      }
    }
    // Primary Analysis Tree.  Pixel branches only.
    if(strcmp(aState->GetTitle(),"2")==0){
      fChain->SetBranchAddress("fPixel",&fPixel,&b_fPixel);
      //      sprintf("nHits[%i]",THIRA_NMBS); fChain->SetBranchAddress(bName+name,&fHits,&b_fHits);
    }
  }
}


//______________________________________________________________________________
Bool_t THiRA::LoadAdcMap()
{
  // -- Load the adc mapping file.
  //

  Int_t slot,chan,tele,csiNum;

  ifstream file(csiMapPath,ios::in);
  if(!file.is_open()) return kFALSE;

  Char_t line[200];

  //  file >> adc_fill;
  file.getline(line,200);
  file.getline(line,200);
  file.getline(line,200);

  file.getline(line,200);
  sscanf(line,"%i %i %i %i",&slot,&chan,&tele,&csiNum);
  adcTeleMap[slot][chan] = tele;
  CsIMap[slot][chan]     = csiNum;
 
  adcList.push_back(slot);
 
  while (!file.eof()){
    file.getline(line,200);
    sscanf(line,"%i %i %i %i",&slot,&chan,&tele,&csiNum);
    adcTeleMap[slot][chan] = tele;
    CsIMap[slot][chan]     = csiNum;

    if (slot != adcList[adcList.size()-1]) adcList.push_back(slot);
  }
  printf("There are %i ADC's for CsI\n",adcList.size());
  for (UInt_t x=0; x<adcList.size(); x++) printf("   slot # %i\t\t\t\t\t\t\t\t\t [OK]\n",adcList[x]);

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadBadElements(Char_t *filePath)
{
  // --
  //
  
  ifstream file(filePath,ios::in);
  if(!file.is_open()){
    printf("* HiRA Bad Elements file could not be opened!!!\t\t\t\t[!!]\n");
    return kFALSE;
  }
  Char_t line[1000];
  Int_t  tID;
  Char_t tType[20];

  file.getline(line,1000);
  while(line[0]=='#') file.getline(line,1000);

  Int_t bdE[THIRA_SI_CHANNELS];
  Int_t bEF[THIRA_SI_CHANNELS];
  Int_t bEB[THIRA_SI_CHANNELS];
  Int_t bCsI[THIRA_CSI_CHANNELS];

  // Reset all badstrips
  for(Int_t i=0; i<THIRA_NTELES; i++){
    tele[i].DE.fBadStrips.reset();
    tele[i].EF.fBadStrips.reset();
    tele[i].EB.fBadStrips.reset();
    tele[i].CsI.fBadSegs.reset();
  }

  while(!file.eof()){
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++) bdE[j]=99;
    sscanf(line,"%i %s %i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i",
	   &tID,tType,&bdE[0],&bdE[1],&bdE[2],&bdE[3],&bdE[4],&bdE[5],&bdE[6],&bdE[7],&bdE[8],&bdE[9],&bdE[10],&bdE[11],
	   &bdE[12],&bdE[13],&bdE[14],&bdE[15],&bdE[16],&bdE[17],&bdE[18],&bdE[19],&bdE[20],&bdE[21],&bdE[22],&bdE[23],
	   &bdE[24],&bdE[25],&bdE[26],&bdE[27],&bdE[28],&bdE[29],&bdE[30],&bdE[31]);
    for(Int_t i=0; i<THIRA_SI_CHANNELS; i++) tele[GetTeleIndex(tID)].DE.SetBad(bdE[i]);
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++) bEF[j]=99;
    file.getline(line,1000);
    sscanf(line,"%s %i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i",
	   tType,&bEF[0],&bEF[1],&bEF[2],&bEF[3],&bEF[4],&bEF[5],&bEF[6],&bEF[7],&bEF[8],&bEF[9],&bEF[10],&bEF[11],
	   &bEF[12],&bEF[13],&bEF[14],&bEF[15],&bEF[16],&bEF[17],&bEF[18],&bEF[19],&bEF[20],&bEF[21],&bEF[22],&bEF[23],
	   &bEF[24],&bEF[25],&bEF[26],&bEF[27],&bEF[28],&bEF[29],&bEF[30],&bEF[31]);
    for(Int_t i=0; i<THIRA_SI_CHANNELS; i++) tele[GetTeleIndex(tID)].EF.SetBad(bEF[i]);
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++) bEB[j]=99;
    file.getline(line,1000);
    sscanf(line,"%s %i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i,%i",
	   tType,&bEB[0],&bEB[1],&bEB[2],&bEB[3],&bEB[4],&bEB[5],&bEB[6],&bEB[7],&bEB[8],&bEB[9],&bEB[10],&bEB[11],
	   &bEB[12],&bEB[13],&bEB[14],&bEB[15],&bEB[16],&bEB[17],&bEB[18],&bEB[19],&bEB[20],&bEB[21],&bEB[22],&bEB[23],
	   &bEB[24],&bEB[25],&bEB[26],&bEB[27],&bEB[28],&bEB[29],&bEB[30],&bEB[31]);
    for(Int_t i=0; i<THIRA_SI_CHANNELS; i++) tele[GetTeleIndex(tID)].EB.SetBad(bEB[i]);
    for(Int_t j=0; j<THIRA_CSI_CHANNELS; j++) bCsI[j]=99;
    file.getline(line,1000);
    sscanf(line,"%s %i,%i,%i,%i",tType,&bCsI[0],&bCsI[1],&bCsI[2],&bCsI[3]);
    for(Int_t i=0; i<THIRA_CSI_CHANNELS; i++) tele[GetTeleIndex(tID)].CsI.SetBad(bCsI[i]);
    file.getline(line,1000);
    file.getline(line,1000);
  }
  file.close();

  // Print bad elements.
  PrintBad("short");

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadCsICalibration(Char_t *filePath,Int_t nCal)
{
  // --
  // Loads and parses the .vdef file containing the HiRA CsI slopes and offset
  // values.  The filepath argument is the location of the calibration file.


  ifstream file(filePath,ios::in);
  if(!file.is_open()){
    printf("* CsI calibration %i file could not be opened!!!\t\t\t\t[!!]\n",nCal); 
    return kFALSE;
  }
  Char_t   name[200];
  Char_t   line[200];
  string   name_str,tele_str,channel_str;
  Int_t    channel,telenum,teleIndex;
  Double_t value;

  file.getline(line,200);
  while(line[0]=='#') file.getline(line,200);

  while (!file.eof()){
    sscanf(line,"%*s %*s %s %lf %*s",name,&value);
    
    name_str    = name;
    tele_str    = name_str.substr(9,2);
    telenum     = atoi (tele_str.c_str());
    teleIndex   = GetTeleIndex(telenum);
    channel_str = name_str.substr((name_str.length()-1),2);
    channel     = atoi (channel_str.c_str());

    if (name_str.find("csi")!=string::npos){
      if(nCal==0){
	if (name_str.find("offset")!=string::npos) tele[teleIndex].CsI_Eoffset[channel]= value;
	if (name_str.find("slope") !=string::npos) tele[teleIndex].CsI_Eslope[channel] = value;
      }else if(nCal==1){
	if (name_str.find("offset")!=string::npos) tele[teleIndex].CsI_E2offset[channel]= value;
	if (name_str.find("slope") !=string::npos) tele[teleIndex].CsI_E2slope[channel] = value;
      }
    }
    file.getline(line,200);
  }
  file.close();
  printf("* CsI calibration %i file loaded.\t\t\t\t\t\t[OK]\n",nCal);

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadCsISegVec(Char_t *filePath)
{
  // -- Load the segment position vectors for all CsI crystals.

  ifstream file(filePath,ios::in);
  if(!file.is_open()){printf("* CsI segment position file could not be opened!!!\t\t\t\t [FAILED]\n"); return kFALSE;}
  
  Char_t line[500];
  Char_t tempS[200];
  Int_t  index;
  Int_t  id;

  file.getline(line,500);
  // Skip all comment lines.
  while(line[0] != 'E') file.getline(line,500);
  printf("* Loading CsI segment positions -- ");

  while(!file.eof()){
    Int_t row,col;
    Double_t x,y,z;

    file.getline(line,500);
    sscanf(line,"%s %s %i %s %s %i",tempS,tempS,&index,tempS,tempS,&id);
    if(GetTeleIndex(id) == -99){
      printf("Failed to load pixel matrix.  Telescope %i does not exist!!!\n",id);
      return kFALSE;
    }

    for(Int_t j=0; j<THIRA_SI_CHANNELS*THIRA_SI_CHANNELS; j++){
      file.getline(line,500);
      sscanf(line,"%u %u %lf %lf %lf",&row,&col,&x,&y,&z);
      tele[GetTeleIndex(id)].pixelMatrix[row][col].SetXYZ(x,y,z);
    }
  }
  file.close();
  printf("* CsI Segment Positions Loaded. \t\t\t\t\t [OK]\n");  

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadCsIVCalibration(Char_t *filePath)
{
  // --
  // Loads and parses the .vdef file containing the HiRA CsI voltage calibration
  // values.  The filepath argument is the location of the calibration file.


  ifstream file(filePath,ios::in);
  if(!file.is_open()){
    printf("* CsI voltage calibration file could not be opened!!!\t\t\t[!!]\n"); 
    return kFALSE;
  }
  Char_t   name[200];
  Char_t   line[200];
  string   name_str,tele_str,channel_str;
  Int_t    channel,telenum,teleIndex;
  Double_t value;

  file.getline(line,200);
  while(line[0]=='#') file.getline(line,200);

  while (!file.eof()){
    sscanf(line,"%*s %*s %s %lf %*s",name,&value);
    
    name_str    = name;
    tele_str    = name_str.substr(9,2);
    telenum     = atoi (tele_str.c_str());
    teleIndex   = GetTeleIndex(telenum);
    channel_str = name_str.substr((name_str.length()-1),2);
    channel     = atoi (channel_str.c_str());

    //    cout << tele_str << " " << telenum << " " << teleIndex << " " << channel_str << " " << channel << " " << value << endl; 
    if (name_str.find("csi")!=string::npos){
      if (name_str.find("A0")!=string::npos) tele[teleIndex].CsI_Vcoeff[channel][0]= value;
      if (name_str.find("A1")!=string::npos) tele[teleIndex].CsI_Vcoeff[channel][1]= value;
      if (name_str.find("A2")!=string::npos) tele[teleIndex].CsI_Vcoeff[channel][2]= value;
      if (name_str.find("A3")!=string::npos) tele[teleIndex].CsI_Vcoeff[channel][3]= value;
      if (name_str.find("A4")!=string::npos) tele[teleIndex].CsI_Vcoeff[channel][4]= value;
    }
    file.getline(line,200);
  }
  file.close();
  printf("* CsI voltage calibration file loaded.\t\t\t[OK]\n");

  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadPixelMatrix(Char_t *filePath)
{
  // -- Load a HiRA pixel file.
  // This file should contain the DESIGN VECTORS of each pixel in the array.
  // A script should be used later for calculations using the MEASURED VECTORS.

  ifstream file(filePath,ios::in);
  if(!file.is_open()){printf("Pixel file could not be opened!!!\n"); return kFALSE;}
  
  Char_t line[500];
  Char_t tempS[200];
  Int_t  index;
  Int_t  id;

  file.getline(line,500);
  // Skip all comment lines.
  while(line[0] != 'E') file.getline(line,500);
  printf("* Loading pixel matrix -- ");

  while(!file.eof()){
    Int_t row,col;
    Double_t x,y,z;

    file.getline(line,500);
    sscanf(line,"%s %s %i %s %s %i",tempS,tempS,&index,tempS,tempS,&id);
    if(GetTeleIndex(id) == -99){
      printf("Failed to load pixel matrix.  Telescope %i does not exist!!!\n",id);
      return kFALSE;
    }

    for(Int_t j=0; j<THIRA_SI_CHANNELS*THIRA_SI_CHANNELS; j++){
      file.getline(line,500);
      sscanf(line,"%u %u %lf %lf %lf",&row,&col,&x,&y,&z);
      tele[GetTeleIndex(id)].pixelMatrix[row][col].SetXYZ(x,y,z);
    }
  }
  file.close();
  
  printf("Pixel Matrix Loaded. \t\t\t\t\t [OK]\n");
  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadReadoutShifts(Char_t *filePath)
{
  // -- Load the HiRA readout order shift parameters.
  //
  // This file contains the readout shifts for each channel in HiRA.
  // 
  // The telescopes are ordered starting at the bottom left, increasing as we 
  // move up a tower.  The number continues at the bottom of the next tower 
  // proceeding until the upper right of the array.
  //
  // This file and the LoadReadoutShift() method should be modified if a 
  // more complicated shift id needed.  This file format should satisfy most 
  // basic needs.
  //
  // The shift can be set per channel, per face, or per chip.  It is assumed 
  // that the highest order if O(5).
  //
  // *** SIMPLE PER CHANNEL SHIFT ***
  // The format is IF Complexity==2,
  //    teleIndex   teleID   face   Complexity
  // 	O(0)	O(1) 	O(2) 	O(3) 	O(4) 	O(5) 	//CH0
  //	O(0) 	O(1) 	O(2) 	O(3) 	O(4) 	O(5)    //CH1
  // 	O(0)	O(1) 	O(2) 	O(3) 	O(4) 	O(5) 	//CH2
  //	O(0) 	O(1) 	O(2) 	O(3) 	O(4) 	O(5)    //CH3
  //	.	.	.	.	.	.	
  //	.	.	.	.	.	.
  //	.	.	.	.	.	.
  //	
  //       and so on . . .
  //
  // *** SIMPLE PER CHIP SHIFT ***
  // The format is IF Complexity==1, (EVEN->Chip 0,ODD->Chip 1)
  //    teleIndex   teleID   face   Complexity
  // 	O(0)	O(1) 	O(2) 	O(3) 	O(4) 	O(5) 	//CHIP0
  //	O(0) 	O(1) 	O(2) 	O(3) 	O(4) 	O(5)    //CHIP1
  //
  //       and so on . . .
  //
  // *** SIMPLE PER FACE SHIFT ***
  // The format is IF Complexity==0,
  //    teleIndex   teleID   face   Complexity
  // 	O(0)	O(1) 	O(2) 	O(3) 	O(4) 	O(5) 	//FACE
  //
  // 	and so on . . . 
  //

  ifstream file(filePath,ios::in);
  if(!file.is_open()) {printf("Readout order shift file could not be opened!!!\n"); return kFALSE;}
  Char_t   line[2000],faceC[200];
  Int_t    teleIndex,teleID,opt,order;
  Int_t    face = -1;
  Double_t valueArray[5];
  memset(valueArray,'\0',sizeof(valueArray));

  // Skip commented lines.
  file.getline(line,2000);
  while(!isdigit(line[0])) file.getline(line,2000);    

  while(!file.eof() || isdigit(line[0])){
    sscanf(line,"%i %i %s %i",&teleIndex,&teleID,faceC,&opt);
    //    printf("%i %i %s %i\n",teleIndex,teleID,faceC,opt);
    if(strcmp(faceC,"de")==0) face = 0;
    if(strcmp(faceC,"ef")==0) face = 1;
    if(strcmp(faceC,"eb")==0) face = 2;
    if(opt==0){
      file.getline(line,2000);
      memset(valueArray,'\0',sizeof(valueArray));
      switch (face){
      case 0:  
	sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	for(Int_t j=0; j<2; j++){
	for(Int_t k=0; k<5; k++){  tele[teleIndex].DE_Shift[i][j][k]=valueArray[k];}}}
	break;
      case 1:  
	sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	for(Int_t j=0; j<2; j++){
	for(Int_t k=0; k<5; k++){  tele[teleIndex].EF_Shift[i][j][k]=valueArray[k];}}}
	break;
      case 2:
	sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	for(Int_t j=0; j<2; j++){
	for(Int_t k=0; k<5; k++){  tele[teleIndex].EB_Shift[i][j][k]=valueArray[k];}}}
	break;
      }
    }else if(opt==1){
      for(Int_t m=0; m<2; m++){
	file.getline(line,2000);
	memset(valueArray,'\0',sizeof(valueArray));
	switch (face){
	case 0:  
	  sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	  for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	  for(Int_t k=0; k<5; k++){  tele[teleIndex].DE_Shift[i][m][k]=valueArray[k];}}
	  break;
	case 1:  
	  order=sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	  //printf("%lf \t%lf \t%lf \t%lf \t%lf \t%i\n",valueArray[0],valueArray[1],valueArray[2],valueArray[3],valueArray[4],order);
	  for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	  for(Int_t k=0; k<5; k++){  tele[teleIndex].EF_Shift[i][m][k]=valueArray[k];}}
	  break;
	case 2:  
	  order=sscanf(line,"%lf %lf %lf %lf %lf",&valueArray[0],&valueArray[1],&valueArray[2],&valueArray[3],&valueArray[4]);
	  //printf("%lf \t%lf \t%lf \t%lf \t%lf \t%i\n",valueArray[0],valueArray[1],valueArray[2],valueArray[3],valueArray[4],order);
	  for(Int_t i=0; i<THIRA_SI_CHANNELS; i++){
	  for(Int_t k=0; k<5; k++){  tele[teleIndex].EB_Shift[i][m][k]=valueArray[k];}}
	  break;
	}
      }
    }else if(opt==2){
      printf("* NOT IMPLEMENTED YET.\n");
    }else{printf("Error reading file!!!\n");}
    file.getline(line,2000);
  }
  printf("* Loaded HiRA readout order correction parameters. \t\t\t\t[OK]\n");
  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::LoadSiCalibration(Char_t *filePath)
{
  // -- Load a HiRA Si calibration file.
  // Loads and parses the .vdef file containing the HiRA Si slopes and offset
  // values.  The filepath argument is the location of the calibration file.


  ifstream file(filePath,ios::in);
  if(!file.is_open()){
    printf("* Si calibration file could not be opened!!!\t\t\t\t[!!]\n"); 
    return kFALSE;
  }
  Char_t   name[200];
  Char_t   line[200];
  string   name_str,tele_str,channel_str;
  Int_t    channel,telenum,teleIndex;
  Double_t value;

  while (!file.eof()){
    file.getline(line,200);
    sscanf(line,"%*s %*s %s %lf %*s",name,&value);
    
    name_str    = name;
    tele_str    = name_str.substr(9,2);
    telenum     = atoi (tele_str.c_str());
    teleIndex   = GetTeleIndex(telenum);
    channel_str = name_str.substr((name_str.length()-2),2);
    channel     = atoi (channel_str.c_str());

    if (name_str.find("ef")!=string::npos){
      if (name_str.find("offset")!=string::npos) tele[teleIndex].EF_Eoffset[channel] = value;
      if (name_str.find("slope") !=string::npos) tele[teleIndex].EF_Eslope[channel]  = value;
    }
    if (name_str.find("eb")!=string::npos){
      if (name_str.find("offset")!=string::npos) tele[teleIndex].EB_Eoffset[channel] = value;
      if (name_str.find("slope") !=string::npos) tele[teleIndex].EB_Eslope[channel]  = value;
    }
    if (name_str.find("de")!=string::npos){
      if (name_str.find("offset")!=string::npos) tele[teleIndex].DE_Eoffset[channel] = value;
      if (name_str.find("slope") !=string::npos) tele[teleIndex].DE_Eslope[channel]  = value;
    }
  }

  file.close();
  printf("* Si calibration file loaded.\t\t\t\t\t\t\t[OK]\n");

  return kTRUE;
}


//______________________________________________________________________________
void THiRA::ResetCalibBranches()
{
  // --
  //
  
  for(Int_t i=0; i<THIRA_NTELES; i++){
    tele[i].DE.b_E->Reset();
    tele[i].DE.b_T->Reset();
    tele[i].EF.b_E->Reset();
    tele[i].EF.b_T->Reset();
    tele[i].EB.b_E->Reset();
    tele[i].EB.b_T->Reset();
    tele[i].CsI.b_E->Reset();
  }
}


//______________________________________________________________________________
void THiRA::PrintBad(Option_t* format)
{
  // -- Print all set bad elements to the stdout.
  //

  printf("**************************************************************************\n");
  printf("*Tele :\tFACE : \tBad strips or segments \t\t\t\t\t*\n");
  printf("**************************************************************************\n");
  if(strcmp(format,"short")==0){
  for(Int_t i=0; i<THIRA_NTELES; i++){
    printf("* %2i\t dE: ",teleID[i]);
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].DE.IsBad(j)) printf("%2i,",j);
    printf("\t EF: ");
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].EF.IsBad(j)) printf("%2i,",j);
    printf("\t EB: ");
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].EB.IsBad(j)) printf("%2i,",j);
    printf("\t CsI: ");
    for(Int_t j=0; j<THIRA_CSI_CHANNELS; j++) if(tele[i].CsI.IsBad(j)) printf("%2i,",j);
    printf("\n");
  }
  printf("**************************************************************************\n");
  }else{
  for(Int_t i=0; i<THIRA_NTELES; i++){
    printf("* %i\t dE  :\t",teleID[i]);
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].DE.IsBad(j)) printf("%i,",j);
    printf("\n");
    printf("* \t EF  :\t");
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].EF.IsBad(j)) printf("%i,",j);
    printf("\n");
    printf("* \t EB  :\t");
    for(Int_t j=0; j<THIRA_SI_CHANNELS; j++)  if(tele[i].EB.IsBad(j)) printf("%i,",j);
    printf("\n");
    printf("* \t CsI :\t");
    for(Int_t j=0; j<THIRA_CSI_CHANNELS; j++) if(tele[i].CsI.IsBad(j)) printf("%i,",j);
    printf("\n");
    printf("**************************************************************************\n");
  }
  }

}


//______________________________________________________________________________
Bool_t THiRA::SetChipMapFile(Char_t *filePath)
{
  // -- Set the path of the ASIC chip mapping file.
  //

  Bool_t fFound = kFALSE;

  asic.chipMapPath = filePath;
  chipMapFile      = asic.chipMapPath;

  fFound = kTRUE;

  return fFound;
}


//______________________________________________________________________________
Bool_t THiRA::SetCsIMapFile(Char_t *filePath)
{
  // -- Set the CsI mapping file.
  //

  Bool_t fFound = kFALSE;

  csiMapPath = filePath;
  fFound = kTRUE;

  return fFound;
}


//______________________________________________________________________________
Bool_t THiRA::SetReadoutShiftPath(Char_t *filePath)
{
  // --

  fReadoutShiftPath = filePath;
  return kTRUE;
}


//______________________________________________________________________________
void THiRA::SetName(const Char_t *name)
{
  // -- Change (i.e. set) the name of the THiRA.
  //

  fName = name;
}


//______________________________________________________________________________
void THiRA::SetNameTitle(const Char_t *name, const Char_t *title)
{
  // -- Change (i.e. set) all the THiRA parameters (name and title).
  //

  fName  = name;
  fTitle = title;
}


//______________________________________________________________________________
void THiRA::SetTitle(const Char_t *title)
{
  // -- Change (i.e. set) the title of the THiRA.
  //

  fTitle = title;
}


//______________________________________________________________________________
Bool_t THiRA::Unpack(UShort_t *pEvent)
{
  // -- Unpack a HiRA event.  
  // This functions unpacks the raw data contained in a HiRA packet within 
  // an event data buffer.  The argument is a pointer which is the location of 
  // beginning of the packet.
  //
  // When unpacking data from an event file this function will typically be 
  // called whenever a HiRA packet tag is found.  (Packet id's are defined in 
  // the file packetID.h).
  //
  // There has been some confusion regarding the mapping of the strips to their 
  // respective chip and channel that is written to the event file.
  // 
  // Typically, 
  // 1.  The dE and EF strips are vertical.  The EB strips are horizontal.
  // 2.  Chip 0 is the first chip readout.
  // 3.  The chipboards are readout from the bottom slot on a MB to the top.
  // 4.  The convention is to have the upper left pixel of the detector to be (0,0) 
  //     so that the detector pixels are elements of a matrix.
  // 5.  The dE detector is connected in the same way as the EF.  HOWEVER, when 
  //     going through external preamps and the pinout is changed!!!  Therefore dE
  //     strips may map differently than below.
  //
  // Chip 0 is PIN #0 on chipboard (even PINS).
  // Chip 1 is PIN #1 on chipboard (odd  PINS).
  //
  // Ch          15 15 14 14 13 13 12 12 11 11 10 10 09 09 08 08 07 07 06 06 05 05 04 04 03 03 02 02 01 01 00 00
  //
  //    PIN      31 30 29 28 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 09 08 07 06 05 04 03 02 01 00
  //
  //       strip 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
  //             __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ __ 
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 00  00  00 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 00  01  01 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 01  02  02 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 01  03  03 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 02  04  04 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 02  05  05 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 03  06  06 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 03  07  07 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 04  08  08 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 04  09  09 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 05  10  10 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 05  11  11 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 06  12  12 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|  ^
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 06  13  13 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|  |
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 07  14  14 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|  EB
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 07  15  15 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|  |
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 08  16  16 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|  V
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 08  17  17 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 09  18  18 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 09  19  19 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 10  20  20 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 10  21  21 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 11  22  22 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 11  23  23 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 12  24  24 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 12  25  25 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 13  26  26 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 13  27  27 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 14  28  28 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 14  29  29 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 15  30  30 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //            |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
  // 15  31  31 |__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|__|
  //       
  //                                                <--- EF --->   
  //
  //

 
  UShort_t *p = pEvent;
  UShort_t packetSize;                 // Number of words in the HiRA packet.

  UShort_t subEventLength,
    subPacketSize,                     // Number of words in a Tower packet.
    XLM0[2000];

  Short_t stripnum = 99,
    towerNum       = 99,
    teleNum        = 99,
    faceNum        = 99;
  Long64_t nCBErrors = 0;                // Number of CB errors encountered.

  // Set sequence numbers for readout order problem.
  Short_t   cbRead[16];                  // If the CB has been read then set the slot element to 1 (true).
  Short_t   cbSeq[THIRA_NMBS];           // Counts the number of CB read in a tower.  Elments are incremented.
  Short_t   chSeq[THIRA_NMBS];           // Counts the number of channels read in a tower.  Elements are incremented.
  Short_t   chSeqPrior;                  // Number of channels read prior to the current CB.
  UInt_t    hitSeq[16];                  // Represents channels in a CB that have fired.  Bits are flipped to 1 if fired.
  memset(cbRead,'\0',sizeof(cbRead));  
  memset(cbSeq,'\0',sizeof(cbSeq));
  memset(chSeq,'\0',sizeof(chSeq));
  chSeqPrior = 0;
  memset(hitSeq,'\0',sizeof(hitSeq));


  *p--;
  *p--;
  packetSize = *p;                   // Number of words in the HiRA packet.
  packetSize = packetSize - 2;
  *p++;                              // Advance the event pointer to beginning of packet.
  *p++;
  fNwords = packetSize;

  while (packetSize > 0) {
    subPacketSize          = *p++;
    UShort_t subPacketTag  = *p++;     // Tower packet tag.
    Bool_t   towerFlag     = kFALSE;
    
    switch(subPacketTag) {
    case TOWER_0:
      towerNum = 0;
      break;
    case TOWER_1:
      towerNum = 1;
      break;
    case TOWER_2:
      towerNum = 2;
      break;
    case TOWER_3:
      towerNum = 3;
      break;
    case TOWER_4:
      towerNum = 4;
      break;
    case TOWER_5:
      towerNum = 5;
      break;
    }
    // Reset the cbRead array of flags because we are now on a new tower packet.
    memset(cbRead,'\0',sizeof(cbRead));
    
    UShort_t nEvents = *p++;                          // Number of events in the HiRA packet.
    if(nEvents < 1000 && nEvents > 0){                // Check for a physically reasonable number of events. 
      UShort_t XLMlength = *p++&0xfff;                // Number of words from the XLM.
      UShort_t XLMID     = *p++;                      // XLM ID number.
      
      for(Int_t j=0; j<nEvents; j++) XLM0[j] = *p++;
      
      UShort_t FADClength = *p++&0xfff;               // Number of words from the FADC.
      UShort_t FADCID     = *p++;	              // FADC ID number.
      
      if((FADClength-2)==(XLMlength-2)*2){            // Check that the number of words match.

	for(Int_t j=0; j<nEvents; j++){ 
	  //	  Short_t chipNumMB  = (Short_t)((XLM0[j]&0x0fe0)>>5);  // Chip number output from MB.  1-32
	  //	  Short_t chanNumMB  = (Short_t)(XLM0[j]&0x0f);         // Channel number output from MB. 0-15
	  // Channel number is lower 5 bits.
	  Short_t chipNumMB  = (Short_t)((XLM0[j]&0x1fe0)>>5);  // Chip number output from MB.  1-32
	  Short_t chanNumMB  = (Short_t)(XLM0[j]&0x1f);         // Channel number output from MB. 0-15
	  
	  // Check for any errors or mismatches first.
	  if (chipNumMB<1 || chipNumMB>=(Short_t)asic.tower[towerNum].slotList.size()*2+1) {
	    nCBErrors++;
	    printf("\t\t\t\t\t\tChipNumMB on Tower %hi=%hi.  Skipping . . .  Total=%lli\r",towerNum,chipNumMB,nCBErrors);
	    cout.flush();
	    fChipErr[towerNum][chipNumMB]++;
	    p++;
	    p++;
	  }else if(chanNumMB>15 || chanNumMB<0){
	    nCBErrors++;
	    printf("\t\t\t\t\t\tChanNumMB on Tower %hi=%hi.  Skipping . . .  Total=%lli\r",towerNum,chanNumMB,nCBErrors);
	    cout.flush();
	    fChanErr[towerNum][chipNumMB][chanNumMB]++;
	    p++; 
	    p++;
	  }else {
	    Short_t slotNum = asic.tower[towerNum].slotNum[chipNumMB]; // Slot number assigned to chipNumMB.
	    Short_t chipNum = asic.tower[towerNum].chipNum[chipNumMB]; // Chip number on a board.  Either 0 or 1.
	    // Chip 0 is the first read out and corresponds to PIN 0 on the connector which is the outer most chip.
	    // Channels are readout starting with 0.
	    
	    #ifdef THIRA_DEBUG_READORDER
	    printf("%x#08 ",hitSeq[slotNum]);
	    #endif
	    
	    //
	    fHits[towerNum]++;

	    // Increment the number of channels fired in the tower.
	    chSeq[towerNum]++;
	    // If the cb has not been readout then increment the number of CB's read (cbSeq array), set the 
	    // number of prior channels fired, and flag the CB as being readout (cbRead = 1 (true)).
	    if(cbRead[slotNum-1]==0){cbSeq[towerNum]++; chSeqPrior=chSeq[towerNum]; cbRead[slotNum-1] = 1;}
	    // For this slot and channel flip the channel bit from 0 to 1 (hitSeq array). 
	    hitSeq[slotNum] = hitSeq[slotNum]|(0x1<<(chanNumMB+chipNum*16));

	    #ifdef THIRA_DEBUG_READORDER
 	    printf("%x#08 %04i %x#08 %02i %02i %02i %02i %02i\n",cbSeq[towerNum],chSeqPrior,hitSeq[slotNum],
		   chanNumMB,chipNum,slotNum,chipNumMB,chanNumMB);
	    #endif

	    Short_t energy = (Short_t)*p++;
	    Short_t time   = (Short_t)*p++;
	    
	    if(energy == 32767) energy = -2;
	    if(time   == 32767) time   = -2;
	    
	    // Do HiRA_filling here.
	    if (slotNum == -1) {
	      printf("There is no slot defined for Tower %hi Chip %hi.  Skipping . . .\n",towerNum,chipNum);
	    }else {
	      teleNum   = asic.tower[towerNum].slot[slotNum].telenum; // 
	      faceNum   = asic.tower[towerNum].slot[slotNum].facenum; // 0=dE, 1=EF, 2=EB.
	      if(teleNum == -1 || faceNum == -1){
		printf("Unknown teleNum=%hi or faceNum=%hi on Tower=%hi Slot=%hi.  Skipping . . .\n",
		       teleNum,faceNum,towerNum,slotNum);
	      }else {
		// Decode the strip number depending on the Si face.
		if (faceNum==0) {
		  if(asic.tower[towerNum].chipNum[chipNumMB]==0) stripnum = chanNumMB;
		  else                                           stripnum = 16+chanNumMB;
		}else if(faceNum==1){
		  if(asic.tower[towerNum].chipNum[chipNumMB]==0) stripnum = 2*(15-chanNumMB)+1;
		  else                                           stripnum = 2*(15-chanNumMB);
		}else if(faceNum==2){
		  if(asic.tower[towerNum].chipNum[chipNumMB]==0) stripnum = 2*chanNumMB;
		  else                                           stripnum = 2*chanNumMB+1;
		  // Flip spectrum to low channels.
		  if (energy>0) {
		    energy -= (16384 - 1);
		    energy *= -1; 
		  }
		}
		
		// If the strip number is out of range output an error.
		if (stripnum < 0 || stripnum > 31) {
		  printf("Strip number is %hi.  Skipping . . .\n",stripnum);
		  faceNum = -1;  // set to avoid filling on these events!
		}
	      }
	    }
	    // Fill the mapped THiRA class.
	    
            #ifdef EXPERIMENT_05133
	    //  **********************
	    // TEMP FIX FOR 05133!!! MUST ADD INPUT FILE CHECKS!!!
	    //  **********************
	    if(teleNum==18 || teleNum==7  || teleNum==8 || teleNum==9 || teleNum==20 || teleNum==21 || teleNum==22 || teleNum==23){
	    }else{
	      Fill(energy,time,GetTeleIndex(teleNum),faceNum,stripnum,cbSeq[towerNum],chSeq[towerNum],hitSeq[slotNum]);
	    }
            #else
	    Fill(energy,time,GetTeleIndex(teleNum),faceNum,stripnum,cbSeq[towerNum],chSeq[towerNum],hitSeq[slotNum]);
            #endif
	  } 
	}
      }else {
	printf( "The FADC and the XLM are mismatched with FADC length=%i and XLM length=%i.  Skipping Tower=%hi. . .\n",
		FADClength,XLMlength,towerNum);
	p+=(FADClength-2);
      }
    }else printf("Bad number of events=%i ",nEvents);
    
    //    cout << "END TOWER PACKET" << endl;
    packetSize -= subPacketSize;
  }
  //  cout << "END HIRA PACKET" << endl;
  return kTRUE;
}


//______________________________________________________________________________
Bool_t THiRA::UnpackCsI(UShort_t *pEvent)
{
  // -- Unpack raw CsI data.
  //

  UShort_t *p = pEvent;
  UShort_t packetSize;                        // Number of words in a CsI packet.

  UShort_t subEventLength,
    subPacketSize;                            // Number of words in a sub-packet.


  *p--;
  *p--;
  packetSize = *p;                   // Number of words in the CsI packet.
  packetSize = packetSize - 2;
  *p++;                              // Advance the event pointer to beginning of packet.
  *p++;


  #ifdef EXPERIMENT_02023
  Bool_t   subFlag      = kFALSE;
  
  UShort_t		ID = *p++;
  
  if (packetSize >3) p = UnpackCsIEnergy(p);
  packetSize -= subPacketSize;
    
  #else
  while (packetSize > 0) {
    subPacketSize         = *p++;
    UShort_t subPacketTag = *p++;
    UShort_t ID; 
    
    switch(subPacketTag) {
    case VMEADC:
       ID = *p++;
       if (packetSize >3) p = UnpackCsIEnergy(p);
       break;
    }
    packetSize -= subPacketSize;
  }
  #endif

  return kTRUE;
}


//______________________________________________________________________________
UShort_t* THiRA::UnpackCsIEnergy(UShort_t *p)
{
  // -- Unpack the HiRA CsI energies.
  //

  Short_t teleNum,csiNum;
  UShort_t m_nCrate,m_nSlot,num_evts=0;
  Short_t channel;
  Short_t value;

  m_nCrate = (*p&0x0700)>>8;
  m_nSlot  = ((*p++)&0xf800)>>11;
  num_evts = ((*p++)&0x3f00)>>8;

  for (Int_t i=0; i<num_evts; i++){
     channel = (Short_t)((*p++)&(0x003f));
     value   = (Short_t)((*p++)&0x0fff);

     #ifdef EXPERIMENT_05133
     if(m_nSlot!=6){
       teleNum = adcTeleMap[m_nSlot][channel];
       csiNum  = CsIMap[m_nSlot][channel];
     }
     #else
     teleNum = adcTeleMap[m_nSlot][channel];
     csiNum  = CsIMap[m_nSlot][channel];
     #endif

     #ifdef EXPERIMENT_05133
     //  **********************
     // TEMP FIX FOR 05133!!! MUST ADD INPUT FILE CHECKS!!!
     //  **********************
     if(teleNum==18 || teleNum==7  || teleNum==8 || teleNum==9 || teleNum==20 || teleNum==21 || teleNum==22 || teleNum==23){
     }else{
       if(teleNum!=-1) FillCsI(value,GetTeleIndex(teleNum),csiNum);
     }
     #else
     if(teleNum!=-1) FillCsI(value,GetTeleIndex(teleNum),csiNum);
     #endif
  }

  p++;
  p++;

  return p;
}