// ----------------------------------------------------------------------------
// $Id: pandpyth_GoLCDVTopl.C,v 1.4 2001/06/19 16:51:57 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: pandpyth_GoLCDVTopl.C,v $
// Revision 1.4  2001/06/19 16:51:57  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.3  2001/05/11 21:52:29  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.2  2001/05/09 21:56:28  toshi
// Add CVS tag $Id: pandpyth_GoLCDVTopl.C,v 1.4 2001/06/19 16:51:57 toshi Exp $.
//
//
// Run Pandora_Pythia and Fast MC interactively with Root
// Sample program using LCDVTopl (translated from SLD ZVTOP) 
//                                              04/26/2001 Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x pandpyth_GoLCDVTopl.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x pandpyth_GoLCDVTopl.C(300)
//


Int_t pandpyth_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("libPandora");
  gSystem->Load("libPandoraProc");
  gSystem->Load("libPythia6");
  gSystem->Load("libEG");
  gSystem->Load("libEGPythia6");
  gSystem->Load("libTauola");
  gSystem->Load("libPandPyth");
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");
  gSystem->Load("libLCDGenUtil");

  //------------< Setup Generator (Pandora-Pythia) >------------------------
  cout<<"Start Generator Setup!"<<endl;

#include "pandpythconsts.h"

  Int_t iseed_pan = 7;  // default pandora initial seed

  /*    define the pandora event selection in this space     */
  
  double ECM  = 91.26;
  double PolE =  0.0; // Polarization Electron
  double PolP =  0.0; // Polarization Positron

  ebeam b1(ECM/2.0, PolE, electron, electron);
  b1.setup(NLC500);
  b1.ISRoff();            // Turn off initial-state radiation
  b1.beamstrahlungoff();  // Turn off beamstrahlung


  ebeam b2(ECM/2.0, PolP, positron, positron);
  b2.setup(NLC500);
  b2.ISRoff();
  b2.beamstrahlungoff();
  
  double ECMmin = 20.; // Higgs Mass
  eetoqqbar pr(ECMmin);

  //For example:: Z0 decays to anything
  //const Int_t invisibleOnly = 557;
  //const Int_t bOnly = 553;
  //const Int_t cOnly = 552;
  //const Int_t allStates = 546;
  //pr.onlyDecay(allStates); 

  pandora P(b1,b2,pr);
  
  pandorarun  PR(&P, epluseminus, ECM, nEvent, iseed_pan);
 
  //----- For PYTHIA ------
  //Or you can also change PYTHIA seed here.
  int iseed_pyth = 22846; PR.newseed(iseed_pyth); // Change Seed

  //PR.partonshoweroff();   // Turn off parton shower
  //PR.fragmentationoff();  // Turn off fragmentation
  //PR.FSRoff();            // Turn off final-state QCD QED radiation

  /*       end of definition                      */

  PR.initialize(0);
  
  /*
  // If you want/don't want to decay some specific particle(s),
  // you can call PartDecayOn/PartDecayOff here.   

  // For example: Do not allow to decay K0S and Lambda
  int ipart1 = 310;  // K0S
  int ipart2 = 3122; // Lambda
  PR.PartDecayOff(ipart1);
  PR.PartDecayOff(ipart2);
  */

  cout<<"End of Generator Setup!"<<endl;

  //Pyjet -> McPart converter
  LCDPyj2McPart py2mc;

  //------------< End Generator (Pandora-Pythia) setup>-------------------

  //Physics info..
  LCDEvent event;

  //-------<Fast Simulation Set up>---------

  //---> Detector Geometry 
  char* geomFile = "Silicon.par"; // Silicon 
  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";
  TString s_smearFile(parfile_dir);
  s_smearFile += smearFile;

  //---> FastMC
  LCDFastMC fmc(&gp,s_smearFile.Data(),&event); 

  // 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,50.);

  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 iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {

    PR.getevent();//Event Generate
    event.MCparticles()->Delete();//prepare for making McPart
    py2mc.MakeMcPart(event.MCparticles());//Make McPart
 
    fmc.Doit();//do fast detector simulation...

    // 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.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;
    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();
      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) {
      cout<<"iEvent = "<<iEvent+1<<" was done..."<<endl;
    }
  }
  forjet.Delete();
  topl.Clean();
  event.Clean();

  return iEvent;
}