TPC/AliTPCExBEffectiveSector.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 // AliTPCExBEffectiveSector class                                                   //
  19 // Correct for the rest of ExB effect which are not covered yet by physical models
  20 //
  21 // Motivation:
  22 //   ExB correction:
  23 //      dr    =  c0* integral(Er/Ez) + c1* integral(Erphi/Ez)
  24 //      drphi = -c1* integral(Er/Ez) + c0* integral(Erphi/Ez)
  25 //   Where:
  26 //   wt = Bz*(k*vdrift/E)           ~ 0.3 at B=0.5 T
  27 //   c0 = 1/(1+T2*T2*wt*wt)
  28 //   c1 = T1*wt/(1+T1*T1*wt*wt)
  29 //
  30 //
  31 //  3 correction maps 0 implemented as histogram used
  32 //  R-Phi correction map obtained minimizing residuals betwee the track
  33 //        and space points (AliTPCcalibAlign class). Track is defined using
  34 //        the points from the refernce plain at the middle of the TPC
  35 //        and vertex
  36 //        Corrected primar tracks straight pointing to the primary vertex
  37 //
  38 //  R distortion - obtained using the cluster residuals in the setup with
  39 //                 plus and minus field
  40 //                 Only high momenta tracks used for this calibration (1 GeV threshold)
  41 //     drphi_plus-drphi_minus=-2*c1 integral(Er/Ez)
  42 //               - Erphi/Ez cancels
  43 ////////////////////////////////////////////////////////////////////////////
  44 #include "AliMagF.h"
  45 #include "TGeoGlobalMagField.h"
  46 #include "AliTPCcalibDB.h"
  47 #include "AliTPCParam.h"
  48 #include "AliLog.h"
  49
  50 #include "TMath.h"
  51 #include "AliTPCROC.h"
  52 #include "TFile.h"
  53 #include "TAxis.h"
  54 #include "TTree.h"
  55 #include "TTreeStream.h"
  56 #include "THnSparse.h"
  57 #include "THnBase.h"
  58 #include "TProfile.h"
  59 #include "TH2F.h"
  60 #include "TH3F.h"
  61 #include "TROOT.h"
  62 #include "AliTPCExBEffectiveSector.h"
  63 ClassImp(AliTPCExBEffectiveSector)
  64
  65 AliTPCExBEffectiveSector::AliTPCExBEffectiveSector()
  66   : AliTPCCorrection("ExB_effectiveSector","ExB effective sector"),
  67     fC0(1.),fC1(0.),
  68     fCorrectionR(0),        // radial correction
  69     fCorrectionRPhi(0),     // r-phi correction
  70     fCorrectionZ(0)        // z correction
  71 {
  72   //
  73   // default constructor
  74   //
  75 }
  76
  77 AliTPCExBEffectiveSector::~AliTPCExBEffectiveSector() {
  78   //
  79   // default destructor
  80   //
  81   delete fCorrectionR;        // radial correction
  82   delete fCorrectionRPhi;     // r-phi correction
  83   delete fCorrectionZ;        // z correction
  84 }
  85
  86
  87
  88 void AliTPCExBEffectiveSector::Init() {
  89   //
  90   // Initialization funtion
  91   //
  92
  93   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  94   if (!magF) AliError("Magneticd field - not initialized");
  95   Double_t bzField = magF->SolenoidField()/10.; //field in T
  96   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  97   if (!param) AliError("Parameters - not initialized");
  98   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  99   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 100   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 101   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 102   SetOmegaTauT1T2(wt,fT1,fT2);
 103 }
 104
 105 void AliTPCExBEffectiveSector::Update(const TTimeStamp &/*timeStamp*/) {
 106   //
 107   // Update function
 108   //
 109   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 110   if (!magF) AliError("Magneticd field - not initialized");
 111   Double_t bzField = magF->SolenoidField()/10.; //field in T
 112   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 113   if (!param) AliError("Parameters - not initialized");
 114   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 115   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 116   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 117   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 118   SetOmegaTauT1T2(wt,fT1,fT2);
 119 }
 120
 121
 122
 123 void AliTPCExBEffectiveSector::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 124   //
 125   // Calculates the correction using the lookup table (histogram) of distortion
 126   // The histogram is created as poscl - postrack
 127   //
 128   dx[0]=0;
 129   dx[1]=0;
 130   dx[2]=0;
 131   if (!fCorrectionRPhi) return;
 132   Double_t phi      = TMath::ATan2(x[1],x[0]);
 133   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
 134   Double_t sector   = 9.*phi/TMath::Pi();
 135   if (sector<0) sector+=18.;
 136   Double_t        kZ=x[2]/r;
 137   //
 138   if (kZ>1.2)        kZ= 1.2;
 139   if (kZ<-1.2)       kZ= -1.2;
 140   if (roc%36<18)  kZ= TMath::Abs(kZ);
 141   if (roc%36>=18) kZ=-TMath::Abs(kZ);
 142   if (TMath::Abs(kZ)<0.15){
 143     kZ = (roc%36<18) ? 0.15:-0.15;
 144   }
 145   //
 146   Double_t dlR=0;
 147   Double_t dlRPhi=0;
 148   Double_t dlZ=0;
 149   Double_t rr=TMath::Max(r,fCorrectionRPhi->GetYaxis()->GetXmin()+0.01);
 150   rr=TMath::Min(rr,fCorrectionRPhi->GetYaxis()->GetXmax()-0.01);
 151   Double_t kZZ=TMath::Max(kZ,fCorrectionRPhi->GetZaxis()->GetXmin()+0.001);
 152   kZZ=TMath::Min(kZZ,fCorrectionRPhi->GetZaxis()->GetXmax()-0.001);
 153
 154   if (fCorrectionRPhi) {
 155     //    dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
 156     dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
 157   }
 158   if (fCorrectionR)    {
 159     //    dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
 160     dlR= -fCorrectionR->GetBinContent(fCorrectionR->FindBin(sector,rr,kZZ));
 161   }
 162   if (fCorrectionZ)    {
 163     //    dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
 164     dlZ= -fCorrectionZ->GetBinContent(fCorrectionZ->FindBin(sector,rr,kZZ));
 165   }
 166   Double_t dr    = fC0*dlR  + fC1*dlRPhi;
 167   Double_t drphi = -fC1*dlR + fC0*dlRPhi;
 168    // Calculate distorted position
 169   if ( r > 0.0 ) {
 170     r   =  r   + dr;
 171     phi =  phi + drphi/r;
 172   }
 173   // Calculate correction in cartesian coordinates
 174   dx[0] = r * TMath::Cos(phi) - x[0];
 175   dx[1] = r * TMath::Sin(phi) - x[1];
 176   dx[2] = dlZ;
 177
 178 }
 179
 180 void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
 181   //
 182   // Print function to check the settings (e.g. the twist in the X direction)
 183   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 184   //
 185
 186   TString opt = option; opt.ToLower();
 187   printf("%s\t%s\n",GetName(),GetTitle());
 188   if (opt.Contains("a")) { // Print all details
 189     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 190     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
 191   }
 192 }
 193
 194 void  AliTPCExBEffectiveSector::MakeResidualMap(THnBase * hisInput, const char *sname, Int_t ptype, Int_t dtype){
 195   //
 196   // Make cluster residual map from the n-dimensional histogram
 197   // hisInput supposed to have given format:
 198   //          - 4 Dim  - delta, sector, localX, kZ
 199   // Vertex position assumed to be at (0,0,0)
 200   TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
 201   //
 202   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
 203   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
 204   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
 205   TF1 *fgaus=0;
 206   TH3F * hisResMap3D =
 207     new TH3F("his3D","his3D",
 208              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 209              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
 210              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
 211   hisResMap3D->GetXaxis()->SetTitle("sector");
 212   hisResMap3D->GetYaxis()->SetTitle("localX");
 213   hisResMap3D->GetZaxis()->SetTitle("kZ");
 214
 215   TH2F * hisResMap2D[4] ={0,0,0,0};
 216   for (Int_t i=0; i<4; i++){
 217     hisResMap2D[i]=
 218       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
 219                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 220                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
 221     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
 222     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
 223   }
 224   //
 225   //
 226   //
 227   TF1 * f1= 0;
 228   Int_t axis0[4]={0,1,2,3};
 229   Int_t axis1[4]={0,1,2,3};
 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       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
 276
 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 }