]>
Commit | Line | Data |
---|---|---|
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" | |
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 | } |