// Author: Andrew M. Rogers and Mark Wallace, NSCL 07/01/2006
//* Copyright (C) 2006-2008 Andrew M. Rogers and Mark Wallace

#include "TMcp.h"

#include <fstream>
#include <math.h>
#include <stdio.h>

#include <TChain.h>
#include <TFile.h>
#include <TMath.h>
#include <TObject.h>
#include <TRandom3.h>
#include <TStopwatch.h>
#include <TTree.h>

ClassImp(TMcp);

using namespace TMath;

////////////////////////////////////////////////////////////////////////////////
/* BEGIN_HTML
<p> The TMcp class inherits from TObject.  This Class provides the structure 
for the MCP (MicroChannel Plate) tracking system.  The tracking system is composed of 
2 MCP devices descibed by the TMcpDet class.</p>  

<p>In a fully functional setup these detectors provide, after proper calibration, x-y 
position information produced from secondary electron emission.  With both MCP 
devices an incoming angle can be calculated yeilding full tracking of incoming
beam particles.  Additionally, fast timing signals from the detectors allow for 
additional and in many cases superior ToF information. 
</p>

<p>During certain experiments additional gain stages have been employed.  To account 
for these additional QDC channels there are 2 global variables declared in the 
rootDefine.h file which specify the number of stages.  Mapping of QDC (high/low/...) 
channels to corner branches is currently only done for high and low gains.
</p>

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

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

<h3><a name="proc">III. Basic Analysis Proceedure</a></h3>

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

<h3>Related References</h3>

<h4>[1] D Shapira, TA Lewis, LD Hulett - Nuclear Inst. and Methods in Physics Research, A, 
2000 - Elsevier.  A fast and accurate position-sensitive timing detector based on secondary 
electron emission <a href="http://dx.doi.org/10.1016/S0168-9002(00)00499-X">[Full Text]</a></h4>

<h4>[2] D Shapira, TA Lewis, LD Hulett, Z Ciao - Nuclear Inst. and Methods in Physics Research, A, 
2000 - Elsevier.  Factors afecting the performance of detectors that use secondary electron 
emission from a thin foil to determine ion impact position 
<a href="http://dx.doi.org/10.1016/S0168-9002(99)01444-8">[Full Text]</a></h4>

<h4>[3] G. Montagnolia, A.M. Stefaninib, M. Trottac, S. Beghinia, M. Bettinia, F. Scarlassaraa, V. Schiavona, L. Corradib, B.R. Beherab, E. Fiorettob, A. Gadeab, A. Latinab, S. Szilnerb, d, L. Donàb, M. Rigatob, N.A. Kondratieve, A. Yu. Chizhove, G. Kniajevae, E.M. Kozuline, I.V. Pokrovskiye, V.M. Voskressenskye and D. Ackermannf 
The large-area micro-channel plate entrance detector of the heavy-ion magnetic spectrometer PRISMA <a href="http://dx.doi.org/10.1016/j.nima.2005.03.158">[Full Text]</a></h4>


END_HTML */
////////////////////////////////////////////////////////////////////////////////


//______________________________________________________________________________
TMcp::TMcp(const TMcp &mcp) : TObject(mcp)
{
  // -- Copy constructor.

  ((TMcp&)mcp).Copy(*this);
}


//______________________________________________________________________________
void TMcp::Copy(TObject &mcp)const
{
  // -- Copy this method.
  //

  ((TMcp&)mcp).fRandom.SetSeed(0);
  ((TMcp&)mcp).fName      = fName;
  ((TMcp&)mcp).fTitle     = fTitle;
  ((TMcp&)mcp).mcpMapFile = mcpMapFile ;
  ((TMcp&)mcp).fAddRndm   = fAddRndm;
  for(Int_t i=0; i<TMCP_NQDCS; i++) ((TMcp&)mcp).fQdcSlot[i] = fQdcSlot[i];
  for(Int_t i=0; i<4; i++) ((TMcp&)mcp).tShift[i]   = tShift[i];
  for(Int_t i=0; i<2; i++) ((TMcp&)mcp).thetaMCP[i] = thetaMCP[i];
  ((TMcp&)mcp).tarMcpPlasticGap = tarMcpPlasticGap;
  ((TMcp&)mcp).mcpMcpGap        = mcpMcpGap;
  ((TMcp&)mcp).upMcpPlasticGap  = upMcpPlasticGap;
  ((TMcp&)mcp).maskTarGap       = maskTarGap;

  mcp0.Copy(((TMcp&)mcp).mcp0);
  mcp1.Copy(((TMcp&)mcp).mcp1);
  ((TMcp&)mcp).mcp0.parent = &((TMcp&)mcp);
  ((TMcp&)mcp).mcp1.parent = &((TMcp&)mcp);
}


