void plot(TString treefile)
{
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  rootlogon();

  Bool_t chain=kTRUE;

  if(chain) {

  TChain mych("cbmsim");

  /*
  mych.Add("/data/panda/pbarCu62_1.root");
  mych.Add("/data/panda/pbarCu62_2.root");
  mych.Add("/data/panda/pbarCu62_3.root");
  mych.Add("/data/panda/pbarCu62_4.root");
  mych.Add("/data/panda/pbarCu62_5.root");
  mych.Add("/data/panda/pbarCu62_6.root");
  mych.Add("/data/panda/pbarCu62_7.root");
  mych.Add("/data/panda/pbarCu62_8.root");
  mych.Add("/data/panda/pbarCu62_9.root");
  mych.Add("/data/panda/pbarCu62_10.root");
  mych.Add("/data/panda/pbarCu62_11.root");
  mych.Add("/data/panda/pbarCu62_12.root");
  mych.Add("/data/panda/pbarCu62_13.root");
  mych.Add("/data/panda/pbarCu62_14.root");
  mych.Add("/data/panda/pbarCu62_15.root");
  mych.Add("/data/panda/pbarCu62_16.root");
  mych.Add("/data/panda/pbarCu62_17.root");
  mych.Add("/data/panda/pbarCu62_18.root");
  mych.Add("/data/panda/pbarCu62_19.root");
  */

  /*
  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("/data/panda/pbarp_150_doseSL_TG3_2.root");
  // mych.Add("/data/panda/pbarp_150_doseSL_TG3_1.root");

  mych.Add("/data/panda/pbarp_150_doseSL_TG4_1.root");

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

  Long64_t nentries = (Long64_t)mych.GetEntries();
  
  } else {

  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 *fRadLenPoint = new TClonesArray("FairRadLenPoint");
  
  if(chain) {
    mych.SetBranchAddress("RadLen",&fRadLenPoint);
  } else {
    tree->SetBranchAddress("RadLen",&fRadLenPoint);
  }
  //  if (nevents>nentries) nevents=nentries;

  //  TGraph2D *graph = new TGraph2D();
  Float_t x,y,z,zSL;

  TH2D *histo = new TH2D("histo","histo",100,-100.,100.,100,0.,100.);
  TProfile2D * histop = 
    new TProfile2D("histop","histop",400,-100.,300.,150,0.,150.,1e-22,1e-11);

  TProfile2D * histop2 = 
    new TProfile2D("histop2","histop2",400,-100.,300.,150,0.,150.,1e-11,1e-5);

  //  TProfile2D * histop3 =
  //    new TProfile2D("histop3","histop3",400,-100.,300.,150,0.,150.,1e-4,10);

  TProfile2D * histoq = 
    new TProfile2D("histoq","histoq",400,-100.,300.,150,0.,150.,1e-22,1e-11);

  TProfile2D * histoq2 = 
    new TProfile2D("histoq2","histoq2",400,-100.,300.,150,0.,150.,1e-11,1e-5);

  TProfile2D * histoq3 =
    new TProfile2D("histoq3","histoq3",400,-100.,300.,150,0.,150.,1e-4,10);


  //  TProfile2D * histop2 = 
  //    new TProfile2D("histop2","histop2",260,-100.,30.,90,0.,45.,1e-13,1e-09);

  histop->Draw("COLZ");

  // loop over events
  Long64_t nevents;

  nevents=nentries;
  //  nevents=100;

  Int_t npoints,sumofpoints;

  sumofpoints=0;

  for (Long64_t i=0; i<nevents; i++) {
    
    x=0.;
    y=0.;
    z=0.;
    zSL=0.;
    npoints=0;

    if(chain) {
      mych.GetEntry(i);
      npoints = fRadLenPoint->GetEntries();
    } else {
      //      tree->GetEntry(i);
      //      npoints =  fRadLenPoint->GetEntriesFast();
    }

    // loop over points
    for (Int_t ii=0; ii < npoints; ii++) {
     
      sumofpoints++;

      FairRadLenPoint *p= (FairRadLenPoint *) fRadLenPoint->At(ii);
      
      //      cout << "Point Number " << ii << "  Zout = "<< p->GetZOut() << endl;

      z= p->GetDose();
      zSL = p->GetDoseSL();
      if(z<=0) continue;

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


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

      //      iii++;
      //      graph->SetPoint(iii,x,y,z);
      histo->Fill(x,y,z);
      histop->Fill(x,y,z);
      //      histop3->Fill(x,y,z);
      histop2->Fill(x,y,z);
      
      histoq->Fill(x,y,zSL);
      histoq3->Fill(x,y,zSL);
      histoq2->Fill(x,y,zSL);



    }
    
  }

  //  graph->Draw("surf1");
  Double_t scale=1./nevents;
  histop->Scale(scale);
  histop2->Scale(scale);
  //  histop3->Scale(scale);
  histoq->Scale(scale);
  histoq2->Scale(scale);
  histoq3->Scale(scale);

  histoq->Draw("COLZ");

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

}