#include "Riostream.h"
#include "TFile.h"
#include "TH1.h"
#include "TNtuple.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TAxis.h"
#include <fstream>

void basic (const Char_t*fdata="DreiD_Data12.3d", const Char_t*froot="b.root") 
{
  
    gStyle->SetOptStat(0);
   //gStyle->SetPalette(60);
   
   TCanvas *c1 = new TCanvas("c1","aa",0,0,600,600);
   TCanvas *c2 = new TCanvas("c2","aa",0,0,600,600);
    c1->SetRightMargin(0.15);
    c1->SetBottomMargin(0.2);
    c1->SetLeftMargin(0.2);
    ifstream ReadData;
   
   ReadData.open(fdata);

   if (!ReadData.is_open()){
	   cout<<"could not open the input file..."<<endl;
	   return -1;
   }
  
   double x,y,z,count;
   int nlines = 0;
   TFile *f = new TFile(froot,"RECREATE");
   TH3D *h3 = new TH3D("h3","h3 title",50,-100,100,50,-100,100,50,-100,100);
   
   
   TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","x:y:z:cont");
   
   
     while (ReadData.good()) 
      {
      ReadData >> x >> y >>z >> count;
      if (ReadData.good())
        {
        	 
          h3->Fill(x,y,z,count);
          if(count>0) // the sky blue part removal
            {
             ntuple->Fill(x,y,z,count);
             }
	   nlines++;
         }
   
       }

   cout<<"number of lines: "<<nlines<<endl;

   //c1->cd();
   h3->GetXaxis()->SetTitle("X-Momentum (au)");
   h3->GetYaxis()->SetTitle("Y-Momentum (au)");
   h3->GetZaxis()->SetTitle("Z-Momentum (au)");
   h3->GetXaxis()->CenterTitle(true);
   h3->GetZaxis()->CenterTitle(true);
   h3->GetXaxis()->CenterTitle(true);
   h3->GetXaxis()->SetTitleFont(62);
   h3->GetYaxis()->SetTitleFont(62);
   h3->GetZaxis()->SetTitleFont(62);
   h3->GetXaxis()->SetTitleSize(0.04);
   h3->GetYaxis()->SetTitleSize(0.04);
   h3->GetZaxis()->SetTitleSize(0.04);
   h3->GetYaxis()->SetTitleOffset(1.5);
   h3->GetXaxis()->SetTitleOffset(1.5);
   h3->GetXaxis()->SetLabelFont(62);
   h3->GetXaxis()->SetLabelOffset(0.005);
   h3->GetXaxis()->SetLabelSize(0.04);
   h3->GetYaxis()->SetLabelFont(62);
   h3->GetYaxis()->SetLabelOffset(0.005);
   h3->GetYaxis()->SetLabelSize(0.04);
   h3->GetXaxis()->SetNdivisions(804,kTRUE);
   h3->GetYaxis()->SetNdivisions(804,kTRUE);
   
   
  

   c2->cd();
   ntuple->Draw("x:y:z:cont>>plot4d","","COLZ");
   plot4d->GetXaxis()->SetTitle("X-Momentum (au)");
   plot4d->GetYaxis()->SetTitle("Y-Momentum (au)");
   plot4d->GetZaxis()->SetTitle("Z-Momentum (au)");
   plot4d->GetXaxis()->SetTitleOffset(1.5);
   plot4d->GetYaxis()->SetTitleOffset(1.8);
   plot4d->GetZaxis()->SetTitleOffset(1.5);
   plot4d->GetXaxis()->CenterTitle(true);
   plot4d->GetYaxis()->CenterTitle(true);
   plot4d->GetZaxis()->CenterTitle(true);
   
   

  // c1->SaveAs("momenta.png");
   c2->SaveAs("momenta2.png");
   ReadData.close();
   f->Write();
   delete f;

}