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
  36   //
  37 */
  38
  39
  40 #include "TLinearFitter.h"
  41
  42 #include "Riostream.h"
  43 #include "TChain.h"
  44 #include "TTree.h"
  45 #include "TH1F.h"
  46 #include "TH2F.h"
  47 #include "TH3F.h"
  48 #include "THnSparse.h"
  49 #include "TList.h"
  50 #include "TMath.h"
  51 #include "TCanvas.h"
  52 #include "TFile.h"
  53 #include "TF1.h"
  54 #include "TVectorD.h"
  55 #include "TProfile.h"
  56 #include "TGraphErrors.h"
  57 #include "TCanvas.h"
  58 #include "TCut.h"
  59
  60
  61 #include "AliTPCclusterMI.h"
  62 #include "AliTPCseed.h"
  63 #include "AliESDVertex.h"
  64 #include "AliESDEvent.h"
  65 #include "AliESDfriend.h"
  66 #include "AliESDInputHandler.h"
  67 #include "AliAnalysisManager.h"
  68 #include "AliTracker.h"
  69 #include "AliMagF.h"
  70 #include "AliTPCCalROC.h"
  71
  72 #include "AliLog.h"
  73
  74
  75 #include "TTreeStream.h"
  76 #include "AliTPCTracklet.h"
  77 #include "TTimeStamp.h"
  78 #include "AliTPCcalibDB.h"
  79 #include "AliTPCcalibLaser.h"
  80 #include "AliDCSSensorArray.h"
  81 #include "AliDCSSensor.h"
  82 #include "TLinearFitter.h"
  83 #include "AliTPCClusterParam.h"
  84 #include "TStatToolkit.h"
  85
  86
  87 #include "AliTPCcalibUnlinearity.h"
  88
  89 AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
  90
  91 ClassImp(AliTPCcalibUnlinearity)
  92
  93 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
  94   AliTPCcalibBase(),
  95   fDiffHistoLine(0),
  96   fDiffHistoPar(0),
  97   fFittersOutR(38),
  98   fFittersOutZ(38),
  99   fParamsOutR(38),
 100   fParamsOutZ(38),
 101   fErrorsOutR(38),
 102   fErrorsOutZ(38)
 103 {
 104   //
 105   // Defualt constructor
 106   //
 107 }
 108
 109 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
 110   AliTPCcalibBase(name,title),
 111   fDiffHistoLine(0),
 112   fDiffHistoPar(0),
 113   fFittersOutR(38),
 114   fFittersOutZ(38),
 115   fParamsOutR(38),
 116   fParamsOutZ(38),
 117   fErrorsOutR(38),
 118   fErrorsOutZ(38)
 119 {
 120   //
 121   // Non default constructor
 122   //
 123   MakeFitters();
 124 }
 125
 126 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
 127   //
 128   //
 129   //
 130   if (fDiffHistoLine) delete fDiffHistoLine;
 131   if (fDiffHistoPar)  delete fDiffHistoPar;
 132 }
 133
 134
 135 AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
 136 {
 137   //
 138   // Singleton implementation
 139   // Returns an instance of this class, it is created if neccessary
 140   //
 141   if (fgInstance == 0){
 142     fgInstance = new AliTPCcalibUnlinearity();
 143   }
 144   return fgInstance;
 145 }
 146
 147
 148
 149
 150 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
 151   //
 152   //
 153   //
 154   const Int_t  kdrow=10;
 155   const Int_t kMinCluster=35;
 156   if (!fDiffHistoLine) MakeHisto();
 157   if (!seed) return;
 158   if (TMath::Abs(fMagF)>0.01) return;   // use only data without field
 159   Int_t counter[72];
 160   for (Int_t i=0;i<72;i++) counter[i]=0;
 161   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 162     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 163     if (!cluster) continue;
 164     if (cluster->GetType()<0) continue;
 165     counter[cluster->GetDetector()]++;
 166   }
 167   //
 168   for (Int_t is=0; is<72;is++){
 169     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
 170   }
 171 }
 172
 173
 174 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
 175   //
 176   // process differnce of the cluster in respect with linear and parabolic fit
 177   //
 178   const  Double_t kXmean=134;
 179   const Int_t     kdrow=10;
 180   static TLinearFitter fy1(2,"pol1");
 181   static TLinearFitter fy2(3,"pol2");
 182   static TLinearFitter fz1(2,"pol1");
 183   static TLinearFitter fz2(3,"pol2");
 184   //
 185   static TVectorD py1(2);
 186   static TVectorD py2(3);
 187   static TVectorD pz1(2);
 188   static TVectorD pz2(3);
 189   //
 190   fy1.ClearPoints();
 191   fy2.ClearPoints();
 192   fz1.ClearPoints();
 193   fz2.ClearPoints();
 194   //
 195   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 196     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 197     if (!c) continue;
 198     if (c->GetDetector()!=isec) continue;
 199     if (c->GetType()<0) continue;
 200     Double_t dx = c->GetX()-kXmean;
 201     Double_t x[2]={dx, dx*dx};
 202     fy1.AddPoint(x,c->GetY());
 203     fy2.AddPoint(x,c->GetY());
 204     fz1.AddPoint(x,c->GetZ());
 205     fz2.AddPoint(x,c->GetZ());
 206   }
 207   fy1.Eval(); fy1.GetParameters(py1);
 208   fy2.Eval(); fy2.GetParameters(py2);
 209   fz1.Eval(); fz1.GetParameters(pz1);
 210   fz2.Eval(); fz2.GetParameters(pz2);
 211   //
 212
 213   for (Int_t irow=0;irow<159;irow++) {
 214     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 215     if (!c) continue;
 216     if (c->GetDetector()!=isec) continue;
 217     if (c->GetType()<0) continue;
 218     Double_t dx = c->GetX()-kXmean;
 219     Double_t y1 = py1[0]+py1[1]*dx;
 220     Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
 221     //
 222     Double_t z1 = pz1[0]+pz1[1]*dx;
 223     Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
 224     //
 225     //
 226     Double_t x[10];
 227     x[0]=isec;
 228     x[1]=irow;
 229     x[2]=c->GetY();
 230     x[3]=c->GetZ();
 231     x[4]=py1[1];
 232     x[5]=pz1[1];
 233     x[6]=py2[2]*150*150/4; // sagita 150 cm
 234     //
 235     x[7]=c->GetY()-y1;
 236     x[8]=c->GetZ()-z1;
 237     fDiffHistoLine->Fill(x);
 238     x[7]=c->GetY()-y2;
 239     x[8]=c->GetZ()-z2;
 240     fDiffHistoPar->Fill(x);
 241
 242     if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
 243       // Apply sagita cut
 244       AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
 245                py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
 246     }
 247
 248     if (fStreamLevel>1){
 249       TTreeSRedirector *cstream = GetDebugStreamer();
 250       if (cstream){
 251         (*cstream)<<"Diff"<<
 252           "run="<<fRun<<              //  run number
 253           "event="<<fEvent<<          //  event number
 254           "time="<<fTime<<            //  time stamp of event
 255           "trigger="<<fTrigger<<      //  trigger
 256           "mag="<<fMagF<<             //  magnetic field
 257           "isec="<<isec<<
 258           "Cl.="<<c<<
 259           "y1="<<y1<<
 260           "y2="<<y2<<
 261           "z1="<<z1<<
 262           "z2="<<z2<<
 263           "py1.="<<&py1<<
 264           "py2.="<<&py2<<
 265           "pz1.="<<&pz1<<
 266           "pz2.="<<&pz2<<
 267           "\n";
 268       }
 269     }
 270   }
 271 }
 272
 273
 274 void AliTPCcalibUnlinearity::MakeHisto(){
 275   //
 276   //
 277   //
 278   TString  axisName[10];
 279   Double_t xmin[10],  xmax[10];
 280   Int_t    nbins[10];
 281   //
 282   //
 283   axisName[0]="sector";
 284   xmin[0]=0; xmax[0]=72; nbins[0]=72;
 285   //
 286   axisName[1]="row";
 287   xmin[1]=0; xmax[1]=159; nbins[1]=159;
 288   //
 289   axisName[2]="ly";
 290   xmin[2]=-50; xmax[2]=50; nbins[2]=10;
 291   //
 292   axisName[3]="lz";
 293   xmin[3]=-250; xmax[3]=250; nbins[3]=50;
 294   //
 295   axisName[4]="p2";
 296   xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
 297   //
 298   axisName[5]="p3";
 299   xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
 300   //
 301   axisName[6]="p4";
 302   xmin[6]=-2; xmax[6]=2; nbins[6]=20;
 303   //
 304   axisName[7]="dy";
 305   xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
 306   //
 307   axisName[8]="dz";
 308   xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
 309   //
 310   fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
 311   fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
 312   for (Int_t i=0; i<9;i++){
 313     fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
 314     fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
 315     fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
 316     fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
 317   }
 318 }
 319
 320
 321 void AliTPCcalibUnlinearity::Terminate(){
 322   //
 323   // Terminate function
 324   // call base terminate + Eval of fitters
 325   //
 326   Info("AliTPCcalibUnlinearities","Terminate");
 327   EvalFitters();
 328   DumpTree();
 329   AliTPCcalibBase::Terminate();
 330 }
 331
 332
 333 void AliTPCcalibUnlinearity::DumpTree(){
 334   //
 335   //
 336   //
 337   if (fStreamLevel==0) return;
 338   TTreeSRedirector *cstream = GetDebugStreamer();
 339   if (!cstream) return;
 340   //
 341   THnSparse *his=0;
 342   Double_t position[10];
 343   Double_t value;
 344   Int_t *bins = new Int_t[10];
 345   //
 346   for (Int_t ihis=0; ihis<2; ihis++){
 347     his =  (ihis==0)? fDiffHistoLine:fDiffHistoPar;
 348     if (!his) continue;
 349     //
 350     Int_t ndim = his->GetNdimensions();
 351     //
 352     for (Long64_t i = 0; i < his->GetNbins(); ++i) {
 353       value = his->GetBinContent(i, bins);
 354       for (Int_t idim = 0; idim < ndim; idim++) {
 355         position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
 356       }
 357       (*cstream)<<"Resol"<<
 358         "ihis="<<ihis<<
 359         "bincont="<<value;
 360       for (Int_t idim=0;idim<ndim;idim++){
 361         (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
 362       }
 363       (*cstream)<<"Resol"<<
 364         "\n";
 365     }
 366   }
 367 }
 368
 369
 370 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
 371   //
 372   //
 373   // Simple distortion fit in outer sectors
 374   //
 375   // Model - 2 exponential decay toward the center of chamber
 376   //       - Decay length 10 and 5 cm
 377   //       - Decay amplitude - Z dependent
 378   //
 379   /*
 380   chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
 381   chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
 382   chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
 383   chainUnlin->SetAlias("dzt","dz*side");
 384   chainUnlin->SetAlias("dr","(1-abs(lz/250))");
 385   chainUnlin->SetAlias("dout","(159-row)*1.5");
 386   chainUnlin->SetAlias("din","row*0.75");
 387   chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
 388   chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");
 389   */
 390
 391   Double_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
 392   Double_t dzt    = dz*side;
 393   Double_t dout   = (159-row)*1.5;  // distance to the last pad row - (valid only for long pads)
 394   if (dout>30) return; // use only edge pad rows
 395   dr-=0.5;             // dr shifted to the middle - reduce numerical instabilities
 396
 397   Double_t eout10 = TMath::Exp(-dout/10.);
 398   Double_t eout5 = TMath::Exp(-dout/5.);
 399   //
 400   Double_t eoutp  = (eout10+eout5)*0.5;
 401   Double_t eoutm  = (eout10-eout5)*0.5;
 402
 403   //
 404   Double_t xxxR[6], xxxZ[6], xxxRZ[6];
 405   //TString fstring="";
 406   xxxZ[0]=eoutp;                //fstring+="eoutp++";
 407   xxxZ[1]=eoutm;               //fstring+="eoutm++";
 408   xxxZ[2]=dr*eoutp;            //fstring+="dr*eoutp++";
 409   xxxZ[3]=dr*eoutm;            //fstring+="dr*eoutm++";
 410   xxxZ[4]=dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 411   xxxZ[5]=dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 412   //
 413   xxxR[0]=p2*eoutp;                //fstring+="eoutp++";
 414   xxxR[1]=p2*eoutm;               //fstring+="eoutm++";
 415   xxxR[2]=p2*dr*eoutp;            //fstring+="dr*eoutp++";
 416   xxxR[3]=p2*dr*eoutm;            //fstring+="dr*eoutm++";
 417   xxxR[4]=p2*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 418   xxxR[5]=p2*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 419   //
 420   xxxRZ[0]=p3*eoutp;                //fstring+="eoutp++";
 421   xxxRZ[1]=p3*eoutm;               //fstring+="eoutm++";
 422   xxxRZ[2]=p3*dr*eoutp;            //fstring+="dr*eoutp++";
 423   xxxRZ[3]=p3*dr*eoutm;            //fstring+="dr*eoutm++";
 424   xxxRZ[4]=p3*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";
 425   xxxRZ[5]=p3*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";
 426   //
 427   TLinearFitter * fitter=0;
 428   Double_t err=0.1;
 429   for (Int_t ip=0; ip<npoints; ip++){
 430     //
 431     fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
 432     fitter->AddPoint(xxxR,dy,err);
 433     //fitter->AddPoint(xxxRZ,dz);
 434     fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
 435     fitter->AddPoint(xxxZ,dzt,err);
 436     //
 437     if (side>0){
 438       fitter = (TLinearFitter*)fFittersOutR.At(36);
 439       fitter->AddPoint(xxxR,dy,err);
 440       //fitter->AddPoint(xxxRZ,dz);
 441       fitter = (TLinearFitter*)fFittersOutZ.At(36);
 442       fitter->AddPoint(xxxZ,dzt,err);
 443     }
 444     if (side<0){
 445       fitter = (TLinearFitter*)fFittersOutR.At(37);
 446       fitter->AddPoint(xxxR,dy,err);
 447       //fitter->AddPoint(xxxRZ,dz);
 448       fitter = (TLinearFitter*)fFittersOutZ.At(37);
 449       fitter->AddPoint(xxxZ,dzt,err);
 450     }
 451   }
 452 }
 453
 454
 455 void AliTPCcalibUnlinearity::MakeFitters(){
 456   //
 457   //
 458   //
 459   // Make outer fitters
 460   TLinearFitter * fitter = 0;
 461   for (Int_t ifit=0; ifit<38; ifit++){
 462     fitter = new TLinearFitter(7,"hyp6");
 463     fitter->StoreData(kFALSE);
 464     fFittersOutR.AddAt(fitter,ifit);
 465     fitter = new TLinearFitter(7,"hyp6");
 466     fitter->StoreData(kFALSE);
 467     fFittersOutZ.AddAt(fitter,ifit);
 468   }
 469 }
 470
 471 void AliTPCcalibUnlinearity::EvalFitters(){
 472   //
 473   //
 474   // Evaluate fitters
 475   //
 476   TLinearFitter * fitter = 0;
 477   TVectorD vec;
 478   for (Int_t ifit=0; ifit<38; ifit++){
 479     fitter=(TLinearFitter*)fFittersOutR.At(ifit);
 480     if (fitter) {
 481       fitter->Eval();
 482       fitter->GetParameters(vec);
 483       fParamsOutR.AddAt(vec.Clone(),ifit);
 484       fitter->GetErrors(vec);
 485       fErrorsOutR.AddAt(vec.Clone(),ifit);
 486     }
 487     fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
 488     if (fitter) {
 489       fitter->Eval();
 490       fitter->GetParameters(vec);
 491       fParamsOutZ.AddAt(vec.Clone(),ifit);
 492       fitter->GetErrors(vec);
 493       fErrorsOutZ.AddAt(vec.Clone(),ifit);
 494     }
 495   }
 496 }
 497
 498 Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
 499   //
 500   //
 501   //
 502   TVectorD * vec = GetParamOutR(sector);
 503   if (!vec) return 0;
 504   dr-=0.5;        // dr=0 at the middle drift length
 505   Double_t eout10 = TMath::Exp(-dout/10.);
 506   Double_t eout5 = TMath::Exp(-dout/5.);
 507   Double_t eoutp  = (eout10+eout5)*0.5;
 508   Double_t eoutm  = (eout10-eout5)*0.5;
 509
 510   Double_t result=0;
 511   result+=(*vec)[1]*eoutp;
 512   result+=(*vec)[2]*eoutm;
 513   result+=(*vec)[3]*eoutp*dr;
 514   result+=(*vec)[4]*eoutm*dr;
 515   result+=(*vec)[5]*eoutp*dr*dr;
 516   result+=(*vec)[6]*eoutm*dr*dr;
 517   return result;
 518 }
 519
 520
 521 Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
 522   //
 523   //
 524   //
 525   Double_t alpha    = TMath::ATan2(gy,gx);
 526   if (alpha<0)  alpha+=TMath::Pi()*2;
 527   Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
 528   alpha = (lsector+0.5)*TMath::Pi()/9.;
 529   //
 530   Double_t r[3];
 531   r[0]=gx; r[1]=gy;
 532   Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
 533   r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
 534   //
 535   AliTPCROC * roc = AliTPCROC::Instance();
 536   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
 537   Double_t dout = xout-r[0];
 538   if (dout<0) return 0;
 539   Double_t dr   = 1-TMath::Abs(gz/250.);
 540   if (gz<0) lsector+=18;
 541   if (stype==0) lsector = (gz>0) ? 36:37;
 542   if (stype<0) return lsector;  //
 543   return GetDr(lsector,dout,dr);
 544 }
 545
 546
 547 Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
 548   //
 549   //
 550   //
 551   TVectorD * vec = GetParamOutZ(sector);
 552   if (!vec) return 0;
 553   dr-=0.5; // 0 at the middle
 554   Double_t eout10 = TMath::Exp(-dout/10.);
 555   Double_t eout5 = TMath::Exp(-dout/5.);
 556   Double_t eoutp  = (eout10+eout5)*0.5;
 557   Double_t eoutm  = (eout10-eout5)*0.5;
 558
 559
 560   Double_t result=0;
 561   result+=(*vec)[1]*eoutp;
 562   result+=(*vec)[2]*eoutm;
 563   result+=(*vec)[3]*eoutp*dr;
 564   result+=(*vec)[4]*eoutm*dr;
 565   result+=(*vec)[5]*eoutp*dr*dr;
 566   result+=(*vec)[6]*eoutm*dr*dr;
 567   return result;
 568 }
 569
 570
 571 Double_t      AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
 572   //
 573   //
 574   //
 575   return fgInstance->GetDr(sector,dout,dr);
 576 }
 577 Double_t      AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
 578   //
 579   //
 580   //
 581   return fgInstance->GetGDr(stype,gx,gy,gz);
 582 }
 583 Double_t      AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
 584   //
 585   //
 586   //
 587   return fgInstance->GetDz(sector,dout,dr);
 588 }
 589
 590
 591
 592
 593 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
 594   //
 595   //
 596   //
 597   //  TTree * tree = chainUnlinP;
 598   AliTPCcalibUnlinearity *calib = this;
 599   //
 600   AliTPCclusterMI * cl=0;
 601   Double_t ty,tz;
 602   TVectorD *vy=0, *vz=0;
 603   TVectorD *vy2=0, *vz2=0;
 604   Long64_t nentries = tree->GetEntries();
 605   if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
 606   //
 607   {
 608   for (Long64_t i=0; i<nentries; i++){
 609     tree->SetBranchAddress("Cl.",&cl);
 610     tree->SetBranchAddress("y1",&ty);
 611     tree->SetBranchAddress("z1",&tz);
 612     tree->SetBranchAddress("py1.",&vy);
 613     tree->SetBranchAddress("pz1.",&vz);
 614     tree->SetBranchAddress("py2.",&vy2);
 615     tree->SetBranchAddress("pz2.",&vz2);
 616     //
 617     tree->GetEntry(i);
 618     if (!cl) continue;
 619     if (i%10000==0) printf("%d\n",(Int_t)i);
 620     Int_t row= cl->GetRow();
 621     if (cl->GetDetector()>36) row+=64;
 622     if (cl->GetType()<0) continue;
 623     Double_t dy = cl->GetY()-ty;
 624     Double_t dz = cl->GetZ()-tz;
 625     Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
 626     //
 627     //
 628     if (TMath::Abs(dy)>0.25) continue;
 629     if (TMath::Abs(dz)>0.25) continue;
 630
 631     if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
 632     if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
 633       // Apply sagita cut
 634     calib->AddPoint(cl->GetDetector(), row, dz, dy,
 635                     (*vy)[1],(*vz)[1], dr, 1);
 636   }
 637   }
 638 }
 639
 640
 641 void  AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
 642   //
 643   // Make position corrections
 644   // for the moment Only using debug streamer
 645   // chainUnlinD  - debug tree
 646   // param     - parameters to be updated
 647   // maxPoints - maximal number of points using for fit
 648   // verbose   - print info flag
 649   //
 650   // Current implementation - only using debug streamers
 651   //
 652
 653   /*
 654   //Defaults
 655   Int_t maxPoints=100000;
 656   */
 657   /*
 658     Usage:
 659     //0. Load libraries
 660     gSystem->Load("libANALYSIS");
 661     gSystem->Load("libSTAT");
 662     gSystem->Load("libTPCcalib");
 663
 664
 665     //1. Load Parameters to be modified:
 666     //e.g:
 667
 668     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 669     AliCDBManager::Instance()->SetRun(0)
 670     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
 671     //
 672     //2. Load the debug streamer
 673     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 674     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 675     AliXRDPROOFtoolkit tool;
 676     TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
 677     chainUnlin->Lookup();
 678     TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
 679     chainUnlinD->Lookup();
 680
 681     //3. Do fits and store results
 682     //
 683     AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
 684     TFile f("paramout.root","recreate");
 685     param->Write("clusterParam");
 686     f.Close();
 687     //4. Verification
 688     TFile f2("paramout.root");
 689     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
 690     param2->SetInstance(param2);
 691     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
 692
 693    */
 694
 695
 696   TStatToolkit toolkit;
 697   Double_t chi2;
 698   TVectorD fitParamY0;
 699   TVectorD fitParamY1;
 700   TVectorD fitParamZ0;
 701   TVectorD fitParamZ1;
 702   TMatrixD covMatrix;
 703   Int_t npoints;
 704
 705   chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
 706   chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
 707   chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
 708   chainUnlinD->SetAlias("st","(sin(dt)-dt)");
 709   chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
 710   //
 711   chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
 712   //
 713   //
 714   //
 715   TCut cutE("Cl.fType>=0");
 716   TCut cutDy("abs(Cl.fY-y1)<0.4");
 717   TCut cutDz("abs(Cl.fZ-z1)<0.4");
 718   TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
 719   TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
 720   TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
 721   //
 722   TCut cutA("Cl.fZ>0");
 723   TCut cutC("Cl.fZ<0");
 724
 725   TString fstringY="";
 726   //
 727   fstringY+="(dp)++";            //1
 728   fstringY+="(dp)*di++";         //2
 729   fstringY+="(sp)++";            //3
 730   fstringY+="(sp)*di++";         //4
 731   fstringY+="(dq)++";            //5
 732   TString fstringZ="";
 733   fstringZ+="(dt)++";            //1
 734   fstringZ+="(dt)*di++";         //2
 735   fstringZ+="(st)++";            //3
 736   fstringZ+="(st)*di++";         //4
 737   fstringZ+="(dq)++";            //5
 738   //
 739   // Z corrections
 740   //
 741   TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
 742   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 743   strZ0->Tokenize("++")->Print();
 744   param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
 745   //
 746   TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
 747   //
 748   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 749   strZ1->Tokenize("++")->Print();
 750   param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
 751   param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
 752   //
 753   // Y corrections
 754   //
 755   TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
 756   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 757   strY0->Tokenize("++")->Print();
 758   param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
 759
 760
 761   TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
 762   //
 763   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
 764   strY1->Tokenize("++")->Print();
 765   param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
 766   param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
 767   //
 768   //
 769   //
 770   chainUnlinD->SetAlias("fitZ0",strZ0->Data());
 771   chainUnlinD->SetAlias("fitZ1",strZ1->Data());
 772   chainUnlinD->SetAlias("fitY0",strY0->Data());
 773   chainUnlinD->SetAlias("fitY1",strY1->Data());
 774 }
 775
 776
 777
 778
 779
 780
 781