#include "Riostream.h"

void QPO() {
 
  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
  dir.ReplaceAll("QPO.C","");
  dir.ReplaceAll("/./","/");
   
  TFile *f = new TFile("QPO.root","RECREATE");        //crea file root da file ascii
  TTree *T = new TTree("ntuple","data from ascii file");   //crea tree
  Long64_t nlines = T->ReadFile(Form("/home/mastercl/10095-01-02-00.dat",dir.Data()),"freq:errore:power:err_pot:media");
  printf("found %lld points\n",nlines);
  //T->Draw("freq","power");
  T->Write();

  TNtuple *ntu= (TNtuple*)f->Get("ntuple");   //creazione ntupla
  Float_t freq,errore,power,err_pot,media; 
  ntu->SetBranchAddress("freq",&freq);  //frequenza
  ntu->SetBranchAddress("errore",&errore);  //errore frequenza
  ntu->SetBranchAddress("power",&power);   //potenza
  ntu->SetBranchAddress("err_pot",&err_pot);     //errore potenza
  ntu->SetBranchAddress("media",&media);     //numero spettri mediati


  int n = 0; // Indice del vettore di arrivo   
  int n_orig = 0; // Indice del vettore di partenza

  Float_t p = 1.02; //esponente rebin
  Float_t num=1; //numero righe da sommare ad ogni iterazione

  Double_t frequenza_somme=0, err_freq_somme=0,potenza_somme=0,errore_potenza=0, wght_sum=0;      //variabili per somme nel rebin
  const int dim=330;    //dimensione array rebinnato
  const int dim_array=32768;         //dim array originale
  Double_t frequenza[dim_array],err_freq[dim_array],potenza[dim_array],errore_pot[dim_array];      //array di dati
  Double_t frequenza_rebin[dim+1],err_freq_rebin[dim+1],potenza_rebin[dim],err_pot_rebin[dim];
  Int_t rebin[dim];

  TH1F *pds = new TH1F("pds","Spettro di potenza",1000,0,2100);  //creazione istogramma - n bin provvisorio
 
  Int_t nev= ntu->GetEntries();

  int n_orig_length = nlines;

  for(Int_t i=0;i<nev;i++){          //riempimento istogramma originale e array
    ntu->GetEntry(i); 
   
    frequenza[i]=freq;
    potenza[i]=power;
    err_freq[i]=errore;
    errore_pot[i]=err_pot;

    pds->SetBinContent(i+1,potenza[i]);
  } 
 
   while(n_orig<nlines){
     num = int(TMath::Power(p, n + 1)); // numero di canali da mediare

     int n_orig_end = n_orig + num <= n_orig_length ? n_orig + num : n_orig_length;
     int num_real = n_orig_end - n_orig;

     /*
      somma  gli elementi del vettore di partenza sugli indici
    */
    for (int j = n_orig; j < n_orig_end; j++)   { 
      ntu->GetEvent(j);
 
      frequenza_somme += frequenza[j];       //somme
      err_freq_somme += err_freq[j];
      potenza_somme +=( potenza[j] / TMath::Power(errore_pot[j], 2)); //somma pesata
      errore_potenza=errore_potenza + TMath::Power(errore_pot[j],2);
      wght_sum += 1 / TMath::Power(errore_pot[j],2);

    }

    n_orig = n_orig_end;

    frequenza_somme = frequenza_somme / num_real;      //medie
    potenza_somme = potenza_somme / wght_sum;

    frequenza_rebin[n] = frequenza_somme;   //carica su array
    err_freq_rebin[n] = err_freq_somme;
    potenza_rebin[n] = potenza_somme;
    err_pot_rebin[n] = TMath::Sqrt(1/wght_sum);
    rebin[n] = num_real;

    // printf("%f\t %f\t %f\t %f\n",frequenza_somme,potenza_somme,err_freq_somme,err_pot_rebin[n]);

    frequenza_somme=0;       //azzera tutto
    err_freq_somme=0;
    potenza_somme=0;
    errore_potenza=0;
    wght_sum=0;
    n++;               //incrementa n
      
  }

   TH1F *pds_rebin= new TH1F("pds_rebin","Rebin logaritmico",dim,0,2048); //istogramma rebin log

   for(int i=0; i<dim;i++){
    
      pds_rebin->SetBinContent(i+1,potenza_rebin[i]);  //riempi istogramma rebinnato
      pds_rebin->SetBinError(i+1,err_pot_rebin[i]);
   }

  //rimozione fondo alte freq

   TH1F *clone_pds= (TH1F*)pds_rebin->Clone();
   clone_pds->SetName("Clone");

   TF1 *funFondo= new TF1("funFondo",fondoCost,0,2048,1);      //funzione per il fondo
   funFondo->SetParName(0,"Fondo");
   funFondo->SetParameter(0,1);
   pds_rebin->Fit("funFondo","","",1200,2048);    //fit nel range
   Double_t fondo= funFondo->GetParameter(0);

   clone_pds->Add(funFondo,-1,"");     //sottrazione fondo

   TCanvas *canvas= new TCanvas("pds","Spettro di potenza",1000,800);    //canvas
   canvas->Divide(1,2);
      
   canvas->cd(1);      //primo pad, istogramma originale
   gPad->SetLogy();  //scala log log
   gPad->SetLogx();
   pds_rebin->SetTitle("Rebin logaritmico");
   pds_rebin->SetXTitle("Frequenza");
   pds_rebin->SetYTitle("Potenza");
   pds_rebin->Draw("E0");
   pds_rebin->Draw("SAME Chist");
   
   canvas->cd(2);    //secondo pad
   gPad->SetLogy();
   gPad->SetLogx();
   gStyle->SetOptStat(111111);       //mostra fino a skewness e curtosi
   clone_pds->SetTitle("Senza fondo");
   clone_pds->SetXTitle("Frequenza");
   clone_pds->SetYTitle("Potenza");
   clone_pds->Draw("E0");
   clone_pds->Draw("SAME Chist");

   //Funzioni per parametri iniziali
   TF1 *lorentz_fondo= new TF1("lorentz_fondo",LorentzFondo,0,100,2);         
  
   lorentz_fondo->SetLineColor(kRed);
   lorentz_fondo->SetParName(0,"Ampiezza1f");
   lorentz_fondo->SetParameter(0,0.01);
   lorentz_fondo->SetParName(1,"Gamma1f");
   lorentz_fondo->SetParameter(1,0.1);
   clone_pds->Fit("lorentz_fondo","R");
   lorentz_fondo->Draw("SAME");

   TF1 *lorentz_fondo2= new TF1("lorentz_fondo2",LorentzFondo,100,900,2);         
  
   lorentz_fondo2->SetLineColor(kGreen);
   lorentz_fondo2->SetParName(0,"Ampiezza2f");
   lorentz_fondo2->SetParameter(0,1);
   lorentz_fondo2->SetParName(1,"Gamma2f");
   lorentz_fondo2->SetParameter(1,10);
   clone_pds->Fit("lorentz_fondo2","R+");
   lorentz_fondo2->Draw("SAME");

   TF1 *lorentz_fondo3 = new TF1("lorentz_fondo3",LorentzFondo,900,2048,2);         
  
   lorentz_fondo3->SetLineColor(kMagenta);
   lorentz_fondo3->SetParName(0,"Ampiezza3f");
   lorentz_fondo3->SetParameter(0,100);
   lorentz_fondo3->SetParName(1,"Gamma3f");
   lorentz_fondo3->SetParameter(1,100);
   clone_pds->Fit("lorentz_fondo3","R+");
   lorentz_fondo3->Draw("SAME");

   TF1 *lorentz_peak= new TF1("lorentz_peak",Lorentziana,1500,2000,3);   //funzione di fit picco
   lorentz_peak->SetParName(0,"Ampiezza");
   lorentz_peak->SetParameter(0,1);
   lorentz_peak->SetParName(1,"Gamma");
   lorentz_peak->SetParameter(1,100);
   lorentz_peak->SetParName(2,"Centroide");
   lorentz_peak->SetParameter(2,1200);
   clone_pds->Fit("lorentz_peak","R+"); 

   Double_t par[9];                        //attribuisci parametri iniziali ad array
   lorentz_fondo->GetParameters(&par[0]);
   lorentz_fondo2->GetParameters(&par[2]);
   lorentz_fondo3->GetParameters(&par[4]);
   lorentz_peak->GetParameters(&par[5]);

   //Funzione totale
   TF1 *lorentz_1= new TF1("lorentz_1",LorentzFondo,0,2048,2);         
   lorentz_1->SetLineColor(kCyan);
   lorentz_1->SetParameters(par);
   clone_pds->Fit("lorentz_1","R+");
   lorentz_1->Draw("SAME");

   TF1 *lorentz_2= new TF1("lorentz_2",LorentzFondo2,0,2048,4);         
   lorentz_2->SetLineColor(kYellow);
   lorentz_2->SetParameters(par);
   clone_pds->Fit("lorentz_2","R+");
   lorentz_2->Draw("SAME");

   TF1 *totale = new TF1("totale",Totale,0,2048,9);
   
   totale->SetLineColor(kBlack);
   // totale->SetParName();
   totale->SetParameters(par);

   clone_pds->Fit(totale,"R+");
   totale->Draw("SAME");
  
   /*  Bool_t test= Ftest(lorentz_1->GetChisquare(),lorentz_1->GetNDF(),lorentz_2->GetChisquare(),lorentz_2->GetNDF());

   if (test==true){
     printf("true\n");
   }
   else{
     printf("false\n");
     }*/



      /* //definizione frequenze
      Float_t v1=fit->GetParameter(?);
      Float_t v2=fit->GetParameter(?);
      Float_t v3=fit->GetParameter(?);
   
      //identificazione
      Float_t max,inter,min;
      max=v3;
      if(v2>max){
      v2=max;
      v3=inter;
      }
      else v2=inter;
      if(v1>max){
      inter=max;
      min=inter;
      max=v1;
      }
      else {
      if(v1>inter){
      min=inter;
      inter=v1;
      }
      } 

      //stampa valori
      printf("Upper QPO: %f\t Lower QPO: %f\t HBO QPO: %f\n",max,inter,min);

      //delta e gamma
      Float_t delta=(1-inter/max)**2;
      Float_t gamma=(1-min/max)**2;

      //modello
      Double_t raggio= raggio(delta,gamma);
      printf("Raggio: %g\n", raggio);
      Double_t spin= spin(delta,gamma,raggio);
      printf("Spin: %g\n", spin);
      Double_t massa= massa(max,raggio,spin);
      printf("Massa: %g\n", massa);



  */



}


   
//Funzioni fit
Double_t Lorentziana(Double_t *x, Double_t *par){
 
  return (0.5*par[0]*par[1]/TMath::Pi())/((x[0]-par[2])*(x[0]-par[2])+.25*par[1]*par[1]);       //picco Lorentziana
}

