TPC/Upgrade/macros/makeResidualSpaceChargeOCDB.C

   1 //
   2 // macro to create a residual OCDB
   3 //
   4
   5 /*
   6 //
   7 //
   8 //
   9 .x $ALICE_ROOT/TPC/Upgrade/macros/ConfigOCDB.C(1)
  10 .L $ALICE_ROOT/TPC/Upgrade/macros/makeResidualSpaceChargeOCDB.C
  11
  12 Example usage rescaling fit:
  13
  14
  15 aliroot -b -q $ALICE_ROOT/TPC/Upgrade/macros/makeResidualSpaceChargeOCDB.C\(0\)
  16
  17
  18 Example usage Jens alignment lookup tables:
  19
  20 ln -sf /hera/alice/wiechula/Upgrade/LUTs_fluctuation_eps20/MeasuredResidual2D/residualMap.1.root current.root
  21 ln -sf /hera/alice/wiechula/Upgrade/LUTs_fluctuation_eps20/ResidualResidual2D/residualMap.1.root mean.root
  22 ln -sf $ALICE_ROOT/OCDB/TPC/Calib/Correction/Run0_999999999_v0_s2.root ocdb.root
  23 aliroot -b -q $ALICE_ROOT/TPC/Upgrade/macros/makeResidualSpaceChargeOCDB.C\(1\)
  24
  25
  26
  27
  28 */
  29
  30 void makeResidualSpaceChargeOCDBLookup(const char *ocdbInName="ocdb.root",const char *scCurrentName="current.root", const char *scMeanName="mean.root"){
  31   //
  32   // Macro to create a clone original  OCDB entry with space point distortion and add there space point
  33   // distortions cause by resifual sapce charge
  34   // Output is stored by default in the current directory.
  35   //
  36
  37   // Parameters to Specify:
  38   //    ocdbInName      : path to the original OCDB entry
  39   //    scCurrentName   : path to the fluctuetaed distortion map
  40   //    scMeanName      : path to the mean distortion map
  41   //
  42
  43   /*
  44     scCurrentName="/hera/alice/wiechula/Upgrade/LUTs_fluctuation/dir1/SpaceChargeFluc10_1.lookup.root";
  45     scMeanName="/hera/alice/wiechula/Upgrade/LUTs_fluctuation/average/SpaceChargeFluc0_1.lookup.root";
  46     ocdbInName="$ALICE_ROOT/OCDB/TPC/Calib/Correction/Run0_999999999_v0_s2.root"
  47   */
  48   TFile * finCurrent = TFile::Open(scCurrentName);
  49   TFile * finMean = TFile::Open(scMeanName);
  50   //
  51   AliTPCCorrectionLookupTable *spaceChargeCurrent=  (AliTPCCorrectionLookupTable *)finCurrent->Get("map");
  52   AliTPCCorrectionLookupTable *spaceChargeMean   =  (AliTPCCorrectionLookupTable *)finMean->Get("map");
  53   AliTPCInverseCorrection * spaceChargeInverseMean = new AliTPCInverseCorrection(spaceChargeMean);
  54   //
  55   TObjArray * arraySC = new TObjArray(2);
  56   arraySC->AddAt(spaceChargeCurrent,0);
  57   arraySC->AddAt(spaceChargeMean,1);
  58   AliTPCComposedCorrection *corrSC = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
  59   AliTPCComposedCorrection *corrSCW = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
  60   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeCurrent,1);
  61   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeMean,2);
  62   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeInverseMean,3);
  63   AliTPCComposedCorrection::AddVisualCorrection(corrSC,4);
  64   AliTPCComposedCorrection::AddVisualCorrection(corrSCW,5);
  65   //
  66   //
  67   TLinearFitter fitterR(2,"pol1");
  68   for (Int_t ipoint =0; ipoint<10000; ipoint++){
  69     Double_t z0 = -250+gRandom->Rndm()*500;
  70     Double_t r0   = 85+gRandom->Rndm()*(245-85.);
  71     Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();
  72     // some problematic parts to be skipped  - investigated later
  73     Double_t xvalue[2]={0};
  74     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1))>50) continue;
  75     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2))>50) continue;
  76     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1))>20) continue;
  77     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2))>20) continue;
  78     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1))>50) continue;
  79     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2))>50) continue;
  80     //
  81     //
  82     xvalue[0]=AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2);
  83     fitterR.AddPoint(xvalue,AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1));
  84   }
  85   fitterR->Eval();
  86   TVectorD weights(2);
  87   weights[0]=1;
  88   weights[1]=-1;
  89   corrSC->SetWeights( &weights);
  90   weights[1]=-fitterR->GetParameter(1);
  91   corrSCW->SetWeights( &weights);
  92   //
  93   //
  94   //
  95   TF1 *fSC = new TF1("fSC","AliTPCCorrection::GetCorrXYZ(x,0,10,0,4)",85,245);
  96   TF1 *fSCW = new TF1("fSCW","AliTPCCorrection::GetCorrXYZ(x,0,10,0,5)",85,245);
  97   fSC->SetLineColor(2);
  98   fSC->SetLineColor(4);
  99   fSC->Draw();
 100   fSCW->Draw("same");
 101   gPad->SaveAs("residualMap.pdf");
 102
 103   //
 104   //
 105   //
 106   TFile *fileOCDBIn=TFile::Open(ocdbInName);
 107   AliCDBEntry   *entry =  ( AliCDBEntry   *)fileOCDBIn->Get("AliCDBEntry");
 108   TObjArray * array = (TObjArray*)entry->GetObject();
 109   for (Int_t i=0;  i<3;  i++){
 110     AliTPCComposedCorrection *corr = ( AliTPCComposedCorrection*)array->At(i);
 111     if (corr){
 112       TObjArray* corrArray = corr->GetCorrections();
 113       corrArray->AddLast(corrSCW);
 114     }
 115   }
 116   AliCDBMetaData *metaData= new AliCDBMetaData();
 117   metaData->SetObjectClassName("TObjArray");
 118   metaData->SetResponsible("Marian Ivanov");
 119   metaData->SetBeamPeriod(1);
 120   metaData->SetAliRootVersion("05-25-01"); //root version
 121   metaData->SetComment("Standard+fluctuation");
 122   AliCDBId* id1=NULL;
 123   id1=new AliCDBId("TPC/Calib/Correction", 0, 9999999999999999);
 124   //
 125   TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
 126   pocdbStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
 127   pocdbStorage->Put(array, (*id1), metaData);
 128
 129
 130 }
 131
 132
 133
 134
 135 void makeResidualSpaceChargeOCDBFit(const char *ocdbInName="ocdb.root",const char *scCurrentName="current.root", const char *scMeanName="mean.root"){
 136   //
 137   // Macro to create a clone original  OCDB entry with space point distortion and add there space point
 138   // distortions cause by resifual sapce charge
 139   // Output is stored by default in the current directory.
 140   //
 141
 142   // Parameters to Specify:
 143   //    ocdbInName      : path to the original OCDB entry
 144   //    scCurrentName   : path to the fluctuetaed distortion map
 145   //    scMeanName      : path to the mean distortion map
 146   //
 147
 148   /*
 149     scCurrentName="/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/dir0/SpaceChargeFluc10_1.root";
 150     scMeanName="/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/dir0//SpaceChargeFluc0_1.root";
 151     ocdbInName="$ALICE_ROOT/OCDB/TPC/Calib/Correction/Run0_999999999_v0_s2.root"
 152   */
 153   TFile * finCurrent = TFile::Open(scCurrentName);
 154   TFile * finMean = TFile::Open(scMeanName);
 155   //
 156   AliTPCCorrectionLookupTable *spaceChargeCurrent=  (AliTPCCorrectionLookupTable *)finCurrent->Get("DistRef");
 157   AliTPCCorrectionLookupTable *spaceChargeMean   =  (AliTPCCorrectionLookupTable *)finMean->Get("DistRef");
 158   AliTPCInverseCorrection * spaceChargeInverseMean = new AliTPCInverseCorrection(spaceChargeMean);
 159   //
 160   TObjArray * arraySC = new TObjArray(2);
 161   arraySC->AddAt(spaceChargeCurrent,0);
 162   arraySC->AddAt(spaceChargeMean,1);
 163   AliTPCComposedCorrection *corrSC = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
 164   AliTPCComposedCorrection *corrSCW = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
 165   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeCurrent,1);
 166   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeMean,2);
 167   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeInverseMean,3);
 168   AliTPCComposedCorrection::AddVisualCorrection(corrSC,4);
 169   AliTPCComposedCorrection::AddVisualCorrection(corrSCW,5);
 170   //
 171   //
 172   TLinearFitter fitterR(2,"pol1");
 173   for (Int_t ipoint =0; ipoint<10000; ipoint++){
 174     Double_t z0 = -250+gRandom->Rndm()*500;
 175     Double_t r0   = 85+gRandom->Rndm()*(245-85.);
 176     Double_t phi0 = gRandom->Rndm()*TMath::TwoPi();
 177     // some problematic parts to be skipped  - investigated later
 178     Double_t xvalue[2]={0};
 179     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,1))>50) continue;
 180     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,0,2))>50) continue;
 181     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1))>20) continue;
 182     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2))>20) continue;
 183     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,1))>50) continue;
 184     if (TMath::Abs(AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,2,2))>50) continue;
 185     //
 186     //
 187     xvalue[0]=AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,2);
 188     fitterR.AddPoint(xvalue,AliTPCCorrection::GetCorrXYZ(r0*TMath::Cos(phi0), r0*TMath::Sin(phi0), z0,1,1));
 189   }
 190   fitterR->Eval();
 191   TVectorD weights(2);
 192   weights[0]=1;
 193   weights[1]=-1;
 194   corrSC->SetWeights( &weights);
 195   weights[1]=-fitterR->GetParameter(1);
 196   corrSCW->SetWeights( &weights);
 197   //
 198   //
 199   //
 200   TF1 *fSC = new TF1("fSC","AliTPCCorrection::GetCorrXYZ(x,0,10,1,4)",85,245);
 201   TF1 *fSCW = new TF1("fSCW","AliTPCCorrection::GetCorrXYZ(x,0,10,1,5)",85,245);
 202   fSC->SetLineColor(2);
 203   fSC->SetLineColor(4);
 204   fSC->GetXaxis()->SetTitle(" R (cm)");
 205   fSC->GetYaxis()->SetTitle(" #Delta_{R#phi} (cm)");
 206   TCanvas * canvasDist = new TCanvas("canvasDist","canvasDist",600,500);
 207   TF1 *f1sc = new TF1("f1sc","[0]+[1]*x+[2]*x**2",85,245);
 208   TF1 *f1scw = new TF1("f1scw","[0]+[1]*x+[2]*x**2",85,245);
 209
 210   {
 211     fSC->Draw();
 212     fSCW->Draw();
 213     fSC->SetMinimum(-1);
 214     fSC->SetMaximum(2.);
 215     fSC->GetHistogram()->SetLineColor(2);
 216     fSCW->GetHistogram()->SetLineColor(4);
 217     //
 218     f1sc->SetLineColor(2);
 219     f1scw->SetLineColor(4);
 220     f1sc->SetLineWidth(0.5);
 221     f1scw->SetLineWidth(0.5);
 222     //
 223     fSC->GetHistogram()->Fit(f1sc);
 224     fSCW->GetHistogram()->Fit(f1scw);
 225     //
 226     fSC->GetHistogram()->Draw("");
 227     fSCW->GetHistogram()->Draw("same");
 228     //
 229     TF1 *f1scp=new TF1(*f1sc);
 230     TF1 *f1scm=new TF1(*f1sc);
 231     f1scp->SetLineStyle(2); f1scm->SetLineStyle(2);
 232     f1scp->SetParameter(0, f1sc->GetParameter(0)+ 0.3);
 233     f1scm->SetParameter(0, f1sc->GetParameter(0)- 0.3);
 234     f1scp->Draw("same");
 235     f1scm->Draw("same");
 236     //
 237     TF1 *f1scwp=new TF1(*f1scw);
 238     TF1 *f1scwm=new TF1(*f1scw);
 239     f1scwp->SetLineStyle(2); f1scwm->SetLineStyle(2);
 240     f1scwp->SetParameter(0, f1scw->GetParameter(0)+ 0.3);
 241     f1scwm->SetParameter(0, f1scw->GetParameter(0)- 0.3);
 242     f1scwp->Draw("same");
 243     f1scwm->Draw("same");
 244   }
 245   TLegend * legend = new TLegend(0.3,0.5,0.89,0.89,"Residual #Delta_{R#phi} distortions and helix fit");
 246   legend->SetBorderSize(0);
 247   legend->AddEntry(fSC->GetHistogram(),"Stage 0: #Delta_{R#phic}-#Delta_{R#phimean}");
 248   legend->AddEntry(fSCW->GetHistogram(),"Stage 1: #Delta_{R#phic}-k_{fit}#Delta_{R#phimean}");
 249   legend->Draw();
 250   gPad->SaveAs("residualMapTrackFit.pdf");
 251
 252   //
 253   //
 254   //
 255   TFile *fileOCDBIn=TFile::Open(ocdbInName);
 256   AliCDBEntry   *entry =  ( AliCDBEntry   *)fileOCDBIn->Get("AliCDBEntry");
 257   TObjArray * array = (TObjArray*)entry->GetObject();
 258   for (Int_t i=0;  i<3;  i++){
 259     AliTPCComposedCorrection *corr = ( AliTPCComposedCorrection*)array->At(i);
 260     if (corr){
 261       TObjArray* corrArray = corr->GetCorrections();
 262       corrArray->AddLast(corrSCW);
 263     }
 264   }
 265   AliCDBMetaData *metaData= new AliCDBMetaData();
 266   metaData->SetObjectClassName("TObjArray");
 267   metaData->SetResponsible("Marian Ivanov");
 268   metaData->SetBeamPeriod(1);
 269   metaData->SetAliRootVersion("05-25-01"); //root version
 270   metaData->SetComment("Standard+fluctuation");
 271   AliCDBId* id1=NULL;
 272   id1=new AliCDBId("TPC/Calib/Correction", 0, 9999999999999999);
 273   //
 274   TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
 275   pocdbStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
 276   pocdbStorage->Put(array, (*id1), metaData);
 277
 278
 279 }
 280
 281
 282
 283
 284 void makeResidualSpaceChargeOCDBAlign(const char *ocdbInName="ocdb.root",const char *scCurrentName="current.root", const char *scMeanName="mean.root"){
 285   //
 286   // Macro to create a clone original  OCDB entry with space point distortion and add there space point
 287   // distortion due residual space charge isalignemnt
 288   // Output is stored by default in the current directory.
 289   //
 290
 291   // Parameters to Specify:
 292   //    ocdbInName      : path to the original OCDB entry
 293   //    scCurrentName   : path to the fluctuetaed distortion map
 294   //    scMeanName      : path to the mean distortion map
 295   //
 296
 297   /*
 298     scCurrentName="/hera/alice/wiechula/Upgrade/LUTs_fluctuation_eps20/MeasuredResidual2D/residualMap.10.root";
 299     scMeanName="/hera/alice/wiechula/Upgrade/LUTs_fluctuation_eps20/ResidualResidual2D/residualMap.10.root";
 300     ocdbInName="$ALICE_ROOT/OCDB/TPC/Calib/Correction/Run0_999999999_v0_s2.root"
 301   */
 302   TFile * finCurrent = TFile::Open(scCurrentName);
 303   TFile * finMean = TFile::Open(scMeanName);
 304   //
 305   AliTPCCorrectionLookupTable *spaceChargeCurrent=  (AliTPCCorrectionLookupTable *)finCurrent->Get("map");
 306   AliTPCCorrectionLookupTable *spaceChargeMean   =  (AliTPCCorrectionLookupTable *)finMean->Get("map");
 307   //
 308   TObjArray * arraySC = new TObjArray(2);
 309   arraySC->AddAt(spaceChargeCurrent,0);
 310   arraySC->AddAt(spaceChargeMean,1);
 311   AliTPCComposedCorrection *corrSC = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
 312   AliTPCComposedCorrection *corrSCW = new AliTPCComposedCorrection(arraySC,AliTPCComposedCorrection::kParallel);
 313   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeCurrent,1);        // reconstructed
 314   AliTPCComposedCorrection::AddVisualCorrection(spaceChargeMean,2);           // ideal
 315   AliTPCComposedCorrection::AddVisualCorrection(corrSCW,5);
 316   //
 317   TVectorD weights(2);
 318   weights[0]=1;
 319   weights[1]=1;
 320   corrSCW->SetWeights( &weights);
 321   TTreeSRedirector * pcstream = new TTreeSRedirector("residualDelta.root","recreate");
 322   //
 323   for (Int_t ipoint=0; ipoint<50000; ipoint++){
 324     Double_t phi = gRandom->Rndm()*TMath::TwoPi();
 325     Double_t r   = 85+gRandom->Rndm()*(245-85);
 326     Double_t theta   = -0.9+gRandom->Rndm()*1.8;
 327     Double_t z=r*theta;
 328     Double_t x = r*TMath::Cos(phi);
 329     Double_t y = r*TMath::Sin(phi);
 330     Double_t drInput = AliTPCCorrection::GetCorrXYZ(x,y,z,0,1);
 331     Double_t drOutput = AliTPCCorrection::GetCorrXYZ(x,y,z,0,2);
 332     Double_t drDiff = AliTPCCorrection::GetCorrXYZ(x,y,z,0,5);
 333     Double_t drphiInput = AliTPCCorrection::GetCorrXYZ(x,y,z,1,1);   //
 334     Double_t drphiOutput = AliTPCCorrection::GetCorrXYZ(x,y,z,1,2);  //
 335     Double_t drphiDiff = AliTPCCorrection::GetCorrXYZ(x,y,z,1,5);  //
 336     //
 337     (*pcstream)<<"diff"<<
 338       "phi="<<phi<<
 339       "r="<<r<<
 340       "x="<<x<<
 341       "y="<<y<<
 342       "z="<<z<<
 343       //
 344       "drInput="<<drInput<<
 345       "drOutput="<<drOutput<<
 346       "drDiff="<<drDiff<<
 347       //
 348       "drphiInput="<<drphiInput<<
 349       "drphiOutput="<<drphiOutput<<
 350       "drphiDiff="<<drphiDiff<<
 351       "\n";
 352   }
 353   TTree * tree = ((*pcstream)<<"diff")->GetTree();
 354   //
 355   TCanvas *canvas = new TCanvas("canvasDiff","canvasDiff",900,800);
 356   canvas->Divide(3,2);
 357   canvas->cd(1);
 358   tree->Draw("drphiOutput:drphiInput","abs(drphiDiff)+abs(drphiOutput)+abs(drphiInput)<3","");
 359   canvas->cd(2);
 360   tree->Draw("drphiDiff:r","abs(drphiDiff)+abs(drphiOutput)+abs(drphiInput)<3","colz");
 361   canvas->cd(3);
 362   tree->Draw("drphiDiff","abs(drphiDiff)+abs(drphiOutput)+abs(drphiInput)<3","");
 363   canvas->cd(4);
 364   tree->Draw("drOutput:drInput","abs(drDiff)+abs(drOutput)+abs(drInput)<3","");
 365   canvas->cd(5);
 366   tree->Draw("drDiff:r","abs(drDiff)+abs(drOutput)+abs(drInput)<3","colz");
 367   canvas->cd(6);
 368   tree->Draw("drDiff","abs(drDiff)+abs(drOutput)+abs(drInput)<3","");
 369   canvas->SaveAs("residualSpaceChargeDistortion.pdf");
 370   //
 371   //
 372   //
 373   TFile *fileOCDBIn=TFile::Open(ocdbInName);
 374   AliCDBEntry   *entry =  ( AliCDBEntry   *)fileOCDBIn->Get("AliCDBEntry");
 375   TObjArray * array = (TObjArray*)entry->GetObject();
 376   for (Int_t i=0;  i<3;  i++){
 377     AliTPCComposedCorrection *corr = ( AliTPCComposedCorrection*)array->At(i);
 378     if (corr){
 379       TObjArray* corrArray = corr->GetCorrections();
 380       corrArray->AddLast(corrSCW);
 381     }
 382   }
 383   AliCDBMetaData *metaData= new AliCDBMetaData();
 384   metaData->SetObjectClassName("TObjArray");
 385   metaData->SetResponsible("Marian Ivanov");
 386   metaData->SetBeamPeriod(1);
 387   metaData->SetAliRootVersion("05-25-01"); //root version
 388   metaData->SetComment("Standard+fluctuation");
 389   AliCDBId* id1=NULL;
 390   id1=new AliCDBId("TPC/Calib/Correction", 0, 9999999999999999);
 391   //
 392   TString localStorage = "local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
 393   pocdbStorage = AliCDBManager::Instance()->GetStorage(localStorage.Data());
 394   pocdbStorage->Put(array, (*id1), metaData);
 395
 396
 397 }
 398
 399
 400
 401 void makeResidualSpaceChargeOCDB(Int_t action=1){
 402   //
 403   //
 404   //
 405   if (action==0) makeResidualSpaceChargeOCDBFit();  // fit mean charge to current charge
 406   if (action==1) makeResidualSpaceChargeOCDBAlign();  // use difference of the map and the reconstructed distorion map as a residual OCDB
 407 }
 408