// root macro to analyze BoxGen output
{

  // -----  Load libraries   ------------------------------------------------
  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
  gROOT->LoadMacro("Tools.C");
  LoadPandaStyle();
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------


//input file
  TChain t("cbmsim");
  t->Add("tmpOutput/Box62_1MEvents/Lumi_1M.root");

//load NTupels
  TClonesArray* hit_array=new TClonesArray("PndMvdMCPoint");
  t->SetBranchAddress("MVDPoint",&hit_array);//Branch names

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

  TH2D* hisxy = new TH2D("hisxy","Lumi MC Points, xy view",400,-20.,20.,400,-20.,20.);
  TH2D* hisrz = new TH2D("hisrz","Lumi MC Points, rz view",400,1000.,1250.,400,-20.,20.);
  TH1D* hisde = new TH1D("hisde","Lumi MC Points, Energyloss",500,0.,0.005);
  TH1D* hisak = new TH1D("hisak","Lumi MC Points, Hits / Event",50,-1.,10.);

  TH1D* his34 = new TH1D("his34","Lumi MC Points, Hits / 3-4mrad Event",50,-1.,10.);
  TH1D* his45 = new TH1D("his45","Lumi MC Points, Hits / 4-5mrad Event",50,-1.,10.);
  TH1D* his56 = new TH1D("his56","Lumi MC Points, Hits / 5-6mrad Event",50,-1.,10.);
  TH1D* his67 = new TH1D("his67","Lumi MC Points, Hits / 6-7mrad Event",50,-1.,10.);
  TH1D* his78 = new TH1D("his78","Lumi MC Points, Hits / 7-8mrad Event",50,-1.,10.);

  const int  nEvents = 1000000;
  const bool verbose = false;

//vars for angle-division
  TVector3 vecs,veco;
  Double_t theta[5], relak[5], absev[5];
  for(Int_t j=0; j<5; j++){
    theta[j] = 3.5+j;
    relak[j] = 0;
    absev[j] = 0;
  }

  for (Int_t j=0; j<nEvents && j<t->GetEntries(); j++)
  {
    t->GetEntry(j);
    if(verbose) cout<<"Event No "<<j<<endl;

//read angle of primary particle !just Box!
    double winkel=0.0;
    for (Int_t i=0; i<mc_array->GetEntriesFast(); i++){
      PndMCTrack *track = (PndMCTrack*) mc_array->At(i);
      if(track->GetStartTime() == 0.0){
        TVector3 mom = track->GetMomentum();
        winkel = mom.Theta();
        if(winkel*1000 <3.0 || winkel*1000>=8.0)cout<<"Winkel: "<<winkel<<endl;
        }
      }

    for (Int_t i=0; i<hit_array->GetEntriesFast(); i++)
    {
      if(verbose) cout<<"Point No "<<i<<endl;
      PndMvdMCPoint *hit=(PndMvdMCPoint*)hit_array->At(i);
      int mcpdg = -1;

      vecs.SetXYZ(hit->GetX(), hit->GetY(), hit->GetZ());
      Int_t layer = Int_t(10.*vecs->Mag());
      if(verbose) cout<<vecs.x()<<"\t"<<vecs.y()<<"\t"<<vecs.z()<<endl;

      hisxy->Fill(vecs.x(),vecs.y());
      hisrz->Fill(vecs.z(),((vecs.y()>0.)?1.:-1.)*vecs.Perp());
      hisde->Fill(hit->GetEnergyLoss());
      
    }//end for i (points in event)
  hisak->Fill(hit_array->GetEntriesFast());

//create histos for each angle
  if(winkel>0){
    //winkel /= hit_array->GetEntriesFast(); 
    winkel*=1000;
    int mrad = floor(winkel);
    switch(mrad){
      case 3:{
	absev[0]++;
	if(hit_array->GetEntriesFast() > 3)
          relak[0]++;
        his34->Fill(hit_array->GetEntriesFast());
        break;
        }
      case 4:{
	absev[1]++;
	if(hit_array->GetEntriesFast() > 3)
          relak[1]++;
        his45->Fill(hit_array->GetEntriesFast());
        break;
        }
      case 5:{
	absev[2]++;
	if(hit_array->GetEntriesFast() > 3)
          relak[2]++;
        his56->Fill(hit_array->GetEntriesFast());
        break;
        }
      case 6:{
	absev[3]++;
	if(hit_array->GetEntriesFast() > 3)
          relak[3]++;
        his67->Fill(hit_array->GetEntriesFast());
        break;
        }
      case 7:{
	absev[4]++;
	if(hit_array->GetEntriesFast() > 3)
          relak[4]++;
        his78->Fill(hit_array->GetEntriesFast());
        break;
        }
      }
    }
  }// end for j (events)

//print graphs
  for(Int_t j=0; j<5; j++)
    relak[j]/=absev[j];
  TGraph* graak = new TGraph(5,theta,relak);
  TCanvas* can2 = new TCanvas("can2","rel. Acceptance",0,0,600,400);
  graak->SetFillColor(38);
  graak->SetMinimum(0);
  graak->SetMaximum(1);
  graak->Draw("AB");
  can2->Print("tmpOutput/Box62_1MEvents/acceptance.ps");


  TCanvas* can1 = new TCanvas("can1","MCHit view in Lumi",0,0,1000,1000);
  can1->Divide(3,3);
  can1->cd(1);  
  DrawNice2DHisto(hisxy);
  can1->cd(2);
  DrawNice2DHisto(hisrz);
  can1->cd(3);
  gPad->SetLogy();
  hisde->DrawCopy("");
  can1->cd(4);
  hisak->DrawCopy("");

  can1->cd(5);
  his34->DrawCopy("");
  can1->cd(6);
  his45->DrawCopy("");
  can1->cd(7);
  his56->DrawCopy("");
  can1->cd(8);
  his67->DrawCopy("");
  can1->cd(9);
  his78->DrawCopy("");
  can1->Print("tmpOutput/Box62_1MEvents/outAnaMvdSim.ps");




  // -----   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;
  // ------------------------------------------------------------------------

}