//______________________________________________________________________________
void TMcp::Clear(Option_t *)
{
  // -- Clear the MCP tracking system data members.
  // 

  memset(ene,'\0',sizeof(ene));
  memset(time,'\0',sizeof(time));
  memset(charge,'\0',sizeof(charge));
  for(Int_t i=0; i<32; i++){
    rTime[i]    = sqrt(-1.0);
    timeCorr[i] = sqrt(-1.0);
  }

  maskX  = sqrt(-1.0);
  maskY  = sqrt(-1.0);
  xAngle = sqrt(-1.0);
  yAngle = sqrt(-1.0);
  iAngle = sqrt(-1.0);

  mcp0.Clear();
  mcp1.Clear();
}


//______________________________________________________________________________
Int_t TMcp::Calibrate(Option_t* prefix,Option_t* path,Option_t* treeName,Option_t* option)
{
  // -- Calibrates TMcp components.
  //

  // Check that fChain has been set.
  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.
  TMcp *newMcp = new TMcp(*this);
  newMcp->mcp0.fCon = this;
  newMcp->mcp1.fCon = this;

  // 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.
    newTree->Branch("mcp",&newMcp);    
    newMcp->InitTree(newTree);

    // Calibrate
    Long64_t nEntries = cTree->GetEntries();
    for(Long64_t ientry=0; ientry<nEntries; ientry++){
      // Update current progress.
      Double_t prog = ((Double_t)ientry/nEntries)*100;
      if(ientry%100==0){printf("Entry: %15lli \t%4.2lf %%\r",ientry,prog); cout.flush();}

      newMcp->Clear();
      newMcp->mcp0.Calibrate(ientry);
      newMcp->mcp1.Calibrate(ientry);
      newMcp->mcp0.Calculate(ientry);
      newMcp->mcp1.Calculate(ientry);

      // Calculate relative times.
      b_time->GetEntry(ientry);
      newMcp->rTime[0] = time[2]-time[1]+tShift[0];
      newMcp->rTime[1] = time[2]-time[0]+tShift[1];
      newMcp->rTime[2] = time[1]-time[3]+tShift[2];
      newMcp->rTime[3] = time[3]-time[4]+tShift[3];

      // Fill the new calibrated tree.      
      newTree->Fill();
    }
    newTree->AutoSave();
    newFile.Close();
    InitTree(inifChain);
  }
  stopWatch.Stop();
  stopWatch.Print();

  return nTrees;
}


