#include "TTreeFormula.h"

void MedFilterTree(char* inFile,const char* selection,char* outFile)
{

  TFile* f = new TFile(inFile,"READ");
  cout << "Read in file: " << inFile << endl;
  TTree* t = (TTree*)f->Get("pan");
  cout << "Got TTree." << endl;

  // ttf implements the event selection
  TTreeFormula ttf("filter_tree",selection,t);

  TFile* g = new TFile(outFile,"recreate");
  g->cd();

  TTree* tc = t->CloneTree(0);

  Long64_t nbr=0;
  std::cout.precision(3);
  Long64_t i=0;
  Long64_t npass=0;
  bool keep_going=true;
  Int_t current_tree=-1;

  while(true) {
    Long64_t local_entry=t->LoadTree(i);
    if(t->GetTreeNumber()!=current_tree) {
      current_tree=t->GetTreeNumber(); 
      ttf.UpdateFormulaLeaves();
    }
    // load branches needed by selection
    
    for(Int_t j=0; j< ttf.GetNcodes(); j++){
      Long64_t nb = ttf.GetLeaf(j)->GetBranch()->GetEntry(local_entry);
      if(nb<=0){
	// done reading
	std::cout<<"Should be done (1)? : [i,j] = "<<i<<" , "<<j<<std::endl;
	keep_going=false;
	break;
      }
      nbr+=nb;
    }
    if(!keep_going) break;
    
    // evaluate the selection, if true, fill output
    Double_t val = ttf.EvalInstance();
    if(val > 0.0){
      npass++;
      Long64_t nb=t->GetEntry(i,0);
      if(nb<=0){
	// done reading
	std::cout<<"Should be done (2)?"<<std::endl;
	keep_going=false;
	continue;
      }
      nbr+=nb;      
      tc->Fill();
    }
    // how's it going
    if(i%100000==0){
      std::cout<<i<<" entries  read : "<<npass<<" filled : "
      <<nbr/(1024.0*1024.0)<<" Mb read"<<std::endl;
    }
    
    i++; // don't forget
  }
  std::cout<<"Done                                              "<<std::endl;
  std::cout<<"filter_tree: output has "
	   <<npass<<" entries."<<std::endl;  

  tc->Write();
  g->Write();
  g->Close();

}