using namespace std;
#include<iostream>
#include<fstream>
#include<cstdlib>
#include<string.h>
#include<cstdio>
#include<assert.h>

void SDEMs_ext(){

  TFile *out = new TFile("sdmes_pars_1.3.root","RECREATE");            //define output file 

  TFile *f = TFile::Open("sdmes_1.3.root");                     //open histograms file              


  /**************************** define some histograms *****************************************/

  TH1F *hazimuth_hel_par = (TH1F*)f->Get("azimuth_hel_par");
  hazimuth_hel_par->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazimuth_hel_par->GetYaxis()->SetTitle("Entries");
  TH1F *hazimuth_hel_per = (TH1F*)f->Get("azimuth_hel_per");
  hazimuth_hel_per->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazimuth_hel_per->GetYaxis()->SetTitle("Entries");
  TH1F *hcosthet_hel_par = (TH1F*)f->Get("costhet_hel_par");
  hcosthet_hel_par->GetXaxis()->SetTitle("cos#theta_{#omega}^{hel}");
  hcosthet_hel_par->GetYaxis()->SetTitle("Entries");
  TH1F *hcosthet_hel_per = (TH1F*)f->Get("costhet_hel_per");
  hcosthet_hel_per->GetXaxis()->SetTitle("cos#theta_{#omega}^{hel}");
  hcosthet_hel_per->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_par = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para");
  hpsipluspip_par->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp");
  hpsipluspip_per->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para");
  hpsipluspim_par->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp");
  hpsipluspim_per->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para");
  hpsiminuspip_par->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp");
  hpsiminuspip_per->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para");
  hpsiminuspim_par->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp");
  hpsiminuspim_per->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per->GetYaxis()->SetTitle("Entries");

  TH1F *hazim_hel_par1 = (TH1F*)f->Get("azim_hel_par1");
  hazim_hel_par1->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_par1->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_per1 = (TH1F*)f->Get("azim_hel_per1");
  hazim_hel_per1->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_per1->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_par2 = (TH1F*)f->Get("azim_hel_par2");
  hazim_hel_par2->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_par2->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_per2 = (TH1F*)f->Get("azim_hel_per2");
  hazim_hel_per2->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_per2->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_par3 = (TH1F*)f->Get("azim_hel_par3");
  hazim_hel_par3->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_par3->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_per3 = (TH1F*)f->Get("azim_hel_per3");
  hazim_hel_per3->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_per3->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_par4 = (TH1F*)f->Get("azim_hel_par4");
  hazim_hel_par4->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_par4->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_per4 = (TH1F*)f->Get("azim_hel_per4");
  hazim_hel_per4->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_per4->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_par5 = (TH1F*)f->Get("azim_hel_par5");
  hazim_hel_par5->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_par5->GetYaxis()->SetTitle("Entries");
  TH1F *hazim_hel_per5 = (TH1F*)f->Get("azim_hel_per5");
  hazim_hel_per5->GetXaxis()->SetTitle("#phi_{#omega}^{hel}");
  hazim_hel_per5->GetYaxis()->SetTitle("Entries");

  TH1F *hpsipa = (TH1F*)f->Get("psi_par");
  hpsipa->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe = (TH1F*)f->Get("psi_per");
  hpsipe->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipa_1 = (TH1F*)f->Get("psi_par1");
  hpsipa_1->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa_1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipa_2 = (TH1F*)f->Get("psi_par2");
  hpsipa_2->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa_2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipa_3 = (TH1F*)f->Get("psi_par3");
  hpsipa_3->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa_3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipa_4 = (TH1F*)f->Get("psi_par4");
  hpsipa_4->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa_4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipa_5 = (TH1F*)f->Get("psi_par5");
  hpsipa_5->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipa_5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe_1 = (TH1F*)f->Get("psi_per1");
  hpsipe_1->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe_1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe_2 = (TH1F*)f->Get("psi_per2");
  hpsipe_2->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe_2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe_3 = (TH1F*)f->Get("psi_per3");
  hpsipe_3->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe_3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe_4 = (TH1F*)f->Get("psi_per4");
  hpsipe_4->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe_4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipe_5 = (TH1F*)f->Get("psi_per5");
  hpsipe_5->GetXaxis()->SetTitle("#psi_{#omega}^{hel}");
  hpsipe_5->GetYaxis()->SetTitle("Entries");

  TH1F *hpsipluspip_par1 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para1");
  hpsipluspip_par1->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per1 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp1");
  hpsipluspip_per1->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par1 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para1");
  hpsipluspim_par1->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per1 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp1");
  hpsipluspim_per1->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par1 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para1");
  hpsiminuspip_par1->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per1 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp1");
  hpsiminuspip_per1->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par1 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para1");
  hpsiminuspim_par1->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per1 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp1");
  hpsiminuspim_per1->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per1->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_par2 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para2");
  hpsipluspip_par2->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per2 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp2");
  hpsipluspip_per2->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par2 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para2");
  hpsipluspim_par2->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per2 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp2");
  hpsipluspim_per2->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par2 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para2");
  hpsiminuspip_par2->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per2 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp2");
  hpsiminuspip_per2->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par2 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para2");
  hpsiminuspim_par2->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per2 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp2");
  hpsiminuspim_per2->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per2->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_par3 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para3");
  hpsipluspip_par3->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per3 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp3");
  hpsipluspip_per3->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par3 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para3");
  hpsipluspim_par3->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per3 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp3");
  hpsipluspim_per3->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par3 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para3");
  hpsiminuspip_par3->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per3 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp3");
  hpsiminuspip_per3->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par3 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para3");
  hpsiminuspim_par3->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per3 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp3");
  hpsiminuspim_per3->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per3->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_par4 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para4");
  hpsipluspip_par4->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per4 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp4");
  hpsipluspip_per4->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par4 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para4");
  hpsipluspim_par4->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per4 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp4");
  hpsipluspim_per4->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par4 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para4");
  hpsiminuspip_par4->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per4 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp4");
  hpsiminuspip_per4->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par4 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para4");
  hpsiminuspim_par4->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per4 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp4");
  hpsiminuspim_per4->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per4->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_par5 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_para5");
  hpsipluspip_par5->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_par5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspip_per5 = (TH1F*)f->Get("dsdwplus_psi_phi_pip_perp5");
  hpsipluspip_per5->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspip_per5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_par5 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_para5");
  hpsipluspim_par5->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_par5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsipluspim_per5 = (TH1F*)f->Get("dsdwplus_psi_phi_pim_perp5");
  hpsipluspim_per5->GetXaxis()->SetTitle("(#phi+#Phi)_{#omega}^{hel}");
  hpsipluspim_per5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_par5 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_para5");
  hpsiminuspip_par5->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_par5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspip_per5 = (TH1F*)f->Get("dsdwminus_psi_phi_pip_perp5");
  hpsiminuspip_per5->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspip_per5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_par5 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_para5");
  hpsiminuspim_par5->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_par5->GetYaxis()->SetTitle("Entries");
  TH1F *hpsiminuspim_per5 = (TH1F*)f->Get("dsdwminus_psi_phi_pim_perp5");
  hpsiminuspim_per5->GetXaxis()->SetTitle("(#phi-#Phi)_{#omega}^{hel}");
  hpsiminuspim_per5->GetYaxis()->SetTitle("Entries");

  /********************* read signal to noise ratios to scale phi distributions *****************************************/  
  /***************************/
  Float_t s2npar[8];
  Float_t s2nper[8];
  int nlines;

  ifstream in; 
  //in.open("signal2noise.txt");
  in.open("signal2noise_voigt_sdmes.txt");
   
  while(in.good() && nlines < 8){
    while (!in.eof()){
      in >> s2npar[nlines] >> s2nper[nlines];
      
      nlines++;
    }
  }
  in.close();
  
  //define scale factors
  
  float c1_pare1 = s2npar[0];
  float c2_pare1 = s2npar[1];
  float c3_pare1 = s2npar[2];
  float c4_pare1 = s2npar[3];
  float c5_pare1 = s2npar[4];
  float c6_pare1 = s2npar[5];
  
  float c1_pere1 = s2nper[0];
  float c2_pere1 = s2nper[1];
  float c3_pere1 = s2nper[2];
  float c4_pere1 = s2nper[3];
  float c5_pere1 = s2nper[4];
  float c6_pere1 = s2nper[5];
  
  /********************* call Sumw2 to calculate errors **************/
  hcosthet_hel_par->Sumw2();
  hcosthet_hel_per->Sumw2();
  hazimuth_hel_par->Sumw2();
  hazimuth_hel_per->Sumw2();
  hpsipluspip_par->Sumw2();
  hpsipluspip_per->Sumw2();
  hpsipluspim_par->Sumw2();
  hpsipluspim_per->Sumw2();
  hpsiminuspip_par->Sumw2();
  hpsiminuspip_per->Sumw2();
  hpsiminuspim_par->Sumw2();
  hpsiminuspim_per->Sumw2();

  hpsipa->Sumw2();
  hpsipe->Sumw2();
  hpsipa_1->Sumw2();
  hpsipa_2->Sumw2();
  hpsipa_3->Sumw2();
  hpsipa_4->Sumw2();
  hpsipa_5->Sumw2();
  hpsipe_1->Sumw2();
  hpsipe_2->Sumw2();
  hpsipe_3->Sumw2();
  hpsipe_4->Sumw2();
  hpsipe_5->Sumw2();

  hazim_hel_par1->Sumw2();
  hazim_hel_par2->Sumw2();
  hazim_hel_par3->Sumw2();
  hazim_hel_par4->Sumw2();
  hazim_hel_par5->Sumw2();
  hazim_hel_per1->Sumw2();
  hazim_hel_per2->Sumw2();
  hazim_hel_per3->Sumw2();
  hazim_hel_per4->Sumw2();
  hazim_hel_per5->Sumw2();

  hpsipluspip_par1->Sumw2();
  hpsipluspip_per1->Sumw2();
  hpsipluspim_par1->Sumw2();
  hpsipluspim_per1->Sumw2();
  hpsiminuspip_par1->Sumw2();
  hpsiminuspip_per1->Sumw2();
  hpsiminuspim_par1->Sumw2();
  hpsiminuspim_per1->Sumw2();
  hpsipluspip_par2->Sumw2();
  hpsipluspip_per2->Sumw2();
  hpsipluspim_par2->Sumw2();
  hpsipluspim_per2->Sumw2();
  hpsiminuspip_par2->Sumw2();
  hpsiminuspip_per2->Sumw2();
  hpsiminuspim_par2->Sumw2();
  hpsiminuspim_per2->Sumw2();
  hpsipluspip_par3->Sumw2();
  hpsipluspip_per3->Sumw2();
  hpsipluspim_par3->Sumw2();
  hpsipluspim_per3->Sumw2();
  hpsiminuspip_par3->Sumw2();
  hpsiminuspip_per3->Sumw2();
  hpsiminuspim_par3->Sumw2();
  hpsiminuspim_per3->Sumw2();
  hpsipluspip_par4->Sumw2();
  hpsipluspip_per4->Sumw2();
  hpsipluspim_par4->Sumw2();
  hpsipluspim_per4->Sumw2();
  hpsiminuspip_par4->Sumw2();
  hpsiminuspip_per4->Sumw2();
  hpsiminuspim_par4->Sumw2();
  hpsiminuspim_per4->Sumw2();
  hpsipluspip_par5->Sumw2();
  hpsipluspip_per5->Sumw2();
  hpsipluspim_par5->Sumw2();
  hpsipluspim_per5->Sumw2();
  hpsiminuspip_par5->Sumw2();
  hpsiminuspip_per5->Sumw2();
  hpsiminuspim_par5->Sumw2();
  hpsiminuspim_per5->Sumw2();

  /************** Scale phis distributions ********************************/

   //hcosthet_CM_par_e1->Scale(c1_pare1);
   //hcosthet_CM_per_e1->Scale(c1_pere1);
   hazim_hel_par1->Scale(c1_pare1);
   hazim_hel_par2->Scale(c2_pare1);
   hazim_hel_par3->Scale(c3_pare1);
   hazim_hel_par4->Scale(c4_pare1);
   hazim_hel_par5->Scale(c5_pare1);
   hazim_hel_per1->Scale(c1_pere1);
   hazim_hel_per2->Scale(c2_pere1);
   hazim_hel_per3->Scale(c3_pere1);
   hazim_hel_per4->Scale(c4_pere1);
   hazim_hel_per5->Scale(c5_pere1);
   hpsipa_1->Scale(c1_pare1);
   hpsipa_2->Scale(c2_pare1);
   hpsipa_3->Scale(c3_pare1);
   hpsipa_4->Scale(c4_pare1);
   hpsipa_5->Scale(c5_pare1);
   hpsipe_1->Scale(c1_pere1);
   hpsipe_2->Scale(c2_pere1);
   hpsipe_3->Scale(c3_pere1);
   hpsipe_4->Scale(c4_pere1);
   hpsipe_5->Scale(c5_pere1);
  hpsipluspip_par1->Scale(c1_pare1);
  hpsipluspip_per1->Scale(c1_pere1);
  hpsipluspim_par1->Scale(c1_pare1);
  hpsipluspim_per1->Scale(c1_pere1);
  hpsiminuspip_par1->Scale(c1_pare1);
  hpsiminuspip_per1->Scale(c1_pere1);
  hpsiminuspim_par1->Scale(c1_pare1);
  hpsiminuspim_per1->Scale(c1_pere1);
  hpsipluspip_par2->Scale(c2_pare1);
  hpsipluspip_per2->Scale(c2_pere1);
  hpsipluspim_par2->Scale(c2_pare1);
  hpsipluspim_per2->Scale(c2_pere1);
  hpsiminuspip_par2->Scale(c2_pare1);
  hpsiminuspip_per2->Scale(c2_pere1);
  hpsiminuspim_par2->Scale(c2_pare1);
  hpsiminuspim_per2->Scale(c2_pere1);
  hpsipluspip_par3->Scale(c3_pare1);
  hpsipluspip_per3->Scale(c3_pere1);
  hpsipluspim_par3->Scale(c3_pare1);
  hpsipluspim_per3->Scale(c3_pere1);
  hpsiminuspip_par3->Scale(c3_pare1);
  hpsiminuspip_per3->Scale(c3_pere1);
  hpsiminuspim_par3->Scale(c3_pare1);
  hpsiminuspim_per3->Scale(c3_pere1);
  hpsipluspip_par4->Scale(c4_pare1);
  hpsipluspip_per4->Scale(c4_pere1);
  hpsipluspim_par4->Scale(c4_pare1);
  hpsipluspim_per4->Scale(c4_pere1);
  hpsiminuspip_par4->Scale(c4_pare1);
  hpsiminuspip_per4->Scale(c4_pere1);
  hpsiminuspim_par4->Scale(c4_pare1);
  hpsiminuspim_per4->Scale(c4_pere1);
  hpsipluspip_par5->Scale(c5_pare1);
  hpsipluspip_per5->Scale(c5_pere1);
  hpsipluspim_par5->Scale(c5_pare1);
  hpsipluspim_per5->Scale(c5_pere1);
  hpsiminuspip_par5->Scale(c5_pare1);
  hpsiminuspip_per5->Scale(c5_pere1);
  hpsiminuspim_par5->Scale(c5_pare1);
  hpsiminuspim_per5->Scale(c5_pere1);

  /***** fit the angular distributions *****/

  TF1* fitcos = new TF1("fitcos","[0]*((0.5*(1-[1])*(1-cos(TMath::DegToRad()*(x))*cos(TMath::DegToRad()*(x)))) + [1]*cos(TMath::DegToRad()*(x))*cos(TMath::DegToRad()*(x)))");
  //TF1* fitcos = new TF1("fitcos","[0]*([1]*(cos(TMath::DegToRad()*(x))*cos(TMath::DegToRad()*(x))))");
  TF1* fitphi = new TF1("fitphi","[0]*(1-[1]*cos(TMath::DegToRad()*2*(x)))");
  TF1* fitpsi = new TF1("fitpsi","[0]*(1+2*[1]*0.8*cos(TMath::DegToRad()*2*(x)))");
  //TF1* fitro5 = new TF1("fitro5","[0]*(1-[1]*0.8*cos(TMath::DegToRad()*2*(x)))");//ojo con esta
  //TF1* fitro5 = new TF1("fitasym","([0]*cos(TMath::DegToRad()*2*(x)+[1]))/(2.0+[0]*cos(TMath::DegToRad()*2*(x)+[1]))");
  TF1* fitro5 = new TF1("fitro5","(([0]+[1])*[2]*cos(TMath::DegToRad()*2*(x)))/(2.0+([0]-[1])*[2]*cos(TMath::DegToRad()*2*(x)))");
  //TF1* fitro5 = new TF1("fitasym","([0]*(1.6)*cos(TMath::DegToRad()*2*(x)))/(2.0+[0]*(0.02)*cos(TMath::DegToRad()*2*(x)))");

  //fitcos->SetParameter(0,1.5);
  //fitphi->SetParameter(0,0.15);

  //define files for pars output
  ofstream output("params_sdmes.txt");
  ofstream outcos("cos_sdmes_tbins.txt");
  ofstream outphi("phi_sdmes_tbins.txt");
  ofstream outro3pip("ro3pip_sdmes_tbins.txt");
  ofstream outro4pip("ro4pip_sdmes_tbins.txt");
  ofstream outro3pim("ro3pim_sdmes_tbins.txt");
  ofstream outro4pim("ro4pim_sdmes_tbins.txt");
  ofstream outro5("ro5_sdmes_tbins.txt");

  /***************************** Fit integrated distributions **********************************/
  hcosthet_hel_par->Fit("fitcos","EM");
  float Npar = fitcos->GetParameter(0); 
  float ro1par = fitcos->GetParameter(1); 
  fitcos->SetParNames("N","#rho^{1}");
  hcosthet_hel_per->Fit("fitcos","EM");
  float Nper = fitcos->GetParameter(0); 
  float ro1per = fitcos->GetParameter(1); 
  fitcos->SetParNames("N","#rho^{1}");
  hazimuth_hel_par->Fit("fitphi","EM");
  float N_par = fitphi->GetParameter(0); 
  float ro2par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazimuth_hel_per->Fit("fitphi","EM");
  float N_per = fitphi->GetParameter(0); 
  float ro2per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  /*********************** integrated cos and phi dstributions *********************************/
  outcos << Npar <<"  "<< ro1par <<endl;
  outcos << Nper <<"  "<< ro1per <<endl;

  output << N_par <<"  "<< ro2par <<endl;
  output << N_per <<"  "<< ro2per <<endl;
  /*********************************************************************************************/

  /*********************** Fit asymmetry distributions for rho 5 *******************************/ 
  TH1F *hassymetry = (TH1F*)hpsipe->GetAsymmetry(hpsipa);
  hassymetry->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  TH1F *hassymetry_1 = (TH1F*)hpsipe_1->GetAsymmetry(hpsipa_1);
  hassymetry_1->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry_1->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  TH1F *hassymetry_2 = (TH1F*)hpsipe_2->GetAsymmetry(hpsipa_2);
  hassymetry_2->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry_2->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  TH1F *hassymetry_3 = (TH1F*)hpsipe_3->GetAsymmetry(hpsipa_3);
  hassymetry_3->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry_3->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  TH1F *hassymetry_4 = (TH1F*)hpsipe_4->GetAsymmetry(hpsipa_4);
  hassymetry_4->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry_4->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  TH1F *hassymetry_5 = (TH1F*)hpsipe_5->GetAsymmetry(hpsipa_5);
  hassymetry_5->GetXaxis()->SetTitle("#Psi_{#omega}^{Hel}");
  hassymetry_5->GetYaxis()->SetTitle("#frac{W_{#parallel}-W_{#perp}}{W_{#parallel}+W_{#perp}}");
  //fitasym->SetParameter(0, 0.8);
  //fitasym->SetParameter(1, 0.7);
  hassymetry->Fit("fitro5","EM");
  float a0 = fitro5->GetParameter(0); 
  float a0_er = fitro5->GetParError(0); 
  float a1 = fitro5->GetParameter(1); 
  float a1_er = fitro5->GetParError(1); 
  float a2 = fitro5->GetParameter(2); 
  float a2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  hassymetry_1->Fit("fitro5","EM");
  float b0 = fitro5->GetParameter(0); 
  float b0_er = fitro5->GetParError(0); 
  float b1 = fitro5->GetParameter(1); 
  float b1_er = fitro5->GetParError(1); 
  float b2 = fitro5->GetParameter(2); 
  float b2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  hassymetry_2->Fit("fitro5","EM");
  float c0 = fitro5->GetParameter(0); 
  float c0_er = fitro5->GetParError(0); 
  float c1 = fitro5->GetParameter(1); 
  float c1_er = fitro5->GetParError(1); 
  float c2 = fitro5->GetParameter(2); 
  float c2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  hassymetry_3->Fit("fitro5","EM");
  float d0 = fitro5->GetParameter(0); 
  float d0_er = fitro5->GetParError(0); 
  float d1 = fitro5->GetParameter(1); 
  float d1_er = fitro5->GetParError(1); 
  float d2 = fitro5->GetParameter(2); 
  float d2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  hassymetry_4->Fit("fitro5","EM");
  float e0 = fitro5->GetParameter(0); 
  float e0_er = fitro5->GetParError(0); 
  float e1 = fitro5->GetParameter(1); 
  float e1_er = fitro5->GetParError(1); 
  float e2 = fitro5->GetParameter(2); 
  float e2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  hassymetry_5->Fit("fitro5","EM");
  float f0 = fitro5->GetParameter(0); 
  float f0_er = fitro5->GetParError(0); 
  float f1 = fitro5->GetParameter(1); 
  float f1_er = fitro5->GetParError(1); 
  float f2 = fitro5->GetParameter(2); 
  float f2_er = fitro5->GetParError(2); 
  fitro5->SetParNames("#rho^{5}");
  /************************************* out rho 5 pars ************************************/
  outro5 << a0 <<"  "<< a0_er <<"  "<< a1 <<"  "<< a1_er <<"  "<< a2 <<"  "<< a2_er << endl;
  outro5 << b0 <<"  "<< b0_er <<"  "<< b1 <<"  "<< b1_er <<"  "<< b2 <<"  "<< b2_er << endl;
  outro5 << c0 <<"  "<< c0_er <<"  "<< c1 <<"  "<< c1_er <<"  "<< c2 <<"  "<< c2_er << endl;
  outro5 << d0 <<"  "<< d0_er <<"  "<< d1 <<"  "<< d1_er <<"  "<< d2 <<"  "<< d2_er << endl;
  outro5 << e0 <<"  "<< e0_er <<"  "<< e1 <<"  "<< e1_er <<"  "<< e2 <<"  "<< e2_er << endl;
  outro5 << f0 <<"  "<< f0_er <<"  "<< f1 <<"  "<< f1_er <<"  "<< f2 <<"  "<< f2_er << endl;
  /*****************************************************************************************/
  
  /*********************** Fit phi distributions for t bins ********************************/
  hazim_hel_par1->Fit("fitphi","EM");
  float Nphi1 = fitphi->GetParameter(0); 
  float ro2F1par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_par2->Fit("fitphi","EM");
  float Nphi2 = fitphi->GetParameter(0); 
  float ro2F2par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_par3->Fit("fitphi","EM");
  float Nphi3 = fitphi->GetParameter(0); 
  float ro2F3par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_par4->Fit("fitphi","EM");
  float Nphi4 = fitphi->GetParameter(0); 
  float ro2F4par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_par5->Fit("fitphi","EM");
  float Nphi5 = fitphi->GetParameter(0); 
  float ro2F5par = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");

  hazim_hel_per1->Fit("fitphi","EM");
  float Nphi11 = fitphi->GetParameter(0); 
  float ro2F1per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_per2->Fit("fitphi","EM");
  float Nphi22 = fitphi->GetParameter(0); 
  float ro2F2per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_per3->Fit("fitphi","EM");
  float Nphi33 = fitphi->GetParameter(0); 
  float ro2F3per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_per4->Fit("fitphi","EM");
  float Nphi44 = fitphi->GetParameter(0); 
  float ro2F4per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  hazim_hel_per5->Fit("fitphi","EM");
  float Nphi55 = fitphi->GetParameter(0); 
  float ro2F5per = fitphi->GetParameter(1); 
  fitphi->SetParNames("N","#rho^{2}");
  /************************************ out rho 2 pars **************************************************************************************************************/
  outphi << Nphi1 <<"  "<< ro2F1par <<"  "<< Nphi2 <<"  "<< ro2F2par <<"  "<< Nphi3 <<"  "<< ro2F3par <<"  "<< Nphi4 <<"  "<< ro2F4par <<"  "<< Nphi5 <<"  "<< ro2F5par <<endl;
  outphi << Nphi11 <<"  "<< ro2F1per <<"  "<< Nphi22 <<"  "<< ro2F2per <<"  "<< Nphi33 <<"  "<< ro2F3per <<"  "<< Nphi44 <<"  "<< ro2F4per <<"  "<< Nphi55 <<"  "<< ro2F5per <<endl;
  /******************************************************************************************************************************************************************/

  /******* Fit psi distribuitons for phi +/- Phi and for pi+ and pi- *************************/
  hpsipluspip_par->Fit("fitpsi","EM");
  float Npsi1_par = fitpsi->GetParameter(0); 
  float ro4psi1par = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per->Fit("fitpsi","EM");
  float Npsi1_per = fitpsi->GetParameter(0); 
  float ro4psi1per = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par->Fit("fitpsi","EM");
  float Npsi2_par = fitpsi->GetParameter(0); 
  float ro4psi2par = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per->Fit("fitpsi","EM");
  float Npsi2_per = fitpsi->GetParameter(0); 
  float ro4psi2per = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par->Fit("fitpsi","EM");
  float Npsi3_par = fitpsi->GetParameter(0); 
  float ro3psi3par = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per->Fit("fitpsi","EM");
  float Npsi3_per = fitpsi->GetParameter(0); 
  float ro3psi3per = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par->Fit("fitpsi","EM");
  float Npsi4_par = fitpsi->GetParameter(0); 
  float ro3psi4par = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per->Fit("fitpsi","EM");
  float Npsi4_per = fitpsi->GetParameter(0); 
  float ro3psi4per = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");

  // and for t bins
  hpsipluspip_par1->Fit("fitpsi","EM");
  float N1_par1 = fitpsi->GetParameter(0); 
  float r1par1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per1->Fit("fitpsi","EM");
  float N1_per1 = fitpsi->GetParameter(0); 
  float r1per1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par1->Fit("fitpsi","EM");
  float N2_par1 = fitpsi->GetParameter(0); 
  float r2par1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per1->Fit("fitpsi","EM");
  float N2_per1 = fitpsi->GetParameter(0); 
  float r2per1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par1->Fit("fitpsi","EM");
  float N3_par1 = fitpsi->GetParameter(0); 
  float r3par1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per1->Fit("fitpsi","EM");
  float N3_per1 = fitpsi->GetParameter(0); 
  float r3per1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par1->Fit("fitpsi","EM");
  float N4_par1 = fitpsi->GetParameter(0); 
  float r4par1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per1->Fit("fitpsi","EM");
  float N4_per1 = fitpsi->GetParameter(0); 
  float r4per1 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsipluspip_par2->Fit("fitpsi","EM");
  float N5_par2 = fitpsi->GetParameter(0); 
  float r5par2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per2->Fit("fitpsi","EM");
  float N5_per2 = fitpsi->GetParameter(0); 
  float r5per2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par2->Fit("fitpsi","EM");
  float N6_par2 = fitpsi->GetParameter(0); 
  float r6par2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per2->Fit("fitpsi","EM");
  float N6_per2 = fitpsi->GetParameter(0); 
  float r6per2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par2->Fit("fitpsi","EM");
  float N7_par2 = fitpsi->GetParameter(0); 
  float r7par2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per2->Fit("fitpsi","EM");
  float N7_per2 = fitpsi->GetParameter(0); 
  float r7per2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par2->Fit("fitpsi","EM");
  float N8_par2 = fitpsi->GetParameter(0); 
  float r8par2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per2->Fit("fitpsi","EM");
  float N8_per2 = fitpsi->GetParameter(0); 
  float r8per2 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsipluspip_par3->Fit("fitpsi","EM");
  float N9_par3 = fitpsi->GetParameter(0); 
  float r9par3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per3->Fit("fitpsi","EM");
  float N9_per3 = fitpsi->GetParameter(0); 
  float r9per3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par3->Fit("fitpsi","EM");
  float NA_par3 = fitpsi->GetParameter(0); 
  float rApar3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per3->Fit("fitpsi","EM");
  float NA_per3 = fitpsi->GetParameter(0); 
  float rAper3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par3->Fit("fitpsi","EM");
  float NB_par3 = fitpsi->GetParameter(0); 
  float rBpar3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per3->Fit("fitpsi","EM");
  float NB_per3 = fitpsi->GetParameter(0); 
  float rBper3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par3->Fit("fitpsi","EM");
  float NC_par3 = fitpsi->GetParameter(0); 
  float rCpar3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per3->Fit("fitpsi","EM");
  float NC_per3 = fitpsi->GetParameter(0); 
  float rCper3 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsipluspip_par4->Fit("fitpsi","EM");
  float ND_par4 = fitpsi->GetParameter(0); 
  float rDpar4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per4->Fit("fitpsi","EM");
  float ND_per4 = fitpsi->GetParameter(0); 
  float rDper4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par4->Fit("fitpsi","EM");
  float NE_par4 = fitpsi->GetParameter(0); 
  float rEpar4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per4->Fit("fitpsi","EM");
  float NE_per4 = fitpsi->GetParameter(0); 
  float rEper4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par4->Fit("fitpsi","EM");
  float NF_par4 = fitpsi->GetParameter(0); 
  float rFpar4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per4->Fit("fitpsi","EM");
  float NF_per4 = fitpsi->GetParameter(0); 
  float rFper4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par4->Fit("fitpsi","EM");
  float NG_par4 = fitpsi->GetParameter(0); 
  float rGpar4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per4->Fit("fitpsi","EM");
  float NG_per4 = fitpsi->GetParameter(0); 
  float rGper4 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsipluspip_par5->Fit("fitpsi","EM");
  float NH_par5 = fitpsi->GetParameter(0); 
  float rHpar5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspip_per5->Fit("fitpsi","EM");
  float NH_per5 = fitpsi->GetParameter(0); 
  float rHper5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_par5->Fit("fitpsi","EM");
  float NI_par5 = fitpsi->GetParameter(0); 
  float rIpar5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsipluspim_per5->Fit("fitpsi","EM");
  float NI_per5 = fitpsi->GetParameter(0); 
  float rIper5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{4}");
  hpsiminuspip_par5->Fit("fitpsi","EM");
  float NJ_par5 = fitpsi->GetParameter(0); 
  float rJpar5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspip_per5->Fit("fitpsi","EM");
  float NJ_per5 = fitpsi->GetParameter(0); 
  float rJper5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_par5->Fit("fitpsi","EM");
  float NK_par5 = fitpsi->GetParameter(0); 
  float rKpar5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");
  hpsiminuspim_per5->Fit("fitpsi","EM");
  float NK_per5 = fitpsi->GetParameter(0); 
  float rKper5 = fitpsi->GetParameter(1); 
  fitpsi->SetParNames("N","#rho^{3}");

  /*********************************************************************** out rho 3 pars **************************************************************************/
  outro3pip << N3_par1 <<"  "<< r3par1 <<"  "<< N3_per1 <<"  "<< r4per1 <<"  "<< N7_par2  <<"  "<< r7par2 <<"  "<< N7_per2 <<"  "<< r7per2 <<"  "<< NB_par3  <<"  "<< rBpar3 <<"  "<< NB_per3 <<"  "<< rBper3 <<"  "<< NF_par4  <<"  "<< rFpar4 <<"  "<< NF_per4 <<"  "<< rFper4 <<"  "<< NJ_par5  <<"  "<< rJpar5 <<"  "<< NJ_per5 <<"  "<< rJper5<<endl;
  outro3pim << N4_par1 <<"  "<< r4par1 <<"  "<< N4_per1 <<"  "<< r4per1 <<"  "<< N8_par2  <<"  "<< r8par2 <<"  "<< N8_per2 <<"  "<< r8per2 <<"  "<< NC_par3  <<"  "<< rCpar3 <<"  "<< NC_per3 <<"  "<< rCper3 <<"  "<< NG_par4  <<"  "<< rGpar4 <<"  "<< NG_per4 <<"  "<< rGper4 <<"  "<< NK_par5  <<"  "<< rKpar5 <<"  "<< NK_per5 <<"  "<< rKper5<<endl;
  /******************************************************************************************************************************************************************/
  /*********************************************************************** out rho 4 pars **************************************************************************/
  outro4pip << N1_par1 <<"  "<< r1par1 <<"  "<< N1_per1 <<"  "<< r1per1 <<"  "<< N5_par2 <<"  "<< r5par2 <<"  "<< N5_per2 <<"  "<< r5per2 <<"  "<< N9_par3  <<"  "<< r9par3 <<"  "<< N9_per3 <<"  "<< r9per3 <<"  "<< ND_par4  <<"  "<< rDpar4 <<"  "<< ND_per4 <<"  "<< rDper4 <<"  "<< NH_par5  <<"  "<< rHpar5 <<"  "<< NH_per5 <<"  "<< rHper5<<endl;
  outro4pim << N2_par1 <<"  "<< r2par1 <<"  "<< N2_per1 <<"  "<< r2per1 <<"  "<< N6_par2 <<"  "<< r6par2 <<"  "<< N6_per2 <<"  "<< r6per2 <<"  "<< NA_par3  <<"  "<< rApar3 <<"  "<< NA_per3 <<"  "<< rAper3 <<"  "<< NE_par4  <<"  "<< rEpar4 <<"  "<< NE_per4 <<"  "<< rEper4 <<"  "<< NI_par5  <<"  "<< rIpar5 <<"  "<< NI_per5 <<"  "<< rIper5<<endl;
  /******************************************************************************************************************************************************************/

  out->WriteTObject(hazimuth_hel_par);  
  out->WriteTObject(hcosthet_hel_par);  
  out->WriteTObject(hazimuth_hel_per);  
  out->WriteTObject(hcosthet_hel_per); 
  out->WriteTObject(hpsipluspip_par);
  out->WriteTObject(hpsipluspip_per);
  out->WriteTObject(hpsipluspim_par);
  out->WriteTObject(hpsipluspim_per);
  out->WriteTObject(hpsiminuspip_par);
  out->WriteTObject(hpsiminuspip_per);
  out->WriteTObject(hpsiminuspim_par);
  out->WriteTObject(hpsiminuspim_per);

  out->WriteTObject(hazim_hel_par1);  
  out->WriteTObject(hazim_hel_per1);  
  out->WriteTObject(hazim_hel_par2);  
  out->WriteTObject(hazim_hel_per2);  
  out->WriteTObject(hazim_hel_par3);  
  out->WriteTObject(hazim_hel_per3);  
  out->WriteTObject(hazim_hel_par4);  
  out->WriteTObject(hazim_hel_per4);  
  out->WriteTObject(hazim_hel_par5);  
  out->WriteTObject(hazim_hel_per5);  

  out->WriteTObject(hpsipluspip_par1);
  out->WriteTObject(hpsipluspip_per1);
  out->WriteTObject(hpsipluspim_par1);
  out->WriteTObject(hpsipluspim_per1);
  out->WriteTObject(hpsiminuspip_par1);
  out->WriteTObject(hpsiminuspip_per1);
  out->WriteTObject(hpsiminuspim_par1);
  out->WriteTObject(hpsiminuspim_per1);
  out->WriteTObject(hpsipluspip_par2);
  out->WriteTObject(hpsipluspip_per2);
  out->WriteTObject(hpsipluspim_par2);
  out->WriteTObject(hpsipluspim_per2);
  out->WriteTObject(hpsiminuspip_par2);
  out->WriteTObject(hpsiminuspip_per2);
  out->WriteTObject(hpsiminuspim_par2);
  out->WriteTObject(hpsiminuspim_per2);
  out->WriteTObject(hpsipluspip_par3);
  out->WriteTObject(hpsipluspip_per3);
  out->WriteTObject(hpsipluspim_par3);
  out->WriteTObject(hpsipluspim_per3);
  out->WriteTObject(hpsiminuspip_par3);
  out->WriteTObject(hpsiminuspip_per3);
  out->WriteTObject(hpsiminuspim_par3);
  out->WriteTObject(hpsiminuspim_per3);
  out->WriteTObject(hpsipluspip_par4);
  out->WriteTObject(hpsipluspip_per4);
  out->WriteTObject(hpsipluspim_par4);
  out->WriteTObject(hpsipluspim_per4);
  out->WriteTObject(hpsiminuspip_par4);
  out->WriteTObject(hpsiminuspip_per4);
  out->WriteTObject(hpsiminuspim_par4);
  out->WriteTObject(hpsiminuspim_per4);
  out->WriteTObject(hpsipluspip_par5);
  out->WriteTObject(hpsipluspip_per5);
  out->WriteTObject(hpsipluspim_par5);
  out->WriteTObject(hpsipluspim_per5);
  out->WriteTObject(hpsiminuspip_par5);
  out->WriteTObject(hpsiminuspip_per5);
  out->WriteTObject(hpsiminuspim_par5);
  out->WriteTObject(hpsiminuspim_per5);

  out->WriteTObject(hassymetry);
  out->WriteTObject(hassymetry_1);
  out->WriteTObject(hassymetry_2);
  out->WriteTObject(hassymetry_3);
  out->WriteTObject(hassymetry_4);
  out->WriteTObject(hassymetry_5);

  out->Close();

}