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