//////////////////////////////////////////////////////////////////////////
// $Name$ - $Id$
// Purpose: classify QFS experiment events per run

void qfs_classify(char *fn) {
  // Load calibrations
  TFile fcali("calibration.root");
  TFile f(fn,"update");
  TF1 *fline1 = (TF1 *)fcali.Get("fline1");
  TF1 *fline2 = (TF1 *)fcali.Get("fline2");
  TF1 *flinem = (TF1 *)fcali.Get("flinem");
  f.cd();

  Float_t  ph1,tof1,psd1,ph2,tof2,psd2,tof12,phm,tofm,psdm,deltaE,E;
  TNtuple *nt = (TNtuple *)f.Get("kmax_event");
  Int_t N = nt->GetEntries();
  nt->SetBranchAddress("var0",&ph1);
  nt->SetBranchAddress("var1",&tof1);
  nt->SetBranchAddress("var2",&psd1);
  nt->SetBranchAddress("var3",&ph2);
  nt->SetBranchAddress("var4",&tof2);
  nt->SetBranchAddress("var5",&psd2);
  nt->SetBranchAddress("var6",&tof12);
  nt->SetBranchAddress("var7",&phm);
  nt->SetBranchAddress("var8",&tofm);
  nt->SetBranchAddress("var9",&psdm);
  nt->SetBranchAddress("var10",&deltaE);
  nt->SetBranchAddress("var11",&E);

  // QFS events -- TTree *qfs
  TTree qfs("qfs","QFS events");
  qfs.Branch("ph1",&ph1,"ph1/F");
  qfs.Branch("psd1",&psd1,"psd1/F");
  qfs.Branch("ph2",&ph2,"ph2/F");
  qfs.Branch("psd2",&psd2,"psd2/F");
  qfs.Branch("tof1",&tof1,"tof1/F");
  qfs.Branch("tof2",&tof2,"tof2/F");
  qfs.Branch("tof12",&tof12,"tof12/F");
  // Monitor events -- TTree *monitor
  TTree monitor("monitor","monitor events");
  monitor.Branch("phm",&phm,"phm/F");
  monitor.Branch("psdm",&psdm,"psdm/F");
  monitor.Branch("tofm",&tofm,"tofm/F");
  // deltaE-E events -- TTree *dE_E
  TTree dE_E("dE_E","deltaE-E events");
  dE_E.Branch("E",&E,"E/F");
  dE_E.Branch("deltaE",&deltaE,"deltaE/F");

  // Processing events
  printf("[--] Processing events ");
  for(Int_t idx=0;idx<N;idx++) {
    nt->GetEntry(idx);
  
    ph1 /= 8.; tof1 /= 8.; psd1 /= 8.;
    ph2 /= 8.; tof2 /= 8.; psd2 /= 8.;
    tof12 /= 8.;
    if( (ph1>fline1->Eval(472.)) && (ph2>fline2->Eval(472.)) )
      qfs.Fill();
     
    phm /= 8.; tofm /= 8.; psdm /= 8.;
    if(phm>flinem->Eval(472.))  monitor.Fill();
    
    deltaE /= 8.; E /= 8.;
    dE_E.Fill();
    
    if( ((idx+1)%1000)==0 ) printf(".");
  }
  printf("\n");

  qfs.Write();
  qfs.Print();
  monitor.Write();
  monitor.Print();
  dE_E.Write();
  dE_E.Print();

  f.Close();
  fcali.Close();
}