//______________________________________________________________________________
Int_t TMcp::Calculate(Long64_t entry)
{
  // -- Calculate the MCP tracking system quantities.
  //

//   // Check that a TTree has been initialized.
//   if(!fChain){
//     printf("TTree or TChain has not been initialized!!!\n");
//     return kFALSE;
//   }

//   b_xAngle->Reset();
//   b_yAngle->Reset();
//   b_iAngle->Reset();
//   b_maskX->Reset();
//   b_maskY->Reset();
//   b_timeCorr->Reset();

//   Double_t zprime;
//   Double_t gap;

//   // Interaction position vectors for MCP0 and MCP1.
//   TVector3 posV0,posV1,aVec;

//   for(Long64_t i=0; i<fChain->GetEntries(); i++){
//     mcp0.b_x->GetEntry(i);
//     mcp0.b_y->GetEntry(i);
//     mcp0.b_xHG->GetEntry(i);
//     mcp0.b_yHG->GetEntry(i);
//     mcp0.b_z->GetEntry(i);
//     mcp1.b_xHG->GetEntry(i);
//     mcp1.b_yHG->GetEntry(i);
//     mcp1.b_z->GetEntry(i);
//     b_rTime->GetEntry(i);

//     // Set the interaction point for each MCP.
//     // Assumes that mcp0 is upstream.
//     posV0.SetXYZ(mcp0.x,mcp0.y,mcp0.z-mcpMcpGap);
//     posV1.SetXYZ(mcp1.x,mcp1.y,mcp1.z);
  
//     // Calculate 
//     aVec   = posV0-posV1;
//     iAngle = aVec.Theta();
//     xAngle = ATan(aVec.X()/aVec.Z());
//     yAngle = ATan(aVec.Y()/aVec.Z());
    
//     //Double_t z0 = mcp0.x*Sin(thetaMCP[0]*Pi()/180.0);
//     //Double_t z1 = mcp1.x*Sin(thetaMCP[1]*Pi()/180.0);

//     // Calculate xAngle.
// //     //    if (mcp0.x>0 && mcp1.x>0) {
// //    xAngle = atan((mcp0.x - mcp1.x)/(mcpMcpGap + z0 + z1))*180.0/Pi();
// //       //    }
// //     // Calculate yAngle.
// //     if(mcp0.y>0 && mcp1.y>0) {
// //       yAngle = atan((mcp0.y - mcp1.y)/(mcpMcpGap + z0 + z1))*180.0/Pi();
// //     }
  
// //     maskX = mcp1.x + (maskTarGap+z1)*tan(xAngle*Pi()/180.0);
// //     maskY = mcp1.y + (maskTarGap+z1)*tan(yAngle*Pi()/180.0);


// //     // timeCorr[0] = target MCP + Down Stream Plastic
  
// //     gap = tarMcpPlasticGap;
// //     if (rTime[0] > 0) timeCorr[0] = rTime[0]*(gap-z0)/gap;
  
// //     // timeCorr[1] = MCP + MCP
// //     gap = mcpMcpGap;
// //     Double_t dist = sqrt(pow((mcp0.x - mcp1.x),2)+pow((mcp0.y - mcp1.y),2)+pow((z0 - z1),2));
// //     if (rTime[1] > 0) timeCorr[1] = rTime[1]*(gap-dist)/gap;

// //     // timeCorr[2] = upstream MCP + Plastic
// //     gap = upMcpPlasticGap;
// //     if (rTime[2] > 0) timeCorr[0] = rTime[2]*(gap-z1)/gap;
    
//     b_xAngle->Fill();
//     b_yAngle->Fill();
//     b_iAngle->Fill();
//     b_maskX->Fill();
//     b_maskY->Fill();
//     b_timeCorr->Fill();
//   }

//   printf("* MCP tracking system is now calculated.\n");

  return 1;
}


//______________________________________________________________________________
void TMcp::CreateFolders()
{
  // -- Create default folders.
  //
  
  f_mcp0 = gROOT->GetRootFolder()->AddFolder("MCP0","MCP0 Folder");
  f_mcp1 = gROOT->GetRootFolder()->AddFolder("MCP0","MCP1 Folder");

  f_mcp0->Write();
  f_mcp1->Write();

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


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

  // Check that a TTree or TChain has been initialized.
  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.
  TMcp *newMcp = new TMcp(*this);
  newMcp->mcp0.fCon = this;
  newMcp->mcp1.fCon = this;
  newMcp->mcp0.fConD = &mcp0;
  newMcp->mcp1.fConD = &mcp1;

  // 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(),"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*)cTree->GetUserInfo()->FindObject("Run Number");
    TString fPath(path);
    TString fPrefix(prefix);
    fileName = runNumber->GetTitle();
    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.
    if(newFile.FindObject(oldTName+"Pri")!=0){printf("Tree already exists!\n"); return -2;}
    else if(strcmp(newTName,"")==0){
      newTName = oldTName+"Pri";
      newTree->SetName(newTName);
      TNamed aState("Analysis State","2");
      newTree->GetUserInfo()->Add(&aState);
    }else{
      newTree->SetName(newTName+"Pri");
      TNamed aState("Analysis State","2");
      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);
    }
    
    // Create Branches.
    newTree->Branch("mcp",&newMcp);
    newMcp->InitTree(newTree);

    // Generate Primary Data.
    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();}
      
      Clear();
      newMcp->mcp0.GenPrime(ientry);
      newMcp->mcp1.GenPrime(ientry);
      newMcp->GenPrime(ientry);

      newTree->Fill();
    }
    newTree->AutoSave();
    newFile.Close();
    InitTree(inifChain);
  }
  stopWatch.Stop();
  stopWatch.Print();

  delete newMcp;

  
  return nTrees;
}


