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