Modification - default figures for the smoothing studies
[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / makeCurrentCorrection.C
1 /*
2 .x $HOME/rootlogon.C
3 .x $ALICE_ROOT/TPC/Upgrade/macros/NimStyle.C
4
5 .L $ALICE_ROOT/TPC/Upgrade/macros/makeCurrentCorrection.C+
6
7 */
8
9 #include "TMath.h"
10 #include "TRandom.h"
11 #include "TTreeStream.h"
12 #include "TVectorD.h"
13 #include "TCanvas.h"
14 #include "TStopwatch.h"
15 #include "AliTPCParam.h"
16 #include "AliTPCcalibDB.h"
17 #include "AliTPCAltroMapping.h"
18 #include "AliAltroRawStream.h"
19 #include "AliSysInfo.h"
20 #include "AliTPCRawStreamV3.h"
21 #include "AliCDBManager.h"
22 #include "TGeoGlobalMagField.h"
23 #include "AliMagF.h"
24 #include "AliRawReaderRoot.h"
25 #include "AliRawReader.h"
26 #include "TH3.h"
27 #include "TH2.h"
28 #include "AliTPCCalPad.h"
29 #include "AliTPCCalROC.h"
30 #include "TChain.h"
31 #include "AliXRDPROOFtoolkit.h"
32 #include "TLegend.h"
33 #include "TCut.h"
34 #include "TGraphErrors.h"
35 #include "TStatToolkit.h"
36 #include "TF2.h"
37
38 #include "AliDCSSensor.h"
39 #include "AliCDBEntry.h"
40 #include "AliDCSSensorArray.h"
41 #include "TStyle.h"
42 #include "AliTPCSpaceCharge3D.h"
43 #include "AliExternalTrackParam.h"
44 #include "AliTrackerBase.h"
45 #include "TDatabasePDG.h"
46 #include "TROOT.h"
47 #include "AliMathBase.h"
48 #include "TLatex.h"
49 #include "AliTPCCorrectionLookupTable.h"
50 //
51
52 const Int_t kColors[6]={1,2,3,4,6,7};
53 const Int_t kMarkers[6]={20,21,24,25,24,25};
54 TObjArray garrayFit(3);
55
56
57 void MakeLocalDistortionCurrentTree(Int_t npointsZ=20000, Int_t id=0);
58 void MakeSmoothKernelStudy(Int_t npoints, Int_t nkernels, Int_t mapID);
59
60 void makeCurrentCorrection(Int_t action=0, Int_t arg0=5000, Int_t arg1=0, Int_t arg2=0){
61   //
62   //
63   //
64   if (action==0) MakeLocalDistortionCurrentTree(arg0,arg1);
65   if (action==1) MakeSmoothKernelStudy(arg0,arg1,arg2);
66  
67 }
68
69 void SetGraphTDRStyle(TGraph * graph){
70   graph->GetXaxis()->SetLabelSize(0.08);
71   graph->GetXaxis()->SetTitleSize(0.08);
72   graph->GetYaxis()->SetLabelSize(0.08);
73   graph->GetYaxis()->SetTitleSize(0.08);
74   graph->GetXaxis()->SetNdivisions(510);
75   graph->GetYaxis()->SetNdivisions(505);
76 }
77
78 void MakeLocalDistortionCurrentTree(Int_t npointsZ, Int_t id){
79   //
80   // Macro to make trees with local distortions 
81   // Results are later visualized in the function DrawLocalDistortionPlots()
82   //
83   /*
84     Int_t npointsZ=1000;
85     Int_t id=0;
86   */
87   Int_t nitteration=300;
88   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
89   //
90   TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
91   TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
92   //
93   AliTPCSpaceCharge3D* distortion = ( AliTPCSpaceCharge3D*)fCurrent->Get("DistRef"); 
94   AliTPCSpaceCharge3D* distortionRef = ( AliTPCSpaceCharge3D*)fRef->Get("DistRef"); 
95   //
96   Int_t nz= distortion->GetInputSpaceCharge3D()->GetZaxis()->GetNbins();
97   TH3 *hisOrig = (TH3*)distortionRef->GetInputSpaceCharge3D();
98   printf("Make mean histo\n");
99   TStopwatch timer;
100   //
101   for (Int_t iside=0; iside<2; iside++){
102     for (Int_t iphi=1; iphi<distortion->GetInputSpaceCharge3D()->GetXaxis()->GetNbins(); iphi+=3){
103       for (Int_t ir=1; ir<distortion->GetInputSpaceCharge3D()->GetYaxis()->GetNbins(); ir+=3){
104         Double_t sum=0, sumW=0;
105         if (iside==0) for (Int_t iz2=1; iz2<nz/2; iz2++)   {sum+= hisOrig->GetBinContent(iphi,ir,iz2); sumW++;}
106         if (iside==1) for (Int_t iz2=nz/2; iz2<nz;  iz2++) {sum+= hisOrig->GetBinContent(iphi,ir,iz2); sumW++;}
107         //
108         if (iside==0) for (Int_t iz=1; iz<nz/2; iz++) hisOrig->SetBinContent(iphi,ir,iz,sum/sumW);
109         if (iside==1) for (Int_t iz=nz/2; iz<=nz; iz++) hisOrig->SetBinContent(iphi,ir,iz,sum/sumW);
110       }
111     }
112   }
113   timer.Print();
114   printf("Make mean histo\n");
115   //
116   distortion->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
117   distortionRef->InitSpaceCharge3DPoisson(129, 129, 144,nitteration);
118   //
119   distortion->AddVisualCorrection(distortion,1);
120   distortionRef->AddVisualCorrection(distortionRef,2);
121   //
122   TVectorD normZR(125), normZRPhi(125), normZZ(125), normZPos(125);
123   TVectorD normDZR(125), normDZRPhi(125), normDZZ(125);
124   TVectorD normZRChi2(125), normZRPhiChi2(125), normZZChi2(125);
125   TVectorD qCurrent(125), qRef(125);
126   TVectorD qCurrentInt(125), qRefInt(125);
127   
128   //
129   for (Int_t iz =0; iz<125; iz++){
130     for (Int_t iphi=1; iphi<distortion->GetInputSpaceCharge3D()->GetXaxis()->GetNbins(); iphi+=3)
131       for (Int_t ir=1; ir<distortion->GetInputSpaceCharge3D()->GetYaxis()->GetNbins(); ir+=3){
132         qCurrent[iz]+=  distortion->GetInputSpaceCharge3D()->GetBinContent(iphi,ir,iz+1);
133         qRef[iz]+=      distortionRef->GetInputSpaceCharge3D()->GetBinContent(iphi,ir,iz+1);            
134       }   
135   } 
136   //
137   for (Int_t iz =0; iz<125; iz++){
138     for (Int_t jz =0; jz<125; jz++){
139       if (iz<125/2 && jz<=iz){
140         qCurrentInt[iz]+=qCurrent[jz];
141         qRefInt[iz]+=qRef[jz];
142       }
143       if (iz>125/2 && jz>=iz){
144         qCurrentInt[iz]+=qCurrent[jz];
145         qRefInt[iz]+=qRef[jz];
146       }
147     }
148   }
149   //    
150   //
151   for (Int_t iz =0; iz<125; iz++){
152     Double_t z0 = -250+iz*4;
153     TLinearFitter fitterR(2,"pol1");
154     TLinearFitter fitterRPhi(2,"pol1");
155     TLinearFitter fitterZ(2,"pol1");
156     Double_t xvalue[10]={0};
157     for (Int_t ipoint =0; ipoint<npointsZ; ipoint++){      
158       Double_t r0   = 95+gRandom->Rndm()*(245-95.);
159       Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();      
160       if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1))>20) continue;
161       if (TMath::Abs(distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2))>20) continue;
162       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2);
163       fitterR.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1));
164       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2);
165       fitterRPhi.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1));
166       xvalue[0]=distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2);
167       fitterZ.AddPoint(xvalue,distortion->GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1));
168     }    
169     fitterR.Eval();
170     fitterRPhi.Eval();
171     fitterZ.Eval();    
172     normZR[iz]=fitterR.GetParameter(1);
173     normZRPhi[iz]=fitterRPhi.GetParameter(1);
174     normZZ[iz]=fitterZ.GetParameter(1);
175     normZRChi2[iz]=TMath::Sqrt(fitterR.GetChisquare()/fitterR.GetNpoints());    
176     normZRPhiChi2[iz]=TMath::Sqrt(fitterRPhi.GetChisquare()/fitterRPhi.GetNpoints());    
177     normZZChi2[iz]=TMath::Sqrt(fitterZ.GetChisquare()/fitterZ.GetNpoints());    
178     //
179     if (iz>0){
180       normDZR[iz]=(normZR[iz]-normZR[iz-1]);
181       normDZRPhi[iz]=(normZRPhi[iz]-normZRPhi[iz-1]);
182     }
183     normZPos[iz]=z0;    
184   }
185   {    
186     (*pcstream)<<"meanCurrent"<<
187       "id="<<id<<                       // lookup ID
188       "normZPos.="<<&normZPos<<         // zposition
189       //
190       "qCurrent.="<<&qCurrent<<         // current measuremn 
191       "qRef.="<<&qRef<<                 // current in refenece sample
192       "qCurrentInt.="<<&qCurrentInt<<   // integral of current
193       "qRefInt.="<<&qRefInt<<           // integral of current 
194       //
195       //
196       //
197       "normZR.="<<&normZR<<             // mult. scaling to minimize R distortions
198       "normDZR.="<<&normDZR<<           // mult. scaling to minimize R distortions
199       "normZRPhi.="<<&normZRPhi<<       // mult. scaling 
200       "normDZRPhi.="<<&normDZRPhi<<     // mult. scaling 
201       "normZZ.="<<&normZZ<<
202       //
203       "normZRChi2.="<<&normZRChi2<<            // mult. scaling to minimize R distortions
204       "normZRPhiChi2.="<<&normZRPhiChi2<<      // mult. scaling 
205       "normZZChi2.="<<&normZZChi2<<
206       "\n";
207   }
208   delete pcstream;
209   
210   pcstream = new TTreeSRedirector("localCurrent.root","update");
211   TTree * treeNormZ= (TTree*)pcstream->GetFile()->Get("meanCurrent");
212   TGraphErrors * grZRfit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZR.fElements:normZPos.fElements","",25,2,0.5);
213   TGraphErrors * grZRPhifit= TStatToolkit::MakeGraphErrors( treeNormZ, "normZRPhi.fElements:normZPos.fElements","",25,4,0.5);
214   grZRfit->Draw("alp");
215   grZRPhifit->Draw("lp");
216 }
217
218
219 void Fit(){
220   //
221   // Not good  rsponse should be more complicated
222   //
223   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
224   TTree * treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
225   //
226   TVectorD* pnormZR=0, *pnormZRPhi=0, *pnormZZ=0, *pnormZPos=0;
227   TVectorD* pqCurrent=0, *pqRef=0;
228   TVectorD* pqCurrentInt=0, *pqRefInt=0;
229   //
230   treeCurrent->SetBranchAddress("normZR.",&pnormZR);
231   treeCurrent->SetBranchAddress("normZPos.",&pnormZPos);
232   treeCurrent->SetBranchAddress("normZRPhi.",&pnormZRPhi);
233   treeCurrent->SetBranchAddress("qCurrent.",&pqCurrent);
234   treeCurrent->SetBranchAddress("qRef.",&pqRef);
235   treeCurrent->SetBranchAddress("qCurrentInt.",&pqCurrentInt);
236   treeCurrent->SetBranchAddress("qRefInt.",&pqRefInt);
237   Int_t entries = treeCurrent->GetEntries();
238   //
239   //
240   for (Double_t sigma=1; sigma<40; sigma++){    
241     for (Int_t entry=0; entry<entries; entry++){
242       treeCurrent->GetEntry(entry);
243       //
244       TVectorD vecCurrentRefInt(125);
245       for (Int_t i=0; i<125; i++){
246         Double_t sumW=0;
247         for (Int_t j=0; j<125; j++){
248           if (((*pnormZPos)[i]*(*pnormZPos)[j])<0) continue;
249           Double_t weight = 0;
250           if ((*pnormZPos)[i]<0) weight=0.5*(TMath::Erf(Double_t(i-j)/Double_t(sigma))+1);
251           if ((*pnormZPos)[i]>0) weight=0.5*(TMath::Erf(Double_t(j-i)/Double_t(sigma))+1);
252           vecCurrentRefInt[i]+=weight*(*pqCurrent)[j]/(*pqRef)[j];        
253           sumW+=weight;
254         }
255         vecCurrentRefInt[i]/=sumW;
256       }
257       (*pcstream)<<"sigmaScan"<<
258         "entry="<<entry<<
259         "sigma="<<sigma<<
260         "normZPos.="<<pnormZPos<<
261         "normZR.="<<pnormZR<<
262         "normZRPhi.="<<pnormZRPhi<<
263         "vecCurrent.="<<&vecCurrentRefInt<<
264         "\n";
265     }
266   }
267   delete pcstream; 
268   pcstream = new TTreeSRedirector("localCurrent.root","update");
269   TTree * treeScan= (TTree*)pcstream->GetFile()->Get("sigmaScan"); 
270   treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
271   treeScan->SetMarkerStyle(25);
272   treeScan->SetMarkerSize(0.4); 
273   treeScan->Draw("vecCurrent.fElements:normZR.fElements:sigma","abs(normZPos.fElements-100)<20&&sigma>10","colz");
274
275 }
276
277 void FitLinear(){
278   //
279   // deltaX_i = A_ij*deltaI_j
280   // 
281   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
282   TTree * treeCurrent= (TTree*)pcstream->GetFile()->Get("meanCurrent");
283   //
284   TVectorD* pnormZR=0, *pnormZRPhi=0, *pnormZZ=0, *pnormZPos=0;
285   TVectorD* pqCurrent=0, *pqRef=0;
286   TVectorD* pqCurrentInt=0, *pqRefInt=0;
287   //
288   treeCurrent->SetBranchAddress("normZR.",&pnormZR);
289   treeCurrent->SetBranchAddress("normZPos.",&pnormZPos);
290   treeCurrent->SetBranchAddress("normZRPhi.",&pnormZRPhi);
291   treeCurrent->SetBranchAddress("qCurrent.",&pqCurrent);
292   treeCurrent->SetBranchAddress("qRef.",&pqRef);
293   treeCurrent->SetBranchAddress("qCurrentInt.",&pqCurrentInt);
294   treeCurrent->SetBranchAddress("qRefInt.",&pqRefInt);
295   Int_t entries = treeCurrent->GetEntries();
296   //
297   //
298   //
299   Int_t nParamSide=62;               // number of z bins on 1 side
300   Int_t group=3;                     // grouping of bins
301   Int_t nParams=nParamSide/group;    // parameters
302   Int_t nParamsFit=nParams*nParams;  //
303   TLinearFitter fitter(nParamsFit+1,TString::Format("hyp%d",nParamsFit));
304   TVectorD xVector(nParamsFit+1);
305   TVectorD pVector(nParamsFit+1);
306   //
307   for (Int_t ievent=0; ievent<entries; ievent++){
308     treeCurrent->GetEntry(ievent);
309     for (Int_t iz=0; iz<nParamSide; iz++){      
310       Int_t dPar=iz/group;
311       if (dPar>nParams-1) dPar=nParams-1;
312       // 1.) clear X vectors
313       for (Int_t ipar=0; ipar<nParamsFit; ipar++) xVector[ipar]=0;
314       // 2.) set   X vector of interest
315       for (Int_t cPar=0; cPar<nParams; cPar++){
316         xVector[dPar*nParams+cPar] += (*pqCurrent)[cPar*group]/(*pqRef)[cPar*group];
317         xVector[dPar*nParams+cPar] += (*pqCurrent)[TMath::Min(cPar*group+1,nParamSide)]/(*pqRef)[TMath::Min(cPar*group+1,nParamSide)];
318         xVector[dPar*nParams+cPar] += (*pqCurrent)[TMath::Max(cPar*group-1,0)]/(*pqRef)[TMath::Max(cPar*group-1,0)];
319       }
320       Double_t val = 0;
321       val+=(*pnormZR)[dPar*group];
322       val+=(*pnormZR)[TMath::Max(dPar*group-1,0)];
323       val+=(*pnormZR)[TMath::Min(dPar*group,nParamSide)];
324       val/=3.;
325       fitter.AddPoint(xVector.GetMatrixArray(),val,0.0035);
326     }
327   }
328   for (Int_t cPar=1; cPar<nParams-1; cPar++){
329     //
330     for (Int_t ipar=0; ipar<nParamsFit; ipar++) xVector[ipar]=0;
331     xVector[(cPar-1)*nParams+cPar-1]=0.5;
332     xVector[(cPar)*nParams+cPar]    =-1;
333     xVector[(cPar+1)*nParams+cPar+1]=0.5;
334     fitter.AddPoint(xVector.GetMatrixArray(),0,0.05);
335   }
336   fitter.Eval();
337   //
338   //
339   TVectorD fitVector(nParamsFit);
340   TVectorD czVector(nParamsFit);
341   TVectorD fitVectorErr(nParamsFit);
342   fitter.GetParameters(fitVector);
343   for (Int_t iPar=0; iPar<nParamsFit; iPar++)  {
344     pVector[iPar]=iPar;
345     czVector[iPar]=250.*((iPar-1)%nParams)/nParams;
346     fitVectorErr[iPar]=fitter.GetParError(iPar);
347   }
348   Double_t chi2= TMath::Sqrt(fitter.GetChisquare()/fitter.GetNpoints());
349   printf("Chi2=%f\n",chi2);
350   TGraphErrors * gr  = new TGraphErrors(nParamsFit,  pVector.GetMatrixArray(), fitVector.GetMatrixArray(), 0,fitVectorErr.GetMatrixArray());
351   gr->SetMarkerStyle(25); gr->SetMarkerSize(0.3);
352   gr->Draw("alp");
353   
354   TGraphErrors * vgraphs[20]={0};
355   for (Int_t ipar=0; ipar<20; ipar++){
356     vgraphs[ipar]=new TGraphErrors(nParams, &(czVector.GetMatrixArray()[ipar*nParams+1]), &(fitVector.GetMatrixArray()[ipar*nParams+1]), 0,  &(fitVectorErr.GetMatrixArray()[ipar*nParams+1])); 
357     vgraphs[ipar]->GetXaxis()->SetTitle("z_{I} (cm)");
358     vgraphs[ipar]->GetYaxis()->SetTitle("#it{A}_{i_{I},j_{#DeltaR}}");
359     vgraphs[ipar]->SetMinimum(0);
360     vgraphs[ipar]->SetMaximum(0.1);
361     SetGraphTDRStyle(vgraphs[ipar]);
362   }
363   
364   TCanvas * canvasFit = new TCanvas("canvasFit","canvasFit",700,700); 
365   canvasFit->SetRightMargin(0.05);
366   canvasFit->SetLeftMargin(0.15);
367   canvasFit->SetBottomMargin(0.18);
368   canvasFit->Divide(1,2,0,0);
369   TLegend * legend0 = new TLegend(0.4,0.4,0.95,0.95,"Current fluctuation correction matrix. #DeltaR_{zi}=A_{ij}I_{zj}");
370   legend0->SetBorderSize(0);
371   TLegend * legend1 = new TLegend(0.4,0.5,0.95,0.95,"Current fluctuation correction matrix. #DeltaR_{zi}=A_{ij}I_{zj}");
372   legend1->SetBorderSize(0);
373
374   for (Int_t ipar=0; ipar<10; ipar++){
375     canvasFit->cd(ipar/5+1)->SetTicks(3,3);
376     vgraphs[ipar*2]->SetMarkerStyle(kMarkers[ipar%5]);
377     vgraphs[ipar*2]->SetMarkerColor(kColors[ipar%5]);
378     vgraphs[ipar*2]->SetLineColor(kColors[ipar%5]);
379     if (ipar%5==0) vgraphs[ipar*2]->Draw("alp");
380     vgraphs[ipar*2]->Draw("lp");
381     if (ipar<5) legend0->AddEntry(vgraphs[ipar*2],TString::Format("z_{#DeltaR}=%1.0f (cm)",250.*ipar/10.),"p");
382     if (ipar>=5) legend1->AddEntry(vgraphs[ipar*2],TString::Format("z_{#DeltaR}=%1.0f (cm)",250.*ipar/10.),"p");
383   }     
384   legend0->SetTextSize(0.05);
385   legend1->SetTextSize(0.05);
386   //
387   canvasFit->cd(1);
388   legend0->Draw();
389   canvasFit->cd(2);
390   legend1->Draw();
391   canvasFit->SaveAs("canvasAIJ_CurrenToDR.pdf");
392   canvasFit->SaveAs("canvasAIJ_CurrenToDR.png");
393 }
394
395 void MakeSmoothKernelStudy(Int_t npoints, Int_t nkernels, Int_t mapID){
396   //
397   // Compare the input and reconstructed distortion maps
398   // Input files are expected to have predefined names
399   // 
400   // Values and smoothed vaules using differnt smoothing algorithms are used
401   // Results of given studies will be used to define "optimal" smoothing algorithm
402   // and optimal Kernel parameters
403   // 
404   //
405   gRandom->SetSeed(mapID);
406   const Double_t cutMaxDist=3;
407   const Double_t cutMinDist=0.00001;
408   const Int_t kMinPoints=10;
409   TFile *foutput = TFile::Open("MeasureResidual.root");
410   TFile *finput = TFile::Open("RealResidualScaled.root");
411   AliTPCCorrectionLookupTable *mapIn = (AliTPCCorrectionLookupTable *)finput->Get("map");
412   AliTPCCorrectionLookupTable *mapOut = (AliTPCCorrectionLookupTable *)foutput->Get("map");
413   AliTPCCorrection::AddVisualCorrection(mapIn,1);    // input==1
414   AliTPCCorrection::AddVisualCorrection(mapOut,2);   // output (reconstructed==2)
415   TF1 * fin = new TF1("fin","AliTPCCorrection::GetCorrXYZ(x,5,10,1,1)",85,245);
416   TF1 * fout = new TF1("fout","AliTPCCorrection::GetCorrXYZ(x,5,10,1,2)",85,245);  
417   TTreeSRedirector * pcstream = new TTreeSRedirector("smoothLookupStudy.root","recreate");
418   //
419   //
420   // 1. Generate random point of interest
421   //
422   TVectorD  sigmaKernelR(nkernels);
423   TVectorD  sigmaKernelRPhi(nkernels);
424   TVectorD  sigmaKernelZ(nkernels);
425   //
426   TVectorD  fitValuesOut(nkernels);
427   TVectorD  fitValuesOutR(nkernels);
428   TVectorD  fitValuesIn(nkernels);
429   TVectorD  fitValuesInR(nkernels);
430   //
431   TVectorD  fitValuesOutErr(nkernels);
432   TVectorD  fitValuesOutChi2(nkernels);
433   TVectorD  fitSumW(nkernels);
434   TVectorD  fitValuesOutN(nkernels);
435   //
436   TLinearFitter *fittersOut[nkernels];
437   TLinearFitter *fittersOutR[nkernels];
438   TLinearFitter *fittersIn[nkernels];
439   TLinearFitter *fittersInR[nkernels];
440   for (Int_t ipoint=0; ipoint<npoints; ipoint++){ 
441     if (ipoint%10==0) printf("%d\n",ipoint);
442     for (Int_t i=0; i<nkernels; i++) {
443       fittersOut[i]=new TLinearFitter(7,"hyp6");
444       fittersOutR[i]=new TLinearFitter(7,"hyp6");
445       fittersIn[i]=new TLinearFitter(7,"hyp6");
446       fittersInR[i]=new TLinearFitter(7,"hyp6");
447     }
448     Double_t phi = gRandom->Rndm()*TMath::TwoPi();
449     Double_t r   = 85+gRandom->Rndm()*(245-85);
450     Double_t theta   = -0.9+gRandom->Rndm()*1.8;
451     Double_t z=r*theta;     
452     //
453     Double_t x = r*TMath::Cos(phi);
454     Double_t y = r*TMath::Sin(phi);
455     Double_t drphiInput = AliTPCCorrection::GetCorrXYZ(x,y,z,1,1);   // 
456     Double_t drphiOutput = AliTPCCorrection::GetCorrXYZ(x,y,z,1,2);  //
457     Double_t drInput = AliTPCCorrection::GetCorrXYZ(x,y,z,0,1);
458     Double_t drOutput = AliTPCCorrection::GetCorrXYZ(x,y,z,0,2);
459     if (TMath::Abs(drphiOutput)<cutMinDist) continue; // beter condition needed
460     if (TMath::Abs(drphiInput)<cutMinDist) continue; // beter condition needed
461     if (TMath::Abs(drphiOutput)>cutMaxDist) continue; // beter condition needed
462     if (TMath::Abs(drphiInput)>cutMaxDist) continue; // beter condition needed
463     //
464     for (Int_t i=0; i<nkernels; i++) {
465       sigmaKernelR[i]    = 1.+25.*gRandom->Rndm(); 
466       sigmaKernelRPhi[i] = 1.+25.*gRandom->Rndm();
467       sigmaKernelZ[i]    = 1.+25.*gRandom->Rndm();
468       fitSumW[i]=0;
469     }
470     for (Int_t idphi=-12; idphi<=12; idphi++){
471       Double_t dphi=idphi*TMath::Pi()/90.;            // we need to know actual segentation of the histograms - lookup should by multiple
472       for (Double_t dR=-50; dR<=50.; dR+=5){
473         for (Double_t dZ=-50; dZ<=50.; dZ+=5){
474           if (r+dR<85) continue;
475           if (r+dR>245) continue;
476           if (z+dZ<-240) continue;
477           if (z+dZ>240) continue;
478           if (z*(z+dZ)<0) continue;
479           //
480           Double_t x2=(r+dR)*TMath::Cos(phi+dphi);
481           Double_t y2=(r+dR)*TMath::Sin(phi+dphi);
482           Double_t z2=z+dZ;
483           Double_t drphiInput2=AliTPCCorrection::GetCorrXYZ(x2,y2,z2,1,1);
484           Double_t drInput2=AliTPCCorrection::GetCorrXYZ(x2,y2,z2,0,1);
485           Double_t drphiOutput2=AliTPCCorrection::GetCorrXYZ(x2,y2,z2,1,2);
486           Double_t drOutput2=AliTPCCorrection::GetCorrXYZ(x2,y2,z2,0,2);
487           if (TMath::Abs(drphiOutput2)<cutMinDist) continue; // hard cut beter condition needed
488           if (TMath::Abs(drphiOutput2)>cutMaxDist) continue; // beter condition needed
489           if (TMath::Abs(drphiInput2)<cutMinDist) continue; // hard cut beter condition needed
490           if (TMath::Abs(drphiInput2)>cutMaxDist) continue; // beter condition needed
491
492           Double_t xfit[7]={dphi,dphi*dphi,dR,dR*dR,dZ,dZ*dZ};
493           for (Int_t i=0; i<nkernels; i++) {
494             Double_t weight=1;
495             weight*=TMath::Gaus(dphi*r,0,sigmaKernelRPhi[i]);
496             weight*=TMath::Gaus(dR,0,sigmaKernelR[i]);
497             weight*=TMath::Gaus(dZ,0,sigmaKernelZ[i]);
498             weight+=0.00000001;
499             fitSumW[i]+=weight;
500             fittersOut[i]->AddPoint(xfit,drphiOutput2,0.1/weight);
501             fittersOutR[i]->AddPoint(xfit,drOutput2,0.1/weight);
502             //
503             fittersIn[i]->AddPoint(xfit,drphiInput2,0.1/weight);
504             fittersInR[i]->AddPoint(xfit,drInput2,0.1/weight);
505           }
506         }
507       }
508     }   
509     
510     for (Int_t ifix=0; ifix<=1; ifix++){
511       Bool_t isOK=kTRUE;
512       for (Int_t i=0; i<nkernels; i++) {
513         if( fittersOut[i]->GetNpoints() < kMinPoints) {
514           isOK=kFALSE;
515           break;
516         }
517         if (fitSumW[i]<0.01/*kMinWeight*/){
518           isOK=kFALSE;
519           break;
520         }
521         if (ifix==1){
522           fittersOut[i]->FixParameter(4,0);
523           fittersOut[i]->FixParameter(5,0);
524           fittersOut[i]->FixParameter(6,0);
525           fittersOutR[i]->FixParameter(4,0);
526           fittersOutR[i]->FixParameter(5,0);
527           fittersOutR[i]->FixParameter(6,0);
528           fittersIn[i]->FixParameter(4,0);
529           fittersIn[i]->FixParameter(5,0);
530           fittersIn[i]->FixParameter(6,0);
531           fittersInR[i]->FixParameter(4,0);
532           fittersInR[i]->FixParameter(5,0);
533           fittersInR[i]->FixParameter(6,0);
534         }
535         fittersOut[i]->Eval();
536         fittersOutR[i]->Eval();
537         fittersIn[i]->Eval();
538         fittersInR[i]->Eval();
539         //
540         fitValuesOut[i]=fittersOut[i]->GetParameter(0);
541         fitValuesOutR[i]=fittersOutR[i]->GetParameter(0);
542         fitValuesIn[i]=fittersIn[i]->GetParameter(0);
543         fitValuesInR[i]=fittersInR[i]->GetParameter(0);
544         //
545         fitValuesOutErr[i]=fittersOut[i]->GetParError(0);
546         fitValuesOutChi2[i]=TMath::Sqrt(fittersOut[i]->GetChisquare()/fitSumW[i]);
547         fitValuesOutN[i]=fittersOut[i]->GetNpoints();
548       }
549       if (isOK){
550       (*pcstream)<<"smoothLookup"<<
551         "ipoint"<<ipoint<<
552         "mapID="<<mapID<<
553         "ifix="<<ifix<<
554         // coordinates
555         "x="<<x<<
556         "y="<<y<<
557         "z="<<z<<
558         "phi="<<phi<<
559         "r="<<r<<
560         "z="<<z<<
561         // Input output values
562         "drphiInput="<<drphiInput<<       // input lookup tables disotrions
563         "drphiOutput="<<drphiOutput<<     // reconstructed lookup tables distortion
564         "drInput="<<drInput<<       // input lookup tables disotrions
565         "drOutput="<<drOutput<<     // reconstructed lookup tables distortion
566         // Smoothed values
567         "fitValuesIn.="<<&fitValuesIn<<            // smoothed input values - rphi direction
568         "fitValuesInR.="<<&fitValuesInR<<          // smoothed input values - r direction
569         "fitValuesOut.="<<&fitValuesOut<<          // smoothed measuured values - rphi
570         "fitValuesOutR.="<<&fitValuesOutR<<        // smoothed measuured values -r 
571         "fitValuesOutErr.="<<&fitValuesOutErr<<
572         "fitValuesOutChi2.="<<&fitValuesOutChi2<<
573         "fitValuesOutN.="<<&fitValuesOutN<<
574         "fitSumW.="<<&fitSumW<<
575         // Kernel sigma
576         "sigmaKernelR.="<<&sigmaKernelR<<
577         "sigmaKernelRPhi.="<<&sigmaKernelRPhi<<
578         "sigmaKernelZ.="<<&sigmaKernelZ<<
579         "\n";
580       }
581     }
582     for (Int_t i=0; i<nkernels; i++) {
583       delete fittersOut[i];
584       delete fittersOutR[i];
585       delete fittersIn[i];
586       delete fittersInR[i];
587     }
588   }
589   delete pcstream;
590   //
591   //
592   /* 
593      TFile *ff = TFile::Open("smoothLookupStudy.root")
594      TChain * chain = AliXRDPROOFtoolkit::MakeChain("smooth.list", "smoothLookup", 0,100)
595
596    */
597 }
598 void MakeSmoothKernelStudyDraw(Int_t npoints){
599   //
600   // make figure for kernel smoothing Draw
601   //
602   // npoints=20000
603   TFile *finput = TFile::Open("smoothLookupStudy.root");
604   TTree * chain = (TTree*)finput->Get("smoothLookup");  
605   chain->SetCacheSize(100000000);
606   TH2* hisInputsS[10]={0};
607   TH3* hisInputs3D[10]={0};
608   TH1* hisSigmaS[10]={0};
609   //
610   // 1.) Resolution not smoothed make nice plots
611   //
612   TH2 * his2DBase[10]={0};
613   TH1 * hisSigmas[10]={0};
614   chain->Draw("10*drphiInput:r>>hisIn(30,85,245,100,-2,2)","","colz");
615   his2DBase[0]=(TH2*)chain->GetHistogram()->Clone();
616   chain->Draw("10*(drphiOutput-drphiInput):r>>hisOutIn(30,85,245,100,-2,2)","","colz");
617   his2DBase[1]=(TH2*)chain->GetHistogram()->Clone();
618   his2DBase[0]->FitSlicesY(0,0,-1,0,"QNR",&garrayFit);
619   hisSigmas[0] = (TH1*)garrayFit.At(2)->Clone();
620   his2DBase[1]->FitSlicesY(0,0,-1,0,"QNR",&garrayFit);
621   hisSigmas[1] = (TH1*)garrayFit.At(2)->Clone();
622   //
623   TCanvas *canvasInOut = new TCanvas("deltaInOut","deltaInOut",600,500);
624   for (Int_t i=0; i<2; i++) {
625     hisSigmas[i]->SetMinimum(0); 
626     hisSigmas[i]->SetMaximum(1.5); 
627     hisSigmas[i]->SetMarkerStyle(kMarkers[i]);
628     hisSigmas[i]->SetMarkerColor(kColors[i]);
629     hisSigmas[i]->GetXaxis()->SetTitle("R (cm)");
630     hisSigmas[i]->GetYaxis()->SetTitle("#Delta_{R#phi} (mm)");
631     if (i==0) hisSigmas[i]->Draw("");
632     hisSigmas[i]->Draw("same");
633   }
634   TLegend * legend = new TLegend(0.4,0.5,0.89,0.89,"Residual #Delta_{r#phi}");
635   legend->SetBorderSize(0);
636   legend->AddEntry(hisSigmas[0],"#Delta_{current}-k#Delta_{mean}");
637   legend->AddEntry(hisSigmas[1],"(#Delta_{current}-k#Delta_{mean})_{align}-(#Delta_{current}-k#Delta_{mean})");
638   legend->Draw();
639   canvasInOut->SaveAs("canvasDistortionAlignmentInOut.pdf");
640   canvasInOut->SaveAs("canvasDistortionAlignmentInOut.png");
641   //
642   // 1.b)
643   //
644   TCanvas * canvasPhi= new TCanvas("canvasPhi","canvasPhi",800,600);
645   canvasPhi->Divide(1,3,0,0);
646   gStyle->SetOptTitle(1);
647   chain->SetAlias("sector","9*phi/pi");
648   {
649     chain->SetMarkerStyle(25);
650     chain->SetMarkerSize(0.3);
651     canvasPhi->cd(1);
652     chain->Draw("(drphiOutput-drphiInput):sector","abs((drphiOutput-drphiInput))<0.2&&abs(r-110)<30","",npoints);
653     canvasPhi->cd(2);
654     chain->Draw("(drphiOutput-drphiInput):sector","abs((drphiOutput-drphiInput))<0.2&&abs(r-160)<25","",npoints);
655     canvasPhi->cd(3);
656     chain->Draw("(drphiOutput-drphiInput):sector","abs((drphiOutput-drphiInput))<0.2&&abs(r-220)<30","",npoints);
657   }
658   TCanvas * canvasPhi2= new TCanvas("canvasPhi","canvasPhi",800,600);
659   canvasPhi2->Divide(1,3,0,0);
660   gStyle->SetOptTitle(1);
661   chain->SetAlias("sector","9*phi/pi");
662   {
663     chain->SetMarkerStyle(25);
664     chain->SetMarkerSize(0.3);
665     canvasPhi2->cd(1);
666     chain->Draw("(drphiOutput-drphiInput):sector-int(sector)","abs((drphiOutput-drphiInput))<0.2&&abs(r-110)<30","colz");
667     canvasPhi2->cd(2);
668     chain->Draw("(drphiOutput-drphiInput):sector-int(sector)","abs((drphiOutput-drphiInput))<0.2&&abs(r-160)<25","colz");
669     canvasPhi2->cd(3);
670     chain->Draw("(drphiOutput-drphiInput):sector-int(sector)","abs((drphiOutput-drphiInput))<0.2&&abs(r-220)<30","colz");
671   }
672
673
674   //
675   //
676   // 2.) Draw plot resolution as fucntion of the kernel smooth sigma
677   //      a.) Smoothed input- input
678   //      b.) Smoothed
679   TObjArray fitResols(1000);
680   const char* varSmooth[4]={"Smoothed(Input,Pol1)-Input","Smoothed(Output,Pol1)-Smoothed(Input,Pol1)","Smoothed(Output,Pol2)-Input","Smoothed(Output,Pol1)-Input"};
681   //
682   chain->Draw("10*(fitValuesIn.fElements-drphiInput):sqrt(sigmaKernelR.fElements**2+sigmaKernelRPhi.fElements**2+sigmaKernelZ.fElements**2):r>>hisSmoothDiff0(5,83,245, 15,5,40,100,-2.2,2.2)","ifix==1&&fitValuesOutErr.fElements<0.2","goff",npoints);
683   hisInputs3D[0]= (TH3*)chain->GetHistogram()->Clone();
684   chain->Draw("10*(fitValuesOut.fElements-fitValuesIn.fElements):sqrt(sigmaKernelR.fElements**2+sigmaKernelRPhi.fElements**2+sigmaKernelZ.fElements**2):r>>hisSmoothDiff1(5,83,245, 15,5,40,100,-2.2,2.2)","ifix==1&&fitValuesOutErr.fElements<0.2","goff",npoints);
685   hisInputs3D[1]= (TH3*)chain->GetHistogram()->Clone();
686   chain->Draw("10*(fitValuesOut.fElements-drphiInput):sqrt(sigmaKernelR.fElements**2+sigmaKernelRPhi.fElements**2+sigmaKernelZ.fElements**2):r>>hisSmoothDiff2(5,83,245, 15,5,40,100,-2.2,2.2)","ifix==0&&fitValuesOutErr.fElements<0.2","goff",npoints);
687   hisInputs3D[2]= (TH3*)chain->GetHistogram()->Clone();
688   chain->Draw("10*(fitValuesOut.fElements-drphiInput):sqrt(sigmaKernelR.fElements**2+sigmaKernelRPhi.fElements**2+sigmaKernelZ.fElements**2):r>>hisSmoothDiff3(5,83,245, 15,5,40,100,-2.2,2.2)","ifix==1&&fitValuesOutErr.fElements<0.2","goff",npoints);
689   hisInputs3D[3]= (TH3*)chain->GetHistogram()->Clone();
690   //
691   //
692   //
693   TCanvas *canvasSmooth = new TCanvas("DistortionAlignmentSmooth3D","DistortionAlignmentSmooth3D",800,800);
694   canvasSmooth->Divide(2,2,0,0);
695   gStyle->SetOptStat(0);  
696   for (Int_t itype=0; itype<4; itype++){
697     canvasSmooth->cd(itype+1);
698     TLegend * legend = new TLegend(0.2,0.7,0.9,0.99,varSmooth[itype]);
699     legend->SetBorderSize(0);
700     for (Int_t ibin=0; ibin<5; ibin++){
701       hisInputs3D[itype]->GetXaxis()->SetRange(ibin+1,ibin+1);
702       Double_t radius=hisInputs3D[itype]->GetXaxis()->GetBinCenter(ibin+1);
703       TH2 * his2D= (TH2*)hisInputs3D[itype]->Project3D("zy");
704       his2D->FitSlicesY(0,0,-1,0,"QNR",&garrayFit);
705       delete his2D;
706       TH1* his1D = (TH1*)garrayFit.At(2)->Clone();
707       fitResols.AddLast(his1D);
708       his1D->SetMarkerColor(kColors[ibin%6]);
709       his1D->SetMarkerStyle(kMarkers[ibin%6]);
710       his1D->SetMinimum(0);
711       his1D->SetMaximum(0.75);
712       his1D->GetXaxis()->SetTitle("#sigma_{kernel} (cm) in 3D (r,r#phi,z)");
713       his1D->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
714       if (ibin==0) his1D->Draw();
715       his1D->Draw("same");
716       legend->AddEntry(his1D,TString::Format("R=%2.0f (cm)",radius));
717     }
718     legend->Draw();
719   }
720   canvasSmooth->SaveAs("DistortionAlignmentSmooth3D.pdf");
721   canvasSmooth->SaveAs("DistortionAlignmentSmooth3D.pdf");
722
723
724
725 }
726
727 void FFTexample(){
728   //
729   TTreeSRedirector *pcstream = new TTreeSRedirector("localCurrent.root","update");
730   //
731   TFile *fCurrent = TFile::Open("SpaceChargeFluc10_1.root");
732   TFile *fRef = TFile::Open("SpaceChargeFluc0_1.root");
733   //  
734   AliTPCSpaceCharge3D* distortion = ( AliTPCSpaceCharge3D*)fCurrent->Get("DistRef"); 
735   AliTPCSpaceCharge3D* distortionRef = ( AliTPCSpaceCharge3D*)fRef->Get("DistRef"); 
736
737   distortion->AddVisualCorrection(distortion,1);
738   distortionRef->AddVisualCorrection(distortionRef,2);
739   //
740   TF2 f2("f2","AliTPCCorrection::GetCorrXYZ(x,0,y,0,1)-1.07*AliTPCCorrection::GetCorrXYZ(x,0,y,0,2)",85,245,0,250);
741   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);
742   f2.SetNpx(32); f2.SetNpy(32);
743   f2c.SetNpx(32); f2c.SetNpy(32);
744   f2.Draw("colz");
745   f2c.Draw("colz");
746   
747   TH2D * his2D = (TH2D*)f2.GetHistogram();
748   TH2D * his2Dc = (TH2D*)f2c.GetHistogram();
749   TH2D * his2DSmooth=new TH2D(*his2Dc);
750   Double_t sigma = 3;
751   Int_t nx=his2D->GetXaxis()->GetNbins();
752   Int_t ny=his2D->GetYaxis()->GetNbins();
753   //
754   for (Int_t ibin=2; ibin<=nx-1;ibin++){
755     for (Int_t jbin=2; jbin<ny-1;jbin++){
756       TLinearFitter fitter(4,"hyp3");
757       for (Int_t di=-3; di<=3; di++)
758         for (Int_t dj=-3; dj<=3; dj++){
759           if (ibin+di<=1) continue;
760           if (jbin+dj<=1) continue;
761           if (ibin+di>=nx) continue;
762           if (jbin+dj>=ny) continue;
763           Double_t x[3]={di,dj,dj*dj};
764           Double_t w = 1./(TMath::Gaus(di/sigma)*TMath::Gaus(dj/sigma));          
765           Double_t y = his2Dc->GetBinContent(ibin+di,jbin+dj);
766           fitter.AddPoint(x,y,w);
767         }
768       fitter.Eval();
769       printf("%d\t%d\t%3.3f\n",ibin,jbin,10*(fitter.GetParameter(0)-his2Dc->GetBinContent(ibin,jbin)));
770       his2DSmooth->SetBinContent(ibin,jbin,fitter.GetParameter(0));
771     } 
772   }
773   TH2D hisDiff=*his2DSmooth;
774   hisDiff.Add(his2D,-1);
775   
776
777   TH1 * hisMag = his2D->FFT(0,"MAG");
778   TH1 * hisMagc = his2Dc->FFT(0,"MAG");
779   TH1 * hisPhi = his2D->FFT(0,"PH");
780   TH1 * hisPhic = his2Dc->FFT(0,"PH");
781   hisMag->Draw("surf2");
782   hisMagc->Draw("surf2");
783
784     
785 }
786
787 /*
788 void MakeFFT(){
789   //
790   //
791   //
792   Int_t nSize = 125;
793   TVirtualFFT *fft_own = TVirtualFFT::FFT(1, &nSize, "R2C ES K");
794   fft_own->SetPoints(normZR->GetMatrixArray());
795   fft_own->Transform();
796
797   TH1 *hre = 0, *hco=0;  
798   hre = TH1::TransformHisto(fft_own, hre, "RE");
799   hco = TH1::TransformHisto(fft_own, hco, "C");
800
801
802   hr->SetTitle("Real part of the 3rd (array) tranfsorm");
803   hr->Draw();
804   hr->SetStats(kFALSE);
805   hr->GetXaxis()->SetLabelSize(0.05);
806   hr->GetYaxis()->SetLabelSize(0.05);
807   c1_6->cd();
808   TH1 *him = 0;
809   him = TH1::TransformHisto(fft_own, him, "IM");
810   him->SetTitle("Im. part of the 3rd (array) transform");
811   him->Draw();
812   him->SetStats(kFALSE);
813   him->GetXaxis()->SetLabelSize(0.05);
814   him->GetYaxis()->SetLabelSize(0.05);
815   //
816 }
817 */