// Run Fast MC interactively with Root
// read multi files
//
//     Aug 29 20000    M.Iwasaki
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x readMulti.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x readMulti.C(300)
//

#include <stdlib.h>

int readMultiF(int nEvent=200)
{

  gROOT->Reset();

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.
  gSystem->Load("libLCDEvent.so");
  gSystem->Load("libLCDRootAppsUtil.so");
  gSystem->Load("libLCDFastMC.so");

  Int_t   nSrc = 2;  // Number of source files

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

  // Detector parameter file
  char* geomFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/Small.par";

  // Tracker smearing parameter file
  char* smearFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/lookup_small.nobmcon";

  LCDreadStdFile* source = 0;
  LCDFastMC* fmc = new LCDFastMC(source, geomFile, smearFile, 0);
  LCDEvent*  event = 0;

  // 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) {
    source = new LCDreadStdFile(stdFile[iSrc]);
    printf("New source %s opened\n", stdFile[iSrc]);

    fmc->SetSource(source);

    // Event loop
    for (int iEvent = iEvent_tot;iEvent < nEvent; iEvent++) {
      
      event = new LCDEvent();
      event->Create();
      
      Int_t ierr = fmc->FetchEvent(event);
      if (!ierr) break;
      
      fmc->Cleanup();
      fmc->Doit(event); //Through FastMC
      
      // Look at CLusters. Can access them by the [] operator.
      int 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 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);
	}
	
      }

      iEvent_tot++;
      event->Clean(); delete event;

    }
    delete source; source = 0;
    iSrc++;    
  }

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

  TCanvas *c1 = new TCanvas("c1", "c1", 400, 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 1;
}