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 "TProfile.h"
  58 #include "TH2F.h"
  59 #include "TH3F.h"
  60 #include "TROOT.h"
  61 #include "AliTPCExBEffectiveSector.h"
  62 ClassImp(AliTPCExBEffectiveSector)
  63
  64 AliTPCExBEffectiveSector::AliTPCExBEffectiveSector()
  65   : AliTPCCorrection("ExB_effectiveSector","ExB effective sector"),
  66     fC0(1.),fC1(0.),
  67     fCorrectionR(0),        // radial correction
  68     fCorrectionRPhi(0),     // r-phi correction
  69     fCorrectionZ(0)        // z correction
  70 {
  71   //
  72   // default constructor
  73   //
  74 }
  75
  76 AliTPCExBEffectiveSector::~AliTPCExBEffectiveSector() {
  77   //
  78   // default destructor
  79   //
  80   delete fCorrectionR;        // radial correction
  81   delete fCorrectionRPhi;     // r-phi correction
  82   delete fCorrectionZ;        // z correction
  83 }
  84
  85
  86
  87 void AliTPCExBEffectiveSector::Init() {
  88   //
  89   // Initialization funtion
  90   //
  91
  92   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  93   if (!magF) AliError("Magneticd field - not initialized");
  94   Double_t bzField = magF->SolenoidField()/10.; //field in T
  95   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  96   if (!param) AliError("Parameters - not initialized");
  97   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  98   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  99   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 100   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 101   SetOmegaTauT1T2(wt,fT1,fT2);
 102 }
 103
 104 void AliTPCExBEffectiveSector::Update(const TTimeStamp &/*timeStamp*/) {
 105   //
 106   // Update function
 107   //
 108   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 109   if (!magF) AliError("Magneticd field - not initialized");
 110   Double_t bzField = magF->SolenoidField()/10.; //field in T
 111   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
 112   if (!param) AliError("Parameters - not initialized");
 113   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
 114   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
 115   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
 116   // Correction Terms for effective omegaTau; obtained by a laser calibration run
 117   SetOmegaTauT1T2(wt,fT1,fT2);
 118 }
 119
 120
 121
 122 void AliTPCExBEffectiveSector::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 123   //
 124   // Calculates the correction using the lookup table (histogram) of distortion
 125   // The histogram is created as poscl - postrack
 126   //
 127   Double_t phi      = TMath::ATan2(x[1],x[0]);
 128   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
 129   Double_t sector   = 9.*phi/TMath::Pi();
 130   if (sector<0) sector+=18.;
 131   Double_t        kZ=x[2]/r;
 132   //
 133   if (kZ>1.2)        kZ= 1.2;
 134   if (kZ<-1.2)       kZ= -1.2;
 135   if (roc%36<18)  kZ= TMath::Abs(kZ);
 136   if (roc%36>=18) kZ=-TMath::Abs(kZ);
 137   if (TMath::Abs(kZ)<0.15){
 138     kZ = (roc%36<18) ? 0.15:-0.15;
 139   }
 140
 141   //
 142   Double_t dlR=0;
 143   Double_t dlRPhi=0;
 144   Double_t dlZ=0;
 145   if (fCorrectionRPhi) {
 146     Double_t rr=TMath::Max(r,fCorrectionRPhi->GetYaxis()->GetXmin());
 147     rr=TMath::Min(rr,fCorrectionRPhi->GetYaxis()->GetXmax());
 148     Double_t kZZ=TMath::Max(kZ,fCorrectionRPhi->GetZaxis()->GetXmin());
 149     kZZ=TMath::Min(kZZ,fCorrectionRPhi->GetZaxis()->GetXmax());
 150     //    dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
 151     dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
 152   }
 153   if (fCorrectionR)    {
 154     Double_t rr=TMath::Max(r,fCorrectionR->GetYaxis()->GetXmin());
 155     rr=TMath::Min(rr,fCorrectionR->GetYaxis()->GetXmax());
 156     Double_t kZZ=TMath::Max(kZ,fCorrectionR->GetZaxis()->GetXmin());
 157     kZZ=TMath::Min(kZZ,fCorrectionR->GetZaxis()->GetXmax());
 158     dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
 159   }
 160   if (fCorrectionZ)    {
 161     Double_t rr=TMath::Max(r,fCorrectionZ->GetYaxis()->GetXmin());
 162     rr=TMath::Min(rr,fCorrectionZ->GetYaxis()->GetXmax());
 163     Double_t kZZ=TMath::Max(kZ,fCorrectionZ->GetZaxis()->GetXmin());
 164     kZZ=TMath::Min(kZZ,fCorrectionZ->GetZaxis()->GetXmax());
 165     dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
 166   }
 167   Double_t dr    = fC0*dlR  + fC1*dlRPhi;
 168   Double_t drphi = -fC1*dlR + fC0*dlRPhi;
 169    // Calculate distorted position
 170   if ( r > 0.0 ) {
 171     r   =  r   + dr;
 172     phi =  phi + drphi/r;
 173   }
 174   // Calculate correction in cartesian coordinates
 175   dx[0] = r * TMath::Cos(phi) - x[0];
 176   dx[1] = r * TMath::Sin(phi) - x[1];
 177   dx[2] = dlZ;
 178
 179 }
 180
 181 void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
 182   //
 183   // Print function to check the settings (e.g. the twist in the X direction)
 184   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 185   //
 186
 187   TString opt = option; opt.ToLower();
 188   printf("%s\t%s\n",GetName(),GetTitle());
 189   if (opt.Contains("a")) { // Print all details
 190     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 191     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
 192   }
 193 }
 194
 195 void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname){
 196   //
 197   // Make cluster residual map from the n-dimensional histogram
 198   // hisInput supposed to have given format:
 199   //          - 4 Dim  - delta, sector, localX, kZ
 200   // Vertex position assumed to be at (0,0,0)
 201   TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
 202   //
 203   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
 204   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
 205   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
 206
 207   TH3F * hisResMap3D =
 208     new TH3F("his3D","his3D",
 209              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 210              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
 211              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
 212   hisResMap3D->GetXaxis()->SetTitle("sector");
 213   hisResMap3D->GetYaxis()->SetTitle("localX");
 214   hisResMap3D->GetZaxis()->SetTitle("kZ");
 215
 216   TH2F * hisResMap2D[4] ={0,0,0,0};
 217   for (Int_t i=0; i<4; i++){
 218     hisResMap2D[i]=
 219       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
 220                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 221                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
 222     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
 223     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
 224   }
 225   //
 226   //
 227   //
 228   TF1 * f1= 0;
 229   Int_t axis0[4]={0,1,2,3};
 230   Int_t axis1[4]={0,1,2,3};
 231   for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
 232     // phi- sector  range
 233     hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
 234     THnSparse *his1=hisInput->Projection(4,axis0);
 235     Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
 236     //
 237     for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
 238       // local x range
 239       // kz fits
 240       his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
 241       THnSparse *his2=his1->Projection(4,axis1);
 242       Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
 243       //
 244       //A side
 245       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 246       TH1 * hisA = his2->Projection(0);
 247       Double_t meanA= hisA->GetMean();
 248       Double_t rmsA= hisA->GetRMS();
 249       Double_t entriesA= hisA->GetEntries();
 250       delete hisA;
 251       //C side
 252       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 253       TH1 * hisC = his2->Projection(0);
 254       Double_t meanC= hisC->GetMean();
 255       Double_t rmsC= hisC->GetRMS();
 256       Double_t entriesC= hisC->GetEntries();
 257       delete hisC;
 258       his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
 259       TH2 * hisAC       = his2->Projection(0,3);
 260       TProfile *profAC  = hisAC->ProfileX();
 261       delete hisAC;
 262       profAC->Fit("pol1","QNR","QNR",0.05,1);
 263       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
 264       Double_t offsetA=f1->GetParameter(0);
 265       Double_t slopeA=f1->GetParameter(1);
 266       Double_t offsetAE=f1->GetParError(0);
 267       Double_t slopeAE=f1->GetParError(1);
 268       Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
 269       profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
 270       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
 271       Double_t offsetC=f1->GetParameter(0);
 272       Double_t slopeC=f1->GetParameter(1);
 273       Double_t offsetCE=f1->GetParError(0);
 274       Double_t slopeCE=f1->GetParError(1);
 275       Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
 276       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
 277
 278       (*pcstream)<<"deltaFit"<<
 279         "sector="<<sector<<
 280         "localX="<<localX<<
 281         "meanA="<<meanA<<
 282         "rmsA="<<rmsA<<
 283         "entriesA="<<entriesA<<
 284         "meanC="<<meanC<<
 285         "rmsC="<<rmsC<<
 286         "entriesC="<<entriesC<<
 287         "offsetA="<<offsetA<<
 288         "slopeA="<<slopeA<<
 289         "offsetAE="<<offsetAE<<
 290         "slopeAE="<<slopeAE<<
 291         "chi2A="<<chi2A<<
 292         "offsetC="<<offsetC<<
 293         "slopeC="<<slopeC<<
 294         "offsetCE="<<offsetCE<<
 295         "slopeCE="<<slopeCE<<
 296         "chi2C="<<chi2C<<
 297         "\n";
 298       //
 299       hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
 300       hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
 301       hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
 302       hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
 303
 304       for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
 305         Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
 306         if (TMath::Abs(kZ)<0.05) continue;  // crossing
 307         his2->GetAxis(3)->SetRange(ibin3,ibin3);
 308         if (TMath::Abs(kZ)>0.15){
 309           his2->GetAxis(3)->SetRange(ibin3,ibin3);
 310         }
 311         TH1 * his = his2->Projection(0);
 312         Double_t mean= his->GetMean();
 313         Double_t rms= his->GetRMS();
 314         Double_t entries= his->GetEntries();
 315         //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
 316         hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
 317         (*pcstream)<<"delta"<<
 318           "sector="<<sector<<
 319           "localX="<<localX<<
 320           "kZ="<<kZ<<
 321           "mean="<<mean<<
 322           "rms="<<rms<<
 323           "entries="<<entries<<
 324           "meanA="<<meanA<<
 325           "rmsA="<<rmsA<<
 326           "entriesA="<<entriesA<<
 327           "meanC="<<meanC<<
 328           "rmsC="<<rmsC<<
 329           "entriesC="<<entriesC<<
 330           "offsetA="<<offsetA<<
 331           "slopeA="<<slopeA<<
 332           "chi2A="<<chi2A<<
 333           "offsetC="<<offsetC<<
 334           "slopeC="<<slopeC<<
 335           "chi2C="<<chi2C<<
 336           "\n";
 337         delete his;
 338       }
 339       delete his2;
 340     }
 341     delete his1;
 342   }
 343   hisResMap3D->Write();
 344   hisResMap2D[0]->Write();
 345   hisResMap2D[1]->Write();
 346   hisResMap2D[2]->Write();
 347   hisResMap2D[3]->Write();
 348   delete pcstream;
 349 }