// macro reads from the RadMap Branch
// and fills histograms
void plot(Long64_t nevreq, TString outputfile="RadMap_Out.root")
{
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  rootlogon();

  // open an output file
  TFile *fout = new TFile(outputfile,"RECREATE");

  // chain for the cbmsim tree in the simulation files
  TChain mych("cbmsim");

  /*
  mych.Add("/data/panda/pbarH150_1.root");
  mych.Add("/data/panda/pbarH150_2.root");
  mych.Add("/data/panda/pbarH150_3.root");
  mych.Add("/data/panda/pbarH150_4.root");
  mych.Add("/data/panda/pbarH150_5.root");
  */

  /*
  mych.Add("/data/panda/pbarp150_dpm_1.root");
  mych.Add("/data/panda/pbarp150_dpm_2.root");
  mych.Add("/data/panda/pbarp150_dpm_3.root");
  mych.Add("/data/panda/pbarp150_dpm_4.root");
  */
  
  mych.Add("/s/olaf/test.root");

  //  mych.SetBranchStatus("*",0);
  //  mych.SetBranchStatus("Rad*",1);

  Long64_t nentries = (Long64_t)mych.GetEntries();
  
  /*
  this lines can be used for a single tree file instead of a chain  
  TFile *file = new TFile(treefile);
  TTree *tree = (TTree*)file->Get("cbmsim");

  tree->SetBranchStatus("*",0);
  tree->SetBranchStatus("Rad*",1);

  Long64_t nentries = (Long64_t)tree->GetEntriesFast();
  */

  cout << "we have " << nentries << " events" << endl;
  
  TClonesArray *fRadMapPoint = new TClonesArray("FairRadMapPoint");
  
  mych.SetBranchAddress("RadMap",&fRadMapPoint);

  // if a single tree file is used
  //  tree->SetBranchAddress("RadMap",&fRadMapPoint);

  // define some histograms

  TH2D *histo = new TH2D("histo","histo",400,-100.,300.,150,0.,150.);

  TH2D *his2  = new TH2D("his2","his2",1000,-300.,700.,300,0.,300.);

  TH2D *his3 = new TH2D("his3","his3",400,-100.,300.,150,0.,150.);

  TH2D *his4 = new TH2D("his4","his4",2300,-100.,130.,1000,0.,50.);

  TH2D *spezial = new TH2D("spezial","spezial",100,125.,135.,100,45.,55.);

  TH2D *his5 = new TH2D("his5","his5",120,-60.,60.,120,-60.,60.);

  //  TProfile2D * histop = 
  //    new TProfile2D("histop","histop",400,-100.,300.,150,0.,150.,1e-22,1e-11);
  
  // define some variables
  Double_t x,y,z,r,dose,dSL,mass;
  Int_t pid,VolID;

  // loop over events
  Long64_t nevents;

  cout << "you requested " << nevreq << " events" << endl;

  if(nevreq==0){
    nevents=nentries;
  }else{
    nevents=nevreq;
  }

  if(nevents>nentries) nevents=nentries;
  
  Int_t npoints,sumofpoints;

  sumofpoints=0;

  for (Long64_t i=0; i<nevents; i++) {

    cout << "event " << i << endl;
    
    x=0.;
    y=0.;
    z=0.;
    r=0.;
    dose=0.;
    dSL=0.;
    pid=0;
    npoints=0;
    
    
    mych.GetEntry(i);
    npoints = fRadMapPoint->GetEntries();
    
    //      tree->GetEntry(i);
    //      npoints =  fRadMapPoint->GetEntriesFast();
    
    
    // loop over points
    for (Int_t ii=0; ii < npoints; ii++) {
      
      sumofpoints++;
      
      FairRadMapPoint *p= (FairRadMapPoint *) fRadMapPoint->At(ii);
      
      //      cout << "Point Number " << ii << "  Zout = "<< p->GetZOut() << endl;
      
      dose= p->GetDose();
      dSL = p->GetDoseSL();
      if(dose<=0) continue;

      pid=p->GetPdg();
      VolID=p->GetDetectorID(); // volume ID (replaces detector ID here!!)
      mass=p->GetMass();

      z= p->GetZOut();
      y= p->GetYOut();
      x= p->GetXOut();
      r= sqrt(pow(p->GetXOut(),2)+pow(p->GetYOut(),2));

      //      cout << "Dose: " << z << endl;

      histo->Fill(z,r,dose);
      his2->Fill(z,r,dose);
      his4->Fill(z,r,dose);

      // special conditions
      //      if(z>129.5||z<130.5) his5->Fill(x,y,dose);

      //      if(VolID==5509) his5->Fill(x,y,dose);

      spezial->Fill(z,r,dose);

      if(pid==2112)
	his3->Fill(z,r,dose); // neutrons

    }
    
  }


  // scale by the number of events
  Double_t scale=1./nevents;

  histo->Scale(scale);
  his2->Scale(scale);
  his3->Scale(scale);
  his4->Scale(scale);
  //  his5->Scale(scale);
  spezial->Scale(scale);


  // write to output file
  fout->Write("histo");
  fout->Write("his2");
  fout->Write("his3");
  fout->Write("his4");
  fout->Write("his5");
  fout->Write("spezial");

  cout << "total no. of points: " << sumofpoints << endl;

}