//______________________________________________________________________________
Int_t TMcp::GenPrime(Long64_t entry)
{
  // -- Generate primary data for an entry.
  //

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

  Double_t zprime;
  Double_t gap;

  // Interaction position vectors for MCP0 and MCP1.
  TVector3 posV0,posV1,aVec;

//   fCon->mcp0.b_x->GetEntry(entry);
//   fCon->mcp0.b_y->GetEntry(entry);
//   fCon->mcp0.b_xHG->GetEntry(entry);
//   fCon->mcp0.b_yHG->GetEntry(entry);
//   fCon->mcp0.b_z->GetEntry(entry);
//   fCon->mcp1.b_xHG->GetEntry(entry);
//   fCon->mcp1.b_yHG->GetEntry(entry);
//   fCon->mcp1.b_z->GetEntry(entry);
//   fCon->b_rTime->GetEntry(entry);

  // Set the interaction point for each MCP.
  // Assumes that mcp0 is upstream.
  posV0.SetXYZ(mcp0.x,mcp0.y,mcp0.z-mcpMcpGap);
  posV1.SetXYZ(mcp1.x,mcp1.y,mcp1.z);
  
  // Calculate 
//   aVec   = posV0-posV1;
//   iAngle = (Pi()-aVec.Theta())*1000;     // [mrad]
//   xAngle = ATan(aVec.X()/aVec.Z())*1000; // [mrad]
//   yAngle = ATan(aVec.Y()/aVec.Z())*1000; // [mrad]
  aVec   = posV1-posV0;
  iAngle = aVec.Theta()*1000;            // [mrad]
  xAngle = ATan(aVec.X()/aVec.Z())*1000; // [mrad]
  yAngle = ATan(aVec.Y()/aVec.Z())*1000; // [mrad]
  
  //Double_t z0 = mcp0.x*Sin(thetaMCP[0]*Pi()/180.0);
  //Double_t z1 = mcp1.x*Sin(thetaMCP[1]*Pi()/180.0);
  
  // Calculate xAngle.
  //    if (mcp0.x>0 && mcp1.x>0) {
//   xAngle = atan((mcp0.x - mcp1.x)/(mcpMcpGap + z0 + z1))*180.0/Pi();
//   //    }
//   // Calculate yAngle.
//   if(mcp0.y>0 && mcp1.y>0) {
//     yAngle = atan((mcp0.y - mcp1.y)/(mcpMcpGap + z0 + z1))*180.0/Pi();
//   }
  
//   maskX = mcp1.x + (maskTarGap+z1)*tan(xAngle*Pi()/180.0);
//   maskY = mcp1.y + (maskTarGap+z1)*tan(yAngle*Pi()/180.0);
  
//   // timeCorr[0] = target MCP + Down Stream Plastic
  
//   gap = tarMcpPlasticGap;
//   if (rTime[0] > 0) timeCorr[0] = rTime[0]*(gap-z0)/gap;
  
//   // timeCorr[1] = MCP + MCP
//   gap = mcpMcpGap;
//   Double_t dist = sqrt(pow((mcp0.x - mcp1.x),2)+pow((mcp0.y - mcp1.y),2)+pow((z0 - z1),2));
//   if (rTime[1] > 0) timeCorr[1] = rTime[1]*(gap-dist)/gap;
  
//   // timeCorr[2] = upstream MCP + Plastic
//   gap = upMcpPlasticGap;
//   if (rTime[2] > 0) timeCorr[0] = rTime[2]*(gap-z1)/gap;

  return 1;
}


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

  Char_t *stageOp = (Char_t*)stage;

  if(strcmp(stageOp,"Raw")==0){
    printf("* Disabling calibrated branches.\n");
    tree->SetBranchStatus("mcp0*",0);
    tree->SetBranchStatus("mcp1*",0);
    tree->SetBranchStatus("xAngle",0);
    tree->SetBranchStatus("yAngle",0);
    tree->SetBranchStatus("iAngle",0);
    tree->SetBranchStatus("mcp0.corner[12]",1);
    tree->SetBranchStatus("mcp1.corner[12]",1);
    tree->SetBranchStatus("mcp0.cMult",1);
    tree->SetBranchStatus("mcp1.cMult",1);
  }else if(strcmp(stageOp,"Cal")==0){
    printf("* Disabling raw branches.\n");
    tree->SetBranchStatus("mcp*",0);
    tree->SetBranchStatus("mcp0*",1);
    tree->SetBranchStatus("mcp1*",1);
    tree->SetBranchStatus("mcp0.corner[12]",0);
    tree->SetBranchStatus("mcp1.corner[12]",0);
  }else if(strcmp(stageOp,"Pri")==0){
    printf("* Enabling only primary branches.\n");
    tree->SetBranchStatus("mcp*",0);
    tree->SetBranchStatus("mcp0*",0);
    tree->SetBranchStatus("mcp1*",0);
    tree->SetBranchStatus("mcp.xAngle",1);
    tree->SetBranchStatus("mcp.yAngle",1);
    tree->SetBranchStatus("mcp.iAngle",1);
  }  

}


