// ----------------------------------------------------------------------------
// $Id: pythia_GoLCDVTopl.C,v 1.2 2001/06/19 16:52:00 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: pythia_GoLCDVTopl.C,v $
// Revision 1.2  2001/06/19 16:52:00  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.1  2001/05/18 23:53:29  masako
// Create a new example file to run pythia and LCDVTopl
//
//
//  Example file to use LCDVTopl (translated from SLD ZVTOP) 
//  interactively with Root.
//
//  In this file, 1) event generation with Pythia (TPythia6)
//                2) through Fast MC
//                3) jet-finding anc call LCDVtopl
//
//  For example, generate e+e- --> tt-bar events
//
//                               May. 18 2001  M.Iwasaki
//
//   :: How to run ::
// 
//   1) type root
//   2)
//      root [0]  .x pythia_GoLCDVTpol.C
//
// For example, if you like to run 20 events, type like
//      root [0]  .x pythia_GoLCDVTpol.C(20)
//    


Int_t pythia_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("libPythia6");
  gSystem->Load("libEG");
  gSystem->Load("libEGPythia6");
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");
  gSystem->Load("libLCDGenUtil");

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

  double ECM  = 500.;
  int    iseed_pyth =  19518503;// pythia initial seed

  TPythia6 pyth6;

  //...t mass.
  pyth6.SetPMAS(6,1,175.);

  // Select gamma*/Z0 production process
  pyth6.SetMSEL(0);
  pyth6.SetMSUB(1,1);

  // Switch off all Z decay modes except tt-bar.
  int MDCY_23_2 = pyth6.GetMDCY(23,2);
  int MDCY_23_3 = pyth6.GetMDCY(23,3);

  for (int IDCY = MDCY_23_2; IDCY< MDCY_23_2+MDCY_23_3; IDCY++){
    int partID = pyth6.GetKFDP(IDCY,1);
    if (TMath::Abs(partID) != 6 && 
	pyth6.GetMDME(IDCY,1)> 0) pyth6.SetMDME(IDCY,1,0);
  }

  pyth6.Pyinit("CMS","E-","E+",ECM);
 
  // Change pythia seed  
  pyth6.SetMRPY(1,iseed_pyth); 

  //Pyjet -> McPart converter
  LCDPyj2McPart py2mc;

  cout<<"End of Generator Setup!"<<endl;
  //------------< End Generator setup>-------------------

  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* 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++) {

    pyth6.GenerateEvent();                //Event Generation with pythia
    //if (iEvent < 3) pyth6.Pylist(1);      //Print Event list for example

    event.MCparticles()->Delete();        //prepare for making McPart
    py2mc.MakeMcPart(event.MCparticles());//Make McPart
 
    fmc.Doit();                           //Through FastMC

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

  pyth6.Pystat(1); // Pythia summary

  return iEvent;
}