// ----------------------------------------------------------------------------
// $Id: stdhep_GoLCDVToplGhost.C,v 1.1 2001/12/05 18:28:04 toshi Exp $
// ----------------------------------------------------------------------------
//
// Run Fast MC interactively with Root
// Sample program using LCDVToplGhost (translated from SLD ZVTOP3) 
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x root_GoLCDVToplGhost.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x root_GoLCDVToplGhost.C(300)
//


#include <stdlib.h>

Int_t stdhep_GoLCDVToplGhost(Int_t nEvent=300){
  
  gROOT->Reset();
  gStyle->SetOptStat(111111);
  
  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;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");
  gSystem->Load("libLCDStdHEPUtil");

  LCDEvent event;

  // Input generator level root data
  Char_t* stdFile = "test.dat";
  LCDreadStdFile source(stdFile,&event);

  // Open detector parameter file
  Char_t* geomFile = 
    "SDMar01.par"; // SD
  //"LDMar01.par"; // LD
  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
  //gp.SetCoilinHAD();        // Set HAD inside the Coil (This case
                              // you need to change some detector size also..)

  // Tracking smearing parameter file
  Char_t* smearFile = 
    "lookup_silicon.nobmcon";
  //"lookup_large.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();
  // fmc.SetIPSigmaX(0.0);
  // fmc.SetIPSigmaY(0.0);
  // fmc.SetIPSigmaZ(0.0);


  // Define histograms
  TH1F* h_gang   = new TH1F("h_gang","Angle (gaxi and bmeson)",40,0.0, 0.4);
  TH1F* h_jang   = new TH1F("h_jang","Angle (jet and bmeson)" ,40,0.0, 0.4);
  TH1F* h_gwid   = new TH1F("h_gwid","Ghost Axis width"       ,50,0.0, 0.5);
  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_gang  ->SetFillColor(kBlue);
  h_jang  ->SetFillColor(kBlue);
  h_gwid  ->SetFillColor(kBlue);
  h_nvrt  ->SetFillColor(kBlue);
  h_rwm   ->SetFillColor(kBlue);
  h_ptm1  ->SetFillColor(kBlue);
  h_ptm2  ->SetFillColor(kBlue);
  h_mvsp  ->SetFillColor(kBlue);


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

  // Setup of JetFinder
  LCDJetFinder jetfind(0.004);  
  TClonesArray forjet("TLorentzVector",1000);
  
  // Setup of Topological vertex finder.
  LCDVToplGhost topl;
  topl.SetDebugLevel(0);
  // topl.SetDebugLevel(43);

  // Track cut value
  Double_t cosCut=0.9, pCut=0.2;

  // Event loop
  Int_t    ierr;
  Int_t    iEvent;
  // for (iEvent = 0; iEvent < nEvent; iEvent++) {
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    printf("iEvent=%d\n",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;
    TVector3  xtrk;
    Int_t itrk=0;
    for (itrk=0; itrk<nTracks; itrk++) {
      trk = (LCDTrack*)(event->Tracks()->At(itrk));
      ptrk=trk->GetMomentumVector(0.0);
      xtrk=trk->GetPositionVector(0.0);
      if (TMath::Abs(ptrk.CosTheta()) > cosCut ) continue;
      if (ptrk.Mag() < pCut) continue;
      if (xtrk.Perp() > 0.3) continue;
      if (TMath::Abs(xtrk.Z()) > 0.3) continue;
      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(2);        // 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
    for (Int_t iJet=0 ; iJet < nJet ; iJet++) {

      Int_t nTrks = jetfind.nParticlesPerJet(iJet);

      if (nTrks <= 2) continue; 

      jet3vec = jetfind.jet4vec(iJet).Vect();

      topl.SetUp(ipv,bfield,jet3vec);

      //Add track information into LCDVToplGhost and find b-meson
      LCDMcPart* bmcpart=0;
      LCDMcPart* tmcpart=0;
      for (Int_t ipart =0; ipart < nGTracks; ipart++){
	if( jetfind.GetPartIndexAt(ipart) == iJet){
	  LCDTrack *trk=(LCDTrack*)gtrk_list.At(ipart);
	  Double_t *tkpar=trk->GetTrackParameters();
	  Double_t *e_tkpar=trk->GetErrorMatrix();
	  Double_t  chrg=trk->GetCharge();
	  Double_t  mfld=trk->GetMagneticField();
	  Double_t  mcid=trk->GetParticle();
	  Bool_t    f_btrk=kFALSE;
	  if (mcid > -1) { // is this b-track?
	    tmcpart = (LCDMcPart*)(event.MCparticles()->At(mcid));
	    while(tmcpart) {
	      Int_t apid=TMath::Abs(tmcpart->GetParticleID());
	      switch(apid) {
	      case  511: // Bd
	      case  521: // Bu
	      case  531: // Bs
	      case  541: // Bc
		f_btrk=kTRUE;
		bmcpart = tmcpart;
		break;
	      default:
		break;
	      }
	      if (f_btrk) break;
	      Int_t parent_id=tmcpart->GetParentIdx();
	      if (parent_id < 1) break;
	      tmcpart = (LCDMcPart*)event.MCparticles()->At(parent_id);
	    }
	  }
	  // if (f_btrk) {
	    topl.AddTrack(ipart,tkpar,e_tkpar,chrg,mfld,mcid,ipv);
	    // }
	}
      }

      if (bmcpart == 0) continue;

      topl.FindGhostTrack();
      TVector3* gaxi=topl.GetGhostAxis();
      Double_t gang=gaxi->Angle((*bmcpart->Get4MomentumPtr()).Vect());
      Double_t jang=jet3vec.Angle((*bmcpart->Get4MomentumPtr()).Vect());
      h_gang->Fill(gang);
      h_jang->Fill(jang);
      h_gwid->Fill(topl.GetGhostAxisWidth());

      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());
    }    

  canvasdraw:    
    if (iEvent%10 == 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();
      
      c3->cd(1); h_gang->Draw("B");
      c3->cd(2); h_jang->Draw("B");
      c3->cd(3); h_gwid->Draw("B");
      c3->Modified();c3->Update();
      
    }    

    if (iEvent%1000 == 999 || iEvent == nEvent-1) {
      printf("iEvent = %d  was done...\n",iEvent+1);
    }
    if (!ierr) break;
  }
  forjet.Delete();
  topl.Clean();
  event.Clean();
  printf("topl Debug Level=%d\n",topl.GetDebugLevel());

  return iEvent;
}