TPC/Upgrade/macros/makeCurrentCorrection.C

   1 /*
   2 .x $HOME/rootlogon.C
   3 .L $ALICE_ROOT/TPC/Upgrade/macros/makeCurrentCorrection.C+
   4
   5 */
   6
   7 #include "TMath.h"
   8 #include "TRandom.h"
   9 #include "TTreeStream.h"
  10 #include "TVectorD.h"
  11 #include "TCanvas.h"
  12 #include "TStopwatch.h"
  13 #include "AliTPCParam.h"
  14 #include "AliTPCcalibDB.h"
  15 #include "AliTPCAltroMapping.h"
  16 #include "AliAltroRawStream.h"
  17 #include "AliSysInfo.h"
  18 #include "AliTPCRawStreamV3.h"
  19 #include "AliCDBManager.h"
  20 #include "TGeoGlobalMagField.h"
  21 #include "AliMagF.h"
  22 #include "AliRawReaderRoot.h"
  23 #include "AliRawReader.h"
  24 #include "TH3.h"
  25 #include "TH2.h"
  26 #include "AliTPCCalPad.h"
  27 #include "AliTPCCalROC.h"
  28 #include "TChain.h"
  29 #include "AliXRDPROOFtoolkit.h"
  30 #include "TLegend.h"
  31 #include "TCut.h"
  32 #include "TGraphErrors.h"
  33 #include "TStatToolkit.h"
  34 #include "TF2.h"
  35
  36 #include "AliDCSSensor.h"
  37 #include "AliCDBEntry.h"
  38 #include "AliDCSSensorArray.h"
  39 #include "TStyle.h"
  40 #include "AliTPCSpaceCharge3D.h"
  41 #include "AliExternalTrackParam.h"
  42 #include "AliTrackerBase.h"
  43 #include "TDatabasePDG.h"
  44 #include "TROOT.h"
  45 #include "AliMathBase.h"
  46 #include "TLatex.h"
  47 //
  48
  49 const Int_t kColors[6]={1,2,3,4,6,7};
  50 const Int_t kMarkers[6]={20,21,24,25,24,25};
  51
  52 void MakeLocalDistortionCurrentTree(Int_t npointsZ=20000, Int_t id=0);
  53
  54 void makeCurrentCorrection(Int_t action=0, Int_t arg0=5000, Int_t arg1=0){
  55   //
  56   //
  57   //
  58   if (action==0) MakeLocalDistortionCurrentTree(arg0,arg1);
  59 }
  60
  61 void SetGraphTDRStyle(TGraph * graph){
  62   graph->GetXaxis()->SetLabelSize(0.08);
  63   graph->GetXaxis()->SetTitleSize(0.08);
  64   graph->GetYaxis()->SetLabelSize(0.08);
  65   graph->GetYaxis()->SetTitleSize(0.08);
  66   graph->GetXaxis()->SetNdivisions(510);
  67   graph->GetYaxis()->SetNdivisions(505);
  68 }
  69
  70 void MakeLocalDistortionCurrentTree(Int_t npointsZ, Int_t id){
  71   //
  72   // Macro to make trees with local distortions
  73   // Results are later visualized in the function DrawLocalDistortionPlots()
  74   //
  75   /*
  76     Int_t npointsZ=1000;
  77     Int_t id=0;
  78   */
  79   Int_t nitteration=300;
  80   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
  81   //
  82   TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
  83   TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
  84   //
  85   AliTPCSpaceCharge3D* distortion = ( AliTPCSpaceCharge3D*)fCurrent->Get("DistRef");
  86   AliTPCSpaceCharge3D* distortionRef = ( AliTPCSpaceCharge3D*)fRef->Get("DistRef");
  87   //
  88   Int_t nz= distortion->GetInputSpaceCharge3D()->GetZaxis()->GetNbins();
  89   TH3 *hisOrig = (TH3*)distortionRef->GetInputSpaceCharge3D();
  90   printf("Make mean histo\n");
  91   TStopwatch timer;
  92   //
  93   for (Int_t iside=0; iside<2; iside++){
  94     for (Int_t iphi=1; iphi<distortion->GetInputSpaceCharge3D()->GetXaxis()->GetNbins(); iphi+=3){
  95       for (Int_t ir=1; ir<distortion->GetInputSpaceCharge3D()->GetYaxis()->GetNbins(); ir+=3){
  96         Double_t sum=0, sumW=0;
  97         if (iside==0) for (Int_t iz2=1; iz2<nz/2; iz2++)   {sum+= hisOrig->GetBinContent(iphi,ir,iz2); sumW++;}
  98         if (iside==1) for (Int_t iz2=nz/2; iz2<nz;  iz2++) {sum+= hisOrig->GetBinContent(iphi,ir,iz2); sumW++;}
  99         //
 100         if (iside==0) for (Int_t iz=1; iz<nz/2; iz++) hisOrig->SetBinContent(iphi,ir,iz,sum/sumW);
 101         if (iside==1) for (Int_t iz=nz/2; iz<=nz; iz++) hisOrig->SetBinContent(iphi,ir,iz,sum/sumW);
 102       }
 103     }
 104   }
 105   timer.Print();
 106   printf("Make mean histo\n");
 107   //
 108   distortion->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
 109   distortionRef->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
 110   //
 111   distortion->AddVisualCorrection(distortion,1);
 112   distortionRef->AddVisualCorrection(distortionRef,2);
 113   //
 114   TVectorD normZR(125), normZRPhi(125), normZZ(125), normZPos(125);
 115   TVectorD normDZR(125), normDZRPhi(125), normDZZ(125);
 116   TVectorD normZRChi2(125), normZRPhiChi2(125), normZZChi2(125);
 117   TVectorD qCurrent(125), qRef(125);
 118   TVectorD qCurrentInt(125), qRefInt(125);
 119
 120   //
 121   for (Int_t iz =0; iz<125; iz++){
 122     for (Int_t iphi=1; iphi<distortion->GetInputSpaceCharge3D()->GetXaxis()->GetNbins(); iphi+=3)
 123       for (Int_t ir=1; ir<distortion->GetInputSpaceCharge3D()->GetYaxis()->GetNbins(); ir+=3){
 124         qCurrent[iz]+=  distortion->GetInputSpaceCharge3D()->GetBinContent(iphi,ir,iz+1);
 125         qRef[iz]+=      distortionRef->GetInputSpaceCharge3D()->GetBinContent(iphi,ir,iz+1);
 126       }
 127   }
 128   //
 129   for (Int_t iz =0; iz<125; iz++){
 130     for (Int_t jz =0; jz<125; jz++){
 131       if (iz<125/2 && jz<=iz){
 132         qCurrentInt[iz]+=qCurrent[jz];
 133         qRefInt[iz]+=qRef[jz];
 134       }
 135       if (iz>125/2 && jz>=iz){
 136         qCurrentInt[iz]+=qCurrent[jz];
 137         qRefInt[iz]+=qRef[jz];
 138       }
 139     }
 140   }
 141   //
 142   //
 143   for (Int_t iz =0; iz<125; iz++){
 144     Double_t z0 = -250+iz*4;
 145     TLinearFitter fitterR(2,"pol1");
 146     TLinearFitter fitterRPhi(2,"pol1");
 147     TLinearFitter fitterZ(2,"pol1");
 148     Double_t xvalue[10]={0};
 149     for (Int_t ipoint =0; ipoint<npointsZ; ipoint++){
 150       Double_t r0   = 95+gRandom->Rndm()*(245-95.);
 151       Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();
 152       if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1))>20) continue;
 153       if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2))>20) continue;
 154       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2);
 155       fitterR.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1));
 156       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2);
 157       fitterRPhi.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1));
 158       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2);
 159       fitterZ.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1));
 160     }
 161     fitterR.Eval();
 162     fitterRPhi.Eval();
 163     fitterZ.Eval();
 164     normZR[iz]=fitterR.GetParameter(1);
 165     normZRPhi[iz]=fitterRPhi.GetParameter(1);
 166     normZZ[iz]=fitterZ.GetParameter(1);
 167     normZRChi2[iz]=TMath::Sqrt(fitterR.GetChisquare()/fitterR.GetNpoints());
 168     normZRPhiChi2[iz]=TMath::Sqrt(fitterRPhi.GetChisquare()/fitterRPhi.GetNpoints());
 169     normZZChi2[iz]=TMath::Sqrt(fitterZ.GetChisquare()/fitterZ.GetNpoints());
 170     //
 171     if (iz>0){
 172       normDZR[iz]=(normZR[iz]-normZR[iz-1]);
 173       normDZRPhi[iz]=(normZRPhi[iz]-normZRPhi[iz-1]);
 174     }
 175     normZPos[iz]=z0;
 176   }
 177   {
 178     (*pcstream)<<"meanCurrent"<<
 179       "id="<<id<<                       // lookup ID
 180       "normZPos.="<<&normZPos<<         // zposition
 181       //
 182       "qCurrent.="<<&qCurrent<<         // current measuremn
 183       "qRef.="<<&qRef<<                 // current in refenece sample
 184       "qCurrentInt.="<<&qCurrentInt<<   // integral of current
 185       "qRefInt.="<<&qRefInt<<           // integral of current
 186       //
 187       //
 188       //
 189       "normZR.="<<&normZR<<             // mult. scaling to minimize R distortions
 190       "normDZR.="<<&normDZR<<           // mult. scaling to minimize R distortions
 191       "normZRPhi.="<<&normZRPhi<<       // mult. scaling
 192       "normDZRPhi.="<<&normDZRPhi<<     // mult. scaling
 193       "normZZ.="<<&normZZ<<
 194       //
 195       "normZRChi2.="<<&normZRChi2<<            // mult. scaling to minimize R distortions
 196       "normZRPhiChi2.="<<&normZRPhiChi2<<      // mult. scaling
 197       "normZZChi2.="<<&normZZChi2<<
 198       "\n";
 199   }
 200   delete pcstream;
 201
 202   pcstream = new TTreeSRedirector("localCurrent.root","update");
 203   TTree * treeNormZ= (TTree*)pcstream->GetFile()->Get("meanCurrent");
 204   TGraphErrors * grZRfit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZR.fElements:normZPos.fElements","",25,2,0.5);
 205   TGraphErrors * grZRPhifit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZRPhi.fElements:normZPos.fElements","",25,4,0.5);
 206   grZRfit->Draw("alp");
 207   grZRPhifit->Draw("lp");
 208 }
 209
 210
 211 void Fit(){
 212   //
 213   // Not good  rsponse should be more complicated
 214   //
 215   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
 216   TTree * treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
 217   //
 218   TVectorD* pnormZR=0, *pnormZRPhi=0, *pnormZZ=0, *pnormZPos=0;
 219   TVectorD* pqCurrent=0, *pqRef=0;
 220   TVectorD* pqCurrentInt=0, *pqRefInt=0;
 221   //
 222   treeCurrent->SetBranchAddress("normZR.",&pnormZR);
 223   treeCurrent->SetBranchAddress("normZPos.",&pnormZPos);
 224   treeCurrent->SetBranchAddress("normZRPhi.",&pnormZRPhi);
 225   treeCurrent->SetBranchAddress("qCurrent.",&pqCurrent);
 226   treeCurrent->SetBranchAddress("qRef.",&pqRef);
 227   treeCurrent->SetBranchAddress("qCurrentInt.",&pqCurrentInt);
 228   treeCurrent->SetBranchAddress("qRefInt.",&pqRefInt);
 229   Int_t entries = treeCurrent->GetEntries();
 230   //
 231   //
 232   for (Double_t sigma=1; sigma<40; sigma++){
 233     for (Int_t entry=0; entry<entries; entry++){
 234       treeCurrent->GetEntry(entry);
 235       //
 236       TVectorD vecCurrentRefInt(125);
 237       for (Int_t i=0; i<125; i++){
 238         Double_t sumW=0;
 239         for (Int_t j=0; j<125; j++){
 240           if (((*pnormZPos)[i]*(*pnormZPos)[j])<0) continue;
 241           Double_t weight = 0;
 242           if ((*pnormZPos)[i]<0) weight=0.5*(TMath::Erf(Double_t(i-j)/Double_t(sigma))+1);
 243           if ((*pnormZPos)[i]>0) weight=0.5*(TMath::Erf(Double_t(j-i)/Double_t(sigma))+1);
 244           vecCurrentRefInt[i]+=weight*(*pqCurrent)[j]/(*pqRef)[j];
 245           sumW+=weight;
 246         }
 247         vecCurrentRefInt[i]/=sumW;
 248       }
 249       (*pcstream)<<"sigmaScan"<<
 250         "entry="<<entry<<
 251         "sigma="<<sigma<<
 252         "normZPos.="<<pnormZPos<<
 253         "normZR.="<<pnormZR<<
 254         "normZRPhi.="<<pnormZRPhi<<
 255         "vecCurrent.="<<&vecCurrentRefInt<<
 256         "\n";
 257     }
 258   }
 259   delete pcstream;
 260   pcstream = new TTreeSRedirector("localCurrent.root","update");
 261   TTree * treeScan= (TTree*)pcstream->GetFile()->Get("sigmaScan");
 262   treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
 263   treeScan->SetMarkerStyle(25);
 264   treeScan->SetMarkerSize(0.4);
 265   treeScan->Draw("vecCurrent.fElements:normZR.fElements:sigma","abs(normZPos.fElements-100)<20&&sigma>10","colz");
 266
 267 }
 268
 269 void FitLinear(){
 270   //
 271   // deltaX_i = A_ij*deltaI_j
 272   //
 273   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
 274   TTree * treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
 275   //
 276   TVectorD* pnormZR=0, *pnormZRPhi=0, *pnormZZ=0, *pnormZPos=0;
 277   TVectorD* pqCurrent=0, *pqRef=0;
 278   TVectorD* pqCurrentInt=0, *pqRefInt=0;
 279   //
 280   treeCurrent->SetBranchAddress("normZR.",&pnormZR);
 281   treeCurrent->SetBranchAddress("normZPos.",&pnormZPos);
 282   treeCurrent->SetBranchAddress("normZRPhi.",&pnormZRPhi);
 283   treeCurrent->SetBranchAddress("qCurrent.",&pqCurrent);
 284   treeCurrent->SetBranchAddress("qRef.",&pqRef);
 285   treeCurrent->SetBranchAddress("qCurrentInt.",&pqCurrentInt);
 286   treeCurrent->SetBranchAddress("qRefInt.",&pqRefInt);
 287   Int_t entries = treeCurrent->GetEntries();
 288   //
 289   //
 290   //
 291   Int_t nParamSide=62;               // number of z bins on 1 side
 292   Int_t group=3;                     // grouping of bins
 293   Int_t nParams=nParamSide/group;    // parameters
 294   Int_t nParamsFit=nParams*nParams;  //
 295   TLinearFitter fitter(nParamsFit+1,TString::Format("hyp%d",nParamsFit));
 296   TVectorD xVector(nParamsFit+1);
 297   TVectorD pVector(nParamsFit+1);
 298   //
 299   for (Int_t ievent=0; ievent<entries; ievent++){
 300     treeCurrent->GetEntry(ievent);
 301     for (Int_t iz=0; iz<nParamSide; iz++){
 302       Int_t dPar=iz/group;
 303       if (dPar>nParams-1) dPar=nParams-1;
 304       // 1.) clear X vectors
 305       for (Int_t ipar=0; ipar<nParamsFit; ipar++) xVector[ipar]=0;
 306       // 2.) set   X vector of interest
 307       for (Int_t cPar=0; cPar<nParams; cPar++){
 308         xVector[dPar*nParams+cPar] += (*pqCurrent)[cPar*group]/(*pqRef)[cPar*group];
 309         xVector[dPar*nParams+cPar] += (*pqCurrent)[TMath::Min(cPar*group+1,nParamSide)]/(*pqRef)[TMath::Min(cPar*group+1,nParamSide)];
 310         xVector[dPar*nParams+cPar] += (*pqCurrent)[TMath::Max(cPar*group-1,0)]/(*pqRef)[TMath::Max(cPar*group-1,0)];
 311       }
 312       Double_t val = 0;
 313       val+=(*pnormZR)[dPar*group];
 314       val+=(*pnormZR)[TMath::Max(dPar*group-1,0)];
 315       val+=(*pnormZR)[TMath::Min(dPar*group,nParamSide)];
 316       val/=3.;
 317       fitter.AddPoint(xVector.GetMatrixArray(),val,0.0035);
 318     }
 319   }
 320   for (Int_t cPar=1; cPar<nParams-1; cPar++){
 321     //
 322     for (Int_t ipar=0; ipar<nParamsFit; ipar++) xVector[ipar]=0;
 323     xVector[(cPar-1)*nParams+cPar-1]=0.5;
 324     xVector[(cPar)*nParams+cPar]    =-1;
 325     xVector[(cPar+1)*nParams+cPar+1]=0.5;
 326     fitter.AddPoint(xVector.GetMatrixArray(),0,0.05);
 327   }
 328   fitter.Eval();
 329   //
 330   //
 331   TVectorD fitVector(nParamsFit);
 332   TVectorD czVector(nParamsFit);
 333   TVectorD fitVectorErr(nParamsFit);
 334   fitter.GetParameters(fitVector);
 335   for (Int_t iPar=0; iPar<nParamsFit; iPar++)  {
 336     pVector[iPar]=iPar;
 337     czVector[iPar]=250.*((iPar-1)%nParams)/nParams;
 338     fitVectorErr[iPar]=fitter.GetParError(iPar);
 339   }
 340   Double_t chi2= TMath::Sqrt(fitter.GetChisquare()/fitter.GetNpoints());
 341   printf("Chi2=%f\n",chi2);
 342   TGraphErrors * gr  = new TGraphErrors(nParamsFit,  pVector.GetMatrixArray(), fitVector.GetMatrixArray(), 0,fitVectorErr.GetMatrixArray());
 343   gr->SetMarkerStyle(25); gr->SetMarkerSize(0.3);
 344   gr->Draw("alp");
 345
 346   TGraphErrors * vgraphs[20]={0};
 347   for (Int_t ipar=0; ipar<20; ipar++){
 348     vgraphs[ipar]=new TGraphErrors(nParams, &(czVector.GetMatrixArray()[ipar*nParams+1]), &(fitVector.GetMatrixArray()[ipar*nParams+1]), 0,  &(fitVectorErr.GetMatrixArray()[ipar*nParams+1]));
 349     vgraphs[ipar]->GetXaxis()->SetTitle("z_{I} (cm)");
 350     vgraphs[ipar]->GetYaxis()->SetTitle("#it{A}_{i_{I},j_{#DeltaR}}");
 351     vgraphs[ipar]->SetMinimum(0);
 352     vgraphs[ipar]->SetMaximum(0.1);
 353     SetGraphTDRStyle(vgraphs[ipar]);
 354   }
 355
 356   TCanvas * canvasFit = new TCanvas("canvasFit","canvasFit",700,700);
 357   canvasFit->SetRightMargin(0.05);
 358   canvasFit->SetLeftMargin(0.15);
 359   canvasFit->SetBottomMargin(0.18);
 360   canvasFit->Divide(1,2,0,0);
 361   TLegend * legend0 = new TLegend(0.4,0.4,0.95,0.95,"Current fluctuation correction matrix. #DeltaR_{zi}=A_{ij}I_{zj}");
 362   legend0->SetBorderSize(0);
 363   TLegend * legend1 = new TLegend(0.4,0.5,0.95,0.95,"Current fluctuation correction matrix. #DeltaR_{zi}=A_{ij}I_{zj}");
 364   legend1->SetBorderSize(0);
 365
 366   for (Int_t ipar=0; ipar<10; ipar++){
 367     canvasFit->cd(ipar/5+1)->SetTicks(3,3);
 368     vgraphs[ipar*2]->SetMarkerStyle(kMarkers[ipar%5]);
 369     vgraphs[ipar*2]->SetMarkerColor(kColors[ipar%5]);
 370     vgraphs[ipar*2]->SetLineColor(kColors[ipar%5]);
 371     if (ipar%5==0) vgraphs[ipar*2]->Draw("alp");
 372     vgraphs[ipar*2]->Draw("lp");
 373     if (ipar<5) legend0->AddEntry(vgraphs[ipar*2],TString::Format("z_{#DeltaR}=%1.0f (cm)",250.*ipar/10.),"p");
 374     if (ipar>=5) legend1->AddEntry(vgraphs[ipar*2],TString::Format("z_{#DeltaR}=%1.0f (cm)",250.*ipar/10.),"p");
 375   }
 376   legend0->SetTextSize(0.05);
 377   legend1->SetTextSize(0.05);
 378   //
 379   canvasFit->cd(1);
 380   legend0->Draw();
 381   canvasFit->cd(2);
 382   legend1->Draw();
 383   canvasFit->SaveAs("canvasAIJ_CurrenToDR.pdf");
 384   canvasFit->SaveAs("canvasAIJ_CurrenToDR.png");
 385 }
 386
 387
 388 void FFTexample(){
 389   //
 390   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
 391   //
 392   TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
 393   TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
 394   //
 395   AliTPCSpaceCharge3D* distortion = ( AliTPCSpaceCharge3D*)fCurrent->Get("DistRef");
 396   AliTPCSpaceCharge3D* distortionRef = ( AliTPCSpaceCharge3D*)fRef->Get("DistRef");
 397
 398   distortion->AddVisualCorrection(distortion,1);
 399   distortionRef->AddVisualCorrection(distortionRef,2);
 400   //
 401   TF2 f2("f2","AliTPCCorrection::GetCorrXYZ(x,0,y,0,1)-1.07*AliTPCCorrection::GetCorrXYZ(x,0,y,0,2)",85,245,0,250);
 402   TF2 f2c("f2c","AliTPCCorrection::GetCorrXYZ(x,0,y,0,1)-1.07*AliTPCCorrection::GetCorrXYZ(x,0,y,0,2)+sqrt(12.)*(rndm-0.5)*0.04",85,245,0,250);
 403   f2.SetNpx(16); f2.SetNpy(32);
 404   f2c.SetNpx(16); f2c.SetNpy(32);
 405   f2.Draw("colz");
 406   f2c.Draw("colz");
 407
 408   TH2D * his2D = (TH2D*)f2.GetHistogram();
 409   TH2D * his2Dc = (TH2D*)f2c.GetHistogram();
 410   TH2D * his2DSmooth=new TH2D(*his2Dc);
 411   Double_t sigma = 3;
 412   Int_t nx=his2D->GetXaxis()->GetNbins();
 413   Int_t ny=his2D->GetYaxis()->GetNbins();
 414   //
 415   for (Int_t ibin=2; ibin<=nx-1;ibin++){
 416     for (Int_t jbin=2; jbin<ny-1;jbin++){
 417       TLinearFitter fitter(4,"hyp3");
 418       for (Int_t di=-3; di<=3; di++)
 419         for (Int_t dj=-3; dj<=3; dj++){
 420           if (ibin+di<=1) continue;
 421           if (jbin+dj<=1) continue;
 422           if (ibin+di>=nx) continue;
 423           if (jbin+dj>=ny) continue;
 424           Double_t x[2]={di,dj,dj*dj};
 425           Double_t w = 1./(TMath::Gaus(di/sigma)*TMath::Gaus(dj/sigma));
 426           Double_t y = his2Dc->GetBinContent(ibin+di,jbin+dj);
 427           fitter.AddPoint(x,y,w);
 428         }
 429       fitter.Eval();
 430       printf("%d\t%d\t%3.3f\n",ibin,jbin,10*(fitter.GetParameter(0)-his2Dc->GetBinContent(ibin,jbin)));
 431       his2DSmooth->SetBinContent(ibin,jbin,fitter.GetParameter(0));
 432     }
 433   }
 434   TH2D hisDiff=*his2DSmooth;
 435   hisDiff.Add(his2D,-1);
 436
 437
 438   TH1 * hisMag = his2D->FFT(0,"MAG");
 439   TH1 * hisMagc = his2Dc->FFT(0,"MAG");
 440   TH1 * hisPhi = his2D->FFT(0,"PH");
 441   TH1 * hisPhic = his2Dc->FFT(0,"PH");
 442   hisMag->Draw("surf2");
 443   hisMagc->Draw("surf2");
 444
 445
 446 }
 447
 448 /*
 449 void MakeFFT(){
 450   //
 451   //
 452   //
 453   Int_t nSize = 125;
 454   TVirtualFFT *fft_own = TVirtualFFT::FFT(1, &nSize, "R2C ES K");
 455   fft_own->SetPoints(normZR->GetMatrixArray());
 456   fft_own->Transform();
 457
 458   TH1 *hre = 0, *hco=0;
 459   hre = TH1::TransformHisto(fft_own, hre, "RE");
 460   hco = TH1::TransformHisto(fft_own, hco, "C");
 461
 462
 463   hr->SetTitle("Real part of the 3rd (array) tranfsorm");
 464   hr->Draw();
 465   hr->SetStats(kFALSE);
 466   hr->GetXaxis()->SetLabelSize(0.05);
 467   hr->GetYaxis()->SetLabelSize(0.05);
 468   c1_6->cd();
 469   TH1 *him = 0;
 470   him = TH1::TransformHisto(fft_own, him, "IM");
 471   him->SetTitle("Im. part of the 3rd (array) transform");
 472   him->Draw();
 473   him->SetStats(kFALSE);
 474   him->GetXaxis()->SetLabelSize(0.05);
 475   him->GetYaxis()->SetLabelSize(0.05);
 476   //
 477 }
 478 */