TPC/AliTPCcalibUnlinearity.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   Class for histogramming of cluster residuals
  18   and fitting of the unlinearities. To be used only for data without
  19   magnetic field. The laser tracks should be also rejected.
  20   //
  21
  22   Track fitting:
  23   The track is fitted using linear and parabolic model
  24   The edge clusters are removed from the fit.
  25   Edge pad-row - first and last 15 padrows removed from the fit
  26
  27   Unlinearities fitting:
  28   Unlinearities at the edge aproximated using two exponential decays.
  29
  30   Model:
  31   dz = dz0(r,z) +dr(r,z)*tan(theta)
  32   dy =          +dr(r,z)*tan(phi)
  33
  34
  35   .x ~/NimStyle.C
  36   gSystem->Load("libANALYSIS");
  37   gSystem->Load("libTPCcalib");
  38   TFile fcalib("CalibObjects.root");
  39   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  40   AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
  41   //
  42
  43 */
  44
  45
  46 #include "TLinearFitter.h"
  47
  48 #include "Riostream.h"
  49 #include "TChain.h"
  50 #include "TTree.h"
  51 #include "TH1F.h"
  52 #include "TH2F.h"
  53 #include "TH3F.h"
  54 #include "THnSparse.h"
  55 #include "TList.h"
  56 #include "TMath.h"
  57 #include "TCanvas.h"
  58 #include "TFile.h"
  59 #include "TF1.h"
  60 #include "TVectorD.h"
  61 #include "TProfile.h"
  62 #include "TGraphErrors.h"
  63 #include "TCanvas.h"
  64 #include "TCut.h"
  65
  66
  67 #include "AliTPCclusterMI.h"
  68 #include "AliTPCseed.h"
  69 #include "AliESDVertex.h"
  70 #include "AliESDEvent.h"
  71 #include "AliESDfriend.h"
  72 #include "AliESDInputHandler.h"
  73 #include "AliAnalysisManager.h"
  74 #include "AliTracker.h"
  75 #include "AliMagF.h"
  76 #include "AliTPCCalROC.h"
  77
  78 #include "AliLog.h"
  79
  80
  81 #include "TTreeStream.h"
  82 #include "AliTPCTracklet.h"
  83 #include "TTimeStamp.h"
  84 #include "AliTPCcalibDB.h"
  85 #include "AliDCSSensorArray.h"
  86 #include "AliDCSSensor.h"
  87 #include "TLinearFitter.h"
  88 #include "AliTPCClusterParam.h"
  89 #include "TStatToolkit.h"
  90
  91
  92 #include "AliTPCcalibUnlinearity.h"
  93 #include "AliTPCPointCorrection.h"
  94
  95
  96 ClassImp(AliTPCcalibUnlinearity)
  97
  98 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
  99   AliTPCcalibBase(),
 100   fInit(kFALSE),
 101   fDiffHistoLineY(0),
 102   fDiffHistoParY(0),
 103   fDiffHistoLineZ(0),
 104   fDiffHistoParZ(0),
 105   fFittersOutR(38),
 106   fFittersOutZ(38),
 107   fParamsOutR(38),
 108   fParamsOutZ(38),
 109   fErrorsOutR(38),
 110   fErrorsOutZ(38),
 111   fDistRPHIPlus(74),
 112   fDistRPHIMinus(74),
 113   fFitterQuadrantY(16*38),
 114   fFitterQuadrantPhi(16*38)
 115 {
 116   //
 117   // Defualt constructor
 118   //
 119 }
 120
 121 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
 122   AliTPCcalibBase(name,title),
 123   fInit(kFALSE),
 124   fDiffHistoLineY(0),
 125   fDiffHistoParY(0),
 126   fDiffHistoLineZ(0),
 127   fDiffHistoParZ(0),
 128   fFittersOutR(38),
 129   fFittersOutZ(38),
 130   fParamsOutR(38),
 131   fParamsOutZ(38),
 132   fErrorsOutR(38),
 133   fErrorsOutZ(38),
 134   fDistRPHIPlus(74),
 135   fDistRPHIMinus(74),
 136   fFitterQuadrantY(16*38),
 137   fFitterQuadrantPhi(16*38)
 138 {
 139   //
 140   // Non default constructor
 141   //
 142   MakeHisto();
 143 }
 144
 145 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
 146   //
 147   //
 148   //
 149   if (fDiffHistoLineZ) delete fDiffHistoLineY;
 150   if (fDiffHistoParZ)  delete fDiffHistoParY;
 151   if (fDiffHistoLineZ) delete fDiffHistoLineZ;
 152   if (fDiffHistoParZ)  delete fDiffHistoParZ;
 153 }
 154
 155
 156
 157
 158
 159 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
 160   //
 161   //
 162   //
 163   if (!fInit) {
 164     Init();   // work around for PROOF - initialize firs time used
 165   }
 166   const Int_t  kdrow=10;
 167   const Int_t kMinCluster=35;
 168   if (!seed) return;
 169   if (TMath::Abs(fMagF)>0.01) return;   // use only data without field
 170   Int_t counter[72];
 171   for (Int_t i=0;i<72;i++) counter[i]=0;
 172   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 173     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 174     if (!cluster) continue;
 175     if (cluster->GetType()<0) continue;
 176     counter[cluster->GetDetector()]++;
 177   }
 178   //
 179   for (Int_t is=0; is<72;is++){
 180     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
 181     if (counter[is]>kMinCluster) {
 182       AlignOROC(seed, is);
 183     }
 184   }
 185 }
 186
 187
 188 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
 189   //
 190   // process differnce of the cluster in respect with linear and parabolic fit
 191   //
 192   const  Double_t kXmean=134;
 193   const Int_t     kdrow=10;
 194   const Float_t kSagitaCut = 1;
 195   static TLinearFitter fy1(2,"pol1");
 196   static TLinearFitter fy2(3,"pol2");
 197   static TLinearFitter fz1(2,"pol1");
 198   static TLinearFitter fz2(3,"pol2");
 199   //
 200   static TVectorD py1(2);
 201   static TVectorD py2(3);
 202   static TVectorD pz1(2);
 203   static TVectorD pz2(3);
 204   //
 205   fy1.ClearPoints();
 206   fy2.ClearPoints();
 207   fz1.ClearPoints();
 208   fz2.ClearPoints();
 209   //
 210   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 211     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 212     if (!c) continue;
 213     if (c->GetDetector()!=isec) continue;
 214     if (c->GetType()<0) continue;
 215     Double_t dx = c->GetX()-kXmean;
 216     Double_t x[2]={dx, dx*dx};
 217     fy1.AddPoint(x,c->GetY());
 218     fy2.AddPoint(x,c->GetY());
 219     fz1.AddPoint(x,c->GetZ());
 220     fz2.AddPoint(x,c->GetZ());
 221   }
 222   fy1.Eval(); fy1.GetParameters(py1);
 223   fy2.Eval(); fy2.GetParameters(py2);
 224   fz1.Eval(); fz1.GetParameters(pz1);
 225   fz2.Eval(); fz2.GetParameters(pz2);
 226   //
 227
 228   for (Int_t irow=0;irow<159;irow++) {
 229     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 230     if (!c) continue;
 231     if (c->GetDetector()!=isec) continue;
 232     Double_t dx = c->GetX()-kXmean;
 233     Double_t y1 = py1[0]+py1[1]*dx;
 234     Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
 235     //
 236     Double_t z1 = pz1[0]+pz1[1]*dx;
 237     Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
 238     //
 239     Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
 240     Double_t edgePlus =  y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
 241     Int_t    npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
 242     Float_t    dpad  = TMath::Min(c->GetPad(), npads-1- c->GetPad());
 243     Float_t dy =c->GetY()-y1;
 244     //
 245     // Corrections
 246     //
 247     AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
 248     Double_t corrclY=0, corrtrY=0, corrR=0;
 249     corrclY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
 250                                        c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5);
 251     corrtrY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
 252                                        c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);
 253     corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
 254     //
 255     //
 256     //
 257     if (fStreamLevel>1){
 258       TTreeSRedirector *cstream = GetDebugStreamer();
 259       if (cstream){
 260         (*cstream)<<"Diff"<<
 261           "run="<<fRun<<              //  run number
 262           "event="<<fEvent<<          //  event number
 263           "time="<<fTime<<            //  time stamp of event
 264           "trigger="<<fTrigger<<      //  trigger
 265           "mag="<<fMagF<<             //  magnetic field
 266           "isec="<<isec<<             // current sector
 267           "npads="<<npads<<           // number of pads at given sector
 268           "dpad="<<dpad<<             // distance to edge pad
 269           //
 270           "crR="<<corrR<<              // radial correction
 271           "cclY="<<corrclY<<            // edge effect correction cl
 272           "ctrY="<<corrtrY<<            // edge effect correction using track
 273           //
 274           "Cl.="<<c<<
 275           "y1="<<y1<<
 276           "y2="<<y2<<
 277           "z1="<<z1<<
 278           "z2="<<z2<<
 279           "py1.="<<&py1<<
 280           "py2.="<<&py2<<
 281           "pz1.="<<&pz1<<
 282           "pz2.="<<&pz2<<
 283           "eP="<<edgePlus<<
 284           "eM="<<edgeMinus<<
 285           "\n";
 286       }
 287     }
 288     if (TMath::Min(edgeMinus,edgePlus)<6){
 289       AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
 290       TTreeSRedirector *cstream = GetDebugStreamer();
 291       if (TMath::Abs(dy)<2 && cstream){
 292         (*cstream)<<"rphi"<<
 293           "isec="<<isec<<             // current sector
 294           "npads="<<npads<<           // number of pads at given sector
 295           "dpad="<<dpad<<             // distance to edge pad
 296           "eP="<<edgePlus<<           // distance to edge
 297           "eM="<<edgeMinus<<          // distance to edge minus
 298           //
 299           "crR="<<corrR<<              // radial correction
 300           "cclY="<<corrclY<<            // edge effect correction cl
 301           "ctrY="<<corrtrY<<            // edge effect correction using track
 302           //
 303           "dy="<<dy<<                 // dy
 304           "Cl.="<<c<<
 305           "y1="<<y1<<
 306           "y2="<<y2<<
 307           "z1="<<z1<<
 308           "z2="<<z2<<
 309           "py1.="<<&py1<<
 310           "pz1.="<<&pz1<<
 311           "\n";
 312       }
 313     }
 314     if (c->GetType()<0) continue;  // don't use edge rphi cluster
 315     //
 316     //
 317     Double_t x[10];
 318     x[0]=c->GetDetector();
 319     x[1]=c->GetRow();
 320     x[2]=c->GetY()/c->GetX();
 321     x[3]=c->GetZ();
 322     //
 323     // apply sagita cut
 324     //
 325     if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut &&
 326         TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
 327       //
 328       x[4]=py1[1];
 329       x[5]=c->GetY()-y1;
 330       fDiffHistoLineY->Fill(x);
 331       x[5]=c->GetY()-y2;
 332       //fDiffHistoParY->Fill(x);
 333       //
 334       x[4]=pz1[1];
 335       x[5]=c->GetY()-z1;
 336       fDiffHistoLineZ->Fill(x);
 337       x[5]=c->GetY()-z2;
 338       //fDiffHistoParZ->Fill(x);
 339       // Apply sagita cut
 340       AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
 341     }
 342
 343   }
 344 }
 345
 346
 347 void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
 348   //
 349   // fit the tracklet in OROC sepatately in Quadrant
 350   //
 351   //
 352   if (isec<36) return;
 353   const Int_t kMinClusterF=35;
 354   const Int_t kMinClusterQ=10;
 355   //
 356   const  Int_t     kdrow1 =8;       // rows to skip at the end
 357   const  Int_t     kdrow0 =2;        // rows to skip at beginning
 358   const  Float_t   kedgey=3;
 359   const  Float_t   kMaxDist=0.5;
 360   const  Float_t   kMaxCorrY=0.1;
 361   const  Float_t   kPRFWidth = 0.6;   //cut 2 sigma of PRF
 362   //
 363   //
 364   AliTPCROC * roc = AliTPCROC::Instance();
 365   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
 366   const  Double_t kXmean=roc->GetPadRowRadii(isec,53);
 367   //
 368   // full track fit parameters
 369   //
 370   TLinearFitter fyf(2,"pol1");
 371   TLinearFitter fzf(2,"pol1");
 372   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
 373   Int_t nf=0;
 374   //
 375   // make full fit as reference
 376   //
 377   for (Int_t iter=0; iter<2; iter++){
 378     fyf.ClearPoints();
 379     for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
 380       AliTPCclusterMI *c=track->GetClusterPointer(irow);
 381       if (!c) continue;
 382       if (c->GetDetector()!=isec) continue;
 383       if (c->GetRow()<kdrow0) continue;
 384       Double_t dx = c->GetX()-kXmean;
 385       Double_t x[2]={dx, dx*dx};
 386       if (iter==0 &&c->GetType()<0) continue;
 387       if (iter==1){
 388         Double_t yfit  =  pyf[0]+pyf[1]*dx;
 389         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
 390         if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
 391         if (dedge<kedgey) continue;
 392         Double_t corrtrY =
 393           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
 394                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
 395         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
 396       }
 397       fyf.AddPoint(x,c->GetY(),0.1);
 398       fzf.AddPoint(x,c->GetZ(),0.1);
 399     }
 400     nf = fyf.GetNpoints();
 401     if (nf<kMinClusterF) return;   // not enough points - skip
 402     fyf.Eval();
 403     fyf.GetParameters(pyf);
 404     fyf.GetErrors(peyf);
 405     fzf.Eval();
 406     fzf.GetParameters(pzf);
 407     fzf.GetErrors(pezf);
 408   }
 409   //
 410   // Make Fitters and params for 3 fitters
 411   //
 412   TLinearFitter *fitters[4];
 413   Int_t npoints[4];
 414   TVectorD params[4];
 415   TVectorD errors[4];
 416   Double_t chi2C[4];
 417   for (Int_t i=0;i<4;i++) {
 418     fitters[i] = new TLinearFitter(2,"pol1");
 419     npoints[i]=0;
 420     params[i].ResizeTo(2);
 421     errors[i].ResizeTo(2);
 422   }
 423   //
 424   // Update fitters
 425   //
 426   for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
 427     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 428     if (!c) continue;
 429     if (c->GetDetector()!=isec) continue;
 430     if (c->GetRow()<kdrow0) continue;
 431     Double_t dx = c->GetX()-kXmean;
 432     Double_t x[2]={dx, dx*dx};
 433     Double_t yfit  =  pyf[0]+pyf[1]*dx;
 434     Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
 435     if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
 436     if (dedge<kedgey) continue;
 437     Double_t corrtrY =
 438       corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
 439                               c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
 440     if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
 441     if (dx<0){
 442       if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
 443       if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
 444     }
 445     if (dx>0){
 446       if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
 447       if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
 448     }
 449   }
 450
 451   for (Int_t i=0;i<4;i++) {
 452     npoints[i] = fitters[i]->GetNpoints();
 453     if (npoints[i]>=kMinClusterQ){
 454       fitters[i]->Eval();
 455       Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
 456       chi2C[i]=chi2Fac;
 457       fitters[i]->GetParameters(params[i]);
 458       fitters[i]->GetErrors(errors[i]);
 459       // renormalize errors
 460       errors[i][0]*=chi2Fac;
 461       errors[i][1]*=chi2Fac;
 462     }
 463   }
 464   //
 465   // Fill fitters
 466   //
 467   for (Int_t i0=0;i0<4;i0++){
 468     for (Int_t i1=0;i1<4;i1++){
 469       if (i0==i1) continue;
 470       if(npoints[i0]<kMinClusterQ) continue;
 471       if(npoints[i1]<kMinClusterQ) continue;
 472       Int_t index0=i0*4+i1;
 473       Double_t xy[1];
 474       Double_t xfi[1];
 475       xy[0] = pyf[1];
 476       xfi[0] = pyf[1];
 477       //
 478       Double_t dy = params[i1][0]-params[i0][0];
 479       Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
 480       Double_t dphi = params[i1][1]-params[i0][1];
 481       Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
 482       //
 483       Int_t sector = isec-36;
 484       //
 485       if (TMath::Abs(dy/(sy+0.1))>5.) continue;          // 5 sigma cut
 486       if (TMath::Abs(dphi/(sphi+0.004))>5.) continue;    // 5 sigma cut
 487       TLinearFitter * fitterY,*fitterPhi;
 488       fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
 489       if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
 490       fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
 491       if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
 492       //
 493       fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
 494       if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
 495       fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
 496       if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
 497     }
 498   }
 499   //
 500   // Dump debug information
 501   //
 502   if (fStreamLevel>0){
 503     TTreeSRedirector *cstream = GetDebugStreamer();
 504     Int_t sector = isec-36;
 505     if (cstream){
 506       for (Int_t i0=0;i0<4;i0++){
 507         for (Int_t i1=i0;i1<4;i1++){
 508           if (i0==i1) continue;
 509           if(npoints[i0]<kMinClusterQ) continue;
 510           if(npoints[i1]<kMinClusterQ) continue;
 511           (*cstream)<<"fitQ"<<
 512             "run="<<fRun<<              //  run number
 513             "event="<<fEvent<<          //  event number
 514             "time="<<fTime<<            //  time stamp of event
 515             "trigger="<<fTrigger<<      //  trigger
 516             "mag="<<fMagF<<             //  magnetic field
 517             "sec="<<sector<<            // current sector from 0 to 36
 518             "isec="<<isec<<             //  current sector
 519             // full fit
 520             "nf="<<nf<<                 //  total number of points
 521             "pyf.="<<&pyf<<             //  full OROC fit y
 522             "pzf.="<<&pzf<<             //  full OROC fit z
 523             // quadrant fit
 524             "i0="<<i0<<                 // quadrant number
 525             "i1="<<i1<<
 526             "n0="<<npoints[i0]<<        // number of points
 527             "n1="<<npoints[i1]<<
 528             "p0.="<<&params[i0]<<       // parameters
 529             "p1.="<<&params[i1]<<
 530             "e0.="<<&errors[i0]<<       // errors
 531             "e1.="<<&errors[i1]<<
 532             "chi0="<<chi2C[i0]<<       // chi2s
 533             "chi1="<<chi2C[i1]<<
 534
 535             "\n";
 536         }
 537       }
 538     }
 539   }
 540 }
 541
 542
 543
 544
 545
 546 void AliTPCcalibUnlinearity::MakeHisto(){
 547   //
 548   //
 549   //
 550   TString  axisName[10];
 551   Double_t xmin[10],  xmax[10];
 552   Int_t    nbins[10];
 553   //
 554   //
 555   axisName[0]="sector";
 556   xmin[0]=0; xmax[0]=72; nbins[0]=72;
 557   //
 558   axisName[1]="row";
 559   xmin[1]=0; xmax[1]=96; nbins[1]=96;
 560   //
 561   axisName[2]="tphi";            // tan phi - ly/lx
 562   xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
 563   //
 564   axisName[3]="lz";              // global z
 565   xmin[3]=-250; xmax[3]=250; nbins[3]=10;
 566   //
 567   axisName[4]="k";              // tangent - in angle
 568   xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
 569   //
 570   //
 571   axisName[5]="delta";          // delta
 572   xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
 573   //
 574   //
 575   fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
 576   fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
 577   fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
 578   fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
 579  //  for (Int_t i=0; i<8;i++){
 580 //     fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
 581 //     fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
 582 //     fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
 583 //     fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
 584 //     fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
 585 //     fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
 586 //     fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
 587 //     fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
 588 //   }
 589
 590   //
 591   //
 592   //
 593   char hname[1000];
 594   TH2F * his=0;
 595   for (Int_t isec=0;isec<74;isec++){
 596     sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
 597     his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
 598     his->SetDirectory(0);
 599     fDistRPHIPlus.AddAt(his,isec);
 600     sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
 601     his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
 602     his->SetDirectory(0);
 603     fDistRPHIMinus.AddAt(his,isec);
 604   }
 605 }
 606
 607
 608 void AliTPCcalibUnlinearity::Terminate(){
 609   //
 610   // Terminate function
 611   // call base terminate + Eval of fitters
 612   //
 613   Info("AliTPCcalibUnlinearities","Terminate");
 614   EvalFitters();
 615   DumpTree();
 616   AliTPCcalibBase::Terminate();
 617 }
 618
 619 Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
 620   //
 621   // merge results
 622   //
 623   TIterator* iter = list->MakeIterator();
 624   AliTPCcalibUnlinearity* cal = 0;
 625   //
 626   while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
 627     if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
 628       return -1;
 629     }
 630     Add(cal);
 631   }
 632   return 0;
 633 }
 634
 635 void    AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
 636   //
 637   //
 638   //
 639   if (!fInit) Init();
 640   if (fDiffHistoLineY && calib->fDiffHistoLineY){
 641     fDiffHistoLineY->Add(calib->fDiffHistoLineY);
 642   }
 643   if (fDiffHistoParY && calib->fDiffHistoParY){
 644     fDiffHistoParY->Add(calib->fDiffHistoParY);
 645   }
 646   if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
 647     fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
 648   }
 649   if (fDiffHistoParZ && calib->fDiffHistoParZ){
 650     fDiffHistoParZ->Add(calib->fDiffHistoParZ);
 651   }
 652   for (Int_t isec=0;isec<38;isec++){
 653     TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
 654     TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
 655     if (f0r&&f1r) f0r->Add(f1r);
 656     TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
 657     TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
 658     if (f0z&&f1z) f0z->Add(f1z);
 659   }
 660
 661   for (Int_t isec=0;isec<16*38;isec++){
 662     TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
 663     TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
 664     if (f0y&&f1y) f0y->Add(f1y);
 665     TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
 666     TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
 667     if (f0phi&&f1phi) f0phi->Add(f1phi);
 668   }
 669
 670
 671   for (Int_t isec=0;isec<74;isec++){
 672     TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
 673     TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
 674     if (h0p&&h1p) h0p->Add(h1p);
 675     TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
 676     TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
 677     if (h0m&&h1m) h0m->Add(h1m);
 678   }
 679
 680
 681 }
 682
 683
 684
 685 void AliTPCcalibUnlinearity::DumpTree(const char *fname){
 686   //
 687   //
 688   //
 689   TTreeSRedirector *cstream =new TTreeSRedirector(fname);
 690   if (!cstream) return;
 691   //
 692   THnSparse *his=0;
 693   Double_t position[10];
 694   Double_t value;
 695   Int_t *bins = new Int_t[10];
 696   //
 697   for (Int_t ihis=0; ihis<2; ihis++){
 698     his =  (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
 699     if (!his) continue;
 700     //
 701     Int_t ndim = his->GetNdimensions();
 702     //
 703     for (Long64_t i = 0; i < his->GetNbins(); ++i) {
 704       value = his->GetBinContent(i, bins);
 705       for (Int_t idim = 0; idim < ndim; idim++) {
 706         position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
 707       }
 708       (*cstream)<<"Resol"<<
 709         "ihis="<<ihis<<
 710         "bincont="<<value;
 711       for (Int_t idim=0;idim<ndim;idim++){
 712         (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
 713       }
 714       (*cstream)<<"Resol"<<
 715         "\n";
 716     }
 717   }
 718   delete cstream;
 719 }
 720
 721
 722 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t ky, Double_t kz, Int_t npoints){
 723   //
 724   //
 725   // Simple distortion fit in outer sectors
 726   //
 727   // Model - 2 exponential decay toward the center of chamber
 728   //       - Decay length 10 and 5 cm
 729   //       - Decay amplitude - Z dependent
 730   //
 731
 732   Double_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
 733   Double_t dzt    = (cz-tz)*side;
 734   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
 735   AliTPCROC * roc = AliTPCROC::Instance();
 736   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
 737   Double_t dout = xout-radius;
 738   if (dout>30) return;
 739   //drift
 740   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
 741   Double_t dy = cy-ty;
 742   Double_t eout10 = TMath::Exp(-dout/10.);
 743   Double_t eout5 = TMath::Exp(-dout/5.);
 744   //
 745   Double_t eoutp  = (eout10+eout5)*0.5;
 746   Double_t eoutm  = (eout10-eout5)*0.5;
 747
 748   //
 749   Double_t xxxR[6], xxxZ[6], xxxRZ[6];
 750   //TString fstring="";
 751   xxxZ[0]=eoutp;                //fstring+="eoutp++";
 752   xxxZ[1]=eoutm;               //fstring+="eoutm++";
 753   xxxZ[2]=dr*eoutp;            //fstring+="dr*eoutp++";
 754   xxxZ[3]=dr*eoutm;            //fstring+="dr*eoutm++";
 755   xxxZ[4]=dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 756   xxxZ[5]=dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 757   //
 758   xxxR[0]=ky*eoutp;                //fstring+="eoutp++";
 759   xxxR[1]=ky*eoutm;               //fstring+="eoutm++";
 760   xxxR[2]=ky*dr*eoutp;            //fstring+="dr*eoutp++";
 761   xxxR[3]=ky*dr*eoutm;            //fstring+="dr*eoutm++";
 762   xxxR[4]=ky*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 763   xxxR[5]=ky*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 764   //
 765   xxxRZ[0]=kz*eoutp;                //fstring+="eoutp++";
 766   xxxRZ[1]=kz*eoutm;               //fstring+="eoutm++";
 767   xxxRZ[2]=kz*dr*eoutp;            //fstring+="dr*eoutp++";
 768   xxxRZ[3]=kz*dr*eoutm;            //fstring+="dr*eoutm++";
 769   xxxRZ[4]=kz*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 770   xxxRZ[5]=kz*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 771   //
 772   TLinearFitter * fitter=0;
 773   Double_t err=0.1;
 774   for (Int_t ip=0; ip<npoints; ip++){
 775     //
 776     fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
 777     fitter->AddPoint(xxxR,dy,err);
 778     //fitter->AddPoint(xxxRZ,dz);
 779     fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
 780     fitter->AddPoint(xxxZ,dzt,err);
 781     //
 782     if (side>0){
 783       fitter = (TLinearFitter*)fFittersOutR.At(36);
 784       fitter->AddPoint(xxxR,dy,err);
 785       //fitter->AddPoint(xxxRZ,dz);
 786       fitter = (TLinearFitter*)fFittersOutZ.At(36);
 787       fitter->AddPoint(xxxZ,dzt,err);
 788     }
 789     if (side<0){
 790       fitter = (TLinearFitter*)fFittersOutR.At(37);
 791       fitter->AddPoint(xxxR,dy,err);
 792       //fitter->AddPoint(xxxRZ,dz);
 793       fitter = (TLinearFitter*)fFittersOutZ.At(37);
 794       fitter->AddPoint(xxxZ,dzt,err);
 795     }
 796   }
 797 }
 798 void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
 799   //
 800   //
 801   //
 802   Float_t kMaxDz = 0.5;  // cut on z in cm
 803   Double_t dy = cy-ty;
 804   Double_t dz = cz-tz;
 805   if (TMath::Abs(dz)>kMaxDz) return;
 806   Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
 807   Double_t edgePlus =  ty+cx*TMath::Tan(TMath::Pi()/18.);
 808   //
 809   TH2F * his =0;
 810   his = (TH2F*)fDistRPHIPlus.At(sector);
 811   his->Fill(edgePlus,dy,npoints);
 812   his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
 813   his->Fill(edgePlus,dy,npoints);
 814   his = (TH2F*)fDistRPHIMinus.At(sector);
 815   his->Fill(edgeMinus,dy,npoints);
 816   his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
 817   his->Fill(edgeMinus,dy,npoints);
 818 }
 819
 820
 821
 822 void AliTPCcalibUnlinearity::Init(){
 823   //
 824   //
 825   //
 826   //
 827   MakeHisto();
 828   // Make outer fitters
 829   TLinearFitter * fitter = 0;
 830   for (Int_t ifit=0; ifit<38; ifit++){
 831     fitter = new TLinearFitter(7,"hyp6");
 832     fitter->StoreData(kFALSE);
 833     fFittersOutR.AddAt(fitter,ifit);
 834     fitter = new TLinearFitter(7,"hyp6");
 835     fitter->StoreData(kFALSE);
 836     fFittersOutZ.AddAt(fitter,ifit);
 837   }
 838   for (Int_t ifit=0; ifit<16*38;ifit++){
 839     fitter = new TLinearFitter(2,"hyp1");
 840     fitter->StoreData(kFALSE);
 841     fFitterQuadrantY.AddAt(fitter,ifit);
 842     fitter = new TLinearFitter(2,"hyp1");
 843     fitter->StoreData(kFALSE);
 844     fFitterQuadrantPhi.AddAt(fitter,ifit);
 845   }
 846   fInit= kTRUE;
 847 }
 848
 849 void AliTPCcalibUnlinearity::EvalFitters(){
 850   //
 851   //
 852   // Evaluate fitters
 853   //
 854   TLinearFitter * fitter = 0;
 855   TVectorD vec;
 856   for (Int_t ifit=0; ifit<38; ifit++){
 857     fitter=(TLinearFitter*)fFittersOutR.At(ifit);
 858     if (fitter) {
 859       fitter->Eval();
 860       fitter->GetParameters(vec);
 861       fParamsOutR.AddAt(vec.Clone(),ifit);
 862       fitter->GetErrors(vec);
 863       fErrorsOutR.AddAt(vec.Clone(),ifit);
 864     }
 865     fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
 866     if (fitter) {
 867       fitter->Eval();
 868       fitter->GetParameters(vec);
 869       fParamsOutZ.AddAt(vec.Clone(),ifit);
 870       fitter->GetErrors(vec);
 871       fErrorsOutZ.AddAt(vec.Clone(),ifit);
 872     }
 873   }
 874 }
 875
 876
 877
 878
 879
 880 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
 881   //
 882   //
 883   //
 884   //  TTree * tree = chainUnlinP;
 885   AliTPCcalibUnlinearity *calib = this;
 886   //
 887   AliTPCclusterMI * cl=0;
 888   Double_t ty,tz;
 889   TVectorD *vy=0, *vz=0;
 890   TVectorD *vy2=0, *vz2=0;
 891   Long64_t nentries = tree->GetEntries();
 892   if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
 893   //
 894   {
 895   for (Long64_t i=0; i<nentries; i++){
 896     tree->SetBranchAddress("Cl.",&cl);
 897     tree->SetBranchAddress("y1",&ty);
 898     tree->SetBranchAddress("z1",&tz);
 899     tree->SetBranchAddress("py1.",&vy);
 900     tree->SetBranchAddress("pz1.",&vz);
 901     tree->SetBranchAddress("py2.",&vy2);
 902     tree->SetBranchAddress("pz2.",&vz2);
 903     //
 904     tree->GetEntry(i);
 905     if (!cl) continue;
 906     if (i%10000==0) printf("%d\n",(Int_t)i);
 907     Int_t row= cl->GetRow();
 908     if (cl->GetDetector()>36) row+=64;
 909     calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
 910                       ty,tz,(*vy)[1],(*vz)[1],1);
 911
 912     if (cl->GetType()<0) continue;
 913     Double_t dy = cl->GetY()-ty;
 914     Double_t dz = cl->GetZ()-tz;
 915     //
 916     //
 917     if (TMath::Abs(dy)>0.25) continue;
 918     if (TMath::Abs(dz)>0.25) continue;
 919
 920     if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
 921     if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
 922       // Apply sagita cut
 923     if (cl->GetType()>=0) {
 924       calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
 925                       ty,tz,(*vy)[1],(*vz)[1],1);
 926     }
 927   }
 928   }
 929 }
 930
 931
 932 void  AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
 933   //
 934   // Make position corrections
 935   // for the moment Only using debug streamer
 936   // chainUnlinD  - debug tree
 937   // param     - parameters to be updated
 938   // maxPoints - maximal number of points using for fit
 939   // verbose   - print info flag
 940   //
 941   // Current implementation - only using debug streamers
 942   //
 943
 944   /*
 945   //Defaults
 946   Int_t maxPoints=100000;
 947   */
 948   /*
 949     Usage:
 950     //0. Load libraries
 951     gSystem->Load("libANALYSIS");
 952     gSystem->Load("libSTAT");
 953     gSystem->Load("libTPCcalib");
 954
 955
 956     //1. Load Parameters to be modified:
 957     //e.g:
 958
 959     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 960     AliCDBManager::Instance()->SetRun(0)
 961     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
 962     //
 963     //2. Load the debug streamer
 964     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 965     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 966     AliXRDPROOFtoolkit tool;
 967     TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
 968     chainUnlin->Lookup();
 969     TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
 970     chainUnlinD->Lookup();
 971
 972     //3. Do fits and store results
 973     //
 974     AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
 975     TFile f("paramout.root","recreate");
 976     param->Write("clusterParam");
 977     f.Close();
 978     //4. Verification
 979     TFile f2("paramout.root");
 980     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
 981     param2->SetInstance(param2);
 982     chainUnlinD->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
 983
 984    */
 985
 986
 987   TStatToolkit toolkit;
 988   Double_t chi2;
 989   TVectorD fitParamY0;
 990   TVectorD fitParamY1;
 991   TVectorD fitParamZ0;
 992   TVectorD fitParamZ1;
 993   TMatrixD covMatrix;
 994   Int_t npoints;
 995
 996   chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
 997   chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
 998   chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
 999   chainUnlinD->SetAlias("st","(sin(dt)-dt)");
