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   Only Linear and parabolic fit used
  19   To be used for laser or for data without field
  20   //
  21 */
  22
  23
  24 #include "TLinearFitter.h"
  25
  26 #include "Riostream.h"
  27 #include "TChain.h"
  28 #include "TTree.h"
  29 #include "TH1F.h"
  30 #include "TH2F.h"
  31 #include "TH3F.h"
  32 #include "THnSparse.h"
  33 #include "TList.h"
  34 #include "TMath.h"
  35 #include "TCanvas.h"
  36 #include "TFile.h"
  37 #include "TF1.h"
  38 #include "TVectorD.h"
  39 #include "TProfile.h"
  40 #include "TGraphErrors.h"
  41 #include "TCanvas.h"
  42
  43 #include "AliTPCclusterMI.h"
  44 #include "AliTPCseed.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDfriend.h"
  48 #include "AliESDInputHandler.h"
  49 #include "AliAnalysisManager.h"
  50 #include "AliTracker.h"
  51 #include "AliMagFMaps.h"
  52 #include "AliTPCCalROC.h"
  53
  54 #include "AliLog.h"
  55
  56
  57 #include "TTreeStream.h"
  58 #include "AliTPCTracklet.h"
  59 #include "TTimeStamp.h"
  60 #include "AliTPCcalibDB.h"
  61 #include "AliTPCcalibLaser.h"
  62 #include "AliDCSSensorArray.h"
  63 #include "AliDCSSensor.h"
  64 #include "TLinearFitter.h"
  65
  66
  67 #include "AliTPCcalibUnlinearity.h"
  68
  69
  70 ClassImp(AliTPCcalibUnlinearity)
  71
  72 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
  73   AliTPCcalibBase(),
  74   fDiffHistoLine(0),
  75   fDiffHistoPar(0),
  76   fFittersOutR(38),
  77   fFittersOutZ(38)
  78 {
  79   //
  80   // Defualt constructor
  81   //
  82 }
  83
  84 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
  85   AliTPCcalibBase(name,title),
  86   fDiffHistoLine(0),
  87   fDiffHistoPar(0),
  88   fFittersOutR(38),
  89   fFittersOutZ(38)
  90 {
  91   //
  92   // Non default constructor
  93   //
  94   MakeFitters();
  95 }
  96
  97 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
  98   //
  99   //
 100   //
 101   if (fDiffHistoLine) delete fDiffHistoLine;
 102   if (fDiffHistoPar)  delete fDiffHistoPar;
 103 }
 104
 105
 106 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
 107   //
 108   //
 109   //
 110   const Int_t  kdrow=10;
 111   const Int_t kMinCluster=35;
 112   if (!fDiffHistoLine) MakeHisto();
 113   if (!seed) return;
 114   Int_t counter[72];
 115   for (Int_t i=0;i<72;i++) counter[i]=0;
 116   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 117     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 118     if (!cluster) continue;
 119     if (cluster->GetType()<0) continue;
 120     counter[cluster->GetDetector()]++;
 121   }
 122   //
 123   for (Int_t is=0; is<72;is++){
 124     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
 125   }
 126 }
 127
 128
 129 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
 130   //
 131   // process differnce of the cluster in respect with linear and parabolic fit
 132   //
 133   const  Double_t kXmean=134;
 134   const Int_t     kdrow=10;
 135   static TLinearFitter fy1(2,"pol1");
 136   static TLinearFitter fy2(3,"pol2");
 137   static TLinearFitter fz1(2,"pol1");
 138   static TLinearFitter fz2(3,"pol2");
 139   //
 140   static TVectorD py1(2);
 141   static TVectorD py2(3);
 142   static TVectorD pz1(2);
 143   static TVectorD pz2(3);
 144   //
 145   fy1.ClearPoints();
 146   fy2.ClearPoints();
 147   fz1.ClearPoints();
 148   fz2.ClearPoints();
 149   //
 150   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
 151     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 152     if (!c) continue;
 153     if (c->GetDetector()!=isec) continue;
 154     if (c->GetType()<0) continue;
 155     Float_t dx = c->GetX()-kXmean;
 156     Double_t x[2]={dx, dx*dx};
 157     fy1.AddPoint(x,c->GetY());
 158     fy2.AddPoint(x,c->GetY());
 159     fz1.AddPoint(x,c->GetZ());
 160     fz2.AddPoint(x,c->GetZ());
 161   }
 162   fy1.Eval(); fy1.GetParameters(py1);
 163   fy2.Eval(); fy2.GetParameters(py2);
 164   fz1.Eval(); fz1.GetParameters(pz1);
 165   fz2.Eval(); fz2.GetParameters(pz2);
 166   //
 167
 168   for (Int_t irow=0;irow<159;irow++) {
 169     AliTPCclusterMI *c=track->GetClusterPointer(irow);
 170     if (!c) continue;
 171     if (c->GetDetector()!=isec) continue;
 172     if (c->GetType()<0) continue;
 173     Float_t dx = c->GetX()-kXmean;
 174     Float_t y1 = py1[0]+py1[1]*dx;
 175     Float_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
 176     //
 177     Float_t z1 = pz1[0]+pz1[1]*dx;
 178     Float_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
 179     //
 180     //
 181     Double_t x[10];
 182     x[0]=isec;
 183     x[1]=irow;
 184     x[2]=c->GetY();
 185     x[3]=c->GetZ();
 186     x[4]=py1[1];
 187     x[5]=pz1[1];
 188     x[6]=py2[2]*150*150/4; // sagita 150 cm
 189     //
 190     x[7]=c->GetY()-y1;
 191     x[8]=c->GetZ()-z1;
 192     fDiffHistoLine->Fill(x);
 193     x[7]=c->GetY()-y2;
 194     x[8]=c->GetZ()-z2;
 195     fDiffHistoPar->Fill(x);
 196
 197     if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(py2[2]*150*150/4)<1){
 198       // Apply sagita cut
 199       AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
 200                py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
 201     }
 202
 203     if (fStreamLevel>1){
 204       TTreeSRedirector *cstream = GetDebugStreamer();
 205       if (cstream){
 206         (*cstream)<<"Diff"<<
 207           "isec="<<isec<<
 208           "Cl.="<<c<<
 209           "y1="<<y1<<
 210           "y2="<<y2<<
 211           "z1="<<z1<<
 212           "z2="<<z2<<
 213           "py1.="<<&py1<<
 214           "py2.="<<&py2<<
 215           "pz1.="<<&pz1<<
 216           "pz2.="<<&pz2<<
 217           "\n";
 218       }
 219     }
 220   }
 221 }
 222
 223
 224 void AliTPCcalibUnlinearity::MakeHisto(){
 225   //
 226   //
 227   //
 228   TString  axisName[10];
 229   Double_t xmin[10],  xmax[10];
 230   Int_t    nbins[10];
 231   //
 232   //
 233   axisName[0]="sector";
 234   xmin[0]=0; xmax[0]=72; nbins[0]=72;
 235   //
 236   axisName[1]="row";
 237   xmin[1]=0; xmax[1]=159; nbins[1]=159;
 238   //
 239   axisName[2]="ly";
 240   xmin[2]=-50; xmax[2]=50; nbins[2]=10;
 241   //
 242   axisName[3]="lz";
 243   xmin[3]=-250; xmax[3]=250; nbins[3]=50;
 244   //
 245   axisName[4]="p2";
 246   xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
 247   //
 248   axisName[5]="p3";
 249   xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
 250   //
 251   axisName[6]="p4";
 252   xmin[6]=-2; xmax[6]=2; nbins[6]=20;
 253   //
 254   axisName[7]="dy";
 255   xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
 256   //
 257   axisName[8]="dz";
 258   xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
 259   //
 260   fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
 261   fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
 262   for (Int_t i=0; i<9;i++){
 263     fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
 264     fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
 265     fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
 266     fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
 267   }
 268 }
 269
 270
 271 void AliTPCcalibUnlinearity::Terminate(){
 272   //
 273   // Terminate function
 274   // call base terminate + Eval of fitters
 275   //
 276   Info("AliTPCcalibUnlinearities","Terminate");
 277   EvalFitters();
 278   DumpTree();
 279   AliTPCcalibBase::Terminate();
 280 }
 281
 282
 283 void AliTPCcalibUnlinearity::DumpTree(){
 284   //
 285   //
 286   //
 287   if (fStreamLevel==0) return;
 288   TTreeSRedirector *cstream = GetDebugStreamer();
 289   if (!cstream) return;
 290   //
 291   THnSparse *his=0;
 292   Double_t position[10];
 293   Double_t value;
 294   Int_t *bins = new Int_t[10];
 295   //
 296   for (Int_t ihis=0; ihis<2; ihis++){
 297     his =  (ihis==0)? fDiffHistoLine:fDiffHistoPar;
 298     if (!his) continue;
 299     //
 300     Int_t ndim = his->GetNdimensions();
 301     //
 302     for (Long64_t i = 0; i < his->GetNbins(); ++i) {
 303       value = his->GetBinContent(i, bins);
 304       for (Int_t idim = 0; idim < ndim; idim++) {
 305         position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
 306       }
 307       (*cstream)<<"Resol"<<
 308         "ihis="<<ihis<<
 309         "bincont="<<value;
 310       for (Int_t idim=0;idim<ndim;idim++){
 311         (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
 312       }
 313       (*cstream)<<"Resol"<<
 314         "\n";
 315     }
 316   }
 317 }
 318
 319
 320 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Float_t dz, Float_t dy, Float_t p2, Float_t p3, Float_t dr, Int_t npoints){
 321   //
 322   //
 323   // Simple distortion fit in outer sectors
 324   //
 325   // Model - 2 exponential decay toward the center of chamber
 326   //       - Decay length 10 and 20 cm
 327   //       - Decay amplitude - Z dependent
 328   //
 329   /*
 330   chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
 331   chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
 332   chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
 333   chainUnlin->SetAlias("dzt","dz*side");
 334   chainUnlin->SetAlias("dr","(1-abs(lz/250))");
 335   chainUnlin->SetAlias("dout","(159-row)*1.5");
 336   chainUnlin->SetAlias("din","row*0.75");
 337   chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
 338   chainUnlin->SetAlias("eout20","exp(-(dout)/20.)");
 339   */
 340   Float_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
 341   Float_t dzt    = dz*side;
 342   Float_t dout   = (159-row)*1.5;  // distance to the last pad row - (valid only for long pads)
 343   Float_t eout10 = TMath::Exp(-dout/10.);
 344   Float_t eout20 = TMath::Exp(-dout/20.);
 345   //
 346   Double_t xxxR[6], xxxZ[6], xxxRZ[6];
 347   //TString fstring="";
 348   xxxZ[0]=eout20;                //fstring+="eout20++";
 349   xxxZ[1]=eout10;               //fstring+="eout10++";
 350   xxxZ[1]=dr*eout20;            //fstring+="dr*eout20++";
 351   xxxZ[1]=dr*eout10;            //fstring+="dr*eout10++";
 352   xxxZ[1]=dr*dr*eout20;         //fstring+="dr*dr*eout20++";
 353   xxxZ[1]=dr*dr*eout10;         //fstring+="dr*dr*eout10++";
 354   //
 355   xxxR[0]=p2*eout20;                //fstring+="eout20++";
 356   xxxR[1]=p2*eout10;               //fstring+="eout10++";
 357   xxxR[1]=p2*dr*eout20;            //fstring+="dr*eout20++";
 358   xxxR[1]=p2*dr*eout10;            //fstring+="dr*eout10++";
 359   xxxR[1]=p2*dr*dr*eout20;         //fstring+="dr*dr*eout20++";
 360   xxxR[1]=p2*dr*dr*eout10;         //fstring+="dr*dr*eout10++";
 361   //
 362   xxxRZ[0]=p3*eout20;                //fstring+="eout20++";
 363   xxxRZ[1]=p3*eout10;               //fstring+="eout10++";
 364   xxxRZ[1]=p3*dr*eout20;            //fstring+="dr*eout20++";
 365   xxxRZ[1]=p3*dr*eout10;            //fstring+="dr*eout10++";
 366   xxxRZ[1]=p3*dr*dr*eout20;         //fstring+="dr*dr*eout20++";
 367   xxxRZ[1]=p3*dr*dr*eout10;         //fstring+="dr*dr*eout10++";
 368   //
 369   TLinearFitter * fitter=0;
 370   for (Int_t ip=0; ip<npoints; ip++){
 371     //
 372     fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
 373     fitter->AddPoint(xxxR,dy);
 374     fitter->AddPoint(xxxRZ,dz);
 375     fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
 376     fitter->AddPoint(xxxZ,dzt);
 377     //
 378     if (side>0){
 379       fitter = (TLinearFitter*)fFittersOutR.At(36);
 380       fitter->AddPoint(xxxR,dy);
 381       fitter->AddPoint(xxxRZ,dz);
 382       fitter = (TLinearFitter*)fFittersOutZ.At(36);
 383       fitter->AddPoint(xxxZ,dzt);
 384     }
 385     if (side<0){
 386       fitter = (TLinearFitter*)fFittersOutR.At(36);
 387       fitter->AddPoint(xxxR,dy);
 388       fitter->AddPoint(xxxRZ,dz);
 389       fitter = (TLinearFitter*)fFittersOutZ.At(36);
 390       fitter->AddPoint(xxxZ,dzt);
 391     }
 392   }
 393 }
 394
 395
 396 void AliTPCcalibUnlinearity::MakeFitters(){
 397   //
 398   //
 399   //
 400   // Make outer fitters
 401   TLinearFitter * fitter = 0;
 402   for (Int_t ifit=0; ifit<38; ifit++){
 403     fitter = new TLinearFitter(7,"hyp6");
 404     fitter->StoreData(kFALSE);
 405     fFittersOutR.AddAt(fitter,ifit);
 406     fitter = new TLinearFitter(7,"hyp6");
 407     fitter->StoreData(kFALSE);
 408     fFittersOutZ.AddAt(fitter,ifit);
 409   }
 410 }
 411
 412 void AliTPCcalibUnlinearity::EvalFitters(){
 413   //
 414   //
 415   // Evaluate fitters
 416   //
 417   TLinearFitter * fitter = 0;
 418   for (Int_t ifit=0; ifit<38; ifit++){
 419     fitter=(TLinearFitter*)fFittersOutR.At(ifit);
 420     if (fitter) fitter->Eval();
 421     fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
 422     if (fitter) fitter->Eval();
 423   }
 424 }
 425
 426
 427 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Int_t nmaxPoints){
 428   //
 429   //
 430   //
 431   //  TTree * tree = chainUnlinP;
 432   AliTPCcalibUnlinearity *calib = this;
 433   //
 434   AliTPCclusterMI * cl=0;
 435   Float_t ty,tz;
 436   TVectorD *vy=0, *vz=0;
 437   Int_t nentries = TMath::Min(Int_t(tree->GetEntries()), nmaxPoints);
 438   //
 439   {
 440   for (Int_t i=0; i<nentries; i++){
 441     tree->SetBranchAddress("Cl.",&cl);
 442     tree->SetBranchAddress("y1",&ty);
 443     tree->SetBranchAddress("z1",&tz);
 444     tree->SetBranchAddress("py1.",&vy);
 445     tree->SetBranchAddress("pz1.",&vz);
 446     //
 447     tree->GetEntry(i);
 448     if (!cl) continue;
 449     calib->AddPoint(cl->GetDetector(), cl->GetRow(), cl->GetZ()-tz, cl->GetY()-ty,
 450                     (*vy)(1),(*vz)(1), 1-TMath::Abs(cl->GetZ()/250),1);
 451   }
 452   }
 453 }