//______________________________________________________________________________
void TMcp::InitClass()
{
  // -- Initialize the MCP tracking system.
  // 1.  All class data members are set to their default values.
  // 2.  The InitClass() method for any subclass is called.
  // 3.  QDC and TDC mapping file is loaded.

  // Set the object identifier and object title.
  SetNameTitle("mcp","MCP Tracking System.");

  // Load the MCP QDC and TDC mapping file.
  printf("* MCP QDC and TDC Map Loading . . .   ");
  if(LoadQdcTdcMap()){
    printf("\t\t\t\t\t\t [OK]\n");
  }else{
    printf("\t\t\t\t\t\t [FAILED!!!]\n");
  }
   
  // Initialize the MCP0 and MCP1 
  mcp0.InitClass("mcp0",0,TMCP0_NSTAGES);
  mcp1.InitClass("mcp1",1,TMCP1_NSTAGES);
  mcp0.parent = this;
  mcp1.parent = this;

  Clear();

  // By default add a random number to Short_t data.
  fAddRndm = kTRUE;

  // Set default parameter values.
  tShift[0] = 0;
  tShift[1] = 0;
  tShift[2] = 0;
  tShift[3] = 0;
  
  thetaMCP[0] = 0;
  thetaMCP[1] = 0;

  tarMcpPlasticGap = 0;
  mcpMcpGap        = 0;
  upMcpPlasticGap  = 0;
  maskTarGap       = 0;
}


//______________________________________________________________________________
void TMcp::InitTree(TTree *tree)
{
  // -- Initialize the MCP fChain.
  // Set fChain to given tree and set the branch addresses.

  fChain   = tree;
  fCurrent = -1;
  //  TMcp  *p  = this;
  //  TMcp **pp = &p;
  //  fChain->SetBranchAddress("mcp",&pp);
  
  Char_t chargeName[200];
  sprintf(chargeName,"charge[%i][32]",TMCP_NQDCS);

  b_time     = fChain->GetBranch("time[32]");
  b_charge   = fChain->GetBranch(chargeName);
  b_ene      = fChain->GetBranch("ene[32]");
  b_rTime    = fChain->GetBranch("rTime[32]");
  b_timeCorr = fChain->GetBranch("timeCorr[32]");
  b_maskX    = fChain->GetBranch("maskX");
  b_maskY    = fChain->GetBranch("maskY");
  b_xAngle   = fChain->GetBranch("xAngle");
  b_yAngle   = fChain->GetBranch("yAngle");
  b_iAngle   = fChain->GetBranch("iAngle");

  //  b_mcp0     = fChain->GetBranch("mcp0");
  //  b_mcp1     = fChain->GetBranch("mcp1");

  // Initialize sub-class fChains.
  mcp0.InitTree(tree);
  mcp1.InitTree(tree);
}


