TPC/AliTPCCorrectionFit.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17
  18 /*
  19   Responsible: marian.ivanov@cern.ch
  20   Tools for fitting of the space point distortion parameters.
  21   Functionality
  22
  23
  24 1. Creation of the distortion maps from the residual histograms
  25
  26 2. Making fit trees
  27
  28 3. Space point distortion not directly observable. Instead a derived variables
  29    like DCA at vertex, local y distortion in the TPC-*TOF,TRD,ITS) matching
  30    in all 5 tracking parameters are obsereved.
  31    In the AliTPCcorrection fir code we calculate the derivative of given variables
  32    dO_{i}/dp_{i}
  33
  34 4. Global fit - later
  35    d0 = sum{ki*dO_{i}/dp_{i}}  - linear fitting of the amplitudes ki
  36
  37  Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
  38  extracted from the old macros
  39
  40
  41 */
  42
  43 #include "Riostream.h"
  44 #include <fstream>
  45 #include "TMap.h"
  46 #include "TGraphErrors.h"
  47 #include "AliExternalTrackParam.h"
  48 #include "TFile.h"
  49 #include "TGraph.h"
  50 #include "TMultiGraph.h"
  51 #include "TCanvas.h"
  52 #include "THnSparse.h"
  53 #include "THnBase.h"
  54 #include "TProfile.h"
  55 #include "TROOT.h"
  56 #include "TLegend.h"
  57 #include "TPad.h"
  58 #include "TH2D.h"
  59 #include "TH3D.h"
  60 #include "AliTPCROC.h"
  61 #include "AliTPCCalROC.h"
  62 #include "AliESDfriend.h"
  63 #include "AliTPCcalibTime.h"
  64 #include "AliSplineFit.h"
  65 #include "AliCDBMetaData.h"
  66 #include "AliCDBId.h"
  67 #include "AliCDBManager.h"
  68 #include "AliCDBStorage.h"
  69 #include "AliTPCcalibBase.h"
  70 #include "AliTPCcalibDB.h"
  71 #include "AliTPCcalibDButil.h"
  72 #include "AliRelAlignerKalman.h"
  73 #include "AliTPCParamSR.h"
  74 #include "AliTPCcalibTimeGain.h"
  75 #include "AliTPCcalibGainMult.h"
  76 #include "AliSplineFit.h"
  77 #include "AliTPCComposedCorrection.h"
  78 #include "AliTPCExBTwist.h"
  79 #include "AliTPCCalibGlobalMisalignment.h"
  80 #include "TStatToolkit.h"
  81 #include "TChain.h"
  82 #include "TCut.h"
  83 #include "AliTrackerBase.h"
  84 #include "AliTPCCorrectionFit.h"
  85 #include "AliTPCLaserTrack.h"
  86 #include "TDatabasePDG.h"
  87 #include "AliTPCcalibAlign.h"
  88 #include "AliLog.h"
  89 #include "AliRieman.h"
  90
  91 ClassImp(AliTPCCorrectionFit)
  92
  93 AliTPCCorrectionFit::AliTPCCorrectionFit():
  94   TNamed("TPCCorrectionFit","TPCCorrectionFit")
  95 {
  96   //
  97   // default constructor
  98   //
  99 }
 100
 101 AliTPCCorrectionFit::~AliTPCCorrectionFit() {
 102   //
 103   // Destructor
 104   //
 105 }
 106
 107
 108 Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){
 109   //
 110   // Evalution at point using the lienar approximation
 111   //
 112   Double_t sector = 9*phi/TMath::Pi();
 113   if (sector<0) sector+=18;
 114   Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
 115   Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
 116   if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
 117   if (ptype==2) return (y245-y85)/(245.-85.);
 118   return 0;
 119 }
 120
 121
 122
 123
 124 Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
 125   //
 126   // Fit the distortion along the line with the parabolic model
 127   // We assume that the track are primaries  - where the vertex is at (0,0,0)
 128   //
 129   // Parameters:
 130   //  phi0  -  phi at the entrance of the TPC
 131   //  snp   -  local inclination angle at the entrance of the TPC
 132   //  refX  -  ref X where the distortion is evanluated
 133   //  theta -  inclination angle
 134   //  corr  -  internal number of the correction and distortion
 135   //  ptype -  0 - local y distortion
 136   //           //1 - local z distortion
 137   //           2 - local derivative distortion
 138   //           //3
 139   //           4 - distortion in the curvature
 140   //  nsteps - number of fit points
 141   //
 142   // return value -  distortion at point refX with type ptype
 143   //
 144   static TLinearFitter fitter(3,"pol2");
 145   fitter.ClearPoints();
 146   if (nsteps<3) nsteps=3;
 147   Double_t deltaX=(245-85)/(nsteps);
 148   for (Int_t istep=0; istep<(nsteps+1); istep++){
 149     //
 150     Double_t localX =85.+deltaX*istep;
 151     Double_t localPhi=phi0+deltaX*snp*istep;
 152     Double_t sector = 9*localPhi/TMath::Pi();
 153     if (sector<0) sector+=18;
 154     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
 155     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
 156     Double_t x[1]={localX-dlocalX};
 157     fitter.AddPoint(x,dy);
 158   }
 159   fitter.Eval();
 160   Double_t par[3];
 161   par[0]=fitter.GetParameter(0);
 162   par[1]=fitter.GetParameter(1);
 163   par[2]=fitter.GetParameter(2);
 164
 165   if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
 166   if (ptype==2) return par[1]+2*par[2]*refX;
 167   if (ptype==4) return par[2];
 168   return 0;
 169 }
 170
 171
 172 Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
 173   //
 174   // Fit the distortion along the line with the helix model
 175   // FIXME - original trajectory to be changed - AliHelix to be used
 176   // We assume that the track are primaries  - where the vertex is at (0,0,0)
 177   //
 178   // Parameters:
 179   //  phi0  -  phi at the entrance of the TPC
 180   //  snp   -  local inclination angle at the entrance of the TPC
 181   //  refX  -  ref X where the distortion is evanluated
 182   //  theta -  inclination angle
 183   //  corr  -  internal number of the correction and distortion
 184   //  ptype -  0 - local y distortion
 185   //           //1 - local z distortion
 186   //           2 - local derivative distortion
 187   //           //3
 188   //           4 - distortion in the curvature
 189   //  nsteps - number of fit points
 190   //
 191   // return value -  distortion at point refX with type ptype
 192   //
 193   if (nsteps<3) nsteps=3;
 194   Double_t deltaX=(245-85)/(nsteps);
 195   AliRieman rieman(nsteps);
 196
 197   for (Int_t istep=0; istep<(nsteps+1); istep++){
 198     //
 199     Double_t localX =85.+deltaX*istep;
 200     Double_t localPhi=phi0+deltaX*snp*istep;
 201     Double_t sector = 9*localPhi/TMath::Pi();
 202     if (sector<0) sector+=18;
 203     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
 204     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
 205     Double_t x[1]={localX-dlocalX};
 206     Double_t z=theta*x[0];
 207     rieman.AddPoint(x[0],dy,z,0.1,0.1);
 208   }
 209   rieman.Update();
 210   //
 211
 212   if (ptype==0) return rieman.GetYat(refX);
 213   if (ptype==2) return rieman.GetDYat(refX);
 214   if (ptype==4) return rieman.GetC();
 215   return 0;
 216 }
 217
 218
 219
 220
 221 void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
 222   //
 223   // Create cluster distortion map
 224   //
 225   TFile *falign = TFile::Open("CalibObjects.root");
 226   TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
 227   if (!arrayAlign) {
 228     AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
 229     return;
 230   }
 231   AliTPCcalibAlign * align =  (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
 232   if (!align) {
 233       AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
 234     return;
 235   }
 236   TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
 237
 238   THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
 239   //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
 240   AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
 241   //  AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
 242
 243   const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
 244   for (Int_t ihis=0; ihis<4; ihis++){
 245     THnSparse * hisAlign =align->GetTrackletDelta(ihis);
 246     AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
 247   }
 248   delete pcstream;
 249   delete falign;
 250 }
 251
 252
 253
 254
 255
 256 void  AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
 257   //
 258   // Make cluster residual map from the n-dimensional histogram
 259   // hisInput supposed to have given format:
 260   //     4 Dim:
 261   //      delta,
 262   //      sector
 263   //      localX
 264   //      kZ
 265   // Vertex position assumed to be at (0,0,0)
 266   //
 267   //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
 268   //
 269   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
 270   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
 271   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
 272   TF1 *fgaus=0;
 273   TH3F * hisResMap3D =
 274     new TH3F("his3D","his3D",
 275              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 276              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
 277              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
 278   hisResMap3D->GetXaxis()->SetTitle("sector");
 279   hisResMap3D->GetYaxis()->SetTitle("localX");
 280   hisResMap3D->GetZaxis()->SetTitle("kZ");
 281
 282   TH2F * hisResMap2D[4] ={0,0,0,0};
 283   for (Int_t i=0; i<4; i++){
 284     hisResMap2D[i]=
 285       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
 286                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 287                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
 288     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
 289     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
 290   }
 291   //
 292   //
 293   //
 294   TF1 * f1= 0;
 295   Int_t axis0[4]={0,1,2,3};
 296   Int_t axis1[4]={0,1,2,3};
 297   Int_t counter=0;
 298   for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
 299     // phi- sector  range
 300     hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
 301     THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
 302     Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
 303     //
 304     for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
 305       // local x range
 306       // kz fits
 307       his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
 308       THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
 309       Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
 310       //
 311       //A side
 312       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 313       TH1 * hisA = his2->Projection(0);
 314       Double_t meanA= hisA->GetMean();
 315       Double_t rmsA= hisA->GetRMS();
 316       Double_t entriesA= hisA->GetEntries();
 317       delete hisA;
 318       //C side
 319       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 320       TH1 * hisC = his2->Projection(0);
 321       Double_t meanC= hisC->GetMean();
 322       Double_t rmsC= hisC->GetRMS();
 323       Double_t entriesC= hisC->GetEntries();
 324       delete hisC;
 325       his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
 326       TH2 * hisAC       = his2->Projection(0,3);
 327       TProfile *profAC  = hisAC->ProfileX();
 328       delete hisAC;
 329       profAC->Fit("pol1","QNR","QNR",0.05,1);
 330       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
 331       Double_t offsetA=f1->GetParameter(0);
 332       Double_t slopeA=f1->GetParameter(1);
 333       Double_t offsetAE=f1->GetParError(0);
 334       Double_t slopeAE=f1->GetParError(1);
 335       Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
 336       profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
 337       f1=(TF1*)gROOT->FindObject("pol1");
 338       Double_t offsetC=f1->GetParameter(0);
 339       Double_t slopeC=f1->GetParameter(1);
 340       Double_t offsetCE=f1->GetParError(0);
 341       Double_t slopeCE=f1->GetParError(1);
 342       Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
 343       if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
 344       counter++;
 345       (*pcstream)<<"deltaFit"<<
 346         "sector="<<sector<<
 347         "localX="<<localX<<
 348         "meanA="<<meanA<<
 349         "rmsA="<<rmsA<<
 350         "entriesA="<<entriesA<<
 351         "meanC="<<meanC<<
 352         "rmsC="<<rmsC<<
 353         "entriesC="<<entriesC<<
 354         "offsetA="<<offsetA<<
 355         "slopeA="<<slopeA<<
 356         "offsetAE="<<offsetAE<<
 357         "slopeAE="<<slopeAE<<
 358         "chi2A="<<chi2A<<
 359         "offsetC="<<offsetC<<
 360         "slopeC="<<slopeC<<
 361         "offsetCE="<<offsetCE<<
 362         "slopeCE="<<slopeCE<<
 363         "chi2C="<<chi2C<<
 364         "\n";
 365       //
 366       hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
 367       hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
 368       hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
 369       hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
 370
 371       for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
 372         Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
 373         if (TMath::Abs(kZ)<0.05) continue;  // crossing
 374         his2->GetAxis(3)->SetRange(ibin3,ibin3);
 375         if (TMath::Abs(kZ)>0.15){
 376           his2->GetAxis(3)->SetRange(ibin3,ibin3);
 377         }
 378         TH1 * his = his2->Projection(0);
 379         Double_t mean= his->GetMean();
 380         Double_t rms= his->GetRMS();
 381         Double_t entries= his->GetEntries();
 382         //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
 383         hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
 384         Double_t phi=TMath::Pi()*sector/9;
 385         if (phi>TMath::Pi()) phi+=TMath::Pi();
 386         Double_t meanG=0;
 387         Double_t rmsG=0;
 388         if (entries>50){
 389           if (!fgaus) {
 390             his->Fit("gaus","Q","goff");
 391             fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
 392           }
 393           if (fgaus) {
 394             his->Fit(fgaus,"Q","goff");
 395             meanG=fgaus->GetParameter(1);
 396             rmsG=fgaus->GetParameter(2);
 397           }
 398         }
 399         Double_t dsec=sector-Int_t(sector)-0.5;
 400         Double_t snp=dsec*TMath::Pi()/9.;
 401         (*pcstream)<<"delta"<<
 402           "ptype="<<ptype<<
 403           "dtype="<<dtype<<
 404           "sector="<<sector<<
 405           "dsec="<<dsec<<
 406           "snp="<<snp<<
 407           "phi="<<phi<<
 408           "localX="<<localX<<
 409           "kZ="<<kZ<<
 410           "theta="<<kZ<<
 411           "mean="<<mean<<
 412           "rms="<<rms<<
 413           "meanG="<<meanG<<
 414           "rmsG="<<rmsG<<
 415           "entries="<<entries<<
 416           "meanA="<<meanA<<
 417           "rmsA="<<rmsA<<
 418           "entriesA="<<entriesA<<
 419           "meanC="<<meanC<<
 420           "rmsC="<<rmsC<<
 421           "entriesC="<<entriesC<<
 422           "offsetA="<<offsetA<<
 423           "slopeA="<<slopeA<<
 424           "chi2A="<<chi2A<<
 425           "offsetC="<<offsetC<<
 426           "slopeC="<<slopeC<<
 427           "chi2C="<<chi2C<<
 428           "\n";
 429         delete his;
 430       }
 431       delete his2;
 432     }
 433     delete his1;
 434   }
 435   hisResMap3D->Write();
 436   hisResMap2D[0]->Write();
 437   hisResMap2D[1]->Write();
 438   hisResMap2D[2]->Write();
 439   hisResMap2D[3]->Write();
 440   //  delete pcstream;
 441 }
 442
 443
 444
 445 void   AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){
 446   //
 447   // make a distortion map out ou fthe residual histogram
 448   // Results are written to the debug streamer - pcstream
 449   // Parameters:
 450   //   his0       - input (4D) residual histogram
 451   //   pcstream   - file to write the tree
 452   //   run        - run number
 453   //   refX       - track matching reference X
 454   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
 455   // THnSparse axes:
 456   // OBJ: TAxis     #Delta  #Delta
 457   // OBJ: TAxis     tanTheta        tan(#Theta)
 458   // OBJ: TAxis     phi     #phi
 459   // OBJ: TAxis     snp     snp
 460
 461   // marian.ivanov@cern.ch
 462   const Int_t kMinEntries=10;
 463   Int_t idim[4]={0,1,2,3};
 464   //
 465   //
 466   //
 467   Int_t nbins3=his0->GetAxis(3)->GetNbins();
 468   Int_t first3=his0->GetAxis(3)->GetFirst();
 469   Int_t last3 =his0->GetAxis(3)->GetLast();
 470   //
 471   for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){   // axis 3 - local angle
 472     his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
 473     Double_t      x3= his0->GetAxis(3)->GetBinCenter(ibin3);
 474     THnSparse * his3= his0->Projection(3,idim);         //projected histogram according selection 3
 475     //
 476     Int_t nbins2    = his3->GetAxis(2)->GetNbins();
 477     Int_t first2    = his3->GetAxis(2)->GetFirst();
 478     Int_t last2     = his3->GetAxis(2)->GetLast();
 479     //
 480     for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){   // axis 2 - phi
 481       his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
 482       Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
 483       THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
 484       Int_t nbins1     = his2->GetAxis(1)->GetNbins();
 485       Int_t first1     = his2->GetAxis(1)->GetFirst();
 486       Int_t last1      = his2->GetAxis(1)->GetLast();
 487       for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
 488         //
 489         Double_t       x1= his2->GetAxis(1)->GetBinCenter(ibin1);
 490         his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
 491         if (TMath::Abs(x1)<0.1){
 492           if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
 493           if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
 494         }
 495         if (TMath::Abs(x1)<0.06){
 496           his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
 497         }
 498         TH1 * hisDelta = his2->Projection(0);
 499         //
 500         Double_t entries = hisDelta->GetEntries();
 501         Double_t mean=0, rms=0;
 502         if (entries>kMinEntries){
 503           mean    = hisDelta->GetMean();
 504           rms = hisDelta->GetRMS();
 505         }
 506         Double_t sector = 9.*x2/TMath::Pi();
 507         if (sector<0) sector+=18;
 508         Double_t dsec = sector-Int_t(sector)-0.5;
 509         Double_t z=refX*x1;
 510         (*pcstream)<<hname<<
 511           "run="<<run<<
 512           "bz="<<bz<<
 513           "theta="<<x1<<
 514           "phi="<<x2<<
 515           "z="<<z<<            // dummy z
 516           "snp="<<x3<<
 517           "entries="<<entries<<
 518           "mean="<<mean<<
 519           "rms="<<rms<<
 520           "refX="<<refX<<   // track matching refernce plane
 521           "type="<<type<<   //
 522           "sector="<<sector<<
 523           "dsec="<<dsec<<
 524           "\n";
 525         delete hisDelta;
 526         //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
 527       }
 528       delete his2;
 529     }
 530     delete his3;
 531   }
 532 }
 533
 534
 535
 536 void   AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
 537   //
 538   // make a distortion map out ou fthe residual histogram
 539   // Results are written to the debug streamer - pcstream
 540   // Parameters:
 541   //   his0       - input (4D) residual histogram
 542   //   pcstream   - file to write the tree
 543   //   run        - run number
 544   //   refX       - track matching reference X
 545   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
 546   // marian.ivanov@cern.ch
 547   //
 548   //  Histo axeses
 549   //   Collection name='TObjArray', class='TObjArray', size=16
 550   //  0. OBJ: TAxis     #Delta  #Delta
 551   //  1. OBJ: TAxis     N_{cl}  N_{cl}
 552   //  2. OBJ: TAxis     dca_{r} (cm)    dca_{r} (cm)
 553   //  3. OBJ: TAxis     z (cm)  z (cm)
 554   //  4. OBJ: TAxis     sin(#phi)       sin(#phi)
 555   //  5. OBJ: TAxis     tan(#theta)     tan(#theta)
 556   //  6. OBJ: TAxis     1/pt (1/GeV)    1/pt (1/GeV)
 557   //  7. OBJ: TAxis     pt (GeV)        pt (GeV)
 558   //  8. OBJ: TAxis     alpha   alpha
 559   const Int_t kMinEntries=10;
 560   //
 561   //  1. make default selections
 562   //
 563   TH1 * hisDelta=0;
 564   Int_t idim0[4]={0 , 5, 8,  3};   // delta, theta, alpha, z
 565   hisInput->GetAxis(1)->SetRangeUser(110,190);   //long tracks
 566   hisInput->GetAxis(2)->SetRangeUser(-10,35);    //tracks close to beam pipe
 567   hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
 568   hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
 569   hisDelta= hisInput->Projection(0);
 570   hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
 571   delete hisDelta;
 572   THnSparse *his0=  hisInput->Projection(4,idim0);
 573   //
 574   // 2. Get mean in diferent bins
 575   //
 576   Int_t nbins1=his0->GetAxis(1)->GetNbins();
 577   Int_t first1=his0->GetAxis(1)->GetFirst();
 578   Int_t last1 =his0->GetAxis(1)->GetLast();
 579   //
 580   Double_t bz=AliTrackerBase::GetBz();
 581   Int_t idim[4]={0,1, 2,  3};  // delta, theta,alpha,z
 582   //
 583   for (Int_t ibin1=first1; ibin1<=last1; ibin1++){   //axis 1 - theta
 584     //
 585     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);
 586     his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
 587     //
 588     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
 589     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
 590     Int_t first3     = his1->GetAxis(3)->GetFirst();
 591     Int_t last3      = his1->GetAxis(3)->GetLast();
 592     //
 593     for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - z at "vertex"
 594       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
 595       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
 596       if (ibin3<first3) {
 597         his1->GetAxis(3)->SetRangeUser(-1,1);
 598         x3=0;
 599       }
 600       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
 601       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
 602       Int_t first2    = his3->GetAxis(2)->GetFirst();
 603       Int_t last2     = his3->GetAxis(2)->GetLast();
 604       //
 605       for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
 606         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
 607         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
 608         hisDelta = his3->Projection(0);
 609         //
 610         Double_t entries = hisDelta->GetEntries();
 611         Double_t mean=0, rms=0;
 612         if (entries>kMinEntries){
 613           mean    = hisDelta->GetMean();
 614           rms = hisDelta->GetRMS();
 615         }
 616         Double_t sector = 9.*x2/TMath::Pi();
 617         if (sector<0) sector+=18;
 618         Double_t dsec = sector-Int_t(sector)-0.5;
 619         Double_t snp=0;  // dummy snp - equal 0
 620         (*pcstream)<<hname<<
 621           "run="<<run<<
 622           "bz="<<bz<<            // magnetic field
 623           "theta="<<x1<<         // theta
 624           "phi="<<x2<<           // phi (alpha)
 625           "z="<<x3<<             // z at "vertex"
 626           "snp="<<snp<<          // dummy snp
 627           "entries="<<entries<<  // entries in bin
 628           "mean="<<mean<<        // mean
 629           "rms="<<rms<<
 630           "refX="<<refX<<        // track matching refernce plane
 631           "type="<<type<<        // parameter type
 632           "sector="<<sector<<    // sector
 633           "dsec="<<dsec<<        // dummy delta sector
 634           "\n";
 635         delete hisDelta;
 636         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
 637       }
 638       delete his3;
 639     }
 640     delete his1;
 641   }
 642   delete his0;
 643 }
 644
 645
 646
 647 void   AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
 648   //
 649   // make a distortion map out of the residual histogram
 650   // Results are written to the debug streamer - pcstream
 651   // Parameters:
 652   //   his0       - input (4D) residual histogram
 653   //   pcstream   - file to write the tree
 654   //   run        - run number
 655   //   type       - 0- y 1-z,2 -snp, 3-theta
 656   //   bz         - magnetic field
 657   // marian.ivanov@cern.ch
 658
 659   //Collection name='TObjArray', class='TObjArray', size=16
 660   //0  OBJ: TAxis     delta   delta
 661   //1  OBJ: TAxis     phi     phi
 662   //2  OBJ: TAxis     localX  localX
 663   //3  OBJ: TAxis     kY      kY
 664   //4  OBJ: TAxis     kZ      kZ
 665   //5  OBJ: TAxis     is1     is1
 666   //6  OBJ: TAxis     is0     is0
 667   //7. OBJ: TAxis     z       z
 668   //8. OBJ: TAxis     IsPrimary       IsPrimary
 669
 670   const Int_t kMinEntries=10;
 671   THnSparse * hisSector0=0;
 672   TH1 * htemp=0;    // histogram to calculate mean value of parameter
 673   //  Double_t bz=AliTrackerBase::GetBz();
 674
 675   //
 676   // Loop over pair of sector:
 677   // isPrim         - 8  ==> 8
 678   // isec0          - 6  ==> 7
 679   //   isec1        - 5  ==> 6
 680   //     refX       - 2  ==> 5
 681   //
 682   //     phi        - 1  ==> 4
 683   //       z        - 7  ==> 3
 684   //         snp    - 3  ==> 2
 685   //           theta- 4  ==> 1
 686   //                  0  ==> 0;
 687   for (Int_t isec0=0; isec0<72; isec0++){
 688     Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
 689     //
 690     //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4);  // select secondaries only ? - get out later ?
 691     hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
 692     hisSector0=hisInput->Projection(7,index0);
 693     //
 694     //
 695     for (Int_t isec1=isec0+1; isec1<72; isec1++){
 696       //if (isec1!=isec0+36) continue;
 697       if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
 698       printf("Sectors %d\t%d\n",isec1,isec0);
 699       hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
 700       TH1 * hisX=hisSector0->Projection(5);
 701       Double_t refX= hisX->GetMean();
 702       delete hisX;
 703       TH1 *hisDelta=hisSector0->Projection(0);
 704       Double_t dmean = hisDelta->GetMean();
 705       Double_t drms = hisDelta->GetRMS();
 706       hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
 707       delete hisDelta;
 708       //
 709       //  1. make default selections
 710       //
 711       Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
 712       THnBase *hisSector1=  hisSector0->ProjectionND(5,idim0);
 713       //
 714       // 2. Get mean in diferent bins
 715       //
 716       Int_t idim[5]={0, 1, 2,  3, 4};  // {delta, theta-1,snp-2 ,z-3, phi-4}
 717       //
 718       //      Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
 719       Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
 720       Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
 721       //
 722       for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){   //axis 4 - phi
 723         //
 724         // Phi loop
 725         //
 726         Double_t       xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
 727         Double_t psec    = (9*xPhi/TMath::Pi());
 728         if (psec<0) psec+=18;
 729         Bool_t isOK0=kFALSE;
 730         Bool_t isOK1=kFALSE;
 731         if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.)  isOK0=kTRUE;
 732         if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.)  isOK1=kTRUE;
 733         if (!isOK0) continue;
 734         if (!isOK1) continue;
 735         //
 736         hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
 737         if (isec1!=isec0+36) {
 738           hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
 739         }
 740         //
 741         htemp = hisSector1->Projection(4);
 742         xPhi=htemp->GetMean();
 743         delete htemp;
 744         THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
 745         //Int_t nbinsZ     = hisPhi->GetAxis(3)->GetNbins();
 746         Int_t firstZ     = hisPhi->GetAxis(3)->GetFirst();
 747         Int_t lastZ      = hisPhi->GetAxis(3)->GetLast();
 748         //
 749         for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){   // axis 3 - z
 750           //
 751           // Z loop
 752           //
 753           hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
 754           if (isec1!=isec0+36) {
 755             hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
 756           }
 757           htemp = hisPhi->Projection(3);
 758           Double_t      xZ= htemp->GetMean();
 759           delete htemp;
 760           THnBase * hisZ= hisPhi->ProjectionND(3,idim);
 761           //projected histogram according selection 3 -z
 762           //
 763           //
 764           //Int_t nbinsSnp    = hisZ->GetAxis(2)->GetNbins();
 765           Int_t firstSnp    = hisZ->GetAxis(2)->GetFirst();
 766           Int_t lastSnp     = hisZ->GetAxis(2)->GetLast();
 767           for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){   // axis 2 - snp
 768             //
 769             // Snp loop
 770             //
 771             hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
 772             if (isec1!=isec0+36) {
 773               hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
 774             }
 775             htemp = hisZ->Projection(2);
 776             Double_t      xSnp= htemp->GetMean();
 777             delete htemp;
 778             THnBase * hisSnp= hisZ->ProjectionND(2,idim);
 779             //projected histogram according selection 2 - snp
 780
 781             //Int_t nbinsTheta    = hisSnp->GetAxis(1)->GetNbins();
 782             Int_t firstTheta    = hisSnp->GetAxis(1)->GetFirst();
 783             Int_t lastTheta     = hisSnp->GetAxis(1)->GetLast();
 784             //
 785             for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){  // axis1 theta
 786
 787
 788               hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
 789               if (isec1!=isec0+36) {
 790                  hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
 791               }
 792               htemp = hisSnp->Projection(1);
 793               Double_t xTheta=htemp->GetMean();
 794               delete htemp;
 795               hisDelta = hisSnp->Projection(0);
 796               //
 797               Double_t entries = hisDelta->GetEntries();
 798               Double_t mean=0, rms=0;
 799               if (entries>kMinEntries){
 800                 mean    = hisDelta->GetMean();
 801                 rms = hisDelta->GetRMS();
 802               }
 803               Double_t sector = 9.*xPhi/TMath::Pi();
 804               if (sector<0) sector+=18;
 805               Double_t dsec = sector-Int_t(sector)-0.5;
 806               Int_t dtype=1;  // TPC alignment type
 807               (*pcstream)<<hname<<
 808                 "run="<<run<<
 809                 "bz="<<bz<<             // magnetic field
 810                 "ptype="<<type<<         // parameter type
 811                 "dtype="<<dtype<<         // parameter type
 812                 "isec0="<<isec0<<       // sector 0
 813                 "isec1="<<isec1<<       // sector 1
 814                 "sector="<<sector<<     // sector as float
 815                 "dsec="<<dsec<<         // delta sector
 816                 //
 817                 "theta="<<xTheta<<      // theta
 818                 "phi="<<xPhi<<          // phi (alpha)
 819                 "z="<<xZ<<              // z
 820                 "snp="<<xSnp<<          // snp
 821
 822                 //
 823                 "entries="<<entries<<  // entries in bin
 824                 "mean="<<mean<<        // mean
 825                 "rms="<<rms<<          // rms
 826                 "refX="<<refX<<        // track matching reference plane
 827                 "\n";
 828               delete hisDelta;
 829               //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
 830               //
 831             }//ibinTheta
 832             delete hisSnp;
 833           } //ibinSnp
 834           delete hisZ;
 835         }//ibinZ
 836         delete hisPhi;
 837       }//ibinPhi
 838       delete hisSector1;
 839     }//isec1
 840     delete hisSector0;
 841   }//isec0
 842 }
 843
 844
 845
 846
 847
 848 // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
 849 //   //
 850 //   // Make a laser fit tree for global minimization
 851 //   //
 852 //   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
 853 //   AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
 854 //   if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
 855 //   correction->AddVisualCorrection(correction,0);  //register correction
 856
 857 //   //  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 858 //   //AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 859 //   //
 860 //   const Double_t cutErrY=0.05;
 861 //   const Double_t kSigmaCut=4;
 862 //   //  const Double_t cutErrZ=0.03;
 863 //   const Double_t kEpsilon=0.00000001;
 864 //   //  const Double_t kMaxDist=1.;  // max distance - space correction
 865 //   TVectorD *vecdY=0;
 866 //   TVectorD *vecdZ=0;
 867 //   TVectorD *veceY=0;
 868 //   TVectorD *veceZ=0;
 869 //   AliTPCLaserTrack *ltr=0;
 870 //   AliTPCLaserTrack::LoadTracks();
 871 //   tree->SetBranchAddress("dY.",&vecdY);
 872 //   tree->SetBranchAddress("dZ.",&vecdZ);
 873 //   tree->SetBranchAddress("eY.",&veceY);
 874 //   tree->SetBranchAddress("eZ.",&veceZ);
 875 //   tree->SetBranchAddress("LTr.",&ltr);
 876 //   Int_t entries= tree->GetEntries();
 877 //   TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
 878 //   Double_t bz=AliTrackerBase::GetBz();
 879 //   //
 880 //   //  Double_t globalXYZ[3];
 881 //   //Double_t globalXYZCorr[3];
 882 //   for (Int_t ientry=0; ientry<entries; ientry++){
 883 //     tree->GetEntry(ientry);
 884 //     if (!ltr->GetVecGX()){
 885 //       ltr->UpdatePoints();
 886 //     }
 887 //     //
 888 //     TVectorD fit10(5);
 889 //     TVectorD fit5(5);
 890 //     printf("Entry\t%d\n",ientry);
 891 //     for (Int_t irow0=0; irow0<158; irow0+=1){
 892 //       //
 893 //       TLinearFitter fitter10(4,"hyp3");
 894 //       TLinearFitter fitter5(2,"hyp1");
 895 //       Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
 896 //       if (sector<0) continue;
 897 //       //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
 898
 899 //       Double_t refX= (*ltr->GetVecLX())[irow0];
 900 //       Int_t firstRow1 = TMath::Max(irow0-10,0);
 901 //       Int_t lastRow1  = TMath::Min(irow0+10,158);
 902 //       Double_t padWidth=(irow0<64)?0.4:0.6;
 903 //       // make long range fit
 904 //       for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
 905 //      if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
 906 //      if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
 907 //      if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
 908 //      Double_t idealX= (*ltr->GetVecLX())[irow1];
 909 //      Double_t idealY= (*ltr->GetVecLY())[irow1];
 910 //      //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
 911 //      Double_t gx= (*ltr->GetVecGX())[irow1];
 912 //      Double_t gy= (*ltr->GetVecGY())[irow1];
 913 //      Double_t gz= (*ltr->GetVecGZ())[irow1];
 914 //      Double_t measY=(*vecdY)[irow1]+idealY;
 915 //      Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
 916 //      // deltaR = R distorted -R ideal
 917 //      Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
 918 //      fitter10.AddPoint(xxx,measY,1);
 919 //       }
 920 //       Bool_t isOK=kTRUE;
 921 //       Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
 922 //       Double_t mean10  =0;//   fitter10.GetParameter(0);
 923 //       Double_t slope10  =0;//   fitter10.GetParameter(0);
 924 //       Double_t cosPart10  = 0;//  fitter10.GetParameter(2);
 925 //       Double_t sinPart10   =0;//  fitter10.GetParameter(3);
 926
 927 //       if (fitter10.GetNpoints()>10){
 928 //      fitter10.Eval();
 929 //      rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
 930 //      mean10      =   fitter10.GetParameter(0);
 931 //      slope10     =   fitter10.GetParameter(1);
 932 //      cosPart10   =   fitter10.GetParameter(2);
 933 //      sinPart10   =  fitter10.GetParameter(3);
 934 //      //
 935 //      // make short range fit
 936 //      //
 937 //      for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
 938 //        if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
 939 //        if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
 940 //        if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
 941 //        Double_t idealX= (*ltr->GetVecLX())[irow1];
 942 //        Double_t idealY= (*ltr->GetVecLY())[irow1];
 943 //        //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
 944 //        Double_t gx= (*ltr->GetVecGX())[irow1];
 945 //        Double_t gy= (*ltr->GetVecGY())[irow1];
 946 //        Double_t gz= (*ltr->GetVecGZ())[irow1];
 947 //        Double_t measY=(*vecdY)[irow1]+idealY;
 948 //        Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
 949 //        // deltaR = R distorted -R ideal
 950 //        Double_t expY= mean10+slope10*(idealX+deltaR-refX);
 951 //        if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
 952 //        //
 953 //        Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
 954 //        Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
 955 //        fitter5.AddPoint(xxx,measY-corr,1);
 956 //      }
 957 //       }else{
 958 //      isOK=kFALSE;
 959 //       }
 960 //       if (fitter5.GetNpoints()<8) isOK=kFALSE;
 961
 962 //       Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
 963 //       Double_t offset5  =0;//  fitter5.GetParameter(0);
 964 //       Double_t slope5   =0;//  fitter5.GetParameter(0);
 965 //       if (isOK){
 966 //      fitter5.Eval();
 967 //      rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
 968 //      offset5  =  fitter5.GetParameter(0);
 969 //      slope5   =  fitter5.GetParameter(0);
 970 //       }
 971 //       //
 972 //       Double_t dtype=5;
 973 //       Double_t ptype=0;
 974 //       Double_t phi   =(*ltr->GetVecPhi())[irow0];
 975 //       Double_t theta =ltr->GetTgl();
 976 //       Double_t mean=(vecdY)->GetMatrixArray()[irow0];
 977 //       Double_t gx=0,gy=0,gz=0;
 978 //       Double_t snp = (*ltr->GetVecP2())[irow0];
 979 //       Int_t bundle= ltr->GetBundle();
 980 //       Int_t id= ltr->GetId();
 981 //       //      Double_t rms = err->GetMatrixArray()[irow];
 982 //       //
 983 //       gx = (*ltr->GetVecGX())[irow0];
 984 //       gy = (*ltr->GetVecGY())[irow0];
 985 //       gz = (*ltr->GetVecGZ())[irow0];
 986 //       Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
 987 //       fitter10.GetParameters(fit10);
 988 //       fitter5.GetParameters(fit5);
 989 //       Double_t idealY= (*ltr->GetVecLY())[irow0];
 990 //       Double_t measY=(*vecdY)[irow0]+idealY;
 991 //       Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
 992 //       if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
 993 //       //
 994 //       (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
 995 //      "bz="<<bz<<         // magnetic filed used
 996 //      "dtype="<<dtype<<   // detector match type
 997 //      "ptype="<<ptype<<   // parameter type
 998 //      "theta="<<theta<<   // theta
 999 //      "phi="<<phi<<       // phi
