#ifdef __CINT__
{
#endif
//////////////////////////////////////////////////
// tof.C : Root Macro
//
//  Example of the Hall D Time-of-Flight
//
// The reaction gamma p -> p K- K+ pi- pi+ was
// generated using "genr8" and processed through
// HDFast to create the data file test1.rdt 
//  
//*-- Author :    Paul Eugenio 1-Jan-99
//*-- CMZ : PME 14-Jan-99
//////////////////////////////////////////////////

//
//Begin_Html
/*
<a href= "smearedTOF.gif" >
Here is the output using smeared momentum in the calculation</a>
<a href= "tracedTOF.gif" >
Here is the output for using the trace point momentum   </a>

*/
//End_Html


#include <iostream.h>
#include <string.h>

#ifndef __CINT__
#include <stdlib.h>
#include "TROOT.h"
#include "TApplication.h"
#include "TFile.h"
#include "TRandom.h"
#include "TTree.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TMCFastHepEvt.h"
#include "TMCFastTOF.h"
#include "TMCFastOfflineTrack.h"
#include "TMCFastCalorimeter.h"
#include "TLGDsmears.h"
#include "TMCesr.h"

extern void InitGui();
VoidFuncPtr_t initfuncs[] = { InitGui, 0 };
TROOT root("GUI", "GUI test environement", initfuncs);
int main(int argc, char **argv){

#endif
 
#ifndef __CINT__
   TApplication theApp("App", &argc, argv);
#endif

   //gROOT->Reset();

#ifdef __CINT__
   gSystem->Load("./libTMCFast.so"); 
   cout<<"gSystem->Loading\n";
#endif
   //   Create a timer object to benchmark this loop
   TStopwatch timer;
   timer.Start();
   
   // This file was generated by HDFast
   
   TFile f("p_K-pi+pi-K+:B8M1.7W0.25.rdt");
   
   
   TTree *tree = (TTree*)f.Get("T");
   
   TBranch *b1 = tree->GetBranch("hepevt");
   TBranch *b2 = tree->GetBranch("tof_trace");
   TBranch *b3 = tree->GetBranch("offtrk");
   TBranch *b4 = tree->GetBranch("bcal");
   TBranch *b5 = tree->GetBranch("lgdSmears");
   //TBranch *b6 = tree->GetBranch("fHepParticle");
   
   TMCFastHepEvt *hepevt= new TMCFastHepEvt();
   //TClonesArray *hepparts= hepevt->GetHepParticles();
   TMCFastTOF *trace= new TMCFastTOF(); 
   TMCFastOfflineTrack *offtrk= new TMCFastOfflineTrack();
   TMCFastCalorimeter *bcal = new TMCFastCalorimeter();
   TLGDsmears *lgdSmears = new TLGDsmears();
   
   
   b1->SetAddress(&hepevt);
   b2->SetAddress(&trace);
   b3->SetAddress(&offtrk);
   b4->SetAddress(&bcal);
   b5->SetAddress(&lgdSmears); 
   //b6->SetAddress(&hepparts); 
   
   // book histograms for storing the raw information:
   
   TH1F *h1 = new TH1F("h1", "CTOF raw mass pi",     100, 0.0, 0.3);
   TH1F *h2 = new TH1F("h2", "CTOF raw mass K ",     100, 0.3, 0.6);
   TH1F *h3 = new TH1F("h3", "CTOF raw mass proton", 100, 0.8, 1.2);
   TH1F *h4 = new TH1F("h4", "CTOF raw mass neutron",100, 0.8, 1.2);
   
   TH1F *h5 = new TH1F("h5", "FTOF raw mass pi",     100, 0.0, 0.3);
   TH1F *h6 = new TH1F("h6", "FTOF raw mass K ",     100, 0.2, 0.6);
   TH1F *h7 = new TH1F("h7", "FTOF raw mass neutron",100, 0.8, 1.2);
   TH1F *h8 = new TH1F("h8", "FTOF raw mass proton", 100, 0.8, 1.2);
   
   // book histograms for storing the tracked information:
   
   TH1F *h11 = new TH1F("h11", "CTOF raw mass pi",     100, 0.0, 0.3);
   TH1F *h12 = new TH1F("h12", "CTOF raw mass K ",     100, 0.3, 0.6);
   TH1F *h13 = new TH1F("h13", "CTOF raw mass proton", 100, 0.8, 1.2);
   TH1F *h14 = new TH1F("h14", "CTOF raw mass neutron",100, 0.8, 1.2);
   
   TH1F *h15 = new TH1F("h15", "FTOF raw mass pi",     100, 0.0, 0.3);
   TH1F *h16 = new TH1F("h16", "FTOF raw mass K ",     100, 0.2, 0.6);
   TH1F *h17 = new TH1F("h17", "FTOF raw mass neutron",100, 0.8, 1.2);
   TH1F *h18 = new TH1F("h18", "FTOF raw mass proton", 100, 0.8, 1.2);
   
   // Now we want to loop over all the events in the data file:
   //TMCesr *esr;
   
   Int_t nentries = (Int_t)tree->GetEntries();
   
   
   cout<<"nentries = "<< nentries <<endl;
   for (Int_t ev = 1; ev < nentries ; ev++) {
     
     tree->GetEvent(ev);        //read  event in memory
     //b1->GetEvent(ev);
     //b2->GetEvent(ev);
     //b3->GetEvent(ev);
     // Now we want to loop over all the particles inside the 
     // event.
     
     
     // Define an iterator "next" for the TcloneArray of TMCFastHepParticle:
     //hepevt->Print(&cout);
     //esr = new TMCesr(*hepevt,*offtrk,*lgdSmears,*bcal);
     
     TIter next (hepevt->GetHepParticles());
     //TIter next (esr->GetParticles());
     // Loop over all the particles:
     
     TMCFastHepParticle *particle;
     static Int_t pdgID[20],hepindex[20],np;
     
     if(ev==1){
       np=0;
       while (( particle  = ((TMCFastHepParticle *)next()))){
	 // Get the PDG code, and the hepindex of the particles:
	 pdgID[np]    = particle->GetIdhep();
	 hepindex[np++] = particle->GetIndex();
       }
     }
     // Now decide what to do based on the particle type:
     for(Int_t i=0;i<np;i++)
       switch(pdgID[i]){
       case 22:               // 22 = Photon
	 break;
       case 211:              //  211 = pi+
       case -211:             // -211 = pi-
	 h1->Fill(trace->CTOFmass(hepindex[i]));
	 h5->Fill(trace->FTOFmass(hepindex[i]));
	 h11->Fill(trace->CTOFmass(hepindex[i],*offtrk));
	 h15->Fill(trace->FTOFmass(hepindex[i],*offtrk));
	 break;
       case 321:             //  321  = K+
       case -321:            // -321  = K-
	 h2->Fill(trace->CTOFmass(hepindex[i]));
	 h6->Fill(trace->FTOFmass(hepindex[i]));
	 h12->Fill(trace->CTOFmass(hepindex[i],*offtrk));
	 h16->Fill(trace->FTOFmass(hepindex[i],*offtrk));
	 break;
       case 2212:            //  2212  proton
	 h3->Fill(trace->CTOFmass(hepindex[i]));
	 h7->Fill(trace->FTOFmass(hepindex[i]));
	 h13->Fill(trace->CTOFmass(hepindex[i],*offtrk));
	 h17->Fill(trace->FTOFmass(hepindex[i],*offtrk));
	 break;
       case 2112:            // 2112  neutron
	 h4->Fill(trace->CTOFmass(hepindex[i]));
	 h8->Fill(trace->FTOFmass(hepindex[i]));
	 h14->Fill(trace->CTOFmass(hepindex[i],*offtrk));
	 h18->Fill(trace->FTOFmass(hepindex[i],*offtrk));
	 break;
       default:
	 cerr << "I could not identify the particle";
	 break;
	 
       }
     
     
     // Get the next event
     //  esr->Clear();
     //delete esr;
     hepevt->Clear();
     trace->Clear();
     offtrk->Clear();
     bcal->Clear();
     lgdSmears->Clear(); 
   }
   
   // now create some pads to put the histo's
   
   TCanvas *c1 = new TCanvas("c1","Halld Time-of-Flight",100,10,700,900);
   c1->SetFillColor(18);
   
TPad *pad1 = new TPad("pad1","The pad with the function",0.05,0.66,0.45,0.96,21);
TPad *   pad2 = new TPad("pad2","The pad with the function",0.55,0.66,0.95,0.96,21);
TPad *   pad3 = new TPad("pad3","The pad with the histogram",0.05,0.34,0.45,0.64,21);
TPad *   pad4 = new TPad("pad4","The pad with the histogram",0.55,0.34,0.95,0.64,21);
TPad *   pad5 = new TPad("pad5","The pad with the histogram",0.05,0.02,0.45,0.32,21);
TPad *   pad6 = new TPad("pad6","The pad with the histogram",0.55,0.02,0.95,0.32,21);
 
 pad1->Draw();
 pad2->Draw();
 pad3->Draw();
 pad4->Draw();
 pad5->Draw();
 pad6->Draw();
 

 // Now draw the histograms
 
 pad1->cd();
 h1->Draw();
 
 pad2->cd();
 h5->Draw();
 
 pad3->cd();
 h2->Draw();

 pad4->cd();
 h6->Draw();
 
 pad5->cd();
 h3->Draw();
 
 pad6->cd();
 h7->Draw();
 
 //  Stop timer and print results
 timer.Stop();
 Double_t rtime = timer.RealTime();
 Double_t ctime = timer.CpuTime();
 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

#ifndef __CINT__
 // Here we don't return from the eventloop. "Exit ROOT" will quit the app.
 theApp.Run();
#endif

 return 0;
 }