void FFT_Nev(TString path = "/panda/myResults/AssymIP/DPM/mom_1_5/", TString file_name="/HitsAssym_x_y_1_5GeV.root", int x_sh = -5, int y_sh = -5,int angx_sh = -5, int angy_sh = -5){
  gROOT->Macro("$VMCWORKDIR/macro/lmd/Anastasia/tests/Assym_IP/rootlogon.C");
  gStyle->SetOptStat(0); //no stat
  //read Ideal parameters ---------------------------------------------------
  TString file_ideal = path;
  file_ideal +="x_0mm_y_0mm_Xtlit_0mrad_Ytlit_0mrad";
  file_ideal +=file_name;
  TFile *finIDEAL = new TFile(file_ideal,"READ");
//N
  TGraphErrors *modN_ideal = (TGraphErrors*)finIDEAL->Get("modN");
  modN_ideal->SetMarkerSize(1.0);
  Double_t *N_ideal; Double_t *errN_ideal;
  N_ideal = modN_ideal->GetY();
  errN_ideal = modN_ideal->GetEY();
 double *modulesID = modN_ideal->GetX();
 // double modulesID_doub[40];
 // for(int i=0;i<40;i++){
 //   modulesID_doub[i] = i;
 // }
 const int Nmod = (modN_ideal->GetN());
  const int Nmod_2 = 0.5*Nmod;
  const int Nmod_4 = 0.25*Nmod;
 //calculate frequency
  double frequencyID[Nmod_2];
  frequencyID[0] = 0;
 
  for(int i=1; i<Nmod_2;i++){
    frequencyID[i] = frequencyID[i-1]+double((Nmod-1))/(Nmod*Nmod);
    //  cout<<"frequencyID[i] = "<<frequencyID[i]<<endl;
  }
  // frequencyID[Nmod_2] = 0;
  // for(int i=Nmod_2+1; i<Nmod;i++){
  //   int k=i-Nmod_2;
  //   frequencyID[i] = -(Nmod-1)/(Nmod*modulesID[k]);
  // }
 // //TODO: think automaticly how extend number of point here????
 //  for(int i=1; i<Nmod;i++){
 //    frequencyID[i] = (2*Nmod-1)/(2*Nmod*modulesID_doub[i]);
 //  }
 //  frequencyID[Nmod] = 0;
 //  for(int i=Nmod+1; i<2*Nmod;i++){
 //    int k=i-Nmod;
 //    frequencyID[i] = -(2*Nmod-1)/(2*Nmod*modulesID_doub[k]);
 //  }
 //  for(int i=0;i<2*Nmod;i++)
 //    cout<<"frequencyID["<<i<<"] ="<<frequencyID[i]<<endl;

 //[END] read Ideal parameters ------------------------------------------

  // read Shifted parameters -----------------------------------------------
  TString this_name = "x_";
  this_name +=x_sh;
  this_name +="mm_y_";
  this_name +=y_sh;
  this_name +="mm_Xtlit_";
  this_name +=angx_sh;
  this_name +="mrad_Ytlit_";
  this_name +=angy_sh;
  this_name +="mrad";


TString this_name2 = "x_";
 if(x_sh<0)  this_name2 +="mn";
 else this_name2 +="pl";
 this_name2 +=fabs(x_sh);
  this_name2 +="mm_y_";
if(y_sh<0)  this_name2 +="mn";
 else this_name2 +="pl";
 this_name2 +=fabs(y_sh);
 this_name2 +="mm_Xtlit_";
 if(angx_sh<0) this_name2 +="mn";
 else this_name2 +="pl";
 this_name2 +=fabs(angx_sh);
  this_name2 +="mrad_Ytlit_";
  if(angy_sh<0) this_name2 +="mn";
  else this_name2 +="pl";
  this_name2 +=fabs(angy_sh);
  this_name2 +="mrad";

  TString file_shifted = path;
  file_shifted +=this_name;
  file_shifted +=file_name;
  TFile *fin_SHIFT = new TFile(file_shifted,"READ");

  TGraphErrors *modN_SHIFT = (TGraphErrors*)fin_SHIFT->Get("modN");
  modN_SHIFT->SetMarkerColor(4);
  modN_SHIFT->SetMarkerSize(1.0);
  Double_t *N_SHIFT = new Double_t[Nmod];
  Double_t *errN_SHIFT = new Double_t[Nmod];
  N_SHIFT = modN_SHIFT->GetY();
  errN_SHIFT = modN_SHIFT->GetEY();
  Double_t *diffN_SHIFT = new Double_t[Nmod];
  Double_t *errdiffN_SHIFT = new Double_t[Nmod];
  double N_ideal_norm = N_ideal[0];
  double N_shifted_norm = N_SHIFT[0];
  for(int i=0;i<Nmod;i++){
    N_ideal[i] *= 1./N_ideal_norm;
    N_SHIFT[i] *= 1./N_shifted_norm;
    diffN_SHIFT[i] = 100.*(N_SHIFT[i] - N_ideal[i])/N_ideal[i];
    //   errdiffN_SHIFT[i] = 0.5*sqrt(errN_ideal[i]*errN_ideal[i]+errN_SHIFT[i]*errN_SHIFT[i]);
    errdiffN_SHIFT[i] = (1/N_ideal[i])*sqrt(TMath::Power(N_SHIFT[i]*errN_ideal[i]/N_ideal[i],2)+TMath::Power(errN_SHIFT[i],2));
  }
  //[END] read Shifted parameters ---------------------------------------

  //create histograms (FFT in ROOT works only for histograms) ----------------
  TH1D *h_Nideal = new TH1D("h_Nideal",";sector;N/N_{0};",Nmod,0,Nmod);
  h_Nideal->SetFillColor(kRed);
  h_Nideal->SetLineColor(kRed);
  h_Nideal->SetFillStyle(3004);
  TH1D *h_Nshifted = new TH1D("h_Nshifted",";sector;N/N_{0};",Nmod,0,Nmod);
  h_Nshifted->SetFillColor(kBlue);

  for(int i=0;i<Nmod;i++){
    h_Nideal->Fill(i,N_ideal[i]);
    h_Nshifted->Fill(i,N_SHIFT[i]);
  }
  // for(int i=modN_ideal->GetN();i<2*modN_ideal->GetN();i++){
  //   h_Nideal->Fill(i,0);
  //   h_Nshifted->Fill(i,0);
  // }
  //[END] create histograms (FFT in ROOT works only for histograms) -------


  //Do FFT ---------------------------------------------------------------------------------

  TVirtualFFT::SetTransform(0);
  // int ikind = 1;
  // TVirtualFFT* fft_my = TVirtualFFT::SineCosine(1,&Nmod,&ikind,"M");
   //  TVirtualFFT::SetTransform(fft_my);
  //Compute the transform and look at the magnitude of the output
   TH1 *hm =0;
   hm = h_Nshifted->FFT(hm, "MAG R2C M");
   // hm = h_Nshifted->FFT(hm, "MAG DHT M");
   hm->SetTitle("Magnitude of FFT");
   //  hm->Scale(1./sqrt(Nmod));
   //NOTE: for "real" frequencies you have to divide the x-axes range with the range of your function 
   //(in this case 4*Pi); y-axes has to be rescaled by a factor of 1/SQRT(n) to be right: this is not done automatically!
   
   //Look at the phase of the output   
   TH1 *hp = 0;
   hp = h_Nshifted->FFT(hp, "PH R2C M");
   // hp = h_Nshifted->FFT(hp, "PH DHT M");
   hp->SetTitle("Phase of FFT");

   TH1 *him = 0;
   him = h_Nshifted->FFT(him, "IM R2C M");
   him->SetTitle("Imaginary part of FFT");
   TH1 *hre = 0;
   hre = h_Nshifted->FFT(hre, "RE R2C M");
   hre->SetTitle("Real part of FFT");

 // //That's the way to get the current transform object:
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
   Double_t *re_full = new Double_t[Nmod_2];
   Double_t *im_full = new Double_t[Nmod_2];
   fft->GetPointsComplex(im_full,re_full);
   // im_full[0]=0;//TEST
   // re_full[0]=0;//TEST
   // for(int i=0;i<10;i++){
   //   cout<<i<<" "<<im_full[i]<<" "<<re_full[i]<<endl;
   // } 
   Double_t amp_full[Nmod_2];
   Double_t phase_full[Nmod_2];
   for(int i=0;i<Nmod_2;i++){
     amp_full[i] = sqrt(im_full[i]*im_full[i]+re_full[i]*re_full[i]);
     phase_full[i] = hp->GetBinContent(i);
   }

   TGraph *gr_m_ph = new TGraph(Nmod);
   gr_m_ph->SetMarkerStyle(20);
   gr_m_ph->SetTitle("");

   TGraph *gr_im = new TGraph(Nmod_2,frequencyID,im_full);
   gr_im->SetMarkerStyle(20);
   gr_im->SetTitle("");

   TGraph *gr_re = new TGraph(Nmod_2,frequencyID,re_full);
   gr_re->SetMarkerStyle(20);
   gr_re->SetTitle("");

  
   TGraph *gr_amp = new TGraph(Nmod_2,frequencyID,amp_full);
   gr_amp->SetMarkerStyle(20);
   gr_amp->SetMarkerSize(0.5);
   gr_amp->SetMarkerColor(kBlue);
   gr_amp->SetTitle("");
   TGraph *gr_ph = new TGraph(Nmod_2,frequencyID,phase_full);
   gr_ph->SetMarkerStyle(20);
   gr_ph->SetMarkerSize(0.5);
   gr_ph->SetMarkerColor(kBlue);
   gr_ph->SetTitle("");
   double Ampl_max = 0;
   for(int i=1;i<Nmod+1;i++)
     if(Ampl_max<hm->GetBinContent(i)) 
       Ampl_max = hm->GetBinContent(i);
   TH2D *hm_ph = new TH2D("hm_ph","; phase; magnitude",1e2,-3.15,3.15,1e2,0,1.1*Ampl_max);
   hm_ph->SetMarkerStyle(21);
   hm_ph->SetMarkerSize(1.1);

   //   gr_im_re->SetPoint(0,im_full[i],re_full[i]);
   for(int i=0;i<Nmod+1;i++){
     double ampl = hm->GetBinContent(i);
     double phase = hp->GetBinContent(i);
     gr_m_ph->SetPoint(i,phase,ampl);
     hm_ph->Fill(phase,ampl);
     //    gr_im_re->SetPoint(i,im_full[i],re_full[i]);
   }
 
//[END]Do FFT -------------------------------------------------------------------------
  

  TLegend *leg = new TLegend(0.71,0.75,0.99,0.99);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.04);
  leg->SetHeader("(X,Y)_{POS} [mm]; (X,Y)_{TILT} [mrad]");
 TLegend *leg2 = new TLegend(0.6,0.8,0.98,0.99);
  leg2->SetFillColor(0);
  leg2->SetTextFont(42);
  leg2->SetTextSize(0.04);
  leg2->SetHeader("(X,Y)_{POS} [mm]; (X,Y)_{TILT} [mrad]");
  TString text_sh = "(";
  text_sh +=x_sh;
  text_sh +=",";
  text_sh +=y_sh;
  text_sh += "); (";
  text_sh += angx_sh;
  text_sh +=",";
  text_sh += angy_sh;
  text_sh += ")";
  leg->AddEntry(h_Nideal,"(0,0); (0,0)","f");	  
  leg->AddEntry(h_Nshifted,text_sh,"f");	  
  leg2->AddEntry(gr_amp,text_sh,"p");
  TCanvas *c1 = new TCanvas("Assym", "Assym", 800, 600);
  c1->Divide(2,2);
  c1->cd(1);
  h_Nshifted->GetXaxis()->SetRangeUser(0,20);
  h_Nshifted->Draw();
  h_Nideal->Draw("same");
  leg->Draw();
  c1->cd(3);
  //hm->SetBinContent(1,0);
  hm->Draw();
  //  hm_ph->Draw();
  c1->cd(4);
 // hp->SetBinContent(1,0);
  hp->Draw();