Double_t fondoCost(Double_t *x, Double_t *par){  //fondo costante

  return (par[0]);   //cost
}

Double_t LorentzFondo(Double_t *x, Double_t *par){

  Float_t freq_fissa= 0.01;     //centroide fissato a 0

  return  (0.5*par[0]*par[1]/TMath::Pi())/((x[0]-freq_fissa)*(x[0]-freq_fissa)+.25*par[1]*par[1]);       //funzione fondo lorentz
}

Double_t LorentzFondo2(Double_t *x, Double_t *par){

  return LorentzFondo(x,par)+LorentzFondo(x,&par[2]);       //funzione fondo 2 lorentziane
}

Double_t LorentzFondo3(Double_t *x, Double_t *par){

  return LorentzFondo(x,par)+LorentzFondo(x,&par[2])+LorentzFondo(x,&par[4]);       //funzione fondo 3 lorentziane
}

Double_t Totale(Double_t *x, Double_t *par){

  return LorentzFondo(x,par)+LorentzFondo(x,&par[2])+LorentzFondo(x,&par[4])+Lorentziana(x,&par[6]);       //funzione totale
}

// funzione F-test

Bool_t Ftest(Double_t chi1, Int_t ndf1, Double_t chi2, Int_t ndf2){

  Float_t F = (chi1/ndf1)/(chi2/ndf2);     //F dati
  printf("chi1:%f\t ndf1:%d\t chi2: %f\t ndf2:%d\n",chi1,ndf1, chi2, ndf2);
  Float_t P_c = 0.01;                     //Probabilità corrispondente alla confidenza 1%
  
  TF1 *distribuzione = new TF1("distribuzione","TMath::Gamma(([0]+[1])/2)/(TMath::Gamma([0]/2)*TMath::Gamma([1]/2))*TMath::Power([0]/[1],[0]/2)*TMath::Power(x,1/(2*([0]-2))/TMath::Power((1+x*[0]/[1]),1/2*[0]+[1]))",0,10);
  distribuzione->SetParameters(ndf1,ndf2);
 
  Double_t probabilita = distribuzione->Integral(0.01,10);        //calcolo probabilità
 
  TCanvas *distr= new TCanvas("distr","F test",1000,800);
  distribuzione->Draw();

  printf("F:%f \t Prob:%f\n",F, probabilita);
  
  if(probabilita < P_c){      //test sul livello di confidenza

    return true;              //aggiungi al modello, bassa probabilita di falso pos
    }
  else {

    return false;            //interrompi modello, alta probabilita falso pos
  }

}


//equazioni modello
Double_t raggio(Double_t delta, Double_t gamma){

  return 0.6*(6-delta-5*gamma+2*sqrt(2*(delta-gamma)*(3-delta-2*gamma)))/((delta+gamma-2)**2);
}

Double_t spin(Double_t delta, Double_t gamma, Double_t raggio){

  return (TMath::(raggio,1.5))/4*(delta+gamma-2+6/raggio);  //tralasciato +/-
}

Double_t massa(Double_t max,Double_t raggio, Double_t spin){

  return (3.2e4)/(max*(TMath::(raggio,1.5)+spin));  //tralasciato +/-
}