// ----------------------------------------------------------------------------
// $Id: root_GoJetsTrk.C,v 1.7 2001/06/19 19:45:42 masako Exp $
// ----------------------------------------------------------------------------
// $Log: root_GoJetsTrk.C,v $
// Revision 1.7  2001/06/19 19:45:42  masako
// Remove cross section.
//
// Revision 1.6  2001/06/19 16:52:04  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.5  2001/05/16 17:47:42  masako
// Show the unit of cross section (fb)
//
// Revision 1.4  2001/05/11 21:52:38  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.3  2001/05/09 22:01:47  toshi
// Add $Id: root_GoJetsTrk.C,v 1.7 2001/06/19 19:45:42 masako Exp $.
//
//
// Read root data, and Run Fast MC interactively with Root
// Sample program using JetFinder (using Charged tracks)
//           
//                            04/26/2001 Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x root_GoJetsTrk.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x root_GoJetsTrk.C(300)
//


Int_t root_GoJetsTrk(Int_t nEvent=100) {

  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.

  // Now Using dev version Event Classes
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");

  LCDEvent event;

  // Input generator level root data
  Char_t* 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 *nJets   = new TH1F("nJets","# of Jets per Event (y=0.004)",20,0.,20);
  TH1F *nJets_1 = new TH1F("nJets_1","# of Jets per Event(y=0.006)",20,0.,20);
  TH1F *nJets_2 = new TH1F("nJets_2","# of Jets per Event(y=0.006, JADE)",20,0.,20);
  TH1F *JetP    = new TH1F("JetP","Jet momentum",100,0.,200);
  TH1F *JetPhi  = new TH1F("JetPhi","JetPhi Angle",100,0.,3.14);
  TH1F *JetCosTh= new TH1F("JetCosTh","Jet Cos Theta",50,-1.0,1.0);

  // Setup of JetFinder
  LCDJetFinder jetfind(0.004);  // Default.. DURHAM algorithm

  TClonesArray veclist("TLorentzVector",100);

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

    fmc.Doit();                    // Through FastMC
    
    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();
        
    // Event cut # track >= 2
    if (nTracks < 2) continue;

    veclist.Delete(); // Clean up veclist

    // Track loop    
    for (Int_t ntrk=0; ntrk<nTracks; ntrk++) {
      LCDTrack* trk = (LCDTrack*)event->Tracks()->At(ntrk);

      TVector3 tP3 = trk->GetMomentumVector(0.0);
      Double_t E_trk = TMath::Sqrt(tP3.Mag2() + 0.1396*0.1396); 
      // Assume pion mass
      TLorentzVector* vec4 = new(veclist[ntrk]) TLorentzVector();
      vec4.SetXYZT(tP3.X(),tP3.Y(),tP3.Z(),E_trk);
    }
    
//  Do Jet Finding
    jetfind.SetPartList(&veclist); // Input particle 4-vector list
    jetfind.SetDURHAM();          // Use DURHAM algorithm
    jetfind.SetYCut(0.004);       // Set YCut value
    jetfind.doFindJets();         // JetFindig
    
    Int_t numjets = jetfind.njets();
    nJets->Fill(numjets);

    for(Int_t j = 0; j < numjets; j++){

	TLorentzVector* Jet_vec = jetfind.jet4vec(j);
	TVector3        Pjet    = Jet_vec ->Vect();   // Jet 3 vector

	Double_t  Jet_mom = Pjet->Mag();
	Double_t  Jet_cos = (Pjet->Z())/Jet_mom;
	Double_t  Jet_phi = Pjet->Phi();

	JetP    ->Fill(Jet_mom);
	JetPhi  ->Fill(Jet_phi);
	JetCosTh->Fill(Jet_cos);
    }
    

//  Example 1  ... Do JetFindong with different Y-cut
    jetfind.SetYCut(0.0060);  // Set different YCut value
    jetfind.doFindJets();     // JetFinding

    Int_t numjets_1 = jetfind.njets();
    nJets_1->Fill(numjets_1);


//  Example 2  ... Do JetFindong with different Algorithm
    jetfind.SetJADE();      // Use JADE algorithm
    jetfind.doFindJets();   // JetFinding

    Int_t numjets_2 = jetfind.njets();
    nJets_2->Fill(numjets_2);

  }
  veclist.Delete();

  nJets   ->SetFillColor(4);
  nJets_1 ->SetFillColor(4);
  nJets_2 ->SetFillColor(4);
  JetP    ->SetFillColor(4);
  JetPhi  ->SetFillColor(4);
  JetCosTh->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 600,600);
  c1->Divide(2,3);
  c1->cd(1); nJets->Draw("B");
  c1->cd(2); nJets_1->Draw("B");
  c1->cd(3); nJets_2->Draw("B");
  c1->cd(4); JetP->Draw("B");
  c1->cd(5); JetPhi->Draw("B");
  c1->cd(6); JetCosTh->Draw("B");

  return iEvent;
}