TPC/TPCupgrade/macros/makeCurrentCorrection.C

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