//// root macro to analyze the simulation output
//#include <map>
//#include <string>
//#include <iostream>
//#include <vector>
//
//// #include <TStopwatch.h>
//#include <TFile.h>
//#include <TTree.h>
//#include <TGeoManager.h>
//#include <TClonesArray.h>
//#include <TH1D.h>
//#include <TROOT.h> // for global root variables
//#include <TH2D.h>
//#include <TVector3.h>
//#include <TPad.h>
//#include <TCanvas.h>
//#include <TStopwatch.h>
//#include <THStack.h>
//#include <TProfile.h>
//#include <TProfile2D.h>
//#include <TColor.h>
//
////#include "$VMCWORKDIR/base/FairRadLenPoint.h"
//#include "/Users/ralfk/Panda/pandaroot/base/FairRadLenPoint.h"
void materialana(int nEvents = 1000, bool verbose = false)
{
  
  // -----  Load libraries   ------------------------------------------------
  gROOT->Macro("$VMCWORKDIR/example/gconfig/basiclibs.C");  
  gSystem->Load("libProof");
  gSystem->Load("libFairTools");
  gSystem->Load("libFairDB");
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  //gSystem->Load("libGen");
  //gSystem->Load("libfield");
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------
  
  // Fetch data arrays
  std::string inFile = "eic_rad_length.mc.root";
  TFile* f = new TFile(inFile.c_str()); // the sim file you want to analyse
  TTree *t=(TTree *) f->Get("cbmsim");
  
  TClonesArray* mc_array=new TClonesArray("FairMCTrack");
  t->SetBranchAddress("MCTrack",&mc_array);//Branch names
  TClonesArray* rad_array=new TClonesArray("FairRadLenPoint");
  t->SetBranchAddress("RadLen",&rad_array);
  
  // Setup Histograms
  Double_t boxwidth=300;
  TH1D* hisx = new TH1D("hX","X",100,-boxwidth,boxwidth);
  TH1D* hisy = new TH1D("hY","Y",100,-boxwidth,boxwidth);
  TH1D* hisz = new TH1D("hZ","Z",100,-boxwidth,boxwidth);
  TProfile* thetaprofile = new TProfile("thetaprof","Radiaion length (#theta)",180,0.,TMath::Pi());
  TProfile* phiprofile = new TProfile("phiprof","Radiaion length (#phi)",180,-1.*TMath::Pi(),TMath::Pi());
  TProfile2D* histhephi = new TProfile2D("hThePhi","",90,0.2,2.8,90,-1.*TMath::Pi(),TMath::Pi());
  histhephi->SetTitle("Radiation length map #theta & #phi;#theta;#phi");
  
  //TString detname,volname;
  double radlen=0.,effradl=0.,effradlsum=0.,theta=0.,phi=0.;
  TVector3 in, out, dist, point;
  TVector3 mom;
  
  for (Int_t event=0; event<nEvents && event<t->GetEntriesFast(); event++)
  {
    t->GetEntry(event);
    if(verbose) cout<<"Event No "<<event<<endl;
    else if (!(event%100)) cout <<"Event No "<<event<<endl;
    
    for (Int_t trackno=0; trackno<mc_array->GetEntriesFast();trackno++){
      FairMCTrack* aTrack = (FairMCTrack*)mc_array->At(trackno);
      aTrack->GetMomentum(mom); // write momentum to mom
      theta = mom.Theta(); phi = mom.Phi();
      effradl=0.; effradlsum=0.;
      for (Int_t k=0; k<rad_array->GetEntriesFast(); k++){
        FairRadLenPoint* radpoint = (FairRadLenPoint*)rad_array->At(k);
        if (radpoint->GetTrackID() != trackno) continue; 
        radlen = radpoint->GetRadLength();
        in = radpoint->GetPosition();
        out = radpoint->GetPositionOut();
        dist = in - out;
        point.SetXYZ(0.5*(in.x()+out.x()),0.5*(in.y()+out.y()),0.5*(in.z()+out.z()));
        if(in.Mag() < 0.02) continue; // cut target
        effradl = dist.Mag()/radlen;
        effradlsum += effradl;
        if(event < 1) // testing output
          printf("Test: phi(%f), theta(%f), radlen(%g), effradl(%g)\n",phi,theta,radlen,effradl);
        hisx->Fill(in.X(),effradl);
        hisy->Fill(in.Y(),effradl);
        hisz->Fill(in.Z(),effradl);
      }//radpoints 
      thetaprofile->Fill(theta,effradlsum);
      phiprofile->Fill(phi,effradlsum);
      histhephi->Fill(theta,phi,effradlsum);
      
    }//tracks
    
  }// end for event
  
  
  // Draw the Histos
  Int_t a=3,b=2;
  int resol = 250;
  TCanvas* can1 = new TCanvas("can1","MCHit view in MVD",0,0,a*resol,b*resol);
  can1->Divide(a,b);
  
  can1->cd(1);
  hisx->SetFillColor(38);
  hisx->DrawCopy();
  
  can1->cd(2);
  hisy->SetFillColor(45);
  hisy->DrawCopy();
  
  can1->cd(3);
  hisz->SetFillColor(30);
  hisz->DrawCopy();
  can1->cd(4);
  gPad->SetLogy();
  thetaprofile->SetFillColor(30);
  thetaprofile->SetLineColor(30);
  thetaprofile->Draw("hist");
  can1->cd(5);
  gPad->SetLogy();
  phiprofile->SetFillColor(45);
  phiprofile->SetLineColor(45);
  phiprofile->Draw("hist");
  
  can1->cd(6);
  histhephi->SetStats(false);
  histhephi->DrawCopy("colz");
  
  // -----   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;
  // ------------------------------------------------------------------------
  
}