// ----------------------------------------------------------------------------
// $Id: root_GoFastMC.C,v 1.8 2001/09/20 01:45:38 toshi Exp $
// ----------------------------------------------------------------------------
//
// Run Fast MC interactively with Root
//
//     Apr 25 20001    Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x root_GoFastMC.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x root_GoFastMC.C(300)
//


#include <stdlib.h>

Int_t root_GoFastMC(int nEvent=100) {

  gROOT->Reset();

  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");

  LCDEvent event; //Create Physics event container first

  // Open a root file (generator output)
  Char_t* rootFile = "test.root";
  LCDreadRootFile source(rootFile,&event);

  // Open detector parameter file
  Char_t* geomFile = 
    "SDMar01.par"; // SDMar01
  // "LDMar01.par"; // LDMar01
  //"Large.par";   // Large
  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();

  // Define a histogram or two
  TH1F *ClsMult = new TH1F("ClsMult","Cluster Multiplicity",50, 0.0, 200.0);
  TH1F *ClsE    = new TH1F("ClsE","Cluster Energy",50, 0.0, 50.0);
  TH1F *TrkMult = new TH1F("TrkMult","Track Multiplicity",50, 0.0, 100.0);
  TH1F *TrkE    = new TH1F("TrkE","Track Energy",50, 0.0, 50.0);

  // Event loop
  Int_t iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {

    if (!source.GetEvent()) break; // Read a given ROOT file.

    fmc.Doit();                    // Through Fast simulation.

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

    TLorentzVector part_mom(0.,0.,0.,0.);
    // 
    TObjArray* McPart = event.MCparticles();
    Int_t nMCpart = McPart->GetEntries();
    for (Int_t Imc=1;Imc<nMCpart;Imc++) {
      LCDMcPart *mcpart = (LCDMcPart *) McPart->At(Imc);
      Int_t  part_st = mcpart->GetStatus();
      Int_t  part_id = mcpart->GetParticleID();
      
      part_mom = mcpart->Get4Momentum();

      /*
      printf("%i ID = %i status = %i momentum ( %f, %f, %f ) E %f \n",
	     Imc, part_id, part_st, part_mom.Px(),
	     part_mom.Py(),part_mom.Pz(),part_mom.E());
      */
    }

    // Look at Clusters. Can access them by the [] operator.
    int nClusters = event.ClusterLst()->GetEntries();
    ClsMult->Fill(nClusters);

    //printf("nClusters %i\n", nClusters);

    for (Int_t iClus = 0; iClus < nClusters; iClus++) {
      LCDCluster* cls = (LCDCluster*)event.ClusterLst()->At(iClus);

      if (cls) {
	Double_t     eClus     = cls->GetEnergy();	
	// Histogram energies
	//printf("energy = %f\n", eClus);
	ClsE->Fill(eClus);
      }

    }
    

    // Look at tracks. Can access them by the [] operator.
    int nTracks = event.Tracks()->GetEntries();
    TrkMult->Fill(nTracks);

    // Track loop    
    for (Int_t ntrk=0; ntrk<nTracks; ntrk++) {
      LCDTrack* trk = (LCDTrack*)event.Tracks()->At(ntrk);
      
      if(trk) {
	Double_t* tP3 = trk->GetMomentum();
	Double_t  E_trk = 0.1396 * 0.1396; // Assume pion mass
	for (int i=0; i<3; i++) { E_trk += (*(tP3+i))*(*(tP3+i)); }
	E_trk = sqrt(E_trk);

	// Histogram energies
	TrkE->Fill(E_trk);
      }

    }

  }

  ClsMult->SetFillColor(4);
  ClsE->SetFillColor(4);
  TrkMult->SetFillColor(4);
  TrkE->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 600, 600);
  c1->Divide(2,2);
  c1->cd(1); ClsMult->Draw("B");
  c1->cd(2); ClsE->Draw("B");
  c1->cd(3); TrkMult->Draw("B");
  c1->cd(4); TrkE->Draw("B");

  return iEvent;
}