]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCExBEffectiveSector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliTPCRecoParam.cxx - Removing old comment
[u/mrichter/AliRoot.git] / TPC / AliTPCExBEffectiveSector.cxx
CommitLineData
40787a37 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"
7f00a31e 57#include "THnBase.h"
40787a37 58#include "TProfile.h"
59#include "TH2F.h"
60#include "TH3F.h"
61#include "TROOT.h"
62#include "AliTPCExBEffectiveSector.h"
63ClassImp(AliTPCExBEffectiveSector)
64
65AliTPCExBEffectiveSector::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
77AliTPCExBEffectiveSector::~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
88void 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
105void 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
123void 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 //
28178a25 128 dx[0]=0;
129 dx[1]=0;
130 dx[2]=0;
131 if (!fCorrectionRPhi) return;
40787a37 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;
28178a25 144 }
40787a37 145 //
146 Double_t dlR=0;
147 Double_t dlRPhi=0;
148 Double_t dlZ=0;
28178a25 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
40787a37 154 if (fCorrectionRPhi) {
28178a25 155 // dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
156 dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
40787a37 157 }
158 if (fCorrectionR) {
28178a25 159 // dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
160 dlR= -fCorrectionR->GetBinContent(fCorrectionR->FindBin(sector,rr,kZZ));
40787a37 161 }
162 if (fCorrectionZ) {
28178a25 163 // dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
164 dlZ= -fCorrectionZ->GetBinContent(fCorrectionZ->FindBin(sector,rr,kZZ));
40787a37 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
180void 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
7f00a31e 194void AliTPCExBEffectiveSector::MakeResidualMap(THnBase * hisInput, const char *sname, Int_t ptype, Int_t dtype){
40787a37 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();
46e89793 205 TF1 *fgaus=0;
40787a37 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);
7f00a31e 233 THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
40787a37 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);
7f00a31e 240 THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
40787a37 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);
0bc13e06 269 f1=(TF1*)gROOT->FindObject("pol1");
40787a37 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);
28178a25 316 Double_t phi=TMath::Pi()*sector/9;
317 if (phi>TMath::Pi()) phi+=TMath::Pi();
46e89793 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());
46e89793 324 }
0bc13e06 325 if (fgaus) {
326 his->Fit(fgaus,"Q","goff");
46e89793 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.;
40787a37 333 (*pcstream)<<"delta"<<
46e89793 334 "ptype="<<ptype<<
335 "dtype="<<dtype<<
40787a37 336 "sector="<<sector<<
46e89793 337 "dsec="<<dsec<<
338 "snp="<<snp<<
28178a25 339 "phi="<<phi<<
40787a37 340 "localX="<<localX<<
341 "kZ="<<kZ<<
46e89793 342 "theta="<<kZ<<
40787a37 343 "mean="<<mean<<
344 "rms="<<rms<<
46e89793 345 "meanG="<<meanG<<
346 "rmsG="<<rmsG<<
40787a37 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}