// ----------------------------------------------------------------------------
// $Id: root_GoNNLCDVTopl.C,v 1.1 2001/06/23 23:26:24 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: root_GoNNLCDVTopl.C,v $
// Revision 1.1  2001/06/23 23:26:24  toshi
// This is a sample file for topological vertexing and neural network heavy
// flavor tag.
//
//
// Run Fast MC interactively with Root
// Sample program using LCDVTopl with NN H.F. tag (translated from SLD ZVTOP
// and B3MASS) 
//                                              June/23/2001 T.Abe
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x GoLCDNNVTopl.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x GoLCDNNVTopl.C(300)
//
// !!!Notice!!!
// NN was optimized for sqrt(s)=91.26 GeV evnets. Hence, you might have to
// reoptimize the NN for your event analysis.
//
/*
#include "TH1.h"
#include "TCanvas.h"

#include "LCDEvent.h"
#include "LCDreadRootFile.h"
#include "LCDFastMC.h"
#include "LCDJetFinder.h"
#include "LCDVTopl.h"

#include "TSystem.h"

TSystem *gSystem=new TSystem();

#include <stdlib.h>
*/
Int_t root_GoNNLCDVTopl(Int_t nEvent=300){
  
  gROOT->Reset();
  
  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;parfile_dir += "/";
  parfile_dir += gSystem->Getenv("LCDVERSION");parfile_dir += "/";
  parfile_dir += "ParFiles/";

  // In order to refer to shareable libraries in this simple 
  // form, include their paths in a .rootrc file.
  
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");

  LCDEvent event;

  // Input generator level root data
  Char_t* rootFile = (Char_t*)"test.root";
  LCDreadRootFile source(rootFile,&event);

  // Open detector parameter file
  Char_t* geomFile = 
    "Silicon.par"; //Silicon
  //"Large.par";   // Large
  //"Small.par";   // Small
  TString s_geomFile(parfile_dir);
  s_geomFile += geomFile;
  FILE* m_geomfile =  fopen(s_geomFile.Data(),"r");
  LCDGetParameters gp(m_geomfile);
  fclose(m_geomfile);

  // If you want to change some detector parameter(s) for FastMC 
  // do here
  // For example..
  //gp.SetMagneticField(2.);  // Change Magnetif field

  // Tracking smearing parameter file
  Char_t* smearFile = 
    "lookup_silicon.nobmcon";
  //"lookup_large.nobmcon";
  //"lookup_small.nobmcon";
  TString s_smearFile(parfile_dir);
  s_smearFile += smearFile;

  // Initialize FastMC
  LCDFastMC fmc(&gp,s_smearFile.Data(),&event); 
  // If you don't want the Cluster merging
  //fmc.SetNoClusterMerge();

  // Define histograms
  TH1F* h_nvrt = new TH1F("h_nvrt" ,"N of LCDVToplVRT"        ,10,0.0,10.0);
  TH1F* h_nvrtl= new TH1F("h_nvrtl","N of LCDVToplVRT (uds+g)",10,0.0,10.0);
  TH1F* h_nvrtc= new TH1F("h_nvrtc","N of LCDVToplVRT (c)"    ,10,0.0,10.0);
  TH1F* h_nvrtb= new TH1F("h_nvrtb","N of LCDVToplVRT (b)"    ,10,0.0,10.0);

  TH1F* h_nn   = new TH1F("h_nn" ,"N.N. for H.F. tag"        ,20,0.0,1.0);
  TH1F* h_nnl  = new TH1F("h_nnl","N.N. for H.F. tag (uds+g)",20,0.0,1.0);
  TH1F* h_nnc  = new TH1F("h_nnc","N.N. for H.F. tag (c)"    ,20,0.0,1.0);
  TH1F* h_nnb  = new TH1F("h_nnb","N.N. for H.F. tag (b)"    ,20,0.0,1.0);

  h_nvrt ->SetFillColor(kBlue);
  h_nvrtl->SetFillColor(kGreen);
  h_nvrtc->SetFillColor(kYellow);
  h_nvrtb->SetFillColor(kMagenta);

  h_nn ->SetFillColor(kBlue);
  h_nnl->SetFillColor(kGreen);
  h_nnc->SetFillColor(kYellow);
  h_nnb->SetFillColor(kMagenta);

  TCanvas *c1 = new TCanvas("c1", "c1", 600,600);
  c1->Divide(2,2);
  TCanvas *c2 = new TCanvas("c2", "c2", 600,600);
  c2->Divide(2,2);

  // Setup of JetFinder
  LCDJetFinder jetfind(0.004);  
  TClonesArray forjet("TLorentzVector",1000);
  
  // Track cut value
  Double_t cosCut=0.9, pCut=0.2;

  // Event loop
  LCDVTopl topl;
  Int_t    ierr;
  Int_t    iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    // Read an event from the ROOT file.
    if (!(ierr=source.GetEvent())) goto canvasdraw;

    // Through FastMC
    fmc.Doit();

    if (iEvent == 0){
      // List Event Header information, for example
      LCDEventHeader* header = (LCDEventHeader*)event.EventHeader();
      printf("\n");
      printf("Process:: %s\n",header->GetProcessName());

      // List beam parameters 
      LCDBeam* beam = (LCDBeam*)event.BEAM();
      printf("Ecm = %f\n",beam->GetEcm());
      printf("electron polarization : %f  positron polarization : %f\n",
	     beam->GetPolElec(), beam->GetPolPosi());
      printf("==> effective polarization for electron : %f \n\n",
	     beam->GetEffectivePol());
    }

    // Look at tracks. Can access them by the [] operator.
    Int_t nTracks = event.Tracks()->GetEntries();
    
    // Cut 1...
    if (nTracks < 5) goto canvasdraw;
    
    // Check number of good tracks
    TObjArray gtrk_list;
    LCDTrack* trk=0;
    TVector3  ptrk;
    Int_t itrk=0;
    for (itrk=0; itrk<nTracks; itrk++) {
      trk = (LCDTrack*)(event.Tracks()->At(itrk));
      ptrk=trk->GetMomentumVector(0.0);
      if (TMath::Abs(ptrk.CosTheta())<cosCut && ptrk.Mag()>pCut) {
	gtrk_list.Add((TObject*)trk);
      }
    }
    
    //Cut 2...
    Int_t nGTracks=gtrk_list.GetEntries();
    if (nGTracks < 5 ) goto canvasdraw;
    TVector3 vec3; 
    // Do Jet Finding
    forjet.Delete(); // Clean up veclist
    for (itrk=0; itrk<nGTracks; itrk++) {
      trk = (LCDTrack*)(gtrk_list.At(itrk));
      vec3 = trk->GetMomentumVector(0.0);
      Double_t Etrk = vec3.Mag2() + 0.1396 * 0.1396;// Assume pion mass
      new((forjet)[itrk]) TLorentzVector(vec3,Etrk);
    }

    jetfind.SetPartList(&forjet); // Input particle 4-vector list
    jetfind.SetDURHAM();         // Use DURHAM algorithm
    jetfind.doFindJets(4);        // JetFindig (force to find 4 jet)
      

    TVector3 ipv=((LCDBeam*)(event.BEAM()))->GetPos(); // IP position
    //TVector3 ipv(0.,0.,0.);
    TVector3 bfield(0.0,0.0,gp.GetMagneticField());    // B field
    
    TLorentzVector p4vrt;
    Double_t rawMass, ptMass, massTag;
        
    TLorentzVector jet4vec; TVector3 jet3vec;
    Int_t nJet = jetfind.njets();

    //Loop Over # of jets
    Int_t* trkList=new Int_t[nGTracks];
    for (Int_t iJet=0 ; iJet < nJet ; iJet++) {

      Int_t nTrks = jetfind.nParticlesPerJet(iJet);

      if (nTrks <= 2) continue; 
      Int_t  iTrk2=0;
      Bool_t f_b=kFALSE;
      Bool_t f_c=kFALSE;
      for (int ipart =0; ipart < nGTracks; ipart++){
	if( jetfind.GetPartIndexAt(ipart) == iJet){
	  trkList[iTrk2++] = ipart;
	  trk = (LCDTrack*)(gtrk_list.At(ipart));
	  LCDMcPart *tkmc=
	    (LCDMcPart*)event.MCparticles()->At(trk->GetParticle());
	  while(tkmc) {
	    Int_t apid=TMath::Abs(tkmc->GetParticleID());
	    switch(apid) {
	    case  511: // Bd
	    case  521: // Bu
	    case  531: // Bs
	    case  541: // Bc
	    case 5122: // Lmbda_b
	      f_b=kTRUE;
	      break;
	    case   411: // D+
	    case   421: // D0
	    case   431: // Ds+
	    case   433: // Ds*+
	    case   441: // eta_c
	    case   443: // J/psi
	    case  4122: // Lambda_c
	      f_c=kTRUE;
	      break;
	    }
	    if (f_b) break;
	    Int_t iparent=tkmc->GetParentIdx();
	    if (iparent < 0) break;
	    tkmc=(LCDMcPart*)event.MCparticles()->At(iparent);
	  }
	}
      }
      jet3vec = jetfind.jet4vec(iJet)->Vect();

      topl.SetDebugLevel(0);
      topl.SetUpLCDVTopl(&gtrk_list,nTrks,trkList,ipv,bfield,jet3vec);
      topl.SetUseNN(kTRUE); // Use NN method
      topl.FindVertices(1); // Find Vertices (1 means V0 check)

      Int_t nVRT=topl.GetVertexEntriesExceptV0();// # of reconstructed vertices
      h_nvrt->Fill(nVRT);
      if (f_b) {
	h_nvrtb->Fill(nVRT);
      } else if (f_c) {
	h_nvrtc->Fill(nVRT);
      } else {
	h_nvrtl->Fill(nVRT);
      }
      if (nVRT <= 1) continue;

      Double_t nnout=topl.GetHFTagNN(); // Get NN output for HF tag

      h_nn->Fill(nnout);
      if (f_b) {
	h_nnb->Fill(nnout);
      } else if (f_c) {
	h_nnc->Fill(nnout);
      } else {
	h_nnl->Fill(nnout);
      }

    }    

    delete trkList; 

  canvasdraw:    
    if (iEvent%100 == 0 || iEvent == nEvent-1) {
      c1->cd(1); h_nvrt ->Draw("B");
      c1->cd(2); h_nvrtl->Draw("B");
      c1->cd(3); h_nvrtc->Draw("B");
      c1->cd(4); h_nvrtb->Draw("B");
      c1->Modified();
      c1->Update();

      c2->cd(1); h_nn ->Draw("B");
      c2->cd(2); h_nnl->Draw("B");
      c2->cd(3); h_nnc->Draw("B");
      c2->cd(4); h_nnb->Draw("B");
      c2->Modified();
      c2->Update();
    }
    if (iEvent%1000 == 999 || iEvent == nEvent-1) {
      printf("iEvent = %d was done...",iEvent+1);
    }
    if (!ierr) break;
  }
  forjet.Delete();
  topl.Clean();
  event.Clean();

  return iEvent;
}