void plottingVars( TH1D* h, TString xname, TString path );
void plottingVars( TH2F* h, TString xname, TString yname, TString path );
void plotProj(TH2F* dp, TH2F* dpaux, TString name, TString input);
Int_t giveMeBin(Double_t x, Double_t min, Double_t max, Int_t nbins);


void DalitzAmpPub(
		  Int_t nbins = 1000,
		  //TString filename = "DalitzAmplitude_dalitzDstarSigGoodDRSTIPdf_KsKK.dat",
		  //TString filename = "DalitzAmplitude_dalitzDstarContGoodDRSTIPdf_KsKK.dat",
		  //TString filename = "DalitzAmplitude_dalitzDstarContGoodDWSTIPdf_KsKK.dat",
		  TString filename = "DalitzAmplitude_dalitzDstarFakeDTIPdf_KsKK.dat",
		  TString pathForPlots = "plots/")
{ 

  Int_t nbins2 = 200;
  
  double mD0_nom  = 1.8645;
  double mKs_nom  = 0.49767;
  double mPi_nom  = 0.139570;
  double mK_nom   = 0.49368;
  double mTrk_nom = (filename.Contains("KsPiPi"))? mPi_nom : mK_nom ;

  double qAB_min = pow(mKs_nom+mTrk_nom,2);
  double qAC_min = pow(mKs_nom+mTrk_nom,2);
  double qBC_min = pow(2*mTrk_nom,2);
  double qAB_max = pow(mD0_nom-mTrk_nom,2);
  double qAC_max = pow(mD0_nom-mTrk_nom,2);
  double qBC_max = pow(mD0_nom-mKs_nom,2);

  double pair1_min = (filename.Contains("KsPiPi"))? qAB_min : qAB_min;
  double pair1_max = (filename.Contains("KsPiPi"))? qAB_max : qAB_max;
  double pair2_min = (filename.Contains("KsPiPi"))? qAC_min : qBC_min;
  double pair2_max = (filename.Contains("KsPiPi"))? qAC_max : qBC_max;
  double pair3_min = (filename.Contains("KsPiPi"))? qBC_min : qAC_min;
  double pair3_max = (filename.Contains("KsPiPi"))? qBC_max : qAC_max;

  ifstream infile;
  infile.open (filename,ios::in);
  if(infile.fail()) {
    cout<<"Error opening the file "<<filename<<endl;
    return;
  }
  cout <<"Reading: "<< filename << endl;
  
  int count(0);
  double pair1,pair2,pair3,argamp,abs2amp,abs2ampNorm;
 
  TH2F* dp_arg = new TH2F("arg","arg",nbins,pair1_min,pair1_max,nbins,pair2_min,pair2_max);
  TH2F* dp_abs2 = new TH2F("abs2","abs2",nbins,pair1_min,pair1_max,nbins,pair2_min,pair2_max);

  TH2F* dp_arg_aux = new TH2F("arg_aux","arg_aux",nbins2,pair3_min,pair3_max,nbins,pair2_min,pair2_max);
  TH2F* dp_abs2_aux = new TH2F("abs2_aux","abs2_aux",nbins2,pair3_min,pair3_max,nbins,pair2_min,pair2_max);
  
  //initialize 2d histos
  for(Int_t j=1 ; j<=nbins ; j++ )
    for(Int_t i=1 ; i<=nbins ; i++ ){
      dp_arg->SetBinContent(i,j,0);
      dp_abs2->SetBinContent(i,j,0);
    }
    for(Int_t i=1 ; i<=nbins2 ; i++ ){
      dp_arg_aux->SetBinContent(i,j,0);
      dp_abs2_aux->SetBinContent(i,j,0);
    }
      
  while(!infile.eof()){ 
    count++; 
    infile >> pair1 >> pair2 >> pair3 >>  argamp >> abs2amp >> abs2ampNorm;
  
    Int_t pair1bin = giveMeBin(pair1, pair1_min, pair1_max, nbins);
    Int_t pair2bin = giveMeBin(pair2, pair2_min, pair2_max, nbins);

    //fill 2d histos
    dp_arg->SetBinContent(pair1bin,pair2bin,argamp);
    dp_abs2->SetBinContent(pair1bin,pair2bin,abs2ampNorm);

    Int_t pair3bin = giveMeBin(pair3, pair3_min, pair3_max, nbins2);

    dp_arg_aux->SetBinContent(pair3bin,pair2bin,argamp);
    dp_abs2_aux->SetBinContent(pair3bin,pair2bin,abs2ampNorm);
   }

  plotProj(dp_arg,dp_arg_aux,"arg",pathForPlots+filename);
  plotProj(dp_abs2,dp_abs2_aux,"abs2",pathForPlots+filename);

}

  
void plotProj(TH2F* dp, TH2F* dpaux, TString name, TString filename)
{
  //plotting two main variables
  TH1D* dp_pair1 =  dp->ProjectionX("_pair1");
  dp_pair1->SetName(name+"_pair1");
  dp_pair1->SetTitle(name+"_pair1");
  plottingVars(dp_pair1,"pair1",filename+"_"+name+"_pair1");

  TH1D* dp_pair2 =  dp->ProjectionY("_pair2");
  dp_pair2->SetName(name+"_pair2");
  dp_pair2->SetTitle(name+"_pair2");
  plottingVars(dp_pair2,"pair2",filename+"_"+name+"_pair2");

  //plotting third component
  TH1D* dp_pair3 =  dpaux->ProjectionX("_pair3");
  dp_pair3->SetName(name+"_pair3");
  dp_pair3->SetTitle(name+"_pair3");
  plottingVars(dp_pair3,"pair3",filename+"_"+name+"_pair3");


  //2D plot
  plottingVars(dp,"pair1","pair2",filename+"_"+name+"_dp");

}

void plottingVars(TH1D* h,
                  TString xname,
                  TString path
                  ){

  gStyle->SetOptStat(0000);
  gStyle->SetOptFit(0000);

  TCanvas*  c = new TCanvas();
  c->cd();

  h->Draw("C");
  h->GetXaxis()->SetLabelSize(0.045);
  h->GetYaxis()->SetLabelSize(0.045);
  h->GetXaxis()->SetTitleSize(0.045);
  h->GetYaxis()->SetTitleSize(0.045);
  h->GetXaxis()->SetTitle(xname);

  c->cd()->Print(path + ".eps");
  delete c;
}

void plottingVars(TH2F* h,
                  TString xname,
                  TString yname,              
		  TString path
                  ){

  gStyle->SetOptStat(0000);
  gStyle->SetOptFit(0000);
  gStyle->SetPalette(1);

  TCanvas*  c = new TCanvas();
  c->cd();

  h->Draw("COLZ");
  h->GetXaxis()->SetLabelSize(0.045);
  h->GetYaxis()->SetLabelSize(0.045);
  h->GetXaxis()->SetTitleSize(0.045);
  h->GetYaxis()->SetTitleSize(0.045);
  h->GetXaxis()->SetTitle(xname);
  h->GetYaxis()->SetTitle(yname);

  c->cd()->Print( path + ".eps");
  delete c;
}

Int_t giveMeBin(Double_t x, Double_t min, Double_t max, Int_t nbins){
  Int_t out;
  Double_t tmp= (Double_t)(x-min)/(Double_t)(max-min);
  tmp *= (Double_t)nbins;
  out = (Int_t) tmp + 1;
  return out;
}