//  gr_ph->Draw("AP");
// gr_ph->GetXaxis()->SetTitle("frequencyID");
//  gr_ph->GetYaxis()->SetTitle("Phase");
leg2->Draw();
 c1->cd(2); 
 gr_amp->Draw("AP");
 // gr_amp->GetXaxis()->SetRangeUser(0,0.5);
 gr_amp->GetXaxis()->SetTitle("frequencyID");
 gr_amp->GetYaxis()->SetTitle("Magnitude");
 leg2->Draw();
//  // hm_ph->SetBinContent(1,0,0);
//  // hm_ph->SetMarkerColor(kBlue);
//  // hm_ph->Draw();
//  // leg2->Draw();
// // gr_im_re->Draw("AP");//TEST
//  //gr_m_ph->Draw("APC");
//  // gr_m_ph->GetXaxis()->SetTitle("phase");
//  // gr_m_ph->GetYaxis()->SetTitle("ampl");
//  c1->cd(5);
//  gr_im->Draw("AP");
//  gr_im->GetXaxis()->SetTitle("frequencyID");
//  gr_im->GetYaxis()->SetTitle("Im");
//  // him->SetBinContent(1,0);
//  // him->Draw();
//  c1->cd(6);
//  gr_re->Draw("AP");
//  gr_re->GetXaxis()->SetTitle("frequencyID");
//  gr_re->GetYaxis()->SetTitle("Re");
//  // hre->SetBinContent(1,0);
//  // hre->Draw();

 TString fout_name = path;
 fout_name +="/AssymResFFT_";
 fout_name += this_name;
 TString fout_name_pdf = fout_name;
 fout_name_pdf  +=".pdf";
 TString fout_name_root = fout_name;
 fout_name_root  +=".root";
 c1->SaveAs(fout_name_pdf);
 TFile *finOUT = new TFile(fout_name_root,"RECREATE");
 h_Nshifted->Write();
 h_Nideal->Write();
 gr_amp->SetName("gr_amp");
 gr_amp->Write();
 leg->Write();
 finOUT->Write();
}