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