1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ////////////////////////////////////////////////////////////////////////////
18 // AliTPCExBEffectiveSector class //
19 // Correct for the rest of ExB effect which are not covered yet by physical models
23 // dr = c0* integral(Er/Ez) + c1* integral(Erphi/Ez)
24 // drphi = -c1* integral(Er/Ez) + c0* integral(Erphi/Ez)
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)
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
36 // Corrected primar tracks straight pointing to the primary vertex
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)
43 ////////////////////////////////////////////////////////////////////////////
45 #include "TGeoGlobalMagField.h"
46 #include "AliTPCcalibDB.h"
47 #include "AliTPCParam.h"
51 #include "AliTPCROC.h"
55 #include "TTreeStream.h"
56 #include "THnSparse.h"
61 #include "AliTPCExBEffectiveSector.h"
62 ClassImp(AliTPCExBEffectiveSector)
64 AliTPCExBEffectiveSector::AliTPCExBEffectiveSector()
65 : AliTPCCorrection("ExB_effectiveSector","ExB effective sector"),
67 fCorrectionR(0), // radial correction
68 fCorrectionRPhi(0), // r-phi correction
69 fCorrectionZ(0) // z correction
72 // default constructor
76 AliTPCExBEffectiveSector::~AliTPCExBEffectiveSector() {
80 delete fCorrectionR; // radial correction
81 delete fCorrectionRPhi; // r-phi correction
82 delete fCorrectionZ; // z correction
87 void AliTPCExBEffectiveSector::Init() {
89 // Initialization funtion
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);
104 void AliTPCExBEffectiveSector::Update(const TTimeStamp &/*timeStamp*/) {
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);
122 void AliTPCExBEffectiveSector::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
124 // Calculates the correction using the lookup table (histogram) of distortion
125 // The histogram is created as poscl - postrack
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.;
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;
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);
153 if (fCorrectionRPhi) {
154 // dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
155 dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
158 // dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
159 dlR= -fCorrectionR->GetBinContent(fCorrectionR->FindBin(sector,rr,kZZ));
162 // dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
163 dlZ= -fCorrectionZ->GetBinContent(fCorrectionZ->FindBin(sector,rr,kZZ));
165 Double_t dr = fC0*dlR + fC1*dlRPhi;
166 Double_t drphi = -fC1*dlR + fC0*dlRPhi;
167 // Calculate distorted position
172 // Calculate correction in cartesian coordinates
173 dx[0] = r * TMath::Cos(phi) - x[0];
174 dx[1] = r * TMath::Sin(phi) - x[1];
179 void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
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
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);
193 void AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname){
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);
201 Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
202 Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
203 Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
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");
214 TH2F * hisResMap2D[4] ={0,0,0,0};
215 for (Int_t i=0; i<4; 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");
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){
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);
235 for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
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);
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();
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();
256 his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
257 TH2 * hisAC = his2->Projection(0,3);
258 TProfile *profAC = hisAC->ProfileX();
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 if (!f1) 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);
276 (*pcstream)<<"deltaFit"<<
281 "entriesA="<<entriesA<<
284 "entriesC="<<entriesC<<
285 "offsetA="<<offsetA<<
287 "offsetAE="<<offsetAE<<
288 "slopeAE="<<slopeAE<<
290 "offsetC="<<offsetC<<
292 "offsetCE="<<offsetCE<<
293 "slopeCE="<<slopeCE<<
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);
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);
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 (*pcstream)<<"delta"<<
324 "entries="<<entries<<
327 "entriesA="<<entriesA<<
330 "entriesC="<<entriesC<<
331 "offsetA="<<offsetA<<
334 "offsetC="<<offsetC<<
344 hisResMap3D->Write();
345 hisResMap2D[0]->Write();
346 hisResMap2D[1]->Write();
347 hisResMap2D[2]->Write();
348 hisResMap2D[3]->Write();