// root macro to analyze the simulation output
//void convertMCPoints()
{

  // -----  Load libraries   ------------------------------------------------
//``gSystem->Load("fstream.h");
  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
  
  gSystem->Load("libSciT");
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  gStyle->SetPalette(1);
  
  TFile* f = new TFile("hit_scit.root"); // the sim file you want to analyse
  TTree *t=(TTree *) f->Get("cbmsim") ;
  t->AddFriend("cbmsim","test.root");
  
  TClonesArray* hit_array=new TClonesArray("PndSciTHit");
  t->SetBranchAddress("SciTHit",&hit_array);//Branch names

  TClonesArray* point_array=new TClonesArray("PndSciTPoint");
  t->SetBranchAddress("SciTPoint",&point_array);//Branch names


  TClonesArray* mc_array=new TClonesArray("PndMCTrack");
  t->SetBranchAddress("MCTrack",&mc_array);//Branch names

   TFile *out = TFile::Open("pid_HitOpt.root","RECREATE");


  // histos
  //TH2D* hisxy = new TH2D("hisxy","MC Points, xy view",200,-70.,70.,200,-70.,70.);
  //TH2D* hisrz = new TH2D("hisrz"," MC Points, rz view",200,-70.,70.,200,0,70.);
  TH1D* hisdt = new TH1D("hisdt"," time (ns)",100,2.1,2.6);
  hisdt->SetTitle("Smearing of the timing measurement");
  TH1D* hist = new TH1D("hist"," time (ns)",100,2.1,2.6);
  hist->SetTitle("MC timing measurement");
  

  //TH2D* hisdedx =new TH2D("hisdedx","",200,0.,0.04,200,0.,1.0);
  //hisdedx->SetTitle("HYP MC Points, dE/dx(P);p/GeV / cm;dE/dx / (GeV/cm)");



  
  //TH1D* hismom = new TH1D("hisximom"," MC Points, momentum",100,0.,1.5);
 

  int  nEvents = 100;
  bool verbose = false;
  int mcpdg;
  
  TVector3 vecs,veco;
  TVector3 vecFront,vecS,vecP;
  Double_t dx,dt,p,dEdX,time;
  
  TString detname;

  for (Int_t j=0; j<t->GetEntriesFast(); j++)
  {
    t->GetEntry(j);
    
    for (Int_t i=0; i<hit_array->GetEntriesFast(); i++)
    {
      dt= 0.;
      time=0.;
      mcpdg = -1;
      
      
	
      if(verbose) cout<<"Hit No "<<i<<endl;
      PndSciTHit *hit=(PndSciTHit*)hit_array->At(i);
      //int mcpdg = -1;
      
       PndSciTPoint *point=(PndSciTPoint*)point_array->At(hit->GetRefIndex());
       if(point==0)continue;
       if(point->GetTrackID()!=0)continue;
       //looking for primary particles (protons)
      
       PndMCTrack *mctruth = (PndMCTrack*)mc_array->At(hit->GetTrackID());
       mcpdg = mctruth->GetPdgCode();
       if(verbose)cout<<"mcpdg="<<mcpdg<<endl;

      dt=hit->GetTime();

      time = point->GetTime();
      
      
      if(dt!=0.)  hisdt->Fill(dt);
     
     if(time!=0.)  hist->Fill(time);

    }//end for i (points in event)
  }// end for j (events)



out->cd();

// hisxy->Write();
// hisrz->Write();
// hisdedx->Write();
hisdt->Write();
hist->Write();


// hismom->Write();
out->Save();


  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  //cout << "Output file is "    << outFile << endl;
  //cout << "Parameter file is " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout << endl;
  // ------------------------------------------------------------------------

}