1000   chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
1001   //
1002   chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1003   //
1004   //
1005   //
1006   TCut cutE("Cl.fType>=0");
1007   TCut cutDy("abs(Cl.fY-y1)<0.4");
1008   TCut cutDz("abs(Cl.fZ-z1)<0.4");
1009   TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
1010   TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
1011   TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
1012   //
1013   TCut cutA("Cl.fZ>0");
1014   TCut cutC("Cl.fZ<0");
1015
1016   TString fstringY="";
1017   //
1018   fstringY+="(dp)++";            //1
1019   fstringY+="(dp)*di++";         //2
1020   fstringY+="(sp)++";            //3
1021   fstringY+="(sp)*di++";         //4
1022   fstringY+="(dq)++";            //5
1023   TString fstringZ="";
1024   fstringZ+="(dt)++";            //1
1025   fstringZ+="(dt)*di++";         //2
1026   fstringZ+="(st)++";            //3
1027   fstringZ+="(st)*di++";         //4
1028   fstringZ+="(dq)++";            //5
1029   //
1030   // Z corrections
1031   //
1032   TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
1033   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1034   strZ0->Tokenize("++")->Print();
1035   param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
1036   //
1037   TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
1038   //
1039   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1040   strZ1->Tokenize("++")->Print();
1041   param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
1042   param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
1043   //
1044   // Y corrections
1045   //
1046   TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
1047   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1048   strY0->Tokenize("++")->Print();
1049   param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
1050
1051
1052   TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
1053   //
1054   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1055   strY1->Tokenize("++")->Print();
1056   param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
1057   param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
1058   //
1059   //
1060   //
1061   chainUnlinD->SetAlias("fitZ0",strZ0->Data());
1062   chainUnlinD->SetAlias("fitZ1",strZ1->Data());
1063   chainUnlinD->SetAlias("fitY0",strY0->Data());
1064   chainUnlinD->SetAlias("fitY1",strY1->Data());
1065 }
1066
1067
1068
1069
1070
1071
1072