//
// Create trees for energy 
// from adc trees and fit64ch.dat
//

void mkEnergyTree() 
{
  gROOT->Reset();
  
  const Int_t ChipChannel = 32;
  const Int_t ChipNum = 2;
  const Int_t MaxChannel = ChipChannel * ChipNum;
  
  Int_t ch, chip, headerFlg=1;
  Char_t line[128];
  Double_t peak, sigma, eres, counts, chisq;
  Double_t rand, gain[MaxChannel]; // energy/channel[keV]
  
  // // Open input and output files // //
  Char_t treeFile[128]="tree.root";
  Char_t energyTreeFile[128]="energyTree.root";
  Char_t dataFile[128]="fit64ch.dat";
  CLtxtrd("input root tree file : ", treeFile, 128);
  CLtxtrd("output root tree file : ", energyTreeFile, 128);
  CLtxtrd("input gain data file : ", dataFile, 128);
  
  // // Read gain data file // // 
  ifstream data(dataFile);
  if (!data) {
    cout << "Cannot open " << dataFile << "." << endl;
    return;
  }
  
  while (!data.eof()) {
    
    if (headerFlg) { 
      data.getline(line, 128);
      headerFlg = 0; 
    }
    
    else if (data>>ch>>peak>>sigma>>eres>>counts>>chisq) {
      gain[ch] = 59.5/peak;
    }
  }
  data.close();
  
  // // Make Tree for Energy // //
  Int_t i_cycles, numEvents, event_seq_no;
  Short_t adc[MaxChannel];  
  Double_t ped[MaxChannel], cmn[ChipNum];
  Float_t energy[MaxChannel];
  
  TFile *infile = new TFile(treeFile);
  TTree *tree = (TTree*)infile->Get("tree");
  tree->SetBranchAddress("i_cycles",&i_cycles);
  tree->SetBranchAddress("event_seq_no",&event_seq_no);
  tree->SetBranchAddress("adc",adc);
  tree->SetBranchAddress("ped",ped);
  tree->SetBranchAddress("cmn",cmn);
  
  TTree *energyTree = new TTree("energyTree", "tree for energy");
  energyTree->Branch("energy", energy, "energy[64]/F");

  TFile *outfile = new TFile(energyTreeFile, "recreate");

  Long_t nentries = (Long_t)tree->GetEntries();
  cout << nentries << " entries. "<< endl;
  
  for (Long_t i=0; i<nentries; i++) {
    tree->GetEntry(i);
    
    for (chip=0; chip<ChipNum; chip++) {
      for( ch=0; ch<ChipChannel; ch++ ) {
	
	rand = gRandom->Rndm() - 0.5;
	energy[ch+chip*ChipChannel] =
	  ( adc[ch+chip*ChipChannel] - ped[ch+chip*ChipChannel]
	    - cmn[chip] + rand ) * gain[ch+chip*ChipChannel];
	energyTree->Fill();
      }
    }
    
    if (i % 1000000 == 1) {
      energyTree->AutoSave("SaveSelf");
      cout << "Filled " << i << " events." << endl;
    }
  }
  energyTree->Write();

  return;
}