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 */