// ----------------------------------------------------------------------------
// $Id: stdhep_readMultF.C,v 1.5 2001/06/19 16:52:21 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: stdhep_readMultF.C,v $
// Revision 1.5  2001/06/19 16:52:21  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.4  2001/05/11 21:53:05  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.3  2001/05/10 23:04:03  toshi
// Add gSystem->Load("libLCDStdHEPUtil.so") according to changes LCDROOT ver3.2.
//
// Revision 1.2  2001/05/09 22:10:19  toshi
// Add $Id: stdhep_readMultF.C,v 1.5 2001/06/19 16:52:21 toshi Exp $.
//
//
// Run Fast MC interactively with Root
// read multi STDHEP files
//
//     04/26/2001 Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x stdhep_readMultiF.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x stdhep_readMultiF.C(300)
//

Int_t stdhep_readMultF(Int_t nEvent=200) {

  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("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDStdHEPUtil");

  Int_t   nSrc = 2;  // Number of source files

  Char_t* stdFile[2] = 
  {    
    "test1.dat", // File 1
    "test2.dat"  // File 2
  };

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

  Int_t iSrc = 0;  Int_t iEvent_tot = 0.;

  // Source file loop
  while (iSrc < nSrc) {
    LCDreadStdFile* source=new LCDreadStdFile(stdFile[iSrc],&event);
    printf("New source %s opened\n", stdFile[iSrc]);

    // Event loop
    for (Int_t iEvent = iEvent_tot ; iEvent < nEvent; iEvent++) {
      if (!source->GetEvent()) break; //Read stdhep data
      fmc.Doit();                     //Through FastMC
      
      // Look at CLusters. Can access them by the [] operator.
      Int_t nClusters = event.ClusterLst()->GetEntries();
      ClsMult->Fill(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
	  ClsE->Fill(eClus);
	}
      }
    
      // Look at tracks. Can access them by the [] operator.
      Int_t 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) {
	  TVector3 tP3 = trk->GetMomentum();
	  // Histogram energies
	  TrkE->Fill(TMath::Sqrt(tP3.Mag2() + 0.1396*0.1396));
	}
	
      }

      iEvent_tot++;
    }
    delete source;
    iSrc++;    
  }

  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_tot;
}