//______________________________________________________________________________
Bool_t TMcp::LoadQdcTdcMap()
{
  // -- Load the QDC/TDC mapping file.
  //

  Int_t  slotNum,chanNum,mcpNum,sigType;

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

  Char_t line[200];

  file.getline(line,200);
  while(!isdigit(line[0])){
    file.getline(line,200);    
  }

  for(Int_t i=0; i<TMCP_NQDCS; i++) fQdcSlot[i] = -1;
  Int_t curSlot = 0;
  while(!file.eof() || isdigit(line[0])){
    sscanf(line,"%i %i %i %i",&slotNum,&chanNum,&mcpNum,&sigType);
    // Order the QDC's by slot.
    if(fQdcSlot[curSlot]!=slotNum && fQdcSlot[curSlot]!=-1){curSlot++;}
    fQdcSlot[curSlot] = slotNum;

//     if(fQdcSlot[0]==-1)      fQdcSlot[0] = slotNum;
//     if(slotNum>fQdcSlot[0])  fQdcSlot[1] = slotNum;
//     if(slotNum<fQdcSlot[0]) {fQdcSlot[1] = fQdcSlot[0]; fQdcSlot[0] = slotNum;}

    if(mcpNum==0 && sigType<9){
      mcp0.fSigMap[sigType]    = chanNum;
      mcp0.fChargeMap[sigType] = curSlot;
    }else if(mcpNum==1 && sigType<9){
      mcp1.fSigMap[sigType]    = chanNum;
      mcp1.fChargeMap[sigType] = curSlot;
    }else{
      printf("Invalid MCP number or signal type!!!\n");
      return kFALSE;
    }
    file.getline(line,200);
  }

  file.close();
  return kTRUE;
}


//______________________________________________________________________________
Bool_t TMcp::SetMcpMapFile(Char_t *filePath)
{
  // -- Set the path of the MCP mapping file.
  //

  Bool_t fFound = kFALSE;

  mcpMapFile = filePath;
  fFound = kTRUE;

  return fFound;
}


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

  fName = name;
}


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

  fName  = name;
  fTitle = title;
}


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

  fTitle = title;
}


//______________________________________________________________________________
Bool_t TMcp::Unpack(UShort_t *pEvent)
{
  // -- Unpack raw MCP data from an event buffer.
  //

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

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


  *p--;
  *p--;
  packetSize = *p;                   // Number of words in the MCP 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++;
    Bool_t   subFlag      = kFALSE;

    UShort_t		ID;

    switch(subPacketTag) {
    case VMEADC:
      ID = *p++;
      if (subPacketSize >3) p = UnpackEnergy(p);
      break;
    case VMETDC:
      ID = *p++;
      if (subPacketSize >3) p = UnpackTime(p);
      break;
    case VMEQDC:
      ID = *p++;
      if (subPacketSize >3) p = UnpackCharge(p);
      break;
    default:
      printf("Unrecognized sub-packet tag: 0x%x\n", subPacketTag);
      Int_t count = 0;
      for (UShort_t ii = 0; ii < *p; ii++) {
	printf("0x%4x ",*p++);
	if (count == 9) {
	  count = 0;
	  printf("\n");
	}
	count++;
      }
      printf("\n");
      
      p += subPacketSize-2;
      return kFALSE;
      break;
    }	
    packetSize -= subPacketSize;
  }

  return kTRUE;
}


//______________________________________________________________________________
UShort_t* TMcp::UnpackEnergy(UShort_t *p)
{
  // -- Unpack the raw MCP energy.
  //

  Int_t m_nCrate,m_nSlot,num_evts=0;
  Int_t channel;
  Short_t value;

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

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

  p++;
  p++;

  return p;
}


//______________________________________________________________________________
UShort_t* TMcp::UnpackCharge(UShort_t *p)
{
  // -- Unpack the raw MCP charge.
  //

  Int_t m_nCrate,m_nSlot,num_evts=0;
  Int_t channel;
  Short_t value;

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

  for (Int_t i=0;i<num_evts;i++){
    channel = (Int_t)((*p++)&(0x003f));
    value   = (Short_t)((*p++)&0x0fff);
    
    for(Int_t j=0; j<TMCP_NQDCS; j++){
      if(m_nSlot==fQdcSlot[j]) charge[j][channel] = value;
    }
  }
  
  p++;
  p++;

  return p;
}


//______________________________________________________________________________
UShort_t* TMcp::UnpackTime(UShort_t *p)
{
  // -- Unpack the raw MCP time.
  //

  Int_t m_nCrate,m_nSlot,num_evts=0;
  Int_t channel;
  Short_t value;

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

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

  return p;
}