1000 //      "snp="<<snp<<       // snp
1001 //      "sector="<<sector<<
1002 //      "bundle="<<bundle<<
1003 // //   //      "dsec="<<dsec<<
1004 //      "refX="<<refX<<      // reference radius
1005 //      "gx="<<gx<<         // global position
1006 //      "gy="<<gy<<         // global position
1007 //      "gz="<<gz<<         // global position
1008 //      "dRrec="<<dRrec<<      // delta Radius in reconstruction
1009 //      "id="<<id<<     //bundle
1010 //      "rms10="<<rms10<<
1011 //      "rms5="<<rms5<<
1012 //      "fit10.="<<&fit10<<
1013 //      "fit5.="<<&fit5<<
1014 //      "measY="<<measY<<
1015 //      "mean="<<mean<<
1016 //      "idealY="<<idealY<<
1017 //      "corr="<<corr<<
1018 //      "isOK="<<isOK<<
1019 //      "\n";
1020 //     }
1021 //   }
1022 //   delete pcstream;
1023 // }
1024
1025
1026 void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1027   //
1028   // Make a fit tree:
1029   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1030   // calculates partial distortions
1031   // Partial distortion is stored in the resulting tree
1032   // Output is storred in the file distortion_<dettype>_<partype>.root
1033   // Partial  distortion is stored with the name given by correction name
1034   //
1035   //
1036   // Parameters of function:
1037   // input     - input tree
1038   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF,  4 - TPCTPC track crossing
1039   // ppype     - parameter type
1040   // corrArray - array with partial corrections
1041   // step      - skipe entries  - if 1 all entries processed - it is slow
1042   // debug     0 if debug on also space points dumped - it is slow
1043
1044   const Double_t kMaxSnp = 0.85;
1045   const Double_t kcutSnp=0.25;
1046   const Double_t kcutTheta=1.;
1047   const Double_t kRadiusTPC=85;
1048   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();
1049   //
1050   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1051   //  const Double_t kB2C=-0.299792458e-3;
1052   const Int_t kMinEntries=20;
1053   Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
1054   Float_t refX;
1055   Int_t run;
1056   tinput->SetBranchAddress("run",&run);
1057   tinput->SetBranchAddress("theta",&theta);
1058   tinput->SetBranchAddress("phi", &phi);
1059   tinput->SetBranchAddress("snp",&snp);
1060   tinput->SetBranchAddress("mean",&mean);
1061   tinput->SetBranchAddress("rms",&rms);
1062   tinput->SetBranchAddress("entries",&entries);
1063   tinput->SetBranchAddress("sector",&sector);
1064   tinput->SetBranchAddress("dsec",&dsec);
1065   tinput->SetBranchAddress("refX",&refX);
1066   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1067   //
1068   Int_t nentries=tinput->GetEntries();
1069   Int_t ncorr=corrArray->GetEntries();
1070   Double_t corrections[100]={0}; //
1071   Double_t tPar[5];
1072   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1073   Int_t dir=0;
1074   if (dtype==5 || dtype==6) dtype=4;
1075   if (dtype==0) { dir=-1;}
1076   if (dtype==1) { dir=1;}
1077   if (dtype==2) { dir=-1;}
1078   if (dtype==3) { dir=1;}
1079   if (dtype==4) { dir=-1;}
1080   //
1081   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1082     tinput->GetEntry(ientry);
1083     if (TMath::Abs(snp)>kMaxSnp) continue;
1084     tPar[0]=0;
1085     tPar[1]=theta*refX;
1086     if (dtype==2)  tPar[1]=theta*kRadiusTPC;
1087     tPar[2]=snp;
1088     tPar[3]=theta;
1089     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
1090     if (dtype==4){
1091       // tracks crossing CE
1092       tPar[1]=0;   // track at the CE
1093       //if (TMath::Abs(theta) <0.05) continue;  // deep cross
1094     }
1095
1096     if (TMath::Abs(snp) >kcutSnp) continue;
1097     if (TMath::Abs(theta) >kcutTheta) continue;
1098     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1099     Double_t bz=AliTrackerBase::GetBz();
1100     if (dtype !=4) { //exclude TPC  - for TPC mainly non primary tracks
1101       if (dtype!=2 && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
1102
1103       if (dtype==2 && TMath::Abs(bz)>0.1 )  {
1104         tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1105         // snp at the TPC inner radius in case the vertex match used
1106       }
1107     }
1108     //
1109     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1110     AliExternalTrackParam track(refX,phi,tPar,cov);
1111     Double_t xyz[3];
1112     track.GetXYZ(xyz);
1113     Int_t id=0;
1114     Double_t pt=1./tPar[4];
1115     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1116     //if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
1117     Double_t refXD=refX;
1118     (*pcstream)<<"fit"<<
1119       "run="<<run<<       // run number
1120       "bz="<<bz<<         // magnetic filed used
1121       "dtype="<<dtype<<   // detector match type
1122       "ptype="<<ptype<<   // parameter type
1123       "theta="<<theta<<   // theta
1124       "phi="<<phi<<       // phi
1125       "snp="<<snp<<       // snp
1126       "mean="<<mean<<     // mean dist value
1127       "rms="<<rms<<       // rms
1128       "sector="<<sector<<
1129       "dsec="<<dsec<<
1130       "refX="<<refXD<<         // referece X as double
1131       "gx="<<xyz[0]<<         // global position at reference
1132       "gy="<<xyz[1]<<         // global position at reference
1133       "gz="<<xyz[2]<<         // global position at reference
1134       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1135       "pt="<<pt<<            // pt
1136       "id="<<id<<             // track id
1137       "entries="<<entries;// number of entries in bin
1138     //
1139     Bool_t isOK=kTRUE;
1140     if (entries<kMinEntries) isOK=kFALSE;
1141     //
1142     if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1143       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1144       corrections[icorr]=0;
1145       if (entries>kMinEntries){
1146         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1147         AliExternalTrackParam *trackOut = 0;
1148         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1149         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1150         if (dtype==0) {dir= -1;}
1151         if (dtype==1) {dir=  1;}
1152         if (dtype==2) {dir= -1;}
1153         if (dtype==3) {dir=  1;}
1154         //
1155         if (trackOut){
1156           if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1157           if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1158           if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1159           //      trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1160           //
1161           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1162           delete trackOut;
1163         }else{
1164           corrections[icorr]=0;
1165           isOK=kFALSE;
1166         }
1167         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out
1168       }
1169       (*pcstream)<<"fit"<<
1170         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1171     }
1172
1173     if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1174       //
1175       // special case of the TPC tracks crossing the CE
1176       //
1177       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1178       corrections[icorr]=0;
1179       if (entries>kMinEntries){
1180         AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1181         AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
1182         AliExternalTrackParam *trackOut0 = 0;
1183         AliExternalTrackParam *trackOut1 = 0;
1184         //
1185         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1186         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1187         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1188         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1189         //
1190         if (trackOut0 && trackOut1){
1191           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp))  isOK=kFALSE;
1192           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1193           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1194           if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1195           //
1196           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1197           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
1198           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1199           if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1200           if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1201           //
1202           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1203           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1204           if (isOK)
1205             if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1206                 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1207                 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1208                 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1209                 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1210                 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1211                 ){
1212               isOK=kFALSE;
1213             }
1214           delete trackOut0;
1215           delete trackOut1;
1216         }else{
1217           corrections[icorr]=0;
1218           isOK=kFALSE;
1219         }
1220         //
1221         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out no in lookup
1222       }
1223       (*pcstream)<<"fit"<<
1224         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1225     }
1226     //
1227     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1228   }
1229
1230
1231   delete pcstream;
1232 }
1233
1234
1235
1236 void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1237   //
1238   // Make a fit tree:
1239   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1240   // calculates partial distortions
1241   // Partial distortion is stored in the resulting tree
1242   // Output is storred in the file distortion_<dettype>_<partype>.root
1243   // Partial  distortion is stored with the name given by correction name
1244   //
1245   //
1246   // Parameters of function:
1247   // input     - input tree
1248   // dtype     - distortion type 10 - IROC-OROC
1249   // ppype     - parameter type
1250   // corrArray - array with partial corrections
1251   // step      - skipe entries  - if 1 all entries processed - it is slow
1252   // debug     0 if debug on also space points dumped - it is slow
1253
1254   const Double_t kMaxSnp = 0.8;
1255   const Int_t kMinEntries=200;
1256   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();
1257   //
1258   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1259   //  const Double_t kB2C=-0.299792458e-3;
1260   Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1261   Int_t isec1, isec0;
1262   Double_t refXD;
1263   Float_t refX;
1264   Int_t run;
1265   tinput->SetBranchAddress("run",&run);
1266   tinput->SetBranchAddress("theta",&theta);
1267   tinput->SetBranchAddress("phi", &phi);
1268   tinput->SetBranchAddress("snp",&snp);
1269   tinput->SetBranchAddress("mean",&mean);
1270   tinput->SetBranchAddress("rms",&rms);
1271   tinput->SetBranchAddress("entries",&entries);
1272   tinput->SetBranchAddress("sector",&sector);
1273   tinput->SetBranchAddress("dsec",&dsec);
1274   tinput->SetBranchAddress("refX",&refXD);
1275   tinput->SetBranchAddress("z",&globalZ);
1276   tinput->SetBranchAddress("isec0",&isec0);
1277   tinput->SetBranchAddress("isec1",&isec1);
1278   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1279   //
1280   Int_t nentries=tinput->GetEntries();
1281   Int_t ncorr=corrArray->GetEntries();
1282   Double_t corrections[100]={0}; //
1283   Double_t tPar[5];
1284   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1285   Int_t dir=0;
1286   //
1287   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1288     tinput->GetEntry(ientry);
1289     refX=refXD;
1290     Int_t id=-1;
1291     if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1;  // IROC-OROC - opposite side
1292     if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2;  // IROC-OROC - same side
1293     if (dtype==10  && id==-1) continue;
1294     //
1295     dir=-1;
1296     tPar[0]=0;
1297     tPar[1]=globalZ;
1298     tPar[2]=snp;
1299     tPar[3]=theta;
1300     tPar[4]=(gRandom->Rndm()-0.1)*0.2;  //
1301     Double_t pt=1./tPar[4];
1302     //
1303     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1304     Double_t bz=AliTrackerBase::GetBz();
1305     AliExternalTrackParam track(refX,phi,tPar,cov);
1306     Double_t xyz[3],xyzIn[3],xyzOut[3];
1307     track.GetXYZ(xyz);
1308     track.GetXYZAt(85,bz,xyzIn);
1309     track.GetXYZAt(245,bz,xyzOut);
1310     Double_t phiIn  = TMath::ATan2(xyzIn[1],xyzIn[0]);
1311     Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1312     Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1313     Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1314     Int_t sectorIn  = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1315     Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1316     //
1317     Bool_t isOK=kTRUE;
1318     if (sectorIn!=sectorOut) isOK=kFALSE;  // requironment - cluster in the same sector
1319     if (sectorIn!=sectorRef) isOK=kFALSE;  // requironment - cluster in the same sector
1320     if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE;  // requironment - minimal amount of tracks in bin
1321     // Do downscale
1322     if (TMath::Abs(theta)>1) isOK=kFALSE;
1323     //
1324     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1325     //
1326     (*pcstream)<<"fit"<<
1327       "run="<<run<<       //run
1328       "bz="<<bz<<         // magnetic filed used
1329       "dtype="<<dtype<<   // detector match type
1330       "ptype="<<ptype<<   // parameter type
1331       "theta="<<theta<<   // theta
1332       "phi="<<phi<<       // phi
1333       "snp="<<snp<<       // snp
1334       "mean="<<mean<<     // mean dist value
1335       "rms="<<rms<<       // rms
1336       "sector="<<sector<<
1337       "dsec="<<dsec<<
1338       "refX="<<refXD<<         // referece X
1339       "gx="<<xyz[0]<<         // global position at reference
1340       "gy="<<xyz[1]<<         // global position at reference
1341       "gz="<<xyz[2]<<         // global position at reference
1342       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1343       "pt="<<pt<<      //pt
1344       "id="<<id<<             // track id
1345       "entries="<<entries;// number of entries in bin
1346     //
1347     AliExternalTrackParam *trackOut0 = 0;
1348     AliExternalTrackParam *trackOut1 = 0;
1349     AliExternalTrackParam *ptrackIn0 = 0;
1350     AliExternalTrackParam *ptrackIn1 = 0;
1351
1352     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1353       //
1354       // special case of the TPC tracks crossing the CE
1355       //
1356       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1357       corrections[icorr]=0;
1358       if (entries>kMinEntries &&isOK){
1359         AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
1360         AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
1361         ptrackIn1=&trackIn0;
1362         ptrackIn0=&trackIn1;
1363         //
1364         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1365         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1366         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1367         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1368         //
1369         if (trackOut0 && trackOut1){
1370           //
1371           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp))  isOK=kFALSE;
1372           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1373           // rotate all tracks to the same frame
1374           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1375           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
1376           if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1377           //
1378           if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1379           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1380           if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1381           //
1382           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1383           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1384           (*pcstream)<<"fitDebug"<< // just to debug the correction
1385             "mean="<<mean<<
1386             "pIn0.="<<ptrackIn0<<
1387             "pIn1.="<<ptrackIn1<<
1388             "pOut0.="<<trackOut0<<
1389             "pOut1.="<<trackOut1<<
1390             "refX="<<refXD<<
1391             "\n";
1392           delete trackOut0;
1393           delete trackOut1;
1394         }else{
1395           corrections[icorr]=0;
1396           isOK=kFALSE;
1397         }
1398       }
1399       (*pcstream)<<"fit"<<
1400         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1401     }
1402     //
1403     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1404   }
1405   delete pcstream;
1406 }
1407
1408