// Run Fast MC interactively with Root
// Sample program using LCDVTopl (transfered from SLD ZVTOP) 
//                                              Sep/28/2000 M. Iwasaki
//
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x GoLCDVTopl.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x GoLCDVTopl.C(300)
//

#include <stdlib.h>

Int_t GoLCDVTopl(Int_t nEvent=300){
  
  gROOT->Reset();
  
  // In order to refer to shareable libraries in this simple 
  // form, include their paths in a .rootrc file.
  
  gSystem->Load("libLCDEvent.so");
  gSystem->Load("libLCDRootAppsUtil.so");
  gSystem->Load("libLCDFastMC.so");
  gSystem->Load("libLCDPhUtil.so");

  // Input STDHEP file
  char* stdFile = "test.dat";

  // Detector parameter file
  char* geomFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/Small.par";

  // Tracker smearing parameter file
  char* smearFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/lookup_small.nobmcon";
  
  LCDreadStdFile*  source    = new LCDreadStdFile(stdFile);
  LCDFastMC*       fmc       = new LCDFastMC(source,geomFile,smearFile,0);
  LCDEvent*        event     = 0;
  FILE*            fgeomFile = fopen(geomFile,"r");
  LCDGetParameters gp(fgeomFile);

  // Define histograms
  TH1F* h_nvtvrt = new TH1F("h_nvrt","N of LCDVToplVRT"       ,10,0.0,10.0);
  TH1F* h_nvtvrtu= new TH1F("h_nvrtu","N of LCDVToplVRT uds"  ,10,0.0,10.0);
  TH1F* h_nvtvrtc= new TH1F("h_nvrtc","N of LCDVToplVRT c"    ,10,0.0,10.0);
  TH1F* h_nvtvrtb= new TH1F("h_nvrtb","N of LCDVToplVRT b"    ,10,0.0,10.0);

  TH1F* h_rwm    = new TH1F("h_rwm" ,"Raw Vrt Mass"           ,40,0.0,10.0);
  TH1F* h_rwmu   = new TH1F("h_rwmu","Raw Vrt Mass uds"       ,40,0.0,10.0);
  TH1F* h_rwmuc  = new TH1F("h_rwmuc","Raw Vrt Mass udsc"     ,40,0.0,10.0);
  TH1F* h_ptm1   = new TH1F("h_ptm1" ,"Pt corr. Vertex Mass"  ,28,0.0,10.0);
  TH1F* h_ptm1u  = new TH1F("h_ptm1u" ,"Pt corr. Vertex Mass" ,28,0.0,10.0);
  TH1F* h_ptm1uc = new TH1F("h_ptm1uc","Pt corr. Vertex Mass" ,28,0.0,10.0);
  TH1F* h_ptm2   = new TH1F("h_ptm2" ,"Pt corr. Vertex Mass"  ,28,0.0,10.0);
  TH1F* h_ptm2u  = new TH1F("h_ptm2u" ,"Pt corr. Vertex Mass" ,28,0.0,10.0);
  TH1F* h_ptm2uc = new TH1F("h_ptm2uc" ,"Pt corr. Vertex Mass",28,0.0,10.0);

  TH2F* h_mvspu  = new TH2F("h_mvspu" ,"Pt Mass vs Vertex P uds" ,
			   20,0.0,2.0,25,0.0,250.);
  TH2F* h_mvspc  = new TH2F("h_mvspc" ,"Pt Mass vs Vertex P c"   ,
			   20,0.0,2.0,25,0.0,250.);
  TH2F* h_mvspb  = new TH2F("h_mvspb" ,"Pt Mass vs Vertex P b"   ,
			   20,0.0,2.0,25,0.0,250.);


  h_nvrt   ->SetFillColor(4);
  h_nvrtu  ->SetFillColor(4);
  h_nvrtc  ->SetFillColor(4);
  h_nvrtb  ->SetFillColor(4);

  h_rwm   ->SetFillColor(4);
  h_ptm1  ->SetFillColor(4);
  h_ptm2  ->SetFillColor(4);
  h_rwmu  ->SetFillColor(3);
  h_ptm1u ->SetFillColor(3);
  h_ptm2u ->SetFillColor(3);
  h_rwmuc ->SetFillColor(2);
  h_ptm1uc->SetFillColor(2);
  h_ptm2uc->SetFillColor(2);

  h_mvspu ->SetFillColor(4);
  h_mvspc ->SetFillColor(4);
  h_mvspb ->SetFillColor(4);


  TCanvas *c1 = new TCanvas("c1", "c1", 700,700);
  c1->Divide(2,2);
  TCanvas *c2 = new TCanvas("c2", "c2", 700,900);
  c2->Divide(2,2);
  TCanvas *c3 = new TCanvas("c3", "c3", 700,700);
  c3->Divide(2,2);

  // Setup of JetFinder
  LCDJetFinder* jetfind   = new LCDJetFinder(0.004);  
  TObjArray*    forjet    = new TObjArray();
  
  // Track cut value
  Double_t cosCut=0.9, pCut=0.2;

  // Event loop
  for (int iEvent = 0; iEvent < nEvent; iEvent++) {
    
    event = new LCDEvent();
    event->Create();
    Int_t ierr = fmc->FetchEvent(event);
    if (!ierr) break;
    
    fmc->Cleanup();
    fmc->Doit(event);


    // Tag flavor
    TObjArray* McPart = event->MCparticles();
    Int_t nMCpart = McPart->GetEntries();
    
    LCDMcPart *mcpart = (LCDMcPart *) McPart->At(6);
    Int_t flav = abs(mcpart->GetType());

    // Look at tracks. Can access them by the [] operator.
    Int_t nTracks = event->Tracks()->GetEntries();
    
    // Cut 1...
    if (nTracks < 5) { event->Clean(); delete event; continue; }
    
    // Check number of good tracks
    TObjArray gtrk_list;
    LCDTrack* trk=0;
    TVector3  ptrk(0.0,0.0,0.0);
    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 ) {
      event->Clean(); delete event; continue; 
    }	
    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
      TLorentzVector* vec4 = new TLorentzVector(vec3,Etrk);
      forjet->Add(vec4);  // Make particle 3-vector list
    }

    jetfind->SetPartList(forjet); // Input particle 4-vector list
    jetfind->SetDURHAM();         // Use DURHAM algorithm
    jetfind->SetYCut(0.004);      // Set YCut value
    jetfind->doFindJets();        // JetFindig
      

    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;
    TArrayI partindex = jetfind->GetPartIndex();
    Int_t nJet = jetfind->njets();

    //Loop Over # of jets
    for (Int_t iJet=0 ; iJet < nJet ; iJet++) {

      Int_t nTrks = jetfind->nParticlesPerJet(iJet);

      if (nTrks <= 2) continue; 
      Int_t* trkList=new Int_t[nTrks];
      Int_t  iTrk2=0;
      for (int ipart =0; ipart<nGTracks; ipart++){
	if( partindex[ipart] == iJet){
	  trkList[iTrk2] = ipart;
	  iTrk2++;
	}
      }	    
      jet3vec = (jetfind->jet4vec(iJet)).Vect();

      
      LCDVTopl* topl=
	new LCDVTopl(&gtrk_list,nTrks,trkList,ipv,bfield,jet3vec);
      
      topl->FindVertices(); //Find Vertices
      topl->CalcPTMass();   //Calcurate Pt corrected mass
      
      Int_t nVRT=topl->GetVertexEntries(); // # of reconstructed vertices
      h_nvtvrt->Fill(nVRT);
      if (flav <= 3) h_nvtvrtu->Fill(nVRT);
      if (flav == 4) h_nvtvrtc->Fill(nVRT);
      if (flav == 5) h_nvtvrtb->Fill(nVRT);
      
      p4vrt=topl->GetP4Vertex();
      
      rawMass  = topl->GetRawMass();
      if (rawMass <= 0.) { delete topl; delete trkList; continue; }
      
      ptMass   = topl->GetPtMass();
      massTag  = topl->GetMassTag();
      
      h_rwm ->Fill(rawMass);
      h_ptm1->Fill(ptMass);
      h_ptm2->Fill(massTag);
      
      if (flav <= 3) {
	h_rwmu ->Fill(rawMass);
	h_ptm1u->Fill(ptMass);
	h_ptm2u->Fill(massTag);
      }
      if (flav <= 4) {
	h_rwmuc ->Fill(rawMass);
	h_ptm1uc->Fill(ptMass);
	h_ptm2uc->Fill(massTag);
      }
      if (massTag<2.) {
	if (flav <= 3) h_mvspu->Fill(massTag,p4vrt.Rho());
	if (flav == 4) h_mvspc->Fill(massTag,p4vrt.Rho());
	if (flav == 5) h_mvspb->Fill(massTag,p4vrt.Rho());
      }
      
      delete topl;
    }    
    
    event->Clean(); delete event;    
    
    if (iEvent%10 == 0 || iEvent == nEvent-1) {
      c1->cd(1); h_nvrt->Draw("B");
      c1->cd(2); h_nvrtu->Draw("B");
      c1->cd(3); h_nvrtc->Draw("B");
      c1->cd(4); h_nvrtb->Draw("B");
      c1->Modified();c1->Update();
      
      c2->cd(1); h_rwm ->Draw("B");
      c2->cd(1); h_rwmuc ->Draw("same");
      c2->cd(1); h_rwmu  ->Draw("same");
      c2->cd(2); h_ptm1 ->Draw("B");
      c2->cd(2); h_ptm1uc ->Draw("same");
      c2->cd(2); h_ptm1u  ->Draw("same");
      c2->cd(3); h_ptm2 ->Draw("B");
      c2->cd(3); h_ptm2uc ->Draw("same");
      c2->cd(3); h_ptm2u  ->Draw("same");
      c2->Modified();c2->Update();
      
      c3->cd(1); h_mvspu->Draw("BOX");
      c3->cd(2); h_mvspc->Draw("BOX");
      c3->cd(3); h_mvspb->Draw("BOX");
      c3->Modified();c3->Update();
    }    
  }

  forjet->Delete(); delete forjet;

  return 1;
}