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   dx[0]=0;
 128   dx[1]=0;
 129   dx[2]=0;
 130   if (!fCorrectionRPhi) return;
 131   Double_t phi      = TMath::ATan2(x[1],x[0]);
 132   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
 133   Double_t sector   = 9.*phi/TMath::Pi();
 134   if (sector<0) sector+=18.;
 135   Double_t        kZ=x[2]/r;
 136   //
 137   if (kZ>1.2)        kZ= 1.2;
 138   if (kZ<-1.2)       kZ= -1.2;
 139   if (roc%36<18)  kZ= TMath::Abs(kZ);
 140   if (roc%36>=18) kZ=-TMath::Abs(kZ);
 141   if (TMath::Abs(kZ)<0.15){
 142     kZ = (roc%36<18) ? 0.15:-0.15;
 143   }
 144   //
 145   Double_t dlR=0;
 146   Double_t dlRPhi=0;
 147   Double_t dlZ=0;
 148   Double_t rr=TMath::Max(r,fCorrectionRPhi->GetYaxis()->GetXmin()+0.01);
 149   rr=TMath::Min(rr,fCorrectionRPhi->GetYaxis()->GetXmax()-0.01);
 150   Double_t kZZ=TMath::Max(kZ,fCorrectionRPhi->GetZaxis()->GetXmin()+0.001);
 151   kZZ=TMath::Min(kZZ,fCorrectionRPhi->GetZaxis()->GetXmax()-0.001);
 152
 153   if (fCorrectionRPhi) {
 154     //    dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
 155     dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
 156   }
 157   if (fCorrectionR)    {
 158     //    dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
 159     dlR= -fCorrectionR->GetBinContent(fCorrectionR->FindBin(sector,rr,kZZ));
 160   }
 161   if (fCorrectionZ)    {
 162     //    dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
 163     dlZ= -fCorrectionZ->GetBinContent(fCorrectionZ->FindBin(sector,rr,kZZ));
 164   }
 165   Double_t dr    = fC0*dlR  + fC1*dlRPhi;
 166   Double_t drphi = -fC1*dlR + fC0*dlRPhi;
 167    // Calculate distorted position
 168   if ( r > 0.0 ) {
 169     r   =  r   + dr;
 170     phi =  phi + drphi/r;
 171   }
 172   // Calculate correction in cartesian coordinates
 173   dx[0] = r * TMath::Cos(phi) - x[0];
 174   dx[1] = r * TMath::Sin(phi) - x[1];
 175   dx[2] = dlZ;
 176
 177 }
 178
 179 void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
 180   //
 181   // Print function to check the settings (e.g. the twist in the X direction)
 182   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 183   //
 184
 185   TString opt = option; opt.ToLower();
 186   printf("%s\t%s\n",GetName(),GetTitle());
 187   if (opt.Contains("a")) { // Print all details
 188     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 189     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
 190   }
 191 }
 192
 193 void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype){
 194   //
 195   // Make cluster residual map from the n-dimensional histogram
 196   // hisInput supposed to have given format:
 197   //          - 4 Dim  - delta, sector, localX, kZ
 198   // Vertex position assumed to be at (0,0,0)
 199   TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
 200   //
 201   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
 202   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
 203   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
 204   TF1 *fgaus=0;
 205   TH3F * hisResMap3D =
 206     new TH3F("his3D","his3D",
 207              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 208              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
 209              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
 210   hisResMap3D->GetXaxis()->SetTitle("sector");
 211   hisResMap3D->GetYaxis()->SetTitle("localX");
 212   hisResMap3D->GetZaxis()->SetTitle("kZ");
 213
 214   TH2F * hisResMap2D[4] ={0,0,0,0};
 215   for (Int_t i=0; i<4; i++){
 216     hisResMap2D[i]=
 217       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
 218                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
 219                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
 220     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
 221     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
 222   }
 223   //
 224   //
 225   //
 226   TF1 * f1= 0;
 227   Int_t axis0[4]={0,1,2,3};
 228   Int_t axis1[4]={0,1,2,3};
 229   for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
 230     // phi- sector  range
 231     hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
 232     THnSparse *his1=hisInput->Projection(4,axis0);
 233     Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
 234     //
 235     for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
 236       // local x range
 237       // kz fits
 238       his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
 239       THnSparse *his2=his1->Projection(4,axis1);
 240       Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
 241       //
 242       //A side
 243       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 244       TH1 * hisA = his2->Projection(0);
 245       Double_t meanA= hisA->GetMean();
 246       Double_t rmsA= hisA->GetRMS();
 247       Double_t entriesA= hisA->GetEntries();
 248       delete hisA;
 249       //C side
 250       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
 251       TH1 * hisC = his2->Projection(0);
 252       Double_t meanC= hisC->GetMean();
 253       Double_t rmsC= hisC->GetRMS();
 254       Double_t entriesC= hisC->GetEntries();
 255       delete hisC;
 256       his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
 257       TH2 * hisAC       = his2->Projection(0,3);
 258       TProfile *profAC  = hisAC->ProfileX();
 259       delete hisAC;
 260       profAC->Fit("pol1","QNR","QNR",0.05,1);
 261       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
 262       Double_t offsetA=f1->GetParameter(0);
 263       Double_t slopeA=f1->GetParameter(1);
 264       Double_t offsetAE=f1->GetParError(0);
 265       Double_t slopeAE=f1->GetParError(1);
 266       Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
 267       profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
 268       f1=(TF1*)gROOT->FindObject("pol1");
 269       Double_t offsetC=f1->GetParameter(0);
 270       Double_t slopeC=f1->GetParameter(1);
 271       Double_t offsetCE=f1->GetParError(0);
 272       Double_t slopeCE=f1->GetParError(1);
 273       Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
 274       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
 275
 276       (*pcstream)<<"deltaFit"<<
 277         "sector="<<sector<<
 278         "localX="<<localX<<
 279         "meanA="<<meanA<<
 280         "rmsA="<<rmsA<<
 281         "entriesA="<<entriesA<<
 282         "meanC="<<meanC<<
 283         "rmsC="<<rmsC<<
 284         "entriesC="<<entriesC<<
 285         "offsetA="<<offsetA<<
 286         "slopeA="<<slopeA<<
 287         "offsetAE="<<offsetAE<<
 288         "slopeAE="<<slopeAE<<
 289         "chi2A="<<chi2A<<
 290         "offsetC="<<offsetC<<
 291         "slopeC="<<slopeC<<
 292         "offsetCE="<<offsetCE<<
 293         "slopeCE="<<slopeCE<<
 294         "chi2C="<<chi2C<<
 295         "\n";
 296       //
 297       hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
 298       hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
 299       hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
 300       hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
 301
 302       for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
 303         Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
 304         if (TMath::Abs(kZ)<0.05) continue;  // crossing
 305         his2->GetAxis(3)->SetRange(ibin3,ibin3);
 306         if (TMath::Abs(kZ)>0.15){
 307           his2->GetAxis(3)->SetRange(ibin3,ibin3);
 308         }
 309         TH1 * his = his2->Projection(0);
 310         Double_t mean= his->GetMean();
 311         Double_t rms= his->GetRMS();
 312         Double_t entries= his->GetEntries();
 313         //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
 314         hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
 315         Double_t phi=TMath::Pi()*sector/9;
 316         if (phi>TMath::Pi()) phi+=TMath::Pi();
 317         Double_t meanG=0;
 318         Double_t rmsG=0;
 319         if (entries>50){
 320           if (!fgaus) {
 321             his->Fit("gaus","Q","goff");
 322             fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
 323           }
 324           if (fgaus) {
 325             his->Fit(fgaus,"Q","goff");
 326             meanG=fgaus->GetParameter(1);
 327             rmsG=fgaus->GetParameter(2);
 328           }
 329         }
 330         Double_t dsec=sector-Int_t(sector)-0.5;
 331         Double_t snp=dsec*TMath::Pi()/9.;
 332         (*pcstream)<<"delta"<<
 333           "ptype="<<ptype<<
 334           "dtype="<<dtype<<
 335           "sector="<<sector<<
 336           "dsec="<<dsec<<
 337           "snp="<<snp<<
 338           "phi="<<phi<<
 339           "localX="<<localX<<
 340           "kZ="<<kZ<<
 341           "theta="<<kZ<<
 342           "mean="<<mean<<
 343           "rms="<<rms<<
 344           "meanG="<<meanG<<
 345           "rmsG="<<rmsG<<
 346           "entries="<<entries<<
 347           "meanA="<<meanA<<
 348           "rmsA="<<rmsA<<
 349           "entriesA="<<entriesA<<
 350           "meanC="<<meanC<<
 351           "rmsC="<<rmsC<<
 352           "entriesC="<<entriesC<<
 353           "offsetA="<<offsetA<<
 354           "slopeA="<<slopeA<<
 355           "chi2A="<<chi2A<<
 356           "offsetC="<<offsetC<<
 357           "slopeC="<<slopeC<<
 358           "chi2C="<<chi2C<<
 359           "\n";
 360         delete his;
 361       }
 362       delete his2;
 363     }
 364     delete his1;
 365   }
 366   hisResMap3D->Write();
 367   hisResMap2D[0]->Write();
 368   hisResMap2D[1]->Write();
 369   hisResMap2D[2]->Write();
 370   hisResMap2D[3]->Write();
 371   delete pcstream;
 372 }