]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCorrection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bastian's Task for Jet Fragmentation, to be merged with other FF tasks
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrection.cxx
CommitLineData
0116859c 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// AliTPCCorrection class //
19// //
20// This class provides a general framework to deal with space point //
21// distortions. An correction class which inherits from here is for example //
22// AliTPCExBBShape or AliTPCExBTwist //
23// //
24// General functions are (for example): //
25// CorrectPoint(x,roc) where x is the vector of inital positions in //
26// cartesian coordinates and roc represents the Read Out chamber number //
27// according to the offline naming convention. The vector x is overwritten //
28// with the corrected coordinates. //
29// //
30// An alternative usage would be CorrectPoint(x,roc,dx), which leaves the //
31// vector x untouched, put returns the distortions via the vector dx //
32// //
33// The class allows "effective Omega Tau" corrections to be shifted to the //
34// single distortion classes. //
35// //
36// Note: This class is normally used via the class AliTPCComposedCorrection //
37// //
38// date: 27/04/2010 //
39// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas //
40////////////////////////////////////////////////////////////////////////////////
be67055b 41#include "Riostream.h"
0116859c 42
43#include <TH2F.h>
44#include <TMath.h>
45#include <TROOT.h>
cf5b0aa0 46#include <TTreeStream.h>
ffab0c37 47#include <TTree.h>
48#include <TFile.h>
e527a1b9 49#include <TTimeStamp.h>
ffab0c37 50#include <AliCDBStorage.h>
51#include <AliCDBId.h>
52#include <AliCDBMetaData.h>
7f4cb119 53#include "TVectorD.h"
54
55#include "TRandom.h"
56#include "AliTPCTransform.h"
57#include "AliTPCcalibDB.h"
58#include "AliTPCExB.h"
59#include "AliTPCCorrection.h"
60#include "AliTPCRecoParam.h"
cf5b0aa0 61
62#include "AliExternalTrackParam.h"
63#include "AliTrackPointArray.h"
64#include "TDatabasePDG.h"
65#include "AliTrackerBase.h"
66#include "AliTPCROC.h"
8b63d99c 67#include "THnSparse.h"
7f4cb119 68#include "AliTPCLaserTrack.h"
0116859c 69
70#include "AliTPCCorrection.h"
71
cf5b0aa0 72ClassImp(AliTPCCorrection)
73
0116859c 74// FIXME: the following values should come from the database
b1f0a2a5 75const Double_t AliTPCCorrection::fgkTPCZ0 =249.7; // nominal gating grid position
0116859c 76const Double_t AliTPCCorrection::fgkIFCRadius= 83.06; // Mean Radius of the Inner Field Cage ( 82.43 min, 83.70 max) (cm)
77const Double_t AliTPCCorrection::fgkOFCRadius=254.5; // Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)
78const Double_t AliTPCCorrection::fgkZOffSet = 0.2; // Offset from CE: calculate all distortions closer to CE as if at this point
79const Double_t AliTPCCorrection::fgkCathodeV =-100000.0; // Cathode Voltage (volts)
80const Double_t AliTPCCorrection::fgkGG =-70.0; // Gating Grid voltage (volts)
81
82
83// FIXME: List of interpolation points (course grid in the middle, fine grid on the borders)
84const Double_t AliTPCCorrection::fgkRList[AliTPCCorrection::kNR] = {
8584.0, 84.5, 85.0, 85.5, 86.0, 87.0, 88.0,
8690.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0,
87110.0, 112.0, 114.0, 116.0, 118.0, 120.0, 122.0, 124.0, 126.0, 128.0,
88130.0, 132.0, 134.0, 136.0, 138.0, 140.0, 142.0, 144.0, 146.0, 148.0,
89150.0, 152.0, 154.0, 156.0, 158.0, 160.0, 162.0, 164.0, 166.0, 168.0,
90170.0, 172.0, 174.0, 176.0, 178.0, 180.0, 182.0, 184.0, 186.0, 188.0,
91190.0, 192.0, 194.0, 196.0, 198.0, 200.0, 202.0, 204.0, 206.0, 208.0,
92210.0, 212.0, 214.0, 216.0, 218.0, 220.0, 222.0, 224.0, 226.0, 228.0,
93230.0, 232.0, 234.0, 236.0, 238.0, 240.0, 242.0, 244.0, 246.0, 248.0,
94249.0, 249.5, 250.0, 251.5, 252.0 } ;
95
96const Double_t AliTPCCorrection::fgkZList[AliTPCCorrection::kNZ] = {
97-249.5, -249.0, -248.5, -248.0, -247.0, -246.0, -245.0, -243.0, -242.0, -241.0,
98-240.0, -238.0, -236.0, -234.0, -232.0, -230.0, -228.0, -226.0, -224.0, -222.0,
99-220.0, -218.0, -216.0, -214.0, -212.0, -210.0, -208.0, -206.0, -204.0, -202.0,
100-200.0, -198.0, -196.0, -194.0, -192.0, -190.0, -188.0, -186.0, -184.0, -182.0,
101-180.0, -178.0, -176.0, -174.0, -172.0, -170.0, -168.0, -166.0, -164.0, -162.0,
102-160.0, -158.0, -156.0, -154.0, -152.0, -150.0, -148.0, -146.0, -144.0, -142.0,
103-140.0, -138.0, -136.0, -134.0, -132.0, -130.0, -128.0, -126.0, -124.0, -122.0,
104-120.0, -118.0, -116.0, -114.0, -112.0, -110.0, -108.0, -106.0, -104.0, -102.0,
105-100.0, -98.0, -96.0, -94.0, -92.0, -90.0, -88.0, -86.0, -84.0, -82.0,
106-80.0, -78.0, -76.0, -74.0, -72.0, -70.0, -68.0, -66.0, -64.0, -62.0,
107-60.0, -58.0, -56.0, -54.0, -52.0, -50.0, -48.0, -46.0, -44.0, -42.0,
108-40.0, -38.0, -36.0, -34.0, -32.0, -30.0, -28.0, -26.0, -24.0, -22.0,
109-20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0,
110-1.0, -0.5, -0.2, -0.1, -0.05, 0.05, 0.1, 0.2, 0.5, 1.0,
111 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
112 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0,
113 42.0, 44.0, 46.0, 48.0, 50.0, 52.0, 54.0, 56.0, 58.0, 60.0,
114 62.0, 64.0, 66.0, 68.0, 70.0, 72.0, 74.0, 76.0, 78.0, 80.0,
115 82.0, 84.0, 86.0, 88.0, 90.0, 92.0, 94.0, 96.0, 98.0, 100.0,
116102.0, 104.0, 106.0, 108.0, 110.0, 112.0, 114.0, 116.0, 118.0, 120.0,
117122.0, 124.0, 126.0, 128.0, 130.0, 132.0, 134.0, 136.0, 138.0, 140.0,
118142.0, 144.0, 146.0, 148.0, 150.0, 152.0, 154.0, 156.0, 158.0, 160.0,
119162.0, 164.0, 166.0, 168.0, 170.0, 172.0, 174.0, 176.0, 178.0, 180.0,
120182.0, 184.0, 186.0, 188.0, 190.0, 192.0, 194.0, 196.0, 198.0, 200.0,
121202.0, 204.0, 206.0, 208.0, 210.0, 212.0, 214.0, 216.0, 218.0, 220.0,
122222.0, 224.0, 226.0, 228.0, 230.0, 232.0, 234.0, 236.0, 238.0, 240.0,
123242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 248.5, 249.0, 249.5 } ;
124
125
126
127AliTPCCorrection::AliTPCCorrection()
534fd34a 128 : TNamed("correction_unity","unity"),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 129{
130 //
131 // default constructor
132 //
133}
134
135AliTPCCorrection::AliTPCCorrection(const char *name,const char *title)
534fd34a 136: TNamed(name,title),fJLow(0),fKLow(0), fT1(1), fT2(1)
0116859c 137{
138 //
139 // default constructor, that set the name and title
140 //
141}
142
143AliTPCCorrection::~AliTPCCorrection() {
144 //
145 // virtual destructor
146 //
147}
148
149void AliTPCCorrection::CorrectPoint(Float_t x[],const Short_t roc) {
150 //
151 // Corrects the initial coordinates x (cartesian coordinates)
152 // according to the given effect (inherited classes)
153 // roc represents the TPC read out chamber (offline numbering convention)
154 //
155 Float_t dx[3];
156 GetCorrection(x,roc,dx);
157 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
158}
159
160void AliTPCCorrection::CorrectPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
161 //
162 // Corrects the initial coordinates x (cartesian coordinates) and stores the new
163 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
164 // roc represents the TPC read out chamber (offline numbering convention)
165 //
166 Float_t dx[3];
167 GetCorrection(x,roc,dx);
168 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
169}
170
171void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
172 //
173 // Distorts the initial coordinates x (cartesian coordinates)
174 // according to the given effect (inherited classes)
175 // roc represents the TPC read out chamber (offline numbering convention)
176 //
177 Float_t dx[3];
178 GetDistortion(x,roc,dx);
179 for (Int_t j=0;j<3;++j) x[j]+=dx[j];
180}
181
182void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
183 //
184 // Distorts the initial coordinates x (cartesian coordinates) and stores the new
185 // (distorted) coordinates in xp. The distortion is set according to the given effect (inherited classes)
186 // roc represents the TPC read out chamber (offline numbering convention)
187 //
188 Float_t dx[3];
189 GetDistortion(x,roc,dx);
190 for (Int_t j=0;j<3;++j) xp[j]=x[j]+dx[j];
191}
192
193void AliTPCCorrection::GetCorrection(const Float_t /*x*/[],const Short_t /*roc*/,Float_t dx[]) {
194 //
195 // This function delivers the correction values dx in respect to the inital coordinates x
196 // roc represents the TPC read out chamber (offline numbering convention)
197 // Note: The dx is overwritten by the inherited effectice class ...
198 //
199 for (Int_t j=0;j<3;++j) { dx[j]=0.; }
200}
201
202void AliTPCCorrection::GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]) {
203 //
204 // This function delivers the distortion values dx in respect to the inital coordinates x
205 // roc represents the TPC read out chamber (offline numbering convention)
206 //
207 GetCorrection(x,roc,dx);
208 for (Int_t j=0;j<3;++j) dx[j]=-dx[j];
209}
210
211void AliTPCCorrection::Init() {
212 //
213 // Initialization funtion (not used at the moment)
214 //
215}
216
e527a1b9 217void AliTPCCorrection::Update(const TTimeStamp &/*timeStamp*/) {
218 //
219 // Update function
220 //
221}
222
0116859c 223void AliTPCCorrection::Print(Option_t* /*option*/) const {
224 //
225 // Print function to check which correction classes are used
226 // option=="d" prints details regarding the setted magnitude
227 // option=="a" prints the C0 and C1 coefficents for calibration purposes
228 //
229 printf("TPC spacepoint correction: \"%s\"\n",GetTitle());
230}
231
534fd34a 232void AliTPCCorrection:: SetOmegaTauT1T2(Float_t /*omegaTau*/,Float_t t1,Float_t t2) {
0116859c 233 //
234 // Virtual funtion to pass the wt values (might become event dependent) to the inherited classes
235 // t1 and t2 represent the "effective omegaTau" corrections and were measured in a dedicated
236 // calibration run
237 //
534fd34a 238 fT1=t1;
239 fT2=t2;
240 //SetOmegaTauT1T2(omegaTau, t1, t2);
0116859c 241}
242
243TH2F* AliTPCCorrection::CreateHistoDRinXY(Float_t z,Int_t nx,Int_t ny) {
244 //
245 // Simple plot functionality.
246 // Returns a 2d hisogram which represents the corrections in radial direction (dr)
247 // in respect to position z within the XY plane.
248 // The histogramm has nx times ny entries.
249 //
250
251 TH2F *h=CreateTH2F("dr_xy",GetTitle(),"x [cm]","y [cm]","dr [cm]",
252 nx,-250.,250.,ny,-250.,250.);
253 Float_t x[3],dx[3];
254 x[2]=z;
255 Int_t roc=z>0.?0:18; // FIXME
256 for (Int_t iy=1;iy<=ny;++iy) {
257 x[1]=h->GetYaxis()->GetBinCenter(iy);
258 for (Int_t ix=1;ix<=nx;++ix) {
259 x[0]=h->GetXaxis()->GetBinCenter(ix);
260 GetCorrection(x,roc,dx);
261 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
262 if (90.<=r0 && r0<=250.) {
263 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
264 h->SetBinContent(ix,iy,r1-r0);
265 }
266 else
267 h->SetBinContent(ix,iy,0.);
268 }
269 }
270 return h;
271}
272
273TH2F* AliTPCCorrection::CreateHistoDRPhiinXY(Float_t z,Int_t nx,Int_t ny) {
274 //
275 // Simple plot functionality.
276 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
277 // in respect to position z within the XY plane.
278 // The histogramm has nx times ny entries.
279 //
280
281 TH2F *h=CreateTH2F("drphi_xy",GetTitle(),"x [cm]","y [cm]","drphi [cm]",
282 nx,-250.,250.,ny,-250.,250.);
283 Float_t x[3],dx[3];
284 x[2]=z;
285 Int_t roc=z>0.?0:18; // FIXME
286 for (Int_t iy=1;iy<=ny;++iy) {
287 x[1]=h->GetYaxis()->GetBinCenter(iy);
288 for (Int_t ix=1;ix<=nx;++ix) {
289 x[0]=h->GetXaxis()->GetBinCenter(ix);
290 GetCorrection(x,roc,dx);
291 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
292 if (90.<=r0 && r0<=250.) {
293 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
294 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
295
296 Float_t dphi=phi1-phi0;
297 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
298 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
299
300 h->SetBinContent(ix,iy,r0*dphi);
301 }
302 else
303 h->SetBinContent(ix,iy,0.);
304 }
305 }
306 return h;
307}
308
309TH2F* AliTPCCorrection::CreateHistoDRinZR(Float_t phi,Int_t nz,Int_t nr) {
310 //
311 // Simple plot functionality.
312 // Returns a 2d hisogram which represents the corrections in r direction (dr)
313 // in respect to angle phi within the ZR plane.
314 // The histogramm has nx times ny entries.
315 //
316 TH2F *h=CreateTH2F("dr_zr",GetTitle(),"z [cm]","r [cm]","dr [cm]",
317 nz,-250.,250.,nr,85.,250.);
318 Float_t x[3],dx[3];
319 for (Int_t ir=1;ir<=nr;++ir) {
320 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
321 x[0]=radius*TMath::Cos(phi);
322 x[1]=radius*TMath::Sin(phi);
323 for (Int_t iz=1;iz<=nz;++iz) {
324 x[2]=h->GetXaxis()->GetBinCenter(iz);
325 Int_t roc=x[2]>0.?0:18; // FIXME
326 GetCorrection(x,roc,dx);
327 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
328 Float_t r1=TMath::Sqrt((x[0]+dx[0])*(x[0]+dx[0])+(x[1]+dx[1])*(x[1]+dx[1]));
329 h->SetBinContent(iz,ir,r1-r0);
330 }
331 }
332 printf("SDF\n");
333 return h;
334
335}
336
337TH2F* AliTPCCorrection::CreateHistoDRPhiinZR(Float_t phi,Int_t nz,Int_t nr) {
338 //
339 // Simple plot functionality.
340 // Returns a 2d hisogram which represents the corrections in rphi direction (drphi)
341 // in respect to angle phi within the ZR plane.
342 // The histogramm has nx times ny entries.
343 //
344 TH2F *h=CreateTH2F("drphi_zr",GetTitle(),"z [cm]","r [cm]","drphi [cm]",
345 nz,-250.,250.,nr,85.,250.);
346 Float_t x[3],dx[3];
347 for (Int_t iz=1;iz<=nz;++iz) {
348 x[2]=h->GetXaxis()->GetBinCenter(iz);
349 Int_t roc=x[2]>0.?0:18; // FIXME
350 for (Int_t ir=1;ir<=nr;++ir) {
351 Float_t radius=h->GetYaxis()->GetBinCenter(ir);
352 x[0]=radius*TMath::Cos(phi);
353 x[1]=radius*TMath::Sin(phi);
354 GetCorrection(x,roc,dx);
355 Float_t r0=TMath::Sqrt((x[0] )*(x[0] )+(x[1] )*(x[1] ));
356 Float_t phi0=TMath::ATan2(x[1] ,x[0] );
357 Float_t phi1=TMath::ATan2(x[1]+dx[1],x[0]+dx[0]);
358
359 Float_t dphi=phi1-phi0;
360 if (dphi<TMath::Pi()) dphi+=TMath::TwoPi();
361 if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
362
363 h->SetBinContent(iz,ir,r0*dphi);
364 }
365 }
366 return h;
367}
368
369TH2F* AliTPCCorrection::CreateTH2F(const char *name,const char *title,
370 const char *xlabel,const char *ylabel,const char *zlabel,
371 Int_t nbinsx,Double_t xlow,Double_t xup,
372 Int_t nbinsy,Double_t ylow,Double_t yup) {
373 //
374 // Helper function to create a 2d histogramm of given size
375 //
376
377 TString hname=name;
378 Int_t i=0;
379 if (gDirectory) {
380 while (gDirectory->FindObject(hname.Data())) {
381 hname =name;
382 hname+="_";
383 hname+=i;
384 ++i;
385 }
386 }
387 TH2F *h=new TH2F(hname.Data(),title,
388 nbinsx,xlow,xup,
389 nbinsy,ylow,yup);
390 h->GetXaxis()->SetTitle(xlabel);
391 h->GetYaxis()->SetTitle(ylabel);
392 h->GetZaxis()->SetTitle(zlabel);
393 h->SetStats(0);
394 return h;
395}
396
397
398// Simple Interpolation functions: e.g. with bi(tri)cubic interpolations (not yet in TH2 and TH3)
399
400void AliTPCCorrection::Interpolate2DEdistortion( const Int_t order, const Double_t r, const Double_t z,
b1f0a2a5 401 const Double_t er[kNZ][kNR], Double_t &erValue ) {
0116859c 402 //
403 // Interpolate table - 2D interpolation
404 //
b1f0a2a5 405 Double_t saveEr[10] ;
0116859c 406
407 Search( kNZ, fgkZList, z, fJLow ) ;
408 Search( kNR, fgkRList, r, fKLow ) ;
409 if ( fJLow < 0 ) fJLow = 0 ; // check if out of range
410 if ( fKLow < 0 ) fKLow = 0 ;
411 if ( fJLow + order >= kNZ - 1 ) fJLow = kNZ - 1 - order ;
412 if ( fKLow + order >= kNR - 1 ) fKLow = kNR - 1 - order ;
413
414 for ( Int_t j = fJLow ; j < fJLow + order + 1 ; j++ ) {
b1f0a2a5 415 saveEr[j-fJLow] = Interpolate( &fgkRList[fKLow], &er[j][fKLow], order, r ) ;
0116859c 416 }
b1f0a2a5 417 erValue = Interpolate( &fgkZList[fJLow], saveEr, order, z ) ;
0116859c 418
419}
420
421
422Double_t AliTPCCorrection::Interpolate( const Double_t xArray[], const Double_t yArray[],
b1f0a2a5 423 const Int_t order, const Double_t x ) {
0116859c 424 //
425 // Interpolate function Y(x) using linear (order=1) or quadratic (order=2) interpolation.
426 //
427
428 Double_t y ;
429 if ( order == 2 ) { // Quadratic Interpolation = 2
430 y = (x-xArray[1]) * (x-xArray[2]) * yArray[0] / ( (xArray[0]-xArray[1]) * (xArray[0]-xArray[2]) ) ;
431 y += (x-xArray[2]) * (x-xArray[0]) * yArray[1] / ( (xArray[1]-xArray[2]) * (xArray[1]-xArray[0]) ) ;
432 y += (x-xArray[0]) * (x-xArray[1]) * yArray[2] / ( (xArray[2]-xArray[0]) * (xArray[2]-xArray[1]) ) ;
433 } else { // Linear Interpolation = 1
434 y = yArray[0] + ( yArray[1]-yArray[0] ) * ( x-xArray[0] ) / ( xArray[1] - xArray[0] ) ;
435 }
436
437 return (y);
438
439}
440
441
b1f0a2a5 442void AliTPCCorrection::Search( const Int_t n, const Double_t xArray[], const Double_t x, Int_t &low ) {
0116859c 443 //
444 // Search an ordered table by starting at the most recently used point
445 //
446
447 Long_t middle, high ;
448 Int_t ascend = 0, increment = 1 ;
449
450 if ( xArray[n-1] >= xArray[0] ) ascend = 1 ; // Ascending ordered table if true
451
452 if ( low < 0 || low > n-1 ) {
453 low = -1 ; high = n ;
454 } else { // Ordered Search phase
455 if ( (Int_t)( x >= xArray[low] ) == ascend ) {
456 if ( low == n-1 ) return ;
457 high = low + 1 ;
458 while ( (Int_t)( x >= xArray[high] ) == ascend ) {
459 low = high ;
460 increment *= 2 ;
461 high = low + increment ;
462 if ( high > n-1 ) { high = n ; break ; }
463 }
464 } else {
465 if ( low == 0 ) { low = -1 ; return ; }
466 high = low - 1 ;
467 while ( (Int_t)( x < xArray[low] ) == ascend ) {
468 high = low ;
469 increment *= 2 ;
470 if ( increment >= high ) { low = -1 ; break ; }
471 else low = high - increment ;
472 }
473 }
474 }
475
476 while ( (high-low) != 1 ) { // Binary Search Phase
477 middle = ( high + low ) / 2 ;
478 if ( (Int_t)( x >= xArray[middle] ) == ascend )
479 low = middle ;
480 else
481 high = middle ;
482 }
483
484 if ( x == xArray[n-1] ) low = n-2 ;
485 if ( x == xArray[0] ) low = 0 ;
486
487}
488
cf5b0aa0 489
b1f0a2a5 490AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir, TTreeSRedirector * const pcstream){
cf5b0aa0 491 //
492 // Fit the track parameters - without and with distortion
493 // 1. Space points in the TPC are simulated along the trajectory
494 // 2. Space points distorted
495 // 3. Fits the non distorted and distroted track to the reference plane at refX
496 // 4. For visualization and debugging purposes the space points and tracks can be stored in the tree - using the TTreeSRedirector functionality
497 //
498 // trackIn - input track parameters
499 // refX - reference X to fit the track
500 // dir - direction - out=1 or in=-1
501 // pcstream - debug streamer to check the results
502 //
cad404e1 503 // see AliExternalTrackParam.h documentation:
504 // track1.fP[0] - local y (rphi)
505 // track1.fP[1] - z
506 // track1.fP[2] - sinus of local inclination angle
507 // track1.fP[3] - tangent of deep angle
508 // track1.fP[4] - 1/pt
cf5b0aa0 509 AliTPCROC * roc = AliTPCROC::Instance();
510 const Int_t npoints0=roc->GetNRows(0)+roc->GetNRows(36);
511 const Double_t kRTPC0 =roc->GetPadRowRadii(0,0);
512 const Double_t kRTPC1 =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
513
514 const Double_t kMaxSnp = 0.85;
515 const Double_t kSigmaY=0.1;
516 const Double_t kSigmaZ=0.1;
517 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
518
be67055b 519 AliExternalTrackParam track(trackIn); //
cf5b0aa0 520 // generate points
521 AliTrackPointArray pointArray0(npoints0);
522 AliTrackPointArray pointArray1(npoints0);
523 Double_t xyz[3];
8b63d99c 524 AliTrackerBase::PropagateTrackToBxByBz(&track,kRTPC0,kMass,3,kTRUE,kMaxSnp);
cf5b0aa0 525 //
526 // simulate the track
527 Int_t npoints=0;
528 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
529 for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
8b63d99c 530 AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
cf5b0aa0 531 track.GetXYZ(xyz);
7f4cb119 532 xyz[0]+=gRandom->Gaus(0,0.005);
533 xyz[1]+=gRandom->Gaus(0,0.005);
534 xyz[2]+=gRandom->Gaus(0,0.005);
cf5b0aa0 535 AliTrackPoint pIn0; // space point
536 AliTrackPoint pIn1;
ffab0c37 537 Int_t sector= (xyz[2]>0)? 0:18;
cf5b0aa0 538 pointArray0.GetPoint(pIn0,npoints);
539 pointArray1.GetPoint(pIn1,npoints);
540 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
541 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
ffab0c37 542 DistortPoint(distPoint, sector);
cf5b0aa0 543 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
544 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
545 //
546 track.Rotate(alpha);
547 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
548 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
549 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
550 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
551 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
552 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
553 pointArray0.AddPoint(npoints, &pIn0);
554 pointArray1.AddPoint(npoints, &pIn1);
555 npoints++;
556 if (npoints>=npoints0) break;
557 }
7f4cb119 558 if (npoints<npoints0/2) return 0;
cf5b0aa0 559 //
560 // refit track
561 //
562 AliExternalTrackParam *track0=0;
563 AliExternalTrackParam *track1=0;
564 AliTrackPoint point1,point2,point3;
565 if (dir==1) { //make seed inner
566 pointArray0.GetPoint(point1,1);
4486a91f 567 pointArray0.GetPoint(point2,30);
568 pointArray0.GetPoint(point3,60);
cf5b0aa0 569 }
570 if (dir==-1){ //make seed outer
4486a91f 571 pointArray0.GetPoint(point1,npoints-60);
572 pointArray0.GetPoint(point2,npoints-30);
cf5b0aa0 573 pointArray0.GetPoint(point3,npoints-1);
574 }
575 track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
576 track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
577
cf5b0aa0 578 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
8b63d99c 579 Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
cf5b0aa0 580 //
581 AliTrackPoint pIn0;
582 AliTrackPoint pIn1;
583 pointArray0.GetPoint(pIn0,ipoint);
584 pointArray1.GetPoint(pIn1,ipoint);
585 AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha()); // rotate to the local frame - non distoted point
586 AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha()); // rotate to the local frame - distorted point
587 //
b322e06a 588 AliTrackerBase::PropagateTrackToBxByBz(track0,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
589 AliTrackerBase::PropagateTrackToBxByBz(track1,prot0.GetX(),kMass,3,kFALSE,kMaxSnp);
8b63d99c 590 track.GetXYZ(xyz); // distorted track also propagated to the same reference radius
cf5b0aa0 591 //
592 Double_t pointPos[2]={0,0};
593 Double_t pointCov[3]={0,0,0};
594 pointPos[0]=prot0.GetY();//local y
595 pointPos[1]=prot0.GetZ();//local z
596 pointCov[0]=prot0.GetCov()[3];//simay^2
597 pointCov[1]=prot0.GetCov()[4];//sigmayz
598 pointCov[2]=prot0.GetCov()[5];//sigmaz^2
599 track0->Update(pointPos,pointCov);
600 //
8b63d99c 601 Double_t deltaX=prot1.GetX()-prot0.GetX(); // delta X
602 Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
603 Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
604
7f4cb119 605 pointPos[0]=prot1.GetY()+deltaYX;//local y is sign correct?
606 pointPos[1]=prot1.GetZ()+deltaZX;//local z is sign correct?
cf5b0aa0 607 pointCov[0]=prot1.GetCov()[3];//simay^2
608 pointCov[1]=prot1.GetCov()[4];//sigmayz
609 pointCov[2]=prot1.GetCov()[5];//sigmaz^2
610 track1->Update(pointPos,pointCov);
611 }
612
8b63d99c 613 AliTrackerBase::PropagateTrackToBxByBz(track0,refX,kMass,2.,kTRUE,kMaxSnp);
cf5b0aa0 614 track1->Rotate(track0->GetAlpha());
be67055b 615 track1->PropagateTo(track0->GetX(),AliTrackerBase::GetBz());
cf5b0aa0 616
cad404e1 617 if (pcstream) (*pcstream)<<Form("fitDistort%s",GetName())<<
cf5b0aa0 618 "point0.="<<&pointArray0<< // points
619 "point1.="<<&pointArray1<< // distorted points
620 "trackIn.="<<&track<< // original track
621 "track0.="<<track0<< // fitted track
622 "track1.="<<track1<< // fitted distorted track
623 "\n";
be67055b 624 new(&trackIn) AliExternalTrackParam(*track0);
cf5b0aa0 625 delete track0;
626 return track1;
627}
628
629
ffab0c37 630
631
632
633TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
634 //
635 // create the distortion tree on a mesh with granularity given by step
636 // return the tree with distortions at given position
637 // Map is created on the mesh with given step size
638 //
639 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("correction%s.root",GetName()));
640 Float_t xyz[3];
641 for (Double_t x= -250; x<250; x+=step){
642 for (Double_t y= -250; y<250; y+=step){
643 Double_t r = TMath::Sqrt(x*x+y*y);
644 if (r<80) continue;
645 if (r>250) continue;
646 for (Double_t z= -250; z<250; z+=step){
647 Int_t roc=(z>0)?0:18;
648 xyz[0]=x;
649 xyz[1]=y;
650 xyz[2]=z;
651 Double_t phi = TMath::ATan2(y,x);
652 DistortPoint(xyz,roc);
653 Double_t r1 = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
654 Double_t phi1 = TMath::ATan2(xyz[1],xyz[0]);
655 if ((phi1-phi)>TMath::Pi()) phi1-=TMath::Pi();
656 if ((phi1-phi)<-TMath::Pi()) phi1+=TMath::Pi();
657 Double_t dx = xyz[0]-x;
658 Double_t dy = xyz[1]-y;
659 Double_t dz = xyz[2]-z;
660 Double_t dr=r1-r;
661 Double_t drphi=(phi1-phi)*r;
662 (*pcstream)<<"distortion"<<
663 "x="<<x<< // original position
664 "y="<<y<<
665 "z="<<z<<
666 "r="<<r<<
667 "phi="<<phi<<
668 "x1="<<xyz[0]<< // distorted position
669 "y1="<<xyz[1]<<
670 "z1="<<xyz[2]<<
671 "r1="<<r1<<
672 "phi1="<<phi1<<
673 //
674 "dx="<<dx<< // delta position
675 "dy="<<dy<<
676 "dz="<<dz<<
677 "dr="<<dr<<
678 "drphi="<<drphi<<
679 "\n";
680 }
681 }
682 }
683 delete pcstream;
684 TFile f(Form("correction%s.root",GetName()));
685 TTree * tree = (TTree*)f.Get("distortion");
686 TTree * tree2= tree->CopyTree("1");
687 tree2->SetName(Form("dist%s",GetName()));
688 tree2->SetDirectory(0);
689 delete tree;
690 return tree2;
691}
692
693
694
be67055b 695
b1f0a2a5 696void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
be67055b 697 //
698 // Make a fit tree:
699 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
700 // calculates partial distortions
701 // Partial distortion is stored in the resulting tree
702 // Output is storred in the file distortion_<dettype>_<partype>.root
703 // Partial distortion is stored with the name given by correction name
704 //
705 //
706 // Parameters of function:
707 // input - input tree
708 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex
709 // ppype - parameter type
710 // corrArray - array with partial corrections
711 // step - skipe entries - if 1 all entries processed - it is slow
712 // debug 0 if debug on also space points dumped - it is slow
b322e06a 713 const Double_t kMaxSnp = 0.85;
714 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
715 // const Double_t kB2C=-0.299792458e-3;
7f4cb119 716 const Int_t kMinEntries=50;
be67055b 717 Double_t phi,theta, snp, mean,rms, entries;
718 tinput->SetBranchAddress("theta",&theta);
719 tinput->SetBranchAddress("phi", &phi);
720 tinput->SetBranchAddress("snp",&snp);
721 tinput->SetBranchAddress("mean",&mean);
722 tinput->SetBranchAddress("rms",&rms);
723 tinput->SetBranchAddress("entries",&entries);
724 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
725 //
726 Int_t nentries=tinput->GetEntries();
727 Int_t ncorr=corrArray->GetEntries();
7f4cb119 728 Double_t corrections[100]={0}; //
be67055b 729 Double_t tPar[5];
730 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
731 Double_t refX=0;
732 Int_t dir=0;
7f4cb119 733 if (dtype==0) {refX=85.; dir=-1;}
734 if (dtype==1) {refX=275.; dir=1;}
735 if (dtype==2) {refX=85.; dir=-1;}
736 if (dtype==3) {refX=360.; dir=-1;}
be67055b 737 //
738 for (Int_t ientry=0; ientry<nentries; ientry+=step){
739 tinput->GetEntry(ientry);
7f4cb119 740 if (TMath::Abs(snp)>kMaxSnp) continue;
be67055b 741 tPar[0]=0;
742 tPar[1]=theta*refX;
743 tPar[2]=snp;
744 tPar[3]=theta;
4486a91f 745 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
8b63d99c 746 Double_t bz=AliTrackerBase::GetBz();
4486a91f 747 if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
748 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
7f4cb119 749 AliExternalTrackParam track(refX,phi,tPar,cov);
750 Double_t xyz[3];
751 track.GetXYZ(xyz);
752 Int_t id=0;
753 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
754
be67055b 755 (*pcstream)<<"fit"<<
8b63d99c 756 "bz="<<bz<< // magnetic filed used
be67055b 757 "dtype="<<dtype<< // detector match type
758 "ptype="<<ptype<< // parameter type
759 "theta="<<theta<< // theta
760 "phi="<<phi<< // phi
761 "snp="<<snp<< // snp
762 "mean="<<mean<< // mean dist value
763 "rms="<<rms<< // rms
7f4cb119 764 "gx="<<xyz[0]<< // global position at reference
765 "gy="<<xyz[1]<< // global position at reference
766 "gz="<<xyz[2]<< // global position at reference
767 "dRrec="<<dRrec<< // delta Radius in reconstruction
768 "id="<<id<< // track id
be67055b 769 "entries="<<entries;// number of entries in bin
770 //
771 for (Int_t icorr=0; icorr<ncorr; icorr++) {
772 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
773 corrections[icorr]=0;
774 if (entries>kMinEntries){
775 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
776 AliExternalTrackParam *trackOut = 0;
777 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
778 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
7f4cb119 779 if (dtype==0) {refX=85.; dir=-1;}
780 if (dtype==1) {refX=275.; dir=1;}
b1f0a2a5 781 if (dtype==2) {refX=0; dir=-1;}
7f4cb119 782 if (dtype==3) {refX=360.; dir=-1;}
b1f0a2a5 783 //
7f4cb119 784 if (trackOut){
785 AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
786 trackOut->Rotate(trackIn.GetAlpha());
787 trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
788 //
789 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
790 delete trackOut;
791 }else{
792 corrections[icorr]=0;
793 }
be67055b 794 }
7f4cb119 795 Double_t dRdummy=0;
be67055b 796 (*pcstream)<<"fit"<<
7f4cb119 797 Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
798 Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
799 // for points it is neccessary
be67055b 800 }
801 (*pcstream)<<"fit"<<"\n";
802 }
803 delete pcstream;
804}
805
806
807
7f4cb119 808void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
809 //
810 // Make a laser fit tree for global minimization
811 //
812 const Double_t cutErrY=0.1;
813 const Double_t cutErrZ=0.1;
814 const Double_t kEpsilon=0.00000001;
815 TVectorD *vecdY=0;
816 TVectorD *vecdZ=0;
817 TVectorD *veceY=0;
818 TVectorD *veceZ=0;
819 AliTPCLaserTrack *ltr=0;
820 AliTPCLaserTrack::LoadTracks();
821 tree->SetBranchAddress("dY.",&vecdY);
822 tree->SetBranchAddress("dZ.",&vecdZ);
823 tree->SetBranchAddress("eY.",&veceY);
824 tree->SetBranchAddress("eZ.",&veceZ);
825 tree->SetBranchAddress("LTr.",&ltr);
826 Int_t entries= tree->GetEntries();
827 TTreeSRedirector *pcstream= new TTreeSRedirector("distortion4_0.root");
828 Double_t bz=AliTrackerBase::GetBz();
829 //
830
831 for (Int_t ientry=0; ientry<entries; ientry++){
832 tree->GetEntry(ientry);
833 if (!ltr->GetVecGX()){
834 ltr->UpdatePoints();
835 }
836 TVectorD * delta= (itype==0)? vecdY:vecdZ;
837 TVectorD * err= (itype==0)? veceY:veceZ;
838
839 for (Int_t irow=0; irow<159; irow++){
840 Int_t nentries = 1000;
841 if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
842 if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
843 Int_t dtype=4;
844 Double_t phi =(*ltr->GetVecPhi())[irow];
845 Double_t theta =ltr->GetTgl();
846 Double_t mean=delta->GetMatrixArray()[irow];
847 Double_t gx=0,gy=0,gz=0;
848 Double_t snp = (*ltr->GetVecP2())[irow];
849 Double_t rms = 0.1+err->GetMatrixArray()[irow];
850 gx = (*ltr->GetVecGX())[irow];
851 gy = (*ltr->GetVecGY())[irow];
852 gz = (*ltr->GetVecGZ())[irow];
853 Int_t bundle= ltr->GetBundle();
854 Double_t dRrec=0;
855 //
856 // get delta R used in reconstruction
857 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
858 AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
859 const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
860 Double_t xyz0[3]={gx,gy,gz};
861 Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
862 //
863 // old ExB correction
864 //
865 if(recoParam&&recoParam->GetUseExBCorrection()) {
866 Double_t xyz1[3]={gx,gy,gz};
867 calib->GetExB()->Correct(xyz0,xyz1);
868 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
869 dRrec=oldR-newR;
870 }
871 if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
872 Float_t xyz1[3]={gx,gy,gz};
873 Int_t sector=(gz>0)?0:18;
874 correction->CorrectPoint(xyz1, sector);
875 Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
876 dRrec=oldR-newR;
877 }
878
879
880 (*pcstream)<<"fit"<<
881 "bz="<<bz<< // magnetic filed used
882 "dtype="<<dtype<< // detector match type
883 "ptype="<<itype<< // parameter type
884 "theta="<<theta<< // theta
885 "phi="<<phi<< // phi
886 "snp="<<snp<< // snp
887 "mean="<<mean<< // mean dist value
888 "rms="<<rms<< // rms
889 "gx="<<gx<< // global position
890 "gy="<<gy<< // global position
891 "gz="<<gz<< // global position
892 "dRrec="<<dRrec<< // delta Radius in reconstruction
893 "id="<<bundle<< //bundle
894 "entries="<<nentries;// number of entries in bin
895 //
896 //
897 Double_t ky = TMath::Tan(TMath::ASin(snp));
898 Int_t ncorr = corrArray->GetEntries();
899 Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
900 Double_t phi0 = TMath::ATan2(gy,gx);
901 Double_t distortions[1000]={0};
902 Double_t distortionsR[1000]={0};
903 for (Int_t icorr=0; icorr<ncorr; icorr++) {
904 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
905 Float_t distPoint[3]={gx,gy,gz};
906 Int_t sector= (gz>0)? 0:18;
907 if (r0>80){
908 corr->DistortPoint(distPoint, sector);
909 }
910 Double_t value=distPoint[2]-gz;
911 if (itype==0){
912 Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
913 Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
914 Double_t drphi= r0*(phi1-phi0);
915 Double_t dr = r1-r0;
916 distortions[icorr] = drphi-ky*dr;
917 distortionsR[icorr] = dr;
918 }
919 (*pcstream)<<"fit"<<
920 Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
921 Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
922 }
923 (*pcstream)<<"fit"<<"\n";
924 }
925 }
926 delete pcstream;
927}
928
929
be67055b 930
b1f0a2a5 931void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
8b63d99c 932 //
933 // make a distortion map out ou fthe residual histogram
934 // Results are written to the debug streamer - pcstream
935 // Parameters:
936 // his0 - input (4D) residual histogram
937 // pcstream - file to write the tree
938 // run - run number
939 // marian.ivanov@cern.ch
940 const Int_t kMinEntries=50;
941 Int_t nbins1=his0->GetAxis(1)->GetNbins();
942 Int_t first1=his0->GetAxis(1)->GetFirst();
943 Int_t last1 =his0->GetAxis(1)->GetLast();
944 //
945 Double_t bz=AliTrackerBase::GetBz();
946 Int_t idim[4]={0,1,2,3};
947 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
948 //
949 his0->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
950 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
951 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
952 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
953 Int_t first3 = his1->GetAxis(3)->GetFirst();
954 Int_t last3 = his1->GetAxis(3)->GetLast();
955 //
956
957 for (Int_t ibin3=first3-1; ibin3<last3; ibin3+=1){ // axis 3 - local angle
958 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
959 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
960 if (ibin3<first3) {
961 his1->GetAxis(3)->SetRangeUser(-1,1);
962 x3=0;
963 }
964 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
965 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
966 Int_t first2 = his3->GetAxis(2)->GetFirst();
967 Int_t last2 = his3->GetAxis(2)->GetLast();
968 //
969 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){
970 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
971 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
972 TH1 * hisDelta = his3->Projection(0);
973 //
974 Double_t entries = hisDelta->GetEntries();
975 Double_t mean=0, rms=0;
976 if (entries>kMinEntries){
977 mean = hisDelta->GetMean();
978 rms = hisDelta->GetRMS();
979 }
980 (*pcstream)<<hname<<
981 "run="<<run<<
982 "bz="<<bz<<
983 "theta="<<x1<<
984 "phi="<<x2<<
985 "snp="<<x3<<
986 "entries="<<entries<<
987 "mean="<<mean<<
988 "rms="<<rms<<
989 "\n";
990 delete hisDelta;
991 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
992 }
993 delete his3;
994 }
995 delete his1;
996 }
997}
998
999
1000
1001
1002
ffab0c37 1003void AliTPCCorrection::StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment){
1004 //
1005 // Store object in the OCDB
1006 // By default the object is stored in the current directory
1007 // default comment consit of user name and the date
1008 //
1009 TString ocdbStorage="";
1010 ocdbStorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
1011 AliCDBMetaData *metaData= new AliCDBMetaData();
1012 metaData->SetObjectClassName("AliTPCCorrection");
1013 metaData->SetResponsible("Marian Ivanov");
1014 metaData->SetBeamPeriod(1);
1015 metaData->SetAliRootVersion("05-25-01"); //root version
1016 TString userName=gSystem->GetFromPipe("echo $USER");
1017 TString date=gSystem->GetFromPipe("date");
1018
1019 if (!comment) metaData->SetComment(Form("Space point distortion calibration\n User: %s\n Data%s",userName.Data(),date.Data()));
1020 if (comment) metaData->SetComment(comment);
1021 AliCDBId* id1=NULL;
1022 id1=new AliCDBId("TPC/Calib/Correction", startRun, endRun);
1023 AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1024 gStorage->Put(this, (*id1), metaData);
1025}
1026