// ----------------------------------------------------------------------------
// $Id: root_GoLCDVTopl.C,v 1.8 2001/06/22 16:39:58 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: root_GoLCDVTopl.C,v $
// Revision 1.8  2001/06/22 16:39:58  toshi
// Remove cout due to Windows limitation.
//
// Revision 1.7  2001/06/19 19:45:43  masako
// Remove cross section.
//
// Revision 1.6  2001/06/19 16:52:05  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.5  2001/05/16 17:48:06  masako
// Show the unit of cross section (fb)
//
// Revision 1.4  2001/05/11 21:52:40  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.3  2001/05/09 22:02:32  toshi
// Add $Id: root_GoLCDVTopl.C,v 1.8 2001/06/22 16:39:58 toshi Exp $.
//
//
// Run Fast MC interactively with Root
// Sample program using LCDVTopl (translated 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 root_GoLCDVTopl(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* rootFile = "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_nvtvrt = new TH1F("h_nvrt","N of LCDVToplVRT"       ,10,0.0,10.0);
  TH1F* h_rwm    = new TH1F("h_rwm" ,"Raw Vrt Mass"           ,40,0.0,10.0);
  TH1F* h_ptm1   = new TH1F("h_ptm1" ,"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);
  TH2F* h_mvsp   = new TH2F("h_mvsp" ,"Pt Mass vs Vertex P uds" ,
			   20,0.0,2.0,25,0.0,100.);

  h_nvrt  ->SetFillColor(4);
  h_rwm   ->SetFillColor(4);
  h_ptm1  ->SetFillColor(4);
  h_ptm2  ->SetFillColor(4);
  h_mvsp  ->SetFillColor(4);


  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
      

    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;
      for (int ipart =0; ipart<nGTracks; ipart++){
	if( jetfind.GetPartIndexAt(ipart) == iJet){
	  trkList[iTrk2++] = ipart;
	}
      }
      jet3vec = jetfind.jet4vec(iJet).Vect();

      topl.SetDebugLevel(0);
      topl.SetUpLCDVTopl(&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 (nVRT <= 1) continue;

      p4vrt=topl.GetP4Vertex();
      
      rawMass  = topl.GetRawMass();      

      if (rawMass <= 0.) continue;
      
      ptMass   = topl.GetPtMass();
      massTag  = topl.GetMassTag();
      
      h_rwm ->Fill(rawMass);
      h_ptm1->Fill(ptMass);
      h_ptm2->Fill(massTag);
      
      if (massTag<2.) h_mvsp->Fill(massTag,p4vrt.Rho());

    }    

    delete trkList; 

  canvasdraw:    
    if (iEvent%100 == 0 || iEvent == nEvent-1) {
      c1->cd(1); h_nvrt->Draw("B");
      c1->cd(2); h_mvsp->Draw("BOX");
      c1->Modified();c1->Update();
      
      c2->cd(1); h_rwm ->Draw("B");
      c2->cd(2); h_ptm1 ->Draw("B");
      c2->cd(3); h_ptm